diff --git a/sPuReMD/src/analyze.c b/sPuReMD/src/analyze.c
index 3f63f09d9c151c304c57be1aca93ace4d0885147..1d0b9c36185f1b9b98c8a356c0c6eec59f845e61 100644
--- a/sPuReMD/src/analyze.c
+++ b/sPuReMD/src/analyze.c
@@ -319,9 +319,8 @@ static void Analyze_Molecules( reax_system *system, control_params *control,
static void Report_Bond_Change( reax_system *system, control_params *control,
- static_storage *workspace, reax_list *old_bonds,
- reax_list *new_bonds, int a1, int a2, int flag,
- FILE *fout )
+ static_storage *workspace, reax_list *old_bonds,
+ reax_list *new_bonds, int a1, int a2, int flag, FILE *fout )
{
int i;
int rev1, rev2;
@@ -331,22 +330,36 @@ static void Report_Bond_Change( reax_system *system, control_params *control,
rev1 = workspace->orig_id[a1];
rev2 = workspace->orig_id[a2];
- if ( !strcmp( system->atoms[a1].name, " Si" ) ||
- !strcmp( system->atoms[a1].name, " O" ) )
+ if ( !strncmp( system->atoms[a1].name, " Si", 8 ) ||
+ !strncmp( system->atoms[a1].name, " O", 8 ) )
+ {
mol1 = 0;
- else mol1 = 1;
+ }
+ else
+ {
+ mol1 = 1;
+ }
- if ( !strcmp( system->atoms[a2].name, " Si" ) ||
- !strcmp( system->atoms[a2].name, " O" ) )
+ if ( !strncmp( system->atoms[a2].name, " Si", 8 ) ||
+ !strncmp( system->atoms[a2].name, " O", 8 ) )
+ {
mol2 = 0;
- else mol2 = 1;
-
+ }
+ else
+ {
+ mol2 = 1;
+ }
if ( mol1 == 0 && mol2 == 0 ) // silica-silica
{
if ( flag )
+ {
fprintf( fout, "silica bond formation:" );
- else fprintf( fout, "silica bond breakage :" );
+ }
+ else
+ {
+ fprintf( fout, "silica bond breakage :" );
+ }
fprintf( fout, "%5d(%s)-%5d(%s)\n",
rev1, system->atoms[a1].name, rev2, system->atoms[a2].name );
@@ -354,8 +367,13 @@ static void Report_Bond_Change( reax_system *system, control_params *control,
else if ( mol1 == 1 && mol2 == 1 ) // water-water
{
if ( flag )
+ {
fprintf( fout, "water bond formation:" );
- else fprintf( fout, "water bond breakage :" );
+ }
+ else
+ {
+ fprintf( fout, "water bond breakage :" );
+ }
fprintf( fout, "%5d(%s)-%5d(%s)\n",
rev1, system->atoms[a1].name, rev2, system->atoms[a2].name );
@@ -363,26 +381,39 @@ static void Report_Bond_Change( reax_system *system, control_params *control,
else // water-silica!
{
if ( flag )
+ {
fprintf( fout, "SILICA-WATER bond formation:" );
- else fprintf( fout, "SILICA-WATER bond breakage :" );
+ }
+ else
+ {
+ fprintf( fout, "SILICA-WATER bond breakage :" );
+ }
fprintf( fout, "%5d(%s)-%5d(%s)\n",
rev1, system->atoms[a1].name, rev2, system->atoms[a2].name );
fprintf( fout, "%5d(%s) was connected to:", rev1, system->atoms[a1].name );
for ( i = Start_Index(a1, old_bonds); i < End_Index(a1, old_bonds); ++i )
+ {
if ( old_bonds->select.bond_list[i].bo_data.BO >= control->bg_cut )
+ {
fprintf( fout, " %5d(%s)",
workspace->orig_id[ old_bonds->select.bond_list[i].nbr ],
system->atoms[ old_bonds->select.bond_list[i].nbr ].name );
+ }
+ }
fprintf( fout, "\n" );
fprintf( fout, "%5d(%s) was connected to:", rev2, system->atoms[a2].name );
for ( i = Start_Index(a2, old_bonds); i < End_Index(a2, old_bonds); ++i )
+ {
if ( old_bonds->select.bond_list[i].bo_data.BO >= control->bg_cut )
+ {
fprintf( fout, " %5d(%s)",
workspace->orig_id[ old_bonds->select.bond_list[i].nbr ],
system->atoms[ old_bonds->select.bond_list[i].nbr ].name );
+ }
+ }
fprintf( fout, "\n" );
}
}
@@ -390,8 +421,8 @@ static void Report_Bond_Change( reax_system *system, control_params *control,
/* ASSUMPTION: Bond lists are sorted */
static void Compare_Bonding( int atom, reax_system *system, control_params *control,
- static_storage *workspace, reax_list *old_bonds,
- reax_list *new_bonds, FILE *fout )
+ static_storage *workspace, reax_list *old_bonds,
+ reax_list *new_bonds, FILE *fout )
{
int oldp, newp;
@@ -409,25 +440,29 @@ static void Compare_Bonding( int atom, reax_system *system, control_params *cont
fprintf( fout, "\n" ); */
for ( oldp = Start_Index( atom, old_bonds );
- oldp < End_Index( atom, old_bonds ) &&
- old_bonds->select.bond_list[oldp].nbr < atom;
- ++oldp );
+ oldp < End_Index( atom, old_bonds )
+ && old_bonds->select.bond_list[oldp].nbr < atom; ++oldp )
+ ;
for ( newp = Start_Index( atom, new_bonds );
- newp < End_Index( atom, new_bonds ) &&
- new_bonds->select.bond_list[newp].nbr < atom;
- ++newp );
+ newp < End_Index( atom, new_bonds )
+ && new_bonds->select.bond_list[newp].nbr < atom; ++newp )
+ ;
while ( oldp < End_Index( atom, old_bonds ) ||
newp < End_Index( atom, new_bonds ) )
{
while ( oldp < End_Index( atom, old_bonds ) &&
old_bonds->select.bond_list[oldp].bo_data.BO < control->bg_cut )
+ {
++oldp;
+ }
while ( newp < End_Index( atom, new_bonds ) &&
new_bonds->select.bond_list[newp].bo_data.BO < control->bg_cut )
+ {
++newp;
+ }
/*fprintf( fout, "%d, oldp: %d - %d: %f newp: %d - %d: %f",
atom, oldp, old_bonds->select.bond_list[oldp].nbr,
@@ -461,7 +496,10 @@ static void Compare_Bonding( int atom, reax_system *system, control_params *cont
++newp;
}
else
- ++newp, ++oldp;
+ {
+ ++newp;
+ ++oldp;
+ }
}
else
/* there is no other bond in the new list */
@@ -483,6 +521,7 @@ static void Compare_Bonding( int atom, reax_system *system, control_params *cont
{
/* there are no more bonds in old_bond list */
if ( newp < End_Index( atom, new_bonds ) )
+ {
/* there is at least one other bond in the new list */
while ( newp < End_Index( atom, new_bonds ) )
{
@@ -491,12 +530,12 @@ static void Compare_Bonding( int atom, reax_system *system, control_params *cont
// fprintf( fout, "%5d-%5d bond formed\n",
// atom, new_bonds->select.bond_list[newp].nbr );
Report_Bond_Change( system, control, workspace,
- old_bonds, new_bonds, atom,
- new_bonds->select.bond_list[newp].nbr, 1,
- fout );
+ old_bonds, new_bonds, atom,
+ new_bonds->select.bond_list[newp].nbr, 1, fout );
}
++newp;
}
+ }
else
{
/* there is no other bond in the new list, either --
@@ -516,7 +555,9 @@ static void Visit_Bonds( int atom, int *mark, int *type, reax_system *system,
mark[atom] = 1;
t = system->atoms[atom].type;
if ( ignore && control->ignore[t] )
+ {
return;
+ }
type[t]++;
start = Start_Index( atom, bonds );
@@ -552,6 +593,7 @@ static void Analyze_Fragments( reax_system *system, control_params *control,
memset( mark, 0, system->N * sizeof(int) );
for ( atom = 0; atom < system->N; ++atom )
+ {
if ( !mark[atom] )
{
/* discover a new fragment */
@@ -563,21 +605,24 @@ static void Analyze_Fragments( reax_system *system, control_params *control,
/* check if a similar fragment already exists */
flag = 0;
for ( i = 0; i < num_fragment_types; ++i )
- if ( !strcmp( fragments[i], fragment ) )
+ {
+ if ( !strncmp( fragments[i], fragment, MAX_ATOM_TYPES ) )
{
++fragment_count[i];
flag = 1;
break;
}
+ }
if ( flag == 0 )
{
/* it is a new one, add to the fragments list */
- strcpy( fragments[num_fragment_types], fragment );
+ strncpy( fragments[num_fragment_types], fragment, MAX_ATOM_TYPES );
fragment_count[num_fragment_types] = 1;
++num_fragment_types;
}
}
+ }
/* output the results of fragment analysis */
for ( i = 0; i < num_fragment_types; ++i )
@@ -672,8 +717,11 @@ static void Analyze_Silica( reax_system *system, control_params *control,
SI_O_SI_count = 0;
for ( j = 0; j < system->N; ++j )
+ {
if ( system->atoms[j].type == O_ATOM || system->atoms[j].type == SI_ATOM )
+ {
for ( pi = Start_Index(j, new_bonds); pi < End_Index(j, new_bonds); ++pi )
+ {
if ( new_bonds->select.bond_list[pi].bo_data.BO >= control->bg_cut )
{
i = new_bonds->select.bond_list[pi].nbr;
@@ -713,6 +761,9 @@ static void Analyze_Silica( reax_system *system, control_params *control,
}
}
}
+ }
+ }
+ }
fprintf( fout, "\nAverage O-Si-O angle: %8.2f\n",
RAD2DEG(O_SI_O / O_SI_O_count) );
@@ -735,8 +786,8 @@ static int Get_Type_of_Molecule( molecule *m )
static void Calculate_Dipole_Moment( reax_system *system, control_params *control,
- simulation_data *data, static_storage *workspace,
- reax_list *bonds, FILE *fout )
+ simulation_data *data, static_storage *workspace,
+ reax_list *bonds, FILE *fout )
{
int i, atom, count;
molecule m;
@@ -744,12 +795,11 @@ static void Calculate_Dipole_Moment( reax_system *system, control_params *contro
rvec tmpvec, mu;
int *mark = workspace->mark;
- //fprintf( fout, "starting dipole moment calculations...\n" );
-
mu_sum = 0;
count = 0;
for ( atom = 0; atom < system->N; ++atom )
+ {
/* start discovering water molecules from the central O atom */
if ( !mark[atom] && system->atoms[atom].type == 2 )
{
@@ -773,6 +823,7 @@ static void Calculate_Dipole_Moment( reax_system *system, control_params *contro
mu_sum += rvec_Norm( mu );
}
}
+ }
fprintf( fout, "%7d %10d %10.5f\n",
data->step, count, mu_sum / count * ECxA_to_DEBYE );
@@ -785,13 +836,14 @@ static void Copy_Positions( reax_system *system, static_storage *workspace )
int i;
for ( i = 0; i < system->N; ++i )
+ {
rvec_Copy( workspace->x_old[i], system->atoms[i].x );
+ }
}
static void Calculate_Drift( reax_system *system, control_params *control,
- simulation_data *data, static_storage *workspace,
- FILE *fout )
+ simulation_data *data, static_storage *workspace, FILE *fout )
{
int i, type;
int count[MAX_ATOM_TYPES];
@@ -806,6 +858,7 @@ static void Calculate_Drift( reax_system *system, control_params *control,
}
for ( i = 0; i < system->N; ++i )
+ {
//if( control->restrict_type == -1 ||
// system->atoms[i].type == control->restrict_type )
if ( workspace->x_old[i][0] > -system->box.box_norms[0] &&
@@ -816,7 +869,7 @@ static void Calculate_Drift( reax_system *system, control_params *control,
++count[type];
Distance_on_T3_Gen( workspace->x_old[i], system->atoms[i].x,
- &(system->box), driftvec );
+ &(system->box), driftvec );
if ( FABS( driftvec[0] ) >= system->box.box_norms[0] / 2.0 - 2.0 ||
FABS( driftvec[1] ) >= system->box.box_norms[0] / 2.0 - 2.0 ||
@@ -835,6 +888,7 @@ static void Calculate_Drift( reax_system *system, control_params *control,
sum_sqr_drift[type] += drift;
sum_drift[type] += SQRT( drift );
}
+ }
fprintf( fout, "%7d oxy %6d %10.6f\n",
data->step, count[2], sum_sqr_drift[2] / count[2] );
@@ -846,7 +900,7 @@ static void Calculate_Drift( reax_system *system, control_params *control,
static void Calculate_Density_3DMesh( reax_system *system, simulation_data *data,
- FILE *fout )
+ FILE *fout )
{
int i, j, k;
int occupied_cells;
@@ -858,7 +912,9 @@ static void Calculate_Density_3DMesh( reax_system *system, simulation_data *data
/* determine the mesh dimensions based on the current box size */
for ( i = 0; i < 3; ++i )
+ {
mesh_dims[i] = system->box.box_norms[i] / mesh_cell_lens[i] + 0.99;
+ }
fprintf( stderr, "mesh_dims: %3d %3d %3d\n",
mesh_dims[0], mesh_dims[1], mesh_dims[2] );
@@ -873,8 +929,10 @@ static void Calculate_Density_3DMesh( reax_system *system, simulation_data *data
"Calculate_Density_3DMesh::cell_counter[i]" );
for ( j = 0; j < mesh_dims[1]; ++j )
+ {
cell_counter[i][j] = (int *) scalloc( mesh_dims[2], sizeof(int),
"Calculate_Density_3DMesh::cell_counter[i][j]" );
+ }
}
@@ -892,28 +950,35 @@ static void Calculate_Density_3DMesh( reax_system *system, simulation_data *data
/* calculate volume occupied */
occupied_cells = 0;
for ( i = 0; i < mesh_dims[0]; ++i )
+ {
for ( j = 0; j < mesh_dims[1]; ++j )
+ {
for ( k = 0; k < mesh_dims[2]; ++k )
+ {
if ( cell_counter[i][j][k] )
+ {
++occupied_cells;
+ }
+ }
+ }
+ }
- occupied_vol =
- occupied_cells * mesh_cell_lens[0] * mesh_cell_lens[1] * mesh_cell_lens[2];
+ occupied_vol = occupied_cells * mesh_cell_lens[0]
+ * mesh_cell_lens[1] * mesh_cell_lens[2];
fprintf( stderr, "occupied cells: %8d\n", occupied_cells );
fprintf( stderr, "occupied vol : %8.2f\n", occupied_vol );
fprintf( stderr, "system volume : %8.2f\n", system->box.volume );
fprintf( stderr, "system mass : %8.2f\n", data->M );
- density = data->M * AMU_to_GRAM / (occupied_vol * POW( ANG_to_CM, 3 ));
+ density = data->M * AMU_to_GRAM / (occupied_vol * POW( ANG_to_CM, 3.0 ));
fprintf( stderr, "density : %g\n", density );
fprintf( stderr, "AMU_to_GRAM : %g\n", AMU_to_GRAM );
fprintf( stderr, "ANG_to_CM : %g\n", ANG_to_CM );
}
-static void Calculate_Density_Slice( reax_system *system, simulation_data *data,
- FILE *fout )
+static void Calculate_Density_Slice( reax_system *system, simulation_data *data, FILE *fout )
{
real slice_thickness = 0.5;
int *slice_occ;
@@ -936,16 +1001,20 @@ static void Calculate_Density_Slice( reax_system *system, simulation_data *data,
{
fprintf( stderr, "occ[%d]: %d\n", i, slice_occ[i] );
if ( slice_occ[i] > max_occ )
+ {
max_occ = slice_occ[i];
+ }
}
/* find luzzatti-interface slice */
for ( i = 0; i < num_slices; ++i )
+ {
if ( (real)slice_occ[i] / max_occ > 0.5 )
{
fprintf( stderr, "%d - %d is the luzzatti interface\n", i - 1, i );
break;
}
+ }
}
diff --git a/sPuReMD/src/charges.c b/sPuReMD/src/charges.c
index ca98e58596469ceb69e654eb929320e323354122..4ff5a3426ba83949b084c2e078ebe974311fee90 100644
--- a/sPuReMD/src/charges.c
+++ b/sPuReMD/src/charges.c
@@ -1478,9 +1478,9 @@ static void ACKS2( reax_system * const system, control_params * const control,
if ( data->step % 10 == 0 )
{
snprintf( fname, SIZE + 11, "s_%d_%s.out", data->step, control->sim_name );
- fp = fopen( fname, "w" );
+ fp = sfopen( fname, "w" );
Vector_Print( fp, NULL, workspace->s[0], system->N_cm );
- fclose( fp );
+ sfclose( fp, "ACKS2::fp" );
}
#undef SIZE
#endif
@@ -1547,14 +1547,14 @@ void Compute_Charges( reax_system * const system, control_params * const control
// Print_Sparse_Matrix_Binary( workspace->H, fname );
snprintf( fname, SIZE + 11, "b_s_%d_%s.out", data->step, control->sim_name );
- fp = fopen( fname, "w" );
+ fp = sfopen( fname, "w" );
Vector_Print( fp, NULL, workspace->b_s, system->N_cm );
- fclose( fp );
+ sfclose( fp, "Compute_Charges::fp" );
// snprintf( fname, SIZE + 11, "b_t_%d_%s.out", data->step, control->sim_name );
-// fp = fopen( fname, "w" );
+// fp = sfopen( fname, "w" );
// Vector_Print( fp, NULL, workspace->b_t, system->N_cm );
-// fclose( fp );
+// sfclose( fp, "Compute_Charges::fp" );
}
#undef SIZE
#endif
diff --git a/sPuReMD/src/control.c b/sPuReMD/src/control.c
index 5325a929dfa08855e406d1dc6c3dbf5a0b453326..bd95e4e9ed6de618f0a5100f80c6f71b946ec00d 100644
--- a/sPuReMD/src/control.c
+++ b/sPuReMD/src/control.c
@@ -37,12 +37,12 @@ void Read_Control_File( FILE* fp, reax_system *system, control_params* control,
int ival;
/* assign default values */
- strcpy( control->sim_name, "default.sim" );
+ strncpy( control->sim_name, "default.sim", MAX_STR );
control->restart = 0;
out_control->restart_format = WRITE_BINARY;
out_control->restart_freq = 0;
- strcpy( control->restart_from, "default.res" );
+ strncpy( control->restart_from, "default.res", MAX_STR );
out_control->restart_freq = 0;
control->random_vel = 0;
@@ -119,7 +119,7 @@ void Read_Control_File( FILE* fp, reax_system *system, control_params* control,
(int (*)( reax_system*, control_params*, simulation_data*,
static_storage*, reax_list **, void* )) Append_Custom_Frame;
- strcpy( out_control->traj_title, "default_title" );
+ strncpy( out_control->traj_title, "default_title", 81 );
out_control->atom_format = 0;
out_control->bond_info = 0;
out_control->angle_info = 0;
@@ -128,7 +128,7 @@ void Read_Control_File( FILE* fp, reax_system *system, control_params* control,
control->freq_molec_anal = 0;
control->bg_cut = 0.3;
control->num_ignored = 0;
- memset( control->ignore, 0, sizeof(int)*MAX_ATOM_TYPES );
+ memset( control->ignore, 0, sizeof(int) * MAX_ATOM_TYPES );
control->dipole_anal = 0;
control->freq_dipole_anal = 0;
@@ -147,129 +147,129 @@ void Read_Control_File( FILE* fp, reax_system *system, control_params* control,
}
/* read control parameters file */
- while (fgets(s, MAX_LINE, fp))
+ while ( fgets( s, MAX_LINE, fp ) )
{
- Tokenize(s, &tmp);
+ Tokenize( s, &tmp );
- if ( strcmp(tmp[0], "simulation_name") == 0 )
+ if ( strncmp(tmp[0], "simulation_name", MAX_LINE) == 0 )
{
- strcpy( control->sim_name, tmp[1] );
+ strncpy( control->sim_name, tmp[1], MAX_STR );
}
- //else if( strcmp(tmp[0], "restart") == 0 ) {
+ //else if( strncmp(tmp[0], "restart", MAX_LINE) == 0 ) {
// ival = atoi(tmp[1]);
// control->restart = ival;
//}
- else if ( strcmp(tmp[0], "restart_format") == 0 )
+ else if ( strncmp(tmp[0], "restart_format", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
out_control->restart_format = ival;
}
- else if ( strcmp(tmp[0], "restart_freq") == 0 )
+ else if ( strncmp(tmp[0], "restart_freq", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
out_control->restart_freq = ival;
}
- else if ( strcmp(tmp[0], "random_vel") == 0 )
+ else if ( strncmp(tmp[0], "random_vel", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
control->random_vel = ival;
}
- else if ( strcmp(tmp[0], "reposition_atoms") == 0 )
+ else if ( strncmp(tmp[0], "reposition_atoms", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
control->reposition_atoms = ival;
}
- else if ( strcmp(tmp[0], "ensemble_type") == 0 )
+ else if ( strncmp(tmp[0], "ensemble_type", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
control->ensemble = ival;
}
- else if ( strcmp(tmp[0], "nsteps") == 0 )
+ else if ( strncmp(tmp[0], "nsteps", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
control->nsteps = ival;
}
- else if ( strcmp(tmp[0], "dt") == 0 )
+ else if ( strncmp(tmp[0], "dt", MAX_LINE) == 0 )
{
val = atof(tmp[1]);
control->dt = val * 1.e-3; // convert dt from fs to ps!
}
- else if ( strcmp(tmp[0], "periodic_boundaries") == 0 )
+ else if ( strncmp(tmp[0], "periodic_boundaries", MAX_LINE) == 0 )
{
ival = atoi( tmp[1] );
control->periodic_boundaries = ival;
}
- else if ( strcmp(tmp[0], "geo_format") == 0 )
+ else if ( strncmp(tmp[0], "geo_format", MAX_LINE) == 0 )
{
ival = atoi( tmp[1] );
control->geo_format = ival;
}
- else if ( strcmp(tmp[0], "restrict_bonds") == 0 )
+ else if ( strncmp(tmp[0], "restrict_bonds", MAX_LINE) == 0 )
{
ival = atoi( tmp[1] );
control->restrict_bonds = ival;
}
- else if ( strcmp(tmp[0], "tabulate_long_range") == 0 )
+ else if ( strncmp(tmp[0], "tabulate_long_range", MAX_LINE) == 0 )
{
ival = atoi( tmp[1] );
control->tabulate = ival;
}
- else if ( strcmp(tmp[0], "reneighbor") == 0 )
+ else if ( strncmp(tmp[0], "reneighbor", MAX_LINE) == 0 )
{
ival = atoi( tmp[1] );
control->reneighbor = ival;
}
- else if ( strcmp(tmp[0], "vlist_buffer") == 0 )
+ else if ( strncmp(tmp[0], "vlist_buffer", MAX_LINE) == 0 )
{
val = atof(tmp[1]);
control->vlist_cut = val;
}
- else if ( strcmp(tmp[0], "nbrhood_cutoff") == 0 )
+ else if ( strncmp(tmp[0], "nbrhood_cutoff", MAX_LINE) == 0 )
{
val = atof(tmp[1]);
control->nbr_cut = val;
}
- else if ( strcmp(tmp[0], "thb_cutoff") == 0 )
+ else if ( strncmp(tmp[0], "thb_cutoff", MAX_LINE) == 0 )
{
val = atof(tmp[1]);
control->thb_cut = val;
}
- else if ( strcmp(tmp[0], "hbond_cutoff") == 0 )
+ else if ( strncmp(tmp[0], "hbond_cutoff", MAX_LINE) == 0 )
{
val = atof( tmp[1] );
control->hb_cut = val;
}
- else if ( strcmp(tmp[0], "charge_method") == 0 )
+ else if ( strncmp(tmp[0], "charge_method", MAX_LINE) == 0 )
{
ival = atoi( tmp[1] );
control->charge_method = ival;
}
- else if ( strcmp(tmp[0], "cm_q_net") == 0 )
+ else if ( strncmp(tmp[0], "cm_q_net", MAX_LINE) == 0 )
{
val = atof( tmp[1] );
control->cm_q_net = val;
}
- else if ( strcmp(tmp[0], "cm_solver_type") == 0 )
+ else if ( strncmp(tmp[0], "cm_solver_type", MAX_LINE) == 0 )
{
ival = atoi( tmp[1] );
control->cm_solver_type = ival;
}
- else if ( strcmp(tmp[0], "cm_solver_max_iters") == 0 )
+ else if ( strncmp(tmp[0], "cm_solver_max_iters", MAX_LINE) == 0 )
{
ival = atoi( tmp[1] );
control->cm_solver_max_iters = ival;
}
- else if ( strcmp(tmp[0], "cm_solver_restart") == 0 )
+ else if ( strncmp(tmp[0], "cm_solver_restart", MAX_LINE) == 0 )
{
ival = atoi( tmp[1] );
control->cm_solver_restart = ival;
}
- else if ( strcmp(tmp[0], "cm_solver_q_err") == 0 )
+ else if ( strncmp(tmp[0], "cm_solver_q_err", MAX_LINE) == 0 )
{
val = atof( tmp[1] );
control->cm_solver_q_err = val;
}
- else if ( strcmp(tmp[0], "cm_domain_sparsity") == 0 )
+ else if ( strncmp(tmp[0], "cm_domain_sparsity", MAX_LINE) == 0 )
{
val = atof( tmp[1] );
control->cm_domain_sparsity = val;
@@ -278,52 +278,52 @@ void Read_Control_File( FILE* fp, reax_system *system, control_params* control,
control->cm_domain_sparsify_enabled = TRUE;
}
}
- else if ( strcmp(tmp[0], "cm_init_guess_extrap1") == 0 )
+ else if ( strncmp(tmp[0], "cm_init_guess_extrap1", MAX_LINE) == 0 )
{
ival = atoi( tmp[1] );
control->cm_init_guess_extrap1 = ival;
}
- else if ( strcmp(tmp[0], "cm_init_guess_extrap2") == 0 )
+ else if ( strncmp(tmp[0], "cm_init_guess_extrap2", MAX_LINE) == 0 )
{
ival = atoi( tmp[1] );
control->cm_init_guess_extrap2 = ival;
}
- else if ( strcmp(tmp[0], "cm_solver_pre_comp_type") == 0 )
+ else if ( strncmp(tmp[0], "cm_solver_pre_comp_type", MAX_LINE) == 0 )
{
ival = atoi( tmp[1] );
control->cm_solver_pre_comp_type = ival;
}
- else if ( strcmp(tmp[0], "cm_solver_pre_comp_refactor") == 0 )
+ else if ( strncmp(tmp[0], "cm_solver_pre_comp_refactor", MAX_LINE) == 0 )
{
ival = atoi( tmp[1] );
control->cm_solver_pre_comp_refactor = ival;
}
- else if ( strcmp(tmp[0], "cm_solver_pre_comp_droptol") == 0 )
+ else if ( strncmp(tmp[0], "cm_solver_pre_comp_droptol", MAX_LINE) == 0 )
{
val = atof( tmp[1] );
control->cm_solver_pre_comp_droptol = val;
}
- else if ( strcmp(tmp[0], "cm_solver_pre_comp_sweeps") == 0 )
+ else if ( strncmp(tmp[0], "cm_solver_pre_comp_sweeps", MAX_LINE) == 0 )
{
ival = atoi( tmp[1] );
control->cm_solver_pre_comp_sweeps = ival;
}
- else if ( strcmp(tmp[0], "cm_solver_pre_comp_sai_thres") == 0 )
+ else if ( strncmp(tmp[0], "cm_solver_pre_comp_sai_thres", MAX_LINE) == 0 )
{
val = atof( tmp[1] );
control->cm_solver_pre_comp_sai_thres = val;
}
- else if ( strcmp(tmp[0], "cm_solver_pre_app_type") == 0 )
+ else if ( strncmp(tmp[0], "cm_solver_pre_app_type", MAX_LINE) == 0 )
{
ival = atoi( tmp[1] );
control->cm_solver_pre_app_type = ival;
}
- else if ( strcmp(tmp[0], "cm_solver_pre_app_jacobi_iters") == 0 )
+ else if ( strncmp(tmp[0], "cm_solver_pre_app_jacobi_iters", MAX_LINE) == 0 )
{
ival = atoi( tmp[1] );
control->cm_solver_pre_app_jacobi_iters = ival;
}
- else if ( strcmp(tmp[0], "temp_init") == 0 )
+ else if ( strncmp(tmp[0], "temp_init", MAX_LINE) == 0 )
{
val = atof(tmp[1]);
control->T_init = val;
@@ -333,7 +333,7 @@ void Read_Control_File( FILE* fp, reax_system *system, control_params* control,
control->T_init = 0.001;
}
}
- else if ( strcmp(tmp[0], "temp_final") == 0 )
+ else if ( strncmp(tmp[0], "temp_final", MAX_LINE) == 0 )
{
val = atof(tmp[1]);
control->T_final = val;
@@ -343,28 +343,28 @@ void Read_Control_File( FILE* fp, reax_system *system, control_params* control,
control->T_final = 0.1;
}
}
- else if ( strcmp(tmp[0], "t_mass") == 0 )
+ else if ( strncmp(tmp[0], "t_mass", MAX_LINE) == 0 )
{
val = atof(tmp[1]);
/* convert t_mass from fs to ps */
control->Tau_T = val * 1.e-3;
}
- else if ( strcmp(tmp[0], "t_mode") == 0 )
+ else if ( strncmp(tmp[0], "t_mode", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
control->T_mode = ival;
}
- else if ( strcmp(tmp[0], "t_rate") == 0 )
+ else if ( strncmp(tmp[0], "t_rate", MAX_LINE) == 0 )
{
val = atof(tmp[1]);
control->T_rate = val;
}
- else if ( strcmp(tmp[0], "t_freq") == 0 )
+ else if ( strncmp(tmp[0], "t_freq", MAX_LINE) == 0 )
{
val = atof(tmp[1]);
control->T_freq = val;
}
- else if ( strcmp(tmp[0], "pressure") == 0 )
+ else if ( strncmp(tmp[0], "pressure", MAX_LINE) == 0 )
{
if ( control->ensemble == iNPT )
{
@@ -383,7 +383,7 @@ void Read_Control_File( FILE* fp, reax_system *system, control_params* control,
control->P[2] = val;
}
}
- else if ( strcmp(tmp[0], "p_mass") == 0 )
+ else if ( strncmp(tmp[0], "p_mass", MAX_LINE) == 0 )
{
if ( control->ensemble == iNPT )
{
@@ -402,42 +402,42 @@ void Read_Control_File( FILE* fp, reax_system *system, control_params* control,
control->Tau_P[2] = val * 1.e-3; // convert p_mass from fs to ps
}
}
- else if ( strcmp(tmp[0], "pt_mass") == 0 )
+ else if ( strncmp(tmp[0], "pt_mass", MAX_LINE) == 0 )
{
val = atof(tmp[1]);
control->Tau_PT = val * 1.e-3; // convert pt_mass from fs to ps
}
- else if ( strcmp(tmp[0], "compress") == 0 )
+ else if ( strncmp(tmp[0], "compress", MAX_LINE) == 0 )
{
val = atof(tmp[1]);
control->compressibility = val;
}
- else if ( strcmp(tmp[0], "press_mode") == 0 )
+ else if ( strncmp(tmp[0], "press_mode", MAX_LINE) == 0 )
{
val = atoi(tmp[1]);
control->press_mode = val;
}
- else if ( strcmp(tmp[0], "remove_CoM_vel") == 0 )
+ else if ( strncmp(tmp[0], "remove_CoM_vel", MAX_LINE) == 0 )
{
val = atoi(tmp[1]);
control->remove_CoM_vel = val;
}
- else if ( strcmp(tmp[0], "debug_level") == 0 )
+ else if ( strncmp(tmp[0], "debug_level", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
out_control->debug_level = ival;
}
- else if ( strcmp(tmp[0], "energy_update_freq") == 0 )
+ else if ( strncmp(tmp[0], "energy_update_freq", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
out_control->energy_update_freq = ival;
}
- else if ( strcmp(tmp[0], "write_freq") == 0 )
+ else if ( strncmp(tmp[0], "write_freq", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
out_control->write_steps = ival;
}
- else if ( strcmp(tmp[0], "traj_compress") == 0 )
+ else if ( strncmp(tmp[0], "traj_compress", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
out_control->traj_compress = ival;
@@ -446,7 +446,7 @@ void Read_Control_File( FILE* fp, reax_system *system, control_params* control,
out_control->write = (int (*)(FILE *, const char *, ...)) gzprintf;
else out_control->write = fprintf;
}
- else if ( strcmp(tmp[0], "traj_format") == 0 )
+ else if ( strncmp(tmp[0], "traj_format", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
out_control->traj_format = ival;
@@ -470,81 +470,81 @@ void Read_Control_File( FILE* fp, reax_system *system, control_params* control,
static_storage*, reax_list **, void* )) Append_xyz_Frame;
}
}
- else if ( strcmp(tmp[0], "traj_title") == 0 )
+ else if ( strncmp(tmp[0], "traj_title", MAX_LINE) == 0 )
{
- strcpy( out_control->traj_title, tmp[1] );
+ strncpy( out_control->traj_title, tmp[1], 81 );
}
- else if ( strcmp(tmp[0], "atom_info") == 0 )
+ else if ( strncmp(tmp[0], "atom_info", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
out_control->atom_format += ival * 4;
}
- else if ( strcmp(tmp[0], "atom_velocities") == 0 )
+ else if ( strncmp(tmp[0], "atom_velocities", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
out_control->atom_format += ival * 2;
}
- else if ( strcmp(tmp[0], "atom_forces") == 0 )
+ else if ( strncmp(tmp[0], "atom_forces", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
out_control->atom_format += ival * 1;
}
- else if ( strcmp(tmp[0], "bond_info") == 0 )
+ else if ( strncmp(tmp[0], "bond_info", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
out_control->bond_info = ival;
}
- else if ( strcmp(tmp[0], "angle_info") == 0 )
+ else if ( strncmp(tmp[0], "angle_info", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
out_control->angle_info = ival;
}
- else if ( strcmp(tmp[0], "test_forces") == 0 )
+ else if ( strncmp(tmp[0], "test_forces", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
}
- else if ( strcmp(tmp[0], "molec_anal") == 0 )
+ else if ( strncmp(tmp[0], "molec_anal", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
control->molec_anal = ival;
}
- else if ( strcmp(tmp[0], "freq_molec_anal") == 0 )
+ else if ( strncmp(tmp[0], "freq_molec_anal", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
control->freq_molec_anal = ival;
}
- else if ( strcmp(tmp[0], "bond_graph_cutoff") == 0 )
+ else if ( strncmp(tmp[0], "bond_graph_cutoff", MAX_LINE) == 0 )
{
val = atof(tmp[1]);
control->bg_cut = val;
}
- else if ( strcmp(tmp[0], "ignore") == 0 )
+ else if ( strncmp(tmp[0], "ignore", MAX_LINE) == 0 )
{
control->num_ignored = atoi(tmp[1]);
for ( i = 0; i < control->num_ignored; ++i )
control->ignore[atoi(tmp[i + 2])] = 1;
}
- else if ( strcmp(tmp[0], "dipole_anal") == 0 )
+ else if ( strncmp(tmp[0], "dipole_anal", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
control->dipole_anal = ival;
}
- else if ( strcmp(tmp[0], "freq_dipole_anal") == 0 )
+ else if ( strncmp(tmp[0], "freq_dipole_anal", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
control->freq_dipole_anal = ival;
}
- else if ( strcmp(tmp[0], "diffusion_coef") == 0 )
+ else if ( strncmp(tmp[0], "diffusion_coef", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
control->diffusion_coef = ival;
}
- else if ( strcmp(tmp[0], "freq_diffusion_coef") == 0 )
+ else if ( strncmp(tmp[0], "freq_diffusion_coef", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
control->freq_diffusion_coef = ival;
}
- else if ( strcmp(tmp[0], "restrict_type") == 0 )
+ else if ( strncmp(tmp[0], "restrict_type", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
control->restrict_type = ival;
diff --git a/sPuReMD/src/geo_tools.c b/sPuReMD/src/geo_tools.c
index 05ece8c81c4ea2b2543593b98cffb27d0d378f6a..cbf9ca42b1b9b85f95fe0c0712f2d9989ddd1333 100644
--- a/sPuReMD/src/geo_tools.c
+++ b/sPuReMD/src/geo_tools.c
@@ -125,11 +125,7 @@ void Read_Geo( const char * const geo_file, reax_system* system, control_params
reax_atom *atom;
/* open the geometry file */
- if ( (geo = fopen(geo_file, "r")) == NULL )
- {
- fprintf( stderr, "Error opening the geo file! terminating...\n" );
- exit( FILE_NOT_FOUND );
- }
+ geo = sfopen( geo_file, "r" );
/* read box information */
fscanf( geo, CUSTOM_BOXGEO_FORMAT,
@@ -156,7 +152,7 @@ void Read_Geo( const char * const geo_file, reax_system* system, control_params
atom = &(system->atoms[top]);
workspace->orig_id[i] = serial;
atom->type = Get_Atom_Type( &(system->reaxprm), element );
- strcpy( atom->name, name );
+ strncpy( atom->name, name, 8 );
rvec_Copy( atom->x, x );
rvec_MakeZero( atom->v );
rvec_MakeZero( atom->f );
@@ -165,7 +161,7 @@ void Read_Geo( const char * const geo_file, reax_system* system, control_params
top++;
}
- fclose( geo );
+ sfclose( geo, "Read_Geo::geo" );
}
@@ -212,7 +208,6 @@ static void Count_PDB_Atoms( FILE *geo, reax_system *system )
void Read_PDB( const char * const pdb_file, reax_system* system, control_params *control,
simulation_data *data, static_storage *workspace )
{
-
FILE *pdb;
char **tmp;
char *s, *s1;
@@ -228,11 +223,7 @@ void Read_PDB( const char * const pdb_file, reax_system* system, control_params
reax_atom *atom;
/* open pdb file */
- if ( (pdb = fopen(pdb_file, "r")) == NULL )
- {
- fprintf( stderr, "fopen: error opening the pdb file! terminating...\n" );
- exit( FILE_NOT_FOUND );
- }
+ pdb = sfopen( pdb_file, "r" );
/* allocate memory for tokenizing pdb lines */
Allocate_Tokenizer_Space( &s, &s1, &tmp );
@@ -346,7 +337,7 @@ void Read_PDB( const char * const pdb_file, reax_system* system, control_params
Trim_Spaces( element, 9 );
atom->type = Get_Atom_Type( &(system->reaxprm), element );
- strcpy( atom->name, atom_name );
+ strncpy( atom->name, atom_name, 8 );
rvec_Copy( atom->x, x );
rvec_MakeZero( atom->v );
@@ -403,7 +394,7 @@ void Read_PDB( const char * const pdb_file, reax_system* system, control_params
exit( INVALID_INPUT );
}
- fclose( pdb );
+ sfclose( pdb, "Read_PDB::pdb" );
Deallocate_Tokenizer_Space( &s, &s1, &tmp );
}
@@ -416,72 +407,58 @@ void Read_PDB( const char * const pdb_file, reax_system* system, control_params
void Write_PDB( reax_system* system, reax_list* bonds, simulation_data *data,
control_params *control, static_storage *workspace, output_controls *out_control )
{
- int i, buffer_req, buffer_len;
- //int j, connect[4];
+ int i;
char name[8];
- //real bo;
real alpha, beta, gamma;
+ rvec x;
reax_atom *p_atom;
char fname[MAX_STR];
- char *line;
- char *buffer;
+ char buffer[PDB_ATOM_FORMAT_O_LENGTH + 1];
FILE *pdb;
- /* Allocation */
- line = (char*) smalloc( sizeof(char) * PDB_ATOM_FORMAT_O_LENGTH, "geo:line" );
- buffer_req = system->N * PDB_ATOM_FORMAT_O_LENGTH;
-
- buffer = (char*) smalloc( sizeof(char) * buffer_req, "geo:buffer" );
-
- pdb = NULL;
- line[0] = 0;
- buffer[0] = 0;
/* Writing Box information */
gamma = ACOS( (system->box.box[0][0] * system->box.box[1][0] +
- system->box.box[0][1] * system->box.box[1][1] +
- system->box.box[0][2] * system->box.box[1][2]) /
+ system->box.box[0][1] * system->box.box[1][1] +
+ system->box.box[0][2] * system->box.box[1][2]) /
(system->box.box_norms[0] * system->box.box_norms[1]) );
beta = ACOS( (system->box.box[0][0] * system->box.box[2][0] +
- system->box.box[0][1] * system->box.box[2][1] +
- system->box.box[0][2] * system->box.box[2][2]) /
- (system->box.box_norms[0] * system->box.box_norms[2]) );
+ system->box.box[0][1] * system->box.box[2][1] +
+ system->box.box[0][2] * system->box.box[2][2]) /
+ (system->box.box_norms[0] * system->box.box_norms[2]) );
alpha = ACOS( (system->box.box[2][0] * system->box.box[1][0] +
- system->box.box[2][1] * system->box.box[1][1] +
- system->box.box[2][2] * system->box.box[1][2]) /
- (system->box.box_norms[2] * system->box.box_norms[1]) );
+ system->box.box[2][1] * system->box.box[1][1] +
+ system->box.box[2][2] * system->box.box[1][2]) /
+ (system->box.box_norms[2] * system->box.box_norms[1]) );
/* open pdb and write header */
snprintf( fname, MAX_STR + 9, "%s-%d.pdb", control->sim_name, data->step );
- pdb = fopen(fname, "w");
+ pdb = sfopen( fname, "w" );
fprintf( pdb, PDB_CRYST1_FORMAT_O,
"CRYST1",
system->box.box_norms[0], system->box.box_norms[1],
system->box.box_norms[2],
RAD2DEG(alpha), RAD2DEG(beta), RAD2DEG(gamma), " ", 0 );
- fprintf( out_control->log, "Box written\n" );
- fflush( out_control->log );
/* write atom lines to buffer */
for ( i = 0; i < system->N; i++)
{
p_atom = &(system->atoms[i]);
- strncpy(name, p_atom->name, 8);
+
+ strncpy( name, p_atom->name, 8 );
Trim_Spaces( name, 8 );
- snprintf( line, MAX_STR, PDB_ATOM_FORMAT_O,
- "ATOM ", workspace->orig_id[i], p_atom->name, ' ', "REX", ' ', 1, ' ',
- p_atom->x[0], p_atom->x[1], p_atom->x[2],
- 1.0, 0.0, "0", name, " " );
-#if defined(DEBUG)
- fprintf( stderr, "PDB NAME <%s>\n", p_atom->name );
-#endif
+ memcpy( x, p_atom->x, 3 * sizeof(real) );
+ Fit_to_Periodic_Box( &(system->box), &x );
- strncpy( buffer + i * PDB_ATOM_FORMAT_O_LENGTH, line,
- PDB_ATOM_FORMAT_O_LENGTH );
- }
+ snprintf( buffer, PDB_ATOM_FORMAT_O_LENGTH, PDB_ATOM_FORMAT_O,
+ "ATOM ", workspace->orig_id[i], p_atom->name, ' ', "REX", ' ', 1, ' ',
+ x[0], x[1], x[2],
+ 1.0, 0.0, "0", name, " " );
+
+ buffer[PDB_ATOM_FORMAT_O_LENGTH] = '\n';
- buffer_len = system->N * PDB_ATOM_FORMAT_O_LENGTH;
- buffer[buffer_len] = 0;
+ fprintf( pdb, "%s\n", buffer );
+ }
if ( ferror( pdb ) )
{
@@ -489,30 +466,7 @@ void Write_PDB( reax_system* system, reax_list* bonds, simulation_data *data,
exit( INVALID_INPUT );
}
- fprintf( pdb, "%s", buffer );
- fclose( pdb );
-
- /* Writing connect information */
- /*
- for(i=0; i < system->N; i++) {
- count = 0;
- for(j = Start_Index(i, bonds); j < End_Index(i, bonds); ++j) {
- bo = bonds->bond_list[j].bo_data.BO;
- if (bo > 0.3) {
- connect[count] = bonds->bond_list[j].nbr+1;
- count++;
- }
- }
-
- fprintf( out_control->pdb, "%6s%5d", "CONECT", i+1 );
- for( k=0; k < count; k++ )
- fprintf( out_control->pdb, "%5d", connect[k] );
- fprintf( out_control->pdb, "\n" );
- }
- */
-
- sfree( buffer, "Write_PDB::buffer" );
- sfree( line, "Write_PDB::line" );
+ sfclose( pdb, "Write_PDB::pdb" );
}
@@ -530,14 +484,10 @@ void Read_BGF( const char * const bgf_file, reax_system* system, control_params
char chain_id;
char s_a[12], s_b[12], s_c[12], s_alpha[12], s_beta[12], s_gamma[12];
char *endptr = NULL;
- int i, atom_cnt, token_cnt, bgf_serial, ratom = 0;
+ int i, atom_cnt, token_cnt, bgf_serial, ratom = 0;
/* open biograf file */
- if ( (bgf = fopen( bgf_file, "r" )) == NULL )
- {
- fprintf( stderr, "Error opening the bgf file!\n" );
- exit( FILE_NOT_FOUND );
- }
+ bgf = sfopen( bgf_file, "r" );
/* allocate memory for tokenizing biograf file lines */
line = (char*) smalloc( sizeof(char) * MAX_LINE,
@@ -561,7 +511,8 @@ void Read_BGF( const char * const bgf_file, reax_system* system, control_params
tokens[0][0] = 0;
token_cnt = Tokenize( line, &tokens );
- if ( !strcmp( tokens[0], "ATOM" ) || !strcmp( tokens[0], "HETATM" ) )
+ if ( !strncmp( tokens[0], "ATOM", MAX_TOKEN_LEN )
+ || !strncmp( tokens[0], "HETATM", MAX_TOKEN_LEN ) )
{
++system->N;
}
@@ -578,7 +529,7 @@ void Read_BGF( const char * const bgf_file, reax_system* system, control_params
fprintf( stderr, "system->N: %d\n", system->N );
#endif
- fclose( bgf );
+ sfclose( bgf, "Read_BGF::bgf" );
/* memory allocations for atoms, atom maps, bond restrictions */
system->atoms = (reax_atom*) scalloc( system->N, sizeof(reax_atom),
@@ -604,11 +555,7 @@ void Read_BGF( const char * const bgf_file, reax_system* system, control_params
}
/* start reading and processing bgf file */
- if ( (bgf = fopen( bgf_file, "r" )) == NULL )
- {
- fprintf( stderr, "Error opening the bgf file!\n" );
- exit( FILE_NOT_FOUND );
- }
+ bgf = sfopen( bgf_file, "r" );
atom_cnt = 0;
token_cnt = 0;
@@ -693,7 +640,7 @@ void Read_BGF( const char * const bgf_file, reax_system* system, control_params
system->atoms[atom_cnt].x[2] = strtod( &s_z[0], &endptr );
/* atom name and type */
- strcpy( system->atoms[atom_cnt].name, atom_name );
+ strncpy( system->atoms[atom_cnt].name, atom_name, 8 );
Trim_Spaces( element, 10 );
system->atoms[atom_cnt].type =
Get_Atom_Type( &(system->reaxprm), element );
@@ -710,14 +657,9 @@ void Read_BGF( const char * const bgf_file, reax_system* system, control_params
}
else if (!strncmp( tokens[0], "CRYSTX", 6 ))
{
- sscanf( backup, BGF_CRYSTX_FORMAT,
- &descriptor[0],
- &s_a[0],
- &s_b[0],
- &s_c[0],
- &s_alpha[0],
- &s_beta[0],
- &s_gamma[0] );
+ sscanf( backup, BGF_CRYSTX_FORMAT, &descriptor[0],
+ &s_a[0], &s_b[0], &s_c[0],
+ &s_alpha[0], &s_beta[0], &s_gamma[0] );
/* Compute full volume tensor from the angles */
Setup_Box( atof(s_a), atof(s_b), atof(s_c),
@@ -766,5 +708,5 @@ void Read_BGF( const char * const bgf_file, reax_system* system, control_params
exit( INVALID_INPUT );
}
- fclose( bgf );
+ sfclose( bgf, "Read_BGF::bgf" );
}
diff --git a/sPuReMD/src/grid.c b/sPuReMD/src/grid.c
index b0749b06940cea6d5aa50a0493de2aa77d76fc05..e5f1e60eeb369357be9f77a43a6853a29ebd66a2 100644
--- a/sPuReMD/src/grid.c
+++ b/sPuReMD/src/grid.c
@@ -539,7 +539,7 @@ static inline void reax_atom_Copy( reax_atom *dest, reax_atom *src )
dest->type = src->type;
rvec_Copy( dest->x, src->x );
rvec_Copy( dest->v, src->v );
- strcpy( dest->name, src->name );
+ strncpy( dest->name, src->name, 8 );
}
diff --git a/sPuReMD/src/init_md.c b/sPuReMD/src/init_md.c
index 02a45310848fd60f9da831cb72929b96e24990b7..cc0ce302ecb46c17b71cfb87810ecbae5ff2d6c3 100644
--- a/sPuReMD/src/init_md.c
+++ b/sPuReMD/src/init_md.c
@@ -829,27 +829,31 @@ static void Init_Out_Controls( reax_system *system, control_params *control,
/* Init trajectory file */
if ( out_control->write_steps > 0 )
{
- strcpy( temp, control->sim_name );
+ strncpy( temp, control->sim_name, TEMP_SIZE );
strcat( temp, ".trj" );
- out_control->trj = fopen( temp, "w" );
+ out_control->trj = sfopen( temp, "w" );
out_control->write_header( system, control, workspace, out_control );
}
+ else
+ {
+ out_control->trj = NULL;
+ }
if ( out_control->energy_update_freq > 0 )
{
/* Init out file */
- strcpy( temp, control->sim_name );
+ strncpy( temp, control->sim_name, TEMP_SIZE );
strcat( temp, ".out" );
- out_control->out = fopen( temp, "w" );
+ out_control->out = sfopen( temp, "w" );
fprintf( out_control->out, "%-6s%16s%16s%16s%11s%11s%13s%13s%13s\n",
"step", "total energy", "poten. energy", "kin. energy",
"temp.", "target", "volume", "press.", "target" );
fflush( out_control->out );
/* Init potentials file */
- strcpy( temp, control->sim_name );
+ strncpy( temp, control->sim_name, TEMP_SIZE );
strcat( temp, ".pot" );
- out_control->pot = fopen( temp, "w" );
+ out_control->pot = sfopen( temp, "w" );
fprintf( out_control->pot,
"%-6s%13s%13s%13s%13s%13s%13s%13s%13s%13s%13s%13s\n",
"step", "ebond", "eatom", "elp", "eang", "ecoa", "ehb",
@@ -857,181 +861,204 @@ static void Init_Out_Controls( reax_system *system, control_params *control,
fflush( out_control->pot );
/* Init log file */
- strcpy( temp, control->sim_name );
+ strncpy( temp, control->sim_name, TEMP_SIZE );
strcat( temp, ".log" );
- out_control->log = fopen( temp, "w" );
+ out_control->log = sfopen( temp, "w" );
fprintf( out_control->log, "%-6s %10s %10s %10s %10s %10s %10s %10s %10s %10s %10s %10s %10s %10s %10s\n",
"step", "total", "neighbors", "init", "bonded",
"nonbonded", "CM", "CM Sort", "S iters", "Pre Comp", "Pre App",
"S spmv", "S vec ops", "S orthog", "S tsolve" );
}
+ else
+ {
+ out_control->out = NULL;
+ out_control->pot = NULL;
+ out_control->log = NULL;
+ }
/* Init pressure file */
if ( control->ensemble == NPT ||
control->ensemble == iNPT ||
control->ensemble == sNPT )
{
- strcpy( temp, control->sim_name );
+ strncpy( temp, control->sim_name, MAX_STR );
strcat( temp, ".prs" );
- out_control->prs = fopen( temp, "w" );
+ out_control->prs = sfopen( temp, "w" );
fprintf( out_control->prs, "%-6s%13s%13s%13s%13s%13s%13s%13s%13s\n",
"step", "norm_x", "norm_y", "norm_z",
"press_x", "press_y", "press_z", "target_p", "volume" );
fflush( out_control->prs );
}
+ else
+ {
+ out_control->prs = NULL;
+ }
/* Init molecular analysis file */
if ( control->molec_anal )
{
snprintf( temp, TEMP_SIZE, "%.*s.mol", TEMP_SIZE - 5, control->sim_name );
- out_control->mol = fopen( temp, "w" );
+ out_control->mol = sfopen( temp, "w" );
if ( control->num_ignored )
{
snprintf( temp, TEMP_SIZE, "%.*s.ign", TEMP_SIZE - 5, control->sim_name );
- out_control->ign = fopen( temp, "w" );
+ out_control->ign = sfopen( temp, "w" );
}
}
+ else
+ {
+ out_control->mol = NULL;
+ out_control->ign = NULL;
+ }
/* Init electric dipole moment analysis file */
if ( control->dipole_anal )
{
- strcpy( temp, control->sim_name );
+ strncpy( temp, control->sim_name, TEMP_SIZE );
strcat( temp, ".dpl" );
- out_control->dpl = fopen( temp, "w" );
+ out_control->dpl = sfopen( temp, "w" );
fprintf( out_control->dpl,
"Step Molecule Count Avg. Dipole Moment Norm\n" );
fflush( out_control->dpl );
}
+ else
+ {
+ out_control->dpl = NULL;
+ }
/* Init diffusion coef analysis file */
if ( control->diffusion_coef )
{
- strcpy( temp, control->sim_name );
+ strncpy( temp, control->sim_name, TEMP_SIZE );
strcat( temp, ".drft" );
- out_control->drft = fopen( temp, "w" );
+ out_control->drft = sfopen( temp, "w" );
fprintf( out_control->drft, "Step Type Count Avg Squared Disp\n" );
fflush( out_control->drft );
}
+ else
+ {
+ out_control->drft = NULL;
+ }
#ifdef TEST_ENERGY
/* open bond energy file */
- strcpy( temp, control->sim_name );
+ strncpy( temp, control->sim_name, TEMP_SIZE );
strcat( temp, ".ebond" );
- out_control->ebond = fopen( temp, "w" );
+ out_control->ebond = sfopen( temp, "w" );
/* open lone-pair energy file */
- strcpy( temp, control->sim_name );
+ strncpy( temp, control->sim_name, TEMP_SIZE );
strcat( temp, ".elp" );
- out_control->elp = fopen( temp, "w" );
+ out_control->elp = sfopen( temp, "w" );
/* open overcoordination energy file */
- strcpy( temp, control->sim_name );
+ strncpy( temp, control->sim_name, TEMP_SIZE );
strcat( temp, ".eov" );
- out_control->eov = fopen( temp, "w" );
+ out_control->eov = sfopen( temp, "w" );
/* open undercoordination energy file */
- strcpy( temp, control->sim_name );
+ strncpy( temp, control->sim_name, TEMP_SIZE );
strcat( temp, ".eun" );
- out_control->eun = fopen( temp, "w" );
+ out_control->eun = sfopen( temp, "w" );
/* open angle energy file */
- strcpy( temp, control->sim_name );
+ strncpy( temp, control->sim_name, TEMP_SIZE );
strcat( temp, ".eval" );
- out_control->eval = fopen( temp, "w" );
+ out_control->eval = sfopen( temp, "w" );
/* open penalty energy file */
- strcpy( temp, control->sim_name );
+ strncpy( temp, control->sim_name, TEMP_SIZE );
strcat( temp, ".epen" );
- out_control->epen = fopen( temp, "w" );
+ out_control->epen = sfopen( temp, "w" );
/* open coalition energy file */
- strcpy( temp, control->sim_name );
+ strncpy( temp, control->sim_name, TEMP_SIZE );
strcat( temp, ".ecoa" );
- out_control->ecoa = fopen( temp, "w" );
+ out_control->ecoa = sfopen( temp, "w" );
/* open hydrogen bond energy file */
- strcpy( temp, control->sim_name );
+ strncpy( temp, control->sim_name, TEMP_SIZE );
strcat( temp, ".ehb" );
- out_control->ehb = fopen( temp, "w" );
+ out_control->ehb = sfopen( temp, "w" );
/* open torsion energy file */
- strcpy( temp, control->sim_name );
+ strncpy( temp, control->sim_name, TEMP_SIZE );
strcat( temp, ".etor" );
- out_control->etor = fopen( temp, "w" );
+ out_control->etor = sfopen( temp, "w" );
/* open conjugation energy file */
- strcpy( temp, control->sim_name );
+ strncpy( temp, control->sim_name, TEMP_SIZE );
strcat( temp, ".econ" );
- out_control->econ = fopen( temp, "w" );
+ out_control->econ = sfopen( temp, "w" );
/* open vdWaals energy file */
- strcpy( temp, control->sim_name );
+ strncpy( temp, control->sim_name, TEMP_SIZE );
strcat( temp, ".evdw" );
- out_control->evdw = fopen( temp, "w" );
+ out_control->evdw = sfopen( temp, "w" );
/* open coulomb energy file */
- strcpy( temp, control->sim_name );
+ strncpy( temp, control->sim_name, TEMP_SIZE );
strcat( temp, ".ecou" );
- out_control->ecou = fopen( temp, "w" );
+ out_control->ecou = sfopen( temp, "w" );
#endif
#ifdef TEST_FORCES
/* open bond orders file */
- strcpy( temp, control->sim_name );
+ strncpy( temp, control->sim_name, TEMP_SIZE );
strcat( temp, ".fbo" );
- out_control->fbo = fopen( temp, "w" );
+ out_control->fbo = sfopen( temp, "w" );
/* open bond orders derivatives file */
- strcpy( temp, control->sim_name );
+ strncpy( temp, control->sim_name, TEMP_SIZE );
strcat( temp, ".fdbo" );
- out_control->fdbo = fopen( temp, "w" );
+ out_control->fdbo = sfopen( temp, "w" );
/* open bond forces file */
- strcpy( temp, control->sim_name );
+ strncpy( temp, control->sim_name, TEMP_SIZE );
strcat( temp, ".fbond" );
- out_control->fbond = fopen( temp, "w" );
+ out_control->fbond = sfopen( temp, "w" );
/* open lone-pair forces file */
- strcpy( temp, control->sim_name );
+ strncpy( temp, control->sim_name, TEMP_SIZE );
strcat( temp, ".flp" );
- out_control->flp = fopen( temp, "w" );
+ out_control->flp = sfopen( temp, "w" );
/* open overcoordination forces file */
- strcpy( temp, control->sim_name );
+ strncpy( temp, control->sim_name, TEMP_SIZE );
strcat( temp, ".fatom" );
- out_control->fatom = fopen( temp, "w" );
+ out_control->fatom = sfopen( temp, "w" );
/* open angle forces file */
- strcpy( temp, control->sim_name );
+ strncpy( temp, control->sim_name, TEMP_SIZE );
strcat( temp, ".f3body" );
- out_control->f3body = fopen( temp, "w" );
+ out_control->f3body = sfopen( temp, "w" );
/* open hydrogen bond forces file */
- strcpy( temp, control->sim_name );
+ strncpy( temp, control->sim_name, TEMP_SIZE );
strcat( temp, ".fhb" );
- out_control->fhb = fopen( temp, "w" );
+ out_control->fhb = sfopen( temp, "w" );
/* open torsion forces file */
- strcpy( temp, control->sim_name );
+ strncpy( temp, control->sim_name, TEMP_SIZE );
strcat( temp, ".f4body" );
- out_control->f4body = fopen( temp, "w" );
+ out_control->f4body = sfopen( temp, "w" );
/* open nonbonded forces file */
- strcpy( temp, control->sim_name );
+ strncpy( temp, control->sim_name, TEMP_SIZE );
strcat( temp, ".fnonb" );
- out_control->fnonb = fopen( temp, "w" );
+ out_control->fnonb = sfopen( temp, "w" );
/* open total force file */
- strcpy( temp, control->sim_name );
+ strncpy( temp, control->sim_name, TEMP_SIZE );
strcat( temp, ".ftot" );
- out_control->ftot = fopen( temp, "w" );
+ out_control->ftot = sfopen( temp, "w" );
/* open coulomb forces file */
- strcpy( temp, control->sim_name );
+ strncpy( temp, control->sim_name, TEMP_SIZE );
strcat( temp, ".ftot2" );
- out_control->ftot2 = fopen( temp, "w" );
+ out_control->ftot2 = sfopen( temp, "w" );
#endif
@@ -1406,31 +1433,19 @@ static void Finalize_Out_Controls( reax_system *system, control_params *control,
/* close trajectory file */
if ( out_control->write_steps > 0 )
{
- if ( out_control->trj )
- {
- fclose( out_control->trj );
- }
+ sfclose( out_control->trj, "Finalize_Out_Controls::out_control->trj" );
}
if ( out_control->energy_update_freq > 0 )
{
/* close out file */
- if ( out_control->out )
- {
- fclose( out_control->out );
- }
+ sfclose( out_control->out, "Finalize_Out_Controls::out_control->out" );
/* close potentials file */
- if ( out_control->pot )
- {
- fclose( out_control->pot );
- }
+ sfclose( out_control->pot, "Finalize_Out_Controls::out_control->pot" );
/* close log file */
- if ( out_control->log )
- {
- fclose( out_control->log );
- }
+ sfclose( out_control->log, "Finalize_Out_Controls::out_control->log" );
}
/* close pressure file */
@@ -1438,171 +1453,104 @@ static void Finalize_Out_Controls( reax_system *system, control_params *control,
control->ensemble == iNPT ||
control->ensemble == sNPT )
{
- if ( out_control->prs )
- {
- fclose( out_control->prs );
- }
+ sfclose( out_control->prs, "Finalize_Out_Controls::out_control->prs" );
}
- /* close molecular analysis file */
if ( control->molec_anal )
{
- fclose( out_control->mol );
+ /* close molecular analysis file */
+ sfclose( out_control->mol, "Finalize_Out_Controls::out_control->mol" );
+
+ if ( control->num_ignored )
+ {
+ sfclose( out_control->ign, "Finalize_Out_Controls::out_control->ign" );
+ }
}
/* close electric dipole moment analysis file */
if ( control->dipole_anal )
{
- fclose( out_control->dpl );
+ sfclose( out_control->dpl, "Finalize_Out_Controls::out_control->dpl" );
}
/* close diffusion coef analysis file */
if ( control->diffusion_coef )
{
- fclose( out_control->drft );
+ sfclose( out_control->drft, "Finalize_Out_Controls::out_control->drft" );
}
#ifdef TEST_ENERGY
/* close bond energy file */
- if ( out_control->ebond )
- {
- fclose( out_control->ebond );
- }
+ sfclose( out_control->ebond, "Finalize_Out_Controls::out_control->ebond" );
/* close lone-pair energy file */
- if ( out_control->help )
- {
- fclose( out_control->elp );
- }
+ sfclose( out_control->elp, "Finalize_Out_Controls::out_control->elp" );
/* close overcoordination energy file */
- if ( out_control->eov )
- {
- fclose( out_control->eov );
- }
+ sfclose( out_control->eov, "Finalize_Out_Controls::out_control->eov" );
/* close undercoordination energy file */
- if ( out_control->eun )
- {
- fclose( out_control->eun );
- }
+ sfclose( out_control->eun, "Finalize_Out_Controls::out_control->eun" );
/* close angle energy file */
- if ( out_control->eval )
- {
- fclose( out_control->eval );
- }
+ sfclose( out_control->eval, "Finalize_Out_Controls::out_control->eval" );
/* close penalty energy file */
- if ( out_control->epen )
- {
- fclose( out_control->epen );
- }
+ sfclose( out_control->epen, "Finalize_Out_Controls::out_control->epen" );
/* close coalition energy file */
- if ( out_control->ecoa )
- {
- fclose( out_control->ecoa );
- }
+ sfclose( out_control->ecoa, "Finalize_Out_Controls::out_control->ecoa" );
/* close hydrogen bond energy file */
- if ( out_control->ehb )
- {
- fclose( out_control->ehb );
- }
+ sfclose( out_control->ehb, "Finalize_Out_Controls::out_control->ehb" );
/* close torsion energy file */
- if ( out_control->etor )
- {
- fclose( out_control->etor );
- }
+ sfclose( out_control->etor, "Finalize_Out_Controls::out_control->etor" );
/* close conjugation energy file */
- if ( out_control->econ )
- {
- fclose( out_control->econ );
- }
+ sfclose( out_control->econ, "Finalize_Out_Controls::out_control->econ" );
/* close vdWaals energy file */
- if ( out_control->evdw )
- {
- fclose( out_control->evdw );
- }
+ sfclose( out_control->evdw, "Finalize_Out_Controls::out_control->evdw" );
/* close coulomb energy file */
- if ( out_control->ecou )
- {
- fclose( out_control->ecou );
- }
+ sfclose( out_control->ecou, "Finalize_Out_Controls::out_control->ecou" );
#endif
#ifdef TEST_FORCES
/* close bond orders file */
- if ( out_control->fbo )
- {
- fclose( out_control->fbo );
- }
+ sfclose( out_control->fbo, "Finalize_Out_Controls::out_control->fbo" );
/* close bond orders derivatives file */
- if ( out_control->fdbo )
- {
- fclose( out_control->fdbo );
- }
+ sfclose( out_control->fdbo, "Finalize_Out_Controls::out_control->fdbo" );
/* close bond forces file */
- if ( out_control->fbond )
- {
- fclose( out_control->fbond );
- }
+ sfclose( out_control->fbond, "Finalize_Out_Controls::out_control->fbond" );
/* close lone-pair forces file */
- if ( out_control->flp )
- {
- fclose( out_control->flp );
- }
+ sfclose( out_control->flp, "Finalize_Out_Controls::out_control->flp" );
/* close overcoordination forces file */
- if ( out_control->fatom )
- {
- fclose( out_control->fatom );
- }
+ sfclose( out_control->fatom, "Finalize_Out_Controls::out_control->fatom" );
/* close angle forces file */
- if ( out_control->f3body )
- {
- fclose( out_control->f3body );
- }
+ sfclose( out_control->f3body, "Finalize_Out_Controls::out_control->f3body" );
/* close hydrogen bond forces file */
- if ( out_control->fhb )
- {
- fclose( out_control->fhb );
- }
+ sfclose( out_control->fhb, "Finalize_Out_Controls::out_control->fhb" );
/* close torsion forces file */
- if ( out_control->f4body )
- {
- fclose( out_control->f4body );
- }
+ sfclose( out_control->f4body, "Finalize_Out_Controls::out_control->f4body" );
/* close nonbonded forces file */
- if ( out_control->fnonb )
- {
- fclose( out_control->fnonb );
- }
+ sfclose( out_control->fnonb, "Finalize_Out_Controls::out_control->fnonb" );
/* close total force file */
- if ( out_control->ftot )
- {
- fclose( out_control->ftot );
- }
+ sfclose( out_control->ftot, "Finalize_Out_Controls::out_control->ftot" );
/* close coulomb forces file */
- if ( out_control->ftot2 )
- {
- fclose( out_control->ftot2 );
- }
+ sfclose( out_control->ftot2, "Finalize_Out_Controls::out_control->ftot2" );
#endif
}
diff --git a/sPuReMD/src/print_utils.c b/sPuReMD/src/print_utils.c
index b019dfa0fb7514636d7213c7859be800ec790951..d099c0be0d3bff33667042caaea0f27a97e1ac78 100644
--- a/sPuReMD/src/print_utils.c
+++ b/sPuReMD/src/print_utils.c
@@ -407,7 +407,7 @@ void Print_Near_Neighbors( reax_system *system, control_params *control,
reax_list *near_nbrs = &((*lists)[NEAR_NBRS]);
snprintf( fname, MAX_STR, "%.*s.near_nbrs", MAX_STR - 11, control->sim_name );
- fout = fopen( fname, "w" );
+ fout = sfopen( fname, "w" );
for ( i = 0; i < system->N; ++i )
{
@@ -427,7 +427,7 @@ void Print_Near_Neighbors( reax_system *system, control_params *control,
}
}
- fclose( fout );
+ sfclose( fout, "Print_Near_Neighbors::fout" );
}
@@ -441,7 +441,7 @@ void Print_Near_Neighbors2( reax_system *system, control_params *control,
reax_list *near_nbrs = &((*lists)[NEAR_NBRS]);
snprintf( fname, MAX_STR, "%.*s.near_nbrs_lgj", MAX_STR - 15, control->sim_name );
- fout = fopen( fname, "w" );
+ fout = sfopen( fname, "w" );
for ( i = 0; i < system->N; ++i )
{
@@ -462,7 +462,7 @@ void Print_Near_Neighbors2( reax_system *system, control_params *control,
fprintf( fout, "\n");
}
- fclose( fout );
+ sfclose( fout, "Print_Near_Neighbors2::fout" );
}
@@ -476,7 +476,7 @@ void Print_Far_Neighbors( reax_system *system, control_params *control,
reax_list *far_nbrs = &(*lists)[FAR_NBRS];
snprintf( fname, MAX_STR, "%.*s.far_nbrs", MAX_STR - 10, control->sim_name );
- fout = fopen( fname, "w" );
+ fout = sfopen( fname, "w" );
for ( i = 0; i < system->N; ++i )
{
@@ -502,7 +502,7 @@ void Print_Far_Neighbors( reax_system *system, control_params *control,
}
}
- fclose( fout );
+ sfclose( fout, "Print_Far_Neighbors::fout" );
}
@@ -521,7 +521,7 @@ void Print_Far_Neighbors2( reax_system *system, control_params *control,
reax_list *far_nbrs = &(*lists)[FAR_NBRS];
snprintf( fname, MAX_STR, "%.*s.far_nbrs_lgj", MAX_STR - 14, control->sim_name );
- fout = fopen( fname, "w" );
+ fout = sfopen( fname, "w" );
int num = 0;
int temp[500];
@@ -542,7 +542,7 @@ void Print_Far_Neighbors2( reax_system *system, control_params *control,
fprintf( fout, "\n");
}
- fclose( fout );
+ sfclose( fout, "Print_Far_Neighbors2::fout" );
}
@@ -703,7 +703,7 @@ void Print_Linear_System( reax_system *system, control_params *control,
FILE *out;
snprintf( fname, 100, "%.*s.state%10d.out", 79, control->sim_name, step );
- out = fopen( fname, "w" );
+ out = sfopen( fname, "w" );
for ( i = 0; i < system->N_cm; i++ )
fprintf( out, "%6d%2d%24.15e%24.15e%24.15e%24.15e%24.15e%24.15e%24.15e\n",
@@ -712,16 +712,16 @@ void Print_Linear_System( reax_system *system, control_params *control,
system->atoms[i].x[2],
workspace->s[0][i], workspace->b_s[i],
workspace->t[0][i], workspace->b_t[i] );
- fclose( out );
+ sfclose( out, "Print_Linear_System::out" );
// snprintf( fname, 100, "x2_%d", step );
- // out = fopen( fname, "w" );
+ // out = sfopen( fname, "w" );
// for( i = 0; i < system->N; i++ )
// fprintf( out, "%g\n", workspace->s_t[i+system->N] );
- // fclose( out );
+ // sfclose( out, "Print_Linear_System::out" );
snprintf( fname, 100, "%.*s.H%10d.out", 83, control->sim_name, step );
- out = fopen( fname, "w" );
+ out = sfopen( fname, "w" );
H = workspace->H;
for ( i = 0; i < system->N_cm; ++i )
@@ -741,10 +741,10 @@ void Print_Linear_System( reax_system *system, control_params *control,
workspace->orig_id[i], workspace->orig_id[i], H->val[j] );
}
- fclose( out );
+ sfclose( out, "Print_Linear_System::out" );
snprintf( fname, 100, "%.*s.H_sp%10d.out", 80, control->sim_name, step );
- out = fopen( fname, "w" );
+ out = sfopen( fname, "w" );
H = workspace->H_sp;
for ( i = 0; i < system->N_cm; ++i )
@@ -764,19 +764,19 @@ void Print_Linear_System( reax_system *system, control_params *control,
workspace->orig_id[i], workspace->orig_id[i], H->val[j] );
}
- fclose( out );
+ sfclose( out, "Print_Linear_System::out" );
/*snprintf( fname, 100, "%.*s.b_s%10d", 84, control->sim_name, step );
- out = fopen( fname, "w" );
+ out = sfopen( fname, "w" );
for( i = 0; i < system->N; i++ )
fprintf( out, "%12.7f\n", workspace->b_s[i] );
- fclose( out );
+ sfclose( out, "Print_Linear_System::out" );
snprintf( fname, 100, "%.*s.b_t%10d", 84, control->sim_name, step );
- out = fopen( fname, "w" );
+ out = sfopen( fname, "w" );
for( i = 0; i < system->N; i++ )
fprintf( out, "%12.7f\n", workspace->b_t[i] );
- fclose( out );*/
+ sfclose( out, "Print_Linear_System::out" );*/
}
@@ -788,7 +788,7 @@ void Print_Charges( reax_system *system, control_params *control,
FILE *fout;
snprintf( fname, 100, "%.*s.q%10d", 87, control->sim_name, step );
- fout = fopen( fname, "w" );
+ fout = sfopen( fname, "w" );
for ( i = 0; i < system->N; ++i )
{
@@ -797,7 +797,7 @@ void Print_Charges( reax_system *system, control_params *control,
workspace->s[0][i], workspace->t[0][i], system->atoms[i].q );
}
- fclose( fout );
+ sfclose( fout, "Print_Charges::fout" );
}
@@ -841,11 +841,11 @@ void Print_Sparse_Matrix2( sparse_matrix *A, char *fname, char *mode )
if ( mode == NULL )
{
- f = fopen( fname, "w" );
+ f = sfopen( fname, "w" );
}
else
{
- f = fopen( fname, mode );
+ f = sfopen( fname, mode );
}
for ( i = 0; i < A->n; ++i )
@@ -863,7 +863,7 @@ void Print_Sparse_Matrix2( sparse_matrix *A, char *fname, char *mode )
fprintf( f, "%6d %6d %24.15e\n", i + 1, A->j[A->start[i + 1] - 1] + 1, A->val[A->start[i + 1] - 1] );
}
- fclose( f );
+ sfclose( f, "Print_Sparse_Matrix2::f" );
}
@@ -874,7 +874,7 @@ void Print_Sparse_Matrix_Binary( sparse_matrix *A, char *fname )
int i, j, temp;
FILE *f;
- f = fopen( fname, "wb" );
+ f = sfopen( fname, "wb" );
/* header: # rows, # nonzeros */
fwrite( &(A->n), sizeof(unsigned int), 1, f );
@@ -900,7 +900,7 @@ void Print_Sparse_Matrix_Binary( sparse_matrix *A, char *fname )
}
}
- fclose(f);
+ sfclose( f, "Print_Sparse_Matrix_Binary::f" );
}
@@ -909,7 +909,7 @@ void Print_Bonds( reax_system *system, reax_list *bonds, char *fname )
int i, pj;
bond_data *pbond;
bond_order_data *bo_ij;
- FILE *f = fopen( fname, "w" );
+ FILE *f = sfopen( fname, "w" );
for ( i = 0; i < system->N; ++i )
{
@@ -925,14 +925,14 @@ void Print_Bonds( reax_system *system, reax_list *bonds, char *fname )
}
}
- fclose( f );
+ sfclose( f, "Print_Bonds::f" );
}
void Print_Bond_List2( reax_system *system, reax_list *bonds, char *fname )
{
int i, j, id_i, id_j, nbr, pj;
- FILE *f = fopen( fname, "w" );
+ FILE *f = sfopen( fname, "w" );
int temp[500];
int num = 0;
@@ -970,7 +970,7 @@ Print_XYZ_Serial( reax_system* system, static_storage *workspace )
int i;
snprintf( fname, 100, "READ_PDB.0" );
- fout = fopen( fname, "w" );
+ fout = sfopen( fname, "w" );
for ( i = 0; i < system->N; i++ )
{
@@ -981,6 +981,6 @@ Print_XYZ_Serial( reax_system* system, static_storage *workspace )
p[2] = system->atoms[i].x[2] );
}
- fclose( fout );
+ sfclose( fout, "Print_XYZ_Serial::fout" );
}
#endif
diff --git a/sPuReMD/src/restart.c b/sPuReMD/src/restart.c
index 23b00726848007bb0e03193437d9697937dc396e..9beaaa33107b724899d60c4d9ef6ad6707bf0f1b 100644
--- a/sPuReMD/src/restart.c
+++ b/sPuReMD/src/restart.c
@@ -37,7 +37,7 @@ void Write_Binary_Restart( reax_system *system, control_params *control,
restart_atom res_data;
snprintf( fname, MAX_STR + 8, "%s.res%d", control->sim_name, data->step );
- fres = fopen( fname, "wb" );
+ fres = sfopen( fname, "wb" );
res_header.step = data->step;
res_header.N = system->N;
@@ -60,7 +60,7 @@ void Write_Binary_Restart( reax_system *system, control_params *control,
fwrite( &res_data, sizeof(restart_atom), 1, fres );
}
- fclose( fres );
+ sfclose( fres, "Write_Binary_Restart::fres" );
#if defined(DEBUG_FOCUS)
fprintf( stderr, "write restart - " );
@@ -69,8 +69,8 @@ void Write_Binary_Restart( reax_system *system, control_params *control,
void Read_Binary_Restart( const char * const fname, reax_system *system,
- control_params *control, simulation_data *data,
- static_storage *workspace )
+ control_params *control, simulation_data *data,
+ static_storage *workspace )
{
int i;
FILE *fres;
@@ -78,9 +78,9 @@ void Read_Binary_Restart( const char * const fname, reax_system *system,
restart_header res_header;
restart_atom res_data;
- fres = fopen( fname, "rb" );
+ fres = sfopen( fname, "rb" );
- fread(&res_header, sizeof(restart_header), 1, fres);
+ fread( &res_header, sizeof(restart_header), 1, fres );
data->prev_steps = res_header.step;
system->N = res_header.N;
data->therm.T = res_header.T;
@@ -109,7 +109,9 @@ void Read_Binary_Restart( const char * const fname, reax_system *system,
workspace->map_serials = (int*) scalloc( MAX_ATOM_ID, sizeof(int),
"Read_Binary_Restart::workspace->map_serials" );
for ( i = 0; i < MAX_ATOM_ID; ++i )
+ {
workspace->map_serials[i] = -1;
+ }
workspace->orig_id = (int*) scalloc( system->N, sizeof(int),
"Read_Binary_Restart::workspace->orig_id" );
@@ -132,7 +134,7 @@ void Read_Binary_Restart( const char * const fname, reax_system *system,
p_atom = &( system->atoms[i] );
p_atom->type = res_data.type;
- strcpy( p_atom->name, res_data.name );
+ strncpy( p_atom->name, res_data.name, 8 );
rvec_Copy( p_atom->x, res_data.x );
rvec_Copy( p_atom->v, res_data.v );
}
@@ -141,7 +143,7 @@ void Read_Binary_Restart( const char * const fname, reax_system *system,
fprintf( stderr, "system->N: %d, i: %d\n", system->N, i );
#endif
- fclose( fres );
+ sfclose( fres, "Read_Binary_Restart::fres" );
data->step = data->prev_steps;
// nsteps is updated based on the number of steps in the previous run
@@ -158,7 +160,7 @@ void Write_ASCII_Restart( reax_system *system, control_params *control,
reax_atom *p_atom;
snprintf( fname, MAX_STR + 8, "%s.res%d", control->sim_name, data->step );
- fres = fopen( fname, "w" );
+ fres = sfopen( fname, "w" );
fprintf( fres, RESTART_HEADER,
data->step, system->N, data->therm.T, data->therm.xi,
@@ -166,7 +168,7 @@ void Write_ASCII_Restart( reax_system *system, control_params *control,
system->box.box[0][0], system->box.box[0][1], system->box.box[0][2],
system->box.box[1][0], system->box.box[1][1], system->box.box[1][2],
system->box.box[2][0], system->box.box[2][1], system->box.box[2][2]);
- fflush(fres);
+ fflush( fres );
for ( i = 0; i < system->N; ++i )
{
@@ -177,7 +179,7 @@ void Write_ASCII_Restart( reax_system *system, control_params *control,
p_atom->v[0], p_atom->v[1], p_atom->v[2] );
}
- fclose( fres );
+ sfclose( fres, "Write_ASCII_Restart::fres" );
#if defined(DEBUG_FOCUS)
fprintf( stderr, "write restart - " );
@@ -186,14 +188,14 @@ void Write_ASCII_Restart( reax_system *system, control_params *control,
void Read_ASCII_Restart( const char * const fname, reax_system *system,
- control_params *control, simulation_data *data,
- static_storage *workspace )
+ control_params *control, simulation_data *data,
+ static_storage *workspace )
{
int i;
FILE *fres;
reax_atom *p_atom;
- fres = fopen( fname, "r" );
+ fres = sfopen( fname, "r" );
/* header */
fscanf( fres, READ_RESTART_HEADER,
@@ -249,7 +251,7 @@ void Read_ASCII_Restart( const char * const fname, reax_system *system,
workspace->map_serials[workspace->orig_id[i]] = i;
}
- fclose( fres );
+ sfclose( fres, "Read_ASCII_Restart::fres" );
data->step = data->prev_steps;
// nsteps is updated based on the number of steps in the previous run
@@ -258,11 +260,21 @@ void Read_ASCII_Restart( const char * const fname, reax_system *system,
void Write_Restart( reax_system *system, control_params *control,
- simulation_data *data, static_storage *workspace,
- output_controls *out_control )
+ simulation_data *data, static_storage *workspace,
+ output_controls *out_control )
{
- if ( out_control->restart_format == WRITE_ASCII )
+ switch ( out_control->restart_format )
+ {
+ case WRITE_ASCII:
Write_ASCII_Restart( system, control, data, workspace );
- else if ( out_control->restart_format == WRITE_BINARY )
+ break;
+
+ case WRITE_BINARY:
Write_Binary_Restart( system, control, data, workspace );
+ break;
+
+ default:
+ fprintf( stderr, "[ERROR] invalid restart format\n" );
+ exit( INVALID_INPUT );
+ }
}
diff --git a/sPuReMD/src/single_body_interactions.c b/sPuReMD/src/single_body_interactions.c
index cf22fb6a3b4e1bf4d1929babbd48e538213d5f75..42b97041c0a5f8274b479dc0fb78e51d8cb66cb6 100644
--- a/sPuReMD/src/single_body_interactions.c
+++ b/sPuReMD/src/single_body_interactions.c
@@ -91,7 +91,7 @@ void LonePair_OverUnder_Coordination_Energy( reax_system *system, control_params
/* correction for C2 */
if ( system->reaxprm.gp.l[5] > 0.001 &&
- !strcmp( system->reaxprm.sbp[type_i].name, "C" ) )
+ !strncmp( system->reaxprm.sbp[type_i].name, "C", 15 ) )
{
for ( pj = Start_Index(i, bonds); pj < End_Index(i, bonds); ++pj )
{
@@ -100,7 +100,7 @@ void LonePair_OverUnder_Coordination_Energy( reax_system *system, control_params
j = bonds->select.bond_list[pj].nbr;
type_j = system->atoms[j].type;
- if ( !strcmp( system->reaxprm.sbp[type_j].name, "C" ) )
+ if ( !strncmp( system->reaxprm.sbp[type_j].name, "C", 15 ) )
{
twbp = &( system->reaxprm.tbp[type_i][type_j]);
bo_ij = &( bonds->select.bond_list[pj].bo_data );
diff --git a/sPuReMD/src/spuremd.c b/sPuReMD/src/spuremd.c
index 968847ceabafaf8ed35a0de77bf6ce7669404b1e..6340fd8fe1241c26372dc38cbd60fecc976e1ea8 100644
--- a/sPuReMD/src/spuremd.c
+++ b/sPuReMD/src/spuremd.c
@@ -90,18 +90,8 @@ static void Read_System( const char * const geo_file,
{
FILE *ffield, *ctrl;
- if ( (ffield = fopen( ffield_file, "r" )) == NULL )
- {
- fprintf( stderr, "[ERROR] Error opening the ffield file!\n" );
- fprintf( stderr, " [INFO] (%s)\n", ffield_file );
- exit( FILE_NOT_FOUND );
- }
- if ( (ctrl = fopen( control_file, "r" )) == NULL )
- {
- fprintf( stderr, "[ERROR] Error opening the ffield file!\n" );
- fprintf( stderr, " [INFO] (%s)\n", control_file );
- exit( FILE_NOT_FOUND );
- }
+ ffield = sfopen( ffield_file, "r" );
+ ctrl = sfopen( control_file, "r" );
/* ffield file */
Read_Force_Field( ffield, &(system->reaxprm) );
@@ -138,8 +128,8 @@ static void Read_System( const char * const geo_file,
exit( INVALID_GEO );
}
- fclose( ffield );
- fclose( ctrl );
+ sfclose( ffield, "Read_System::ffield" );
+ sfclose( ctrl, "Read_System::ctrl" );
#if defined(DEBUG_FOCUS)
fprintf( stderr, "input files have been read...\n" );
@@ -301,7 +291,6 @@ int simulate( const void * const handle )
if ( spmd_handle->out_control->write_steps > 0 && spmd_handle->output_enabled == TRUE )
{
- fclose( spmd_handle->out_control->trj );
Write_PDB( spmd_handle->system, &(spmd_handle->lists[BONDS]), spmd_handle->data,
spmd_handle->control, spmd_handle->workspace, spmd_handle->out_control );
}
diff --git a/sPuReMD/src/tool_box.c b/sPuReMD/src/tool_box.c
index 1cc6ba72f6fa14fe5d501194932f5f1a8434e633..c46f8615ee0bc63f031011ca92bc726ebb49b849 100644
--- a/sPuReMD/src/tool_box.c
+++ b/sPuReMD/src/tool_box.c
@@ -22,7 +22,7 @@
#include "tool_box.h"
#include <ctype.h>
-#include <sys/time.h>
+#include <time.h>
/************** taken from box.c **************/
@@ -299,33 +299,33 @@ void Trim_Spaces( char * const element, const size_t size )
/************ from system_props.c *************/
real Get_Time( )
{
- struct timeval tim;
+ int ret;
+ struct timespec t;
- gettimeofday(&tim, NULL );
- return ( tim.tv_sec + (tim.tv_usec / 1000000.0) );
+ ret = clock_gettime( CLOCK_MONOTONIC, &t );
+
+ if ( ret != 0 )
+ {
+ fprintf( stderr, "[WARNING] non-zero error in measuring time\n" );
+ }
+
+ return t.tv_sec + t.tv_nsec / 1.0e9;
}
real Get_Timing_Info( real t_start )
{
- struct timeval tim;
- real t_end;
-
- gettimeofday(&tim, NULL );
- t_end = tim.tv_sec + (tim.tv_usec / 1000000.0);
- return (t_end - t_start);
-}
+ int ret;
+ struct timespec t_end;
+ ret = clock_gettime( CLOCK_MONOTONIC, &t_end );
-void Update_Timing_Info( real *t_start, real *timing )
-{
- struct timeval tim;
- real t_end;
+ if ( ret != 0 )
+ {
+ fprintf( stderr, "[WARNING] non-zero error in measuring time\n" );
+ }
- gettimeofday(&tim, NULL );
- t_end = tim.tv_sec + (tim.tv_usec / 1000000.0);
- *timing += (t_end - *t_start);
- *t_start = t_end;
+ return t_end.tv_sec + t_end.tv_nsec / 1.0e9 - t_start;
}
@@ -336,7 +336,7 @@ int Get_Atom_Type( reax_interaction *reax_param, char *s )
for ( i = 0; i < reax_param->num_atom_types; ++i )
{
- if ( !strcmp( reax_param->sbp[i].name, s ) )
+ if ( !strncmp( reax_param->sbp[i].name, s, 15 ) )
{
return i;
}
@@ -422,7 +422,7 @@ int Tokenize( char* s, char*** tok )
*
* returns: ptr to allocated memory
* */
-void *smalloc( size_t n, const char *name )
+void * smalloc( size_t n, const char *name )
{
void *ptr;
@@ -463,7 +463,7 @@ void *smalloc( size_t n, const char *name )
*
* returns: ptr to reallocated memory
* */
-void* srealloc( void *ptr, size_t n, const char *name )
+void * srealloc( void *ptr, size_t n, const char *name )
{
void *new_ptr;
@@ -513,7 +513,7 @@ void* srealloc( void *ptr, size_t n, const char *name )
*
* returns: ptr to allocated memory, all bits initialized to zeros
* */
-void* scalloc( size_t n, size_t size, const char *name )
+void * scalloc( size_t n, size_t size, const char *name )
{
void *ptr;
@@ -547,7 +547,6 @@ void* scalloc( size_t n, size_t size, const char *name )
}
-/* safe free */
/* Safe wrapper around libc free
*
* ptr: pointer to dynamically allocated memory which will be deallocated
@@ -571,3 +570,66 @@ void sfree( void *ptr, const char *name )
free( ptr );
}
+
+
+/* Safe wrapper around libc fopen
+ *
+ * fname: name of file to be opened
+ * mode: mode in which to open file
+ * */
+FILE * sfopen( const char * fname, const char * mode )
+{
+ FILE * ptr;
+
+ if ( fname == NULL )
+ {
+ fprintf( stderr, "[ERROR] trying to open file\n" );
+ fprintf( stderr, " [INFO] NULL file name\n" );
+ exit( INVALID_INPUT );
+ }
+ if ( mode == NULL )
+ {
+ fprintf( stderr, "[ERROR] trying to open file\n" );
+ fprintf( stderr, " [INFO] NULL mode\n" );
+ exit( INVALID_INPUT );
+ }
+
+ ptr = fopen( fname, mode );
+
+ if ( ptr == NULL )
+ {
+ fprintf( stderr, "[ERROR] failed to open file %s with mode %s\n",
+ fname, mode );
+ exit( INVALID_INPUT );
+ }
+
+ return ptr;
+}
+
+
+/* Safe wrapper around libc fclose
+ *
+ * fname: name of file to be opened
+ * mode: mode in which to open file
+ * msg: message to be printed in case of error
+ * */
+void sfclose( FILE * fp, const char * msg )
+{
+ int ret;
+
+ if ( fp == NULL )
+ {
+ fprintf( stderr, "[WARNING] trying to close NULL file pointer. Returning...\n" );
+ fprintf( stderr, " [INFO] %s\n", msg );
+ return;
+ }
+
+ ret = fclose( fp );
+
+ if ( ret != 0 )
+ {
+ fprintf( stderr, "[ERROR] error detected when closing file\n" );
+ fprintf( stderr, " [INFO] %s\n", msg );
+ exit( INVALID_INPUT );
+ }
+}
diff --git a/sPuReMD/src/tool_box.h b/sPuReMD/src/tool_box.h
index 29781e00dc637633ae8775dd7ac1d1237711986e..a5d181e853157387bf74ad7b1564e60d2a6c0d0f 100644
--- a/sPuReMD/src/tool_box.h
+++ b/sPuReMD/src/tool_box.h
@@ -63,8 +63,6 @@ real Get_Time( );
real Get_Timing_Info( real );
-void Update_Timing_Info( real*, real* );
-
/* from io_tools.h */
int Get_Atom_Type( reax_interaction*, char* );
@@ -87,5 +85,8 @@ void* scalloc( size_t, size_t, const char * );
void sfree( void *, const char * );
+FILE * sfopen( const char *, const char * );
+
+void sfclose( FILE *, const char * );
#endif
diff --git a/sPuReMD/src/two_body_interactions.c b/sPuReMD/src/two_body_interactions.c
index d83db947b82936725eb37ad9ac9deb83029abbed..9a04bb71651d64aeafb43e2a6afa1ecafdea6b1c 100644
--- a/sPuReMD/src/two_body_interactions.c
+++ b/sPuReMD/src/two_body_interactions.c
@@ -395,7 +395,7 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
data->E_vdW = e_vdW_total;
data->E_Ele = e_ele_total;
- // fclose( fout );
+ // sfclose( fout, "vdW_Coulomb_Energy::fout" );
// fprintf( stderr, "nonbonded: ext_press (%24.15e %24.15e %24.15e)\n",
// data->ext_press[0], data->ext_press[1], data->ext_press[2] );