From e01ba7a4c54cc1cf094c5af808042b7fdbac3517 Mon Sep 17 00:00:00 2001
From: "Kurt A. O'Hearn" <ohearnk@msu.edu>
Date: Mon, 12 Nov 2018 09:55:38 -0500
Subject: [PATCH] PuReMD-old: backport control file format changes.
---
PuReMD/src/control.c | 777 ++++++++++++++++++++++-----------------
PuReMD/src/forces.c | 4 +-
PuReMD/src/io_tools.c | 4 +-
PuReMD/src/qEq.c | 13 +-
PuReMD/src/reax_types.h | 299 ++++++++++++---
PuReMD/src/reset_tools.c | 15 +-
PuReMD/src/tool_box.c | 62 ++++
PuReMD/src/tool_box.h | 4 +
PuReMD/src/traj.c | 2 +-
9 files changed, 772 insertions(+), 408 deletions(-)
diff --git a/PuReMD/src/control.c b/PuReMD/src/control.c
index cd95c246..ae9bcee8 100644
--- a/PuReMD/src/control.c
+++ b/PuReMD/src/control.c
@@ -20,12 +20,13 @@
----------------------------------------------------------------------*/
#include "reax_types.h"
+
#if defined(PURE_REAX)
-#include "control.h"
-#include "tool_box.h"
+ #include "control.h"
+ #include "tool_box.h"
#elif defined(LAMMPS_REAX)
-#include "reax_control.h"
-#include "reax_tool_box.h"
+ #include "reax_control.h"
+ #include "reax_tool_box.h"
#endif
@@ -34,49 +35,58 @@ char Read_Control_File( char *control_file, control_params* control,
{
FILE *fp;
char *s, **tmp;
- int c, i, ival;
- real val;
+ int c, i, ival;
+ real val;
- /* open control file */
- if ( (fp = fopen( control_file, "r" ) ) == NULL )
- {
- fprintf( stderr, "error opening the control file! terminating...\n" );
- MPI_Abort( MPI_COMM_WORLD, FILE_NOT_FOUND );
- }
+ fp = sfopen( control_file, "r", "Read_Control_File::fp" );
/* assign default values */
- strcpy( control->sim_name, "simulate" );
- control->ensemble = NVE;
- control->nsteps = 0;
- control->dt = 0.25;
- control->nprocs = 1;
+ strcpy( control->sim_name, "default.sim" );
+ control->ensemble = NVE;
+ control->nsteps = 100;
+ control->dt = 0.25;
+ control->nprocs = 1;
control->procs_by_dim[0] = 1;
control->procs_by_dim[1] = 1;
control->procs_by_dim[2] = 1;
control->geo_format = 1;
- control->restart = 0;
+ control->random_vel = 0;
+ control->restart = 0;
out_control->restart_format = WRITE_BINARY;
out_control->restart_freq = 0;
control->reposition_atoms = 0;
- control->restrict_bonds = 0;
- control->remove_CoM_vel = 25;
- out_control->debug_level = 0;
+ control->restrict_bonds = 0;
+ control->remove_CoM_vel = 25;
+ out_control->debug_level = 0;
out_control->energy_update_freq = 0;
control->reneighbor = 1;
- control->vlist_cut = control->nonb_cut;
control->bond_cut = 5.0;
+ control->vlist_cut = control->nonb_cut;
control->bg_cut = 0.3;
control->thb_cut = 0.001;
control->hbond_cut = 0.0;
control->tabulate = 0;
- control->qeq_freq = 1;
- control->q_err = 1e-6;
- control->refactor = 100;
- control->droptol = 1e-2;;
+ control->charge_method = QEQ_CM;
+ control->charge_freq = 1;
+ control->cm_q_net = 0.0;
+ control->cm_solver_type = CG_S;
+ control->cm_solver_max_iters = 1000;
+ control->cm_solver_restart = 100;
+ control->cm_solver_q_err = 0.000001;
+ control->cm_domain_sparsify_enabled = FALSE;
+ control->cm_init_guess_extrap1 = 3;
+ control->cm_init_guess_extrap2 = 2;
+ control->cm_domain_sparsity = 1.0;
+ control->cm_solver_pre_comp_type = JACOBI_PC;
+ control->cm_solver_pre_comp_sweeps = 3;
+ control->cm_solver_pre_comp_refactor = 100;
+ control->cm_solver_pre_comp_droptol = 0.01;
+ control->cm_solver_pre_app_type = TRI_SOLVE_PA;
+ control->cm_solver_pre_app_jacobi_iters = 50;
control->T_init = 0.;
control->T_final = 300.;
@@ -97,8 +107,8 @@ char Read_Control_File( char *control_file, control_params* control,
out_control->traj_method = REG_TRAJ;
strcpy( out_control->traj_title, "default_title" );
out_control->atom_info = 0;
- out_control->bond_info = 0;
- out_control->angle_info = 0;
+ out_control->bond_info = 1;
+ out_control->angle_info = 1;
control->molecular_analysis = 0;
control->dipole_anal = 0;
@@ -108,339 +118,440 @@ char Read_Control_File( char *control_file, control_params* control,
control->restrict_type = 0;
/* memory allocations */
- s = (char*) malloc(sizeof(char) * MAX_LINE);
- tmp = (char**) malloc(sizeof(char*)*MAX_TOKENS);
+ s = (char*) smalloc( sizeof(char) * MAX_LINE, "Read_Control_File::s", MPI_COMM_WORLD );
+ tmp = (char**) smalloc( sizeof(char*) * MAX_TOKENS, "Read_Control_File::tmp", MPI_COMM_WORLD );
for (i = 0; i < MAX_TOKENS; i++)
- tmp[i] = (char*) malloc(sizeof(char) * MAX_LINE);
+ {
+ tmp[i] = (char*) smalloc( sizeof(char) * MAX_LINE, "Read_Control_File::tmp[i]", MPI_COMM_WORLD );
+ }
/* read control parameters file */
while ( fgets( s, MAX_LINE, fp ) )
{
c = Tokenize( s, &tmp );
- //fprintf( stderr, "%s\n", s );
- if ( strcmp(tmp[0], "simulation_name") == 0 )
- {
- strcpy( control->sim_name, tmp[1] );
- }
- else if ( strcmp(tmp[0], "ensemble_type") == 0 )
- {
- ival = atoi(tmp[1]);
- control->ensemble = ival;
- if ( ival == iNPT || ival == sNPT || ival == NPT )
- control->virial = 1;
- }
- else if ( strcmp(tmp[0], "nsteps") == 0 )
- {
- ival = atoi(tmp[1]);
- control->nsteps = ival;
- }
- else if ( strcmp(tmp[0], "dt") == 0)
- {
- val = atof(tmp[1]);
- control->dt = val * 1.e-3; // convert dt from fs to ps!
- }
- else if ( strcmp(tmp[0], "proc_by_dim") == 0 )
+ if ( c > 0 )
{
- ival = atoi(tmp[1]);
- control->procs_by_dim[0] = ival;
- ival = atoi(tmp[2]);
- control->procs_by_dim[1] = ival;
- ival = atoi(tmp[3]);
- control->procs_by_dim[2] = ival;
+ if ( strcmp(tmp[0], "simulation_name") == 0 )
+ {
+ strcpy( control->sim_name, tmp[1] );
+ }
+ else if ( strcmp(tmp[0], "ensemble_type") == 0 )
+ {
+ ival = atoi(tmp[1]);
+ control->ensemble = ival;
+ }
+ else if ( strcmp(tmp[0], "nsteps") == 0 )
+ {
+ ival = atoi(tmp[1]);
+ control->nsteps = ival;
+ }
+ else if ( strcmp(tmp[0], "dt") == 0)
+ {
+ val = atof(tmp[1]);
+ control->dt = val * 1.e-3; // convert dt from fs to ps!
+ }
+ else if ( strncmp(tmp[0], "gpus_per_node", MAX_LINE) == 0 )
+ {
+ // skip since not applicable to non-gpu distributed memory code
+ ;
+ }
+ else if ( strcmp(tmp[0], "proc_by_dim") == 0 )
+ {
+ if ( c < 4 )
+ {
+ fprintf( stderr, "[ERROR] invalid number of control file parameters (procs_by_dim). terminating!\n" );
+ MPI_Abort( MPI_COMM_WORLD, INVALID_INPUT );
+ }
- control->nprocs = control->procs_by_dim[0] * control->procs_by_dim[1] *
- control->procs_by_dim[2];
- }
- //else if( strcmp(tmp[0], "restart") == 0 ) {
- // ival = atoi(tmp[1]);
- // control->restart = ival;
- //}
- //else if( strcmp(tmp[0], "restart_from") == 0 ) {
- // strcpy( control->restart_from, tmp[1] );
- //}
- else if ( strcmp(tmp[0], "random_vel") == 0 )
- {
- ival = atoi(tmp[1]);
- control->random_vel = ival;
- }
- else if ( strcmp(tmp[0], "restart_format") == 0 )
- {
- ival = atoi(tmp[1]);
- out_control->restart_format = ival;
- }
- else if ( strcmp(tmp[0], "restart_freq") == 0 )
- {
- ival = atoi(tmp[1]);
- out_control->restart_freq = ival;
- }
- else if ( strcmp(tmp[0], "reposition_atoms") == 0 )
- {
- ival = atoi(tmp[1]);
- control->reposition_atoms = ival;
- }
- else if ( strcmp(tmp[0], "restrict_bonds") == 0 )
- {
- ival = atoi( tmp[1] );
- control->restrict_bonds = ival;
- }
- else if ( strcmp(tmp[0], "remove_CoM_vel") == 0 )
- {
- ival = atoi(tmp[1]);
- control->remove_CoM_vel = ival;
- }
- else if ( strcmp(tmp[0], "debug_level") == 0 )
- {
- ival = atoi(tmp[1]);
- out_control->debug_level = ival;
- }
- else if ( strcmp(tmp[0], "energy_update_freq") == 0 )
- {
- ival = atoi(tmp[1]);
- out_control->energy_update_freq = ival;
- }
- else if ( strcmp(tmp[0], "reneighbor") == 0 )
- {
- ival = atoi( tmp[1] );
- control->reneighbor = ival;
- }
- else if ( strcmp(tmp[0], "vlist_buffer") == 0 )
- {
- val = atof(tmp[1]);
- control->vlist_cut = val + control->nonb_cut;
- }
- else if ( strcmp(tmp[0], "nbrhood_cutoff") == 0 )
- {
- val = atof(tmp[1]);
- control->bond_cut = val;
- }
- else if ( strcmp(tmp[0], "bond_graph_cutoff") == 0 )
- {
- val = atof(tmp[1]);
- control->bg_cut = val;
- }
- else if ( strcmp(tmp[0], "thb_cutoff") == 0 )
- {
- val = atof(tmp[1]);
- control->thb_cut = val;
- }
- else if ( strcmp(tmp[0], "hbond_cutoff") == 0 )
- {
- val = atof( tmp[1] );
- control->hbond_cut = val;
- }
- else if ( strcmp(tmp[0], "ghost_cutoff") == 0 )
- {
- val = atof(tmp[1]);
- control->user_ghost_cut = val;
- }
- else if ( strcmp(tmp[0], "tabulate_long_range") == 0 )
- {
- ival = atoi( tmp[1] );
- control->tabulate = ival;
- }
- else if ( strcmp(tmp[0], "qeq_freq") == 0 )
- {
- ival = atoi( tmp[1] );
- control->qeq_freq = ival;
- }
- else if ( strcmp(tmp[0], "q_err") == 0 )
- {
- val = atof( tmp[1] );
- control->q_err = val;
- }
- else if ( strcmp(tmp[0], "ilu_refactor") == 0 )
- {
- ival = atoi( tmp[1] );
- control->refactor = ival;
- }
- else if ( strcmp(tmp[0], "ilu_droptol") == 0 )
- {
- val = atof( tmp[1] );
- control->droptol = val;
- }
- else if ( strcmp(tmp[0], "temp_init") == 0 )
- {
- val = atof(tmp[1]);
- control->T_init = val;
+ ival = atoi(tmp[1]);
+ control->procs_by_dim[0] = ival;
+ ival = atoi(tmp[2]);
+ control->procs_by_dim[1] = ival;
+ ival = atoi(tmp[3]);
+ control->procs_by_dim[2] = ival;
- if ( control->T_init < 0.1 )
- control->T_init = 0.1;
- }
- else if ( strcmp(tmp[0], "temp_final") == 0 )
- {
- val = atof(tmp[1]);
- control->T_final = val;
+ control->nprocs = control->procs_by_dim[0] * control->procs_by_dim[1] *
+ control->procs_by_dim[2];
+ }
+ //else if( strcmp(tmp[0], "restart") == 0 ) {
+ // ival = atoi(tmp[1]);
+ // control->restart = ival;
+ //}
+ //else if( strcmp(tmp[0], "restart_from") == 0 ) {
+ // strcpy( control->restart_from, tmp[1] );
+ //}
+ else if ( strcmp(tmp[0], "random_vel") == 0 )
+ {
+ ival = atoi(tmp[1]);
+ control->random_vel = ival;
+ }
+ else if ( strcmp(tmp[0], "restart_format") == 0 )
+ {
+ ival = atoi(tmp[1]);
+ out_control->restart_format = ival;
+ }
+ else if ( strcmp(tmp[0], "restart_freq") == 0 )
+ {
+ ival = atoi(tmp[1]);
+ out_control->restart_freq = ival;
+ }
+ else if ( strcmp(tmp[0], "reposition_atoms") == 0 )
+ {
+ ival = atoi(tmp[1]);
+ control->reposition_atoms = ival;
+ }
+ else if ( strcmp(tmp[0], "restrict_bonds") == 0 )
+ {
+ ival = atoi( tmp[1] );
+ control->restrict_bonds = ival;
+ }
+ else if ( strcmp(tmp[0], "remove_CoM_vel") == 0 )
+ {
+ ival = atoi(tmp[1]);
+ control->remove_CoM_vel = ival;
+ }
+ else if ( strcmp(tmp[0], "debug_level") == 0 )
+ {
+ ival = atoi(tmp[1]);
+ out_control->debug_level = ival;
+ }
+ else if ( strcmp(tmp[0], "energy_update_freq") == 0 )
+ {
+ ival = atoi(tmp[1]);
+ out_control->energy_update_freq = ival;
+ }
+ else if ( strcmp(tmp[0], "reneighbor") == 0 )
+ {
+ ival = atoi( tmp[1] );
+ control->reneighbor = ival;
+ }
+ else if ( strcmp(tmp[0], "vlist_buffer") == 0 )
+ {
+ val = atof(tmp[1]);
+ control->vlist_cut = val + control->nonb_cut;
+ }
+ else if ( strcmp(tmp[0], "nbrhood_cutoff") == 0 )
+ {
+ val = atof(tmp[1]);
+ control->bond_cut = val;
+ }
+ else if ( strcmp(tmp[0], "bond_graph_cutoff") == 0 )
+ {
+ val = atof(tmp[1]);
+ control->bg_cut = val;
+ }
+ else if ( strcmp(tmp[0], "thb_cutoff") == 0 )
+ {
+ val = atof(tmp[1]);
+ control->thb_cut = val;
+ }
+ else if ( strcmp(tmp[0], "hbond_cutoff") == 0 )
+ {
+ val = atof( tmp[1] );
+ control->hbond_cut = val;
+ }
+ else if ( strcmp(tmp[0], "ghost_cutoff") == 0 )
+ {
+ val = atof(tmp[1]);
+ control->user_ghost_cut = val;
+ }
+ else if ( strcmp(tmp[0], "tabulate_long_range") == 0 )
+ {
+ ival = atoi( tmp[1] );
+ control->tabulate = ival;
+ }
+ else if ( strcmp(tmp[0], "charge_method") == 0 )
+ {
+ ival = atoi( tmp[1] );
+ control->charge_method = ival;
+ }
+ else if ( strcmp(tmp[0], "charge_freq") == 0 )
+ {
+ ival = atoi( tmp[1] );
+ control->charge_freq = ival;
+ }
+ else if ( strcmp(tmp[0], "cm_q_net") == 0 )
+ {
+ val = atof( tmp[1] );
+ control->cm_q_net = val;
+ }
+ else if ( strcmp(tmp[0], "cm_solver_type") == 0 )
+ {
+ ival = atoi( tmp[1] );
+ control->cm_solver_type = ival;
+ }
+ else if ( strcmp(tmp[0], "cm_solver_max_iters") == 0 )
+ {
+ ival = atoi( tmp[1] );
+ control->cm_solver_max_iters = ival;
+ }
+ else if ( strcmp(tmp[0], "cm_solver_restart") == 0 )
+ {
+ ival = atoi( tmp[1] );
+ control->cm_solver_restart = ival;
+ }
+ else if ( strcmp(tmp[0], "cm_solver_q_err") == 0 )
+ {
+ val = atof( tmp[1] );
+ control->cm_solver_q_err = val;
+ }
+ else if ( strcmp(tmp[0], "cm_domain_sparsity") == 0 )
+ {
+ val = atof( tmp[1] );
+ control->cm_domain_sparsity = val;
+ if ( val < 1.0 )
+ {
+ control->cm_domain_sparsify_enabled = TRUE;
+ }
+ }
+ else if ( strncmp(tmp[0], "cm_init_guess_extrap1", MAX_LINE) == 0 )
+ {
+ ival = atoi( tmp[1] );
+ control->cm_init_guess_extrap1 = ival;
+ }
+ 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 )
+ {
+ ival = atoi( tmp[1] );
+ control->cm_solver_pre_comp_type = ival;
+ }
+ else if ( strcmp(tmp[0], "cm_solver_pre_comp_refactor") == 0 )
+ {
+ ival = atoi( tmp[1] );
+ control->cm_solver_pre_comp_refactor = ival;
+ }
+ else if ( strcmp(tmp[0], "cm_solver_pre_comp_droptol") == 0 )
+ {
+ val = atof( tmp[1] );
+ control->cm_solver_pre_comp_droptol = val;
+ }
+ else if ( strcmp(tmp[0], "cm_solver_pre_comp_sweeps") == 0 )
+ {
+ ival = atoi( tmp[1] );
+ control->cm_solver_pre_comp_sweeps = ival;
+ }
+ 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 )
+ {
+ ival = atoi( tmp[1] );
+ control->cm_solver_pre_app_type = ival;
+ }
+ else if ( strcmp(tmp[0], "cm_solver_pre_app_jacobi_iters") == 0 )
+ {
+ ival = atoi( tmp[1] );
+ control->cm_solver_pre_app_jacobi_iters = ival;
+ }
+ else if ( strcmp(tmp[0], "temp_init") == 0 )
+ {
+ val = atof(tmp[1]);
+ control->T_init = val;
- if ( control->T_final < 0.1 )
- control->T_final = 0.1;
- }
- else if ( strcmp(tmp[0], "t_mass") == 0 )
- {
- val = atof(tmp[1]);
- control->Tau_T = val * 1.e-3; // convert t_mass from fs to ps
- }
- else if ( strcmp(tmp[0], "t_mode") == 0 )
- {
- ival = atoi(tmp[1]);
- control->T_mode = ival;
- }
- else if ( strcmp(tmp[0], "t_rate") == 0 )
- {
- val = atof(tmp[1]);
- control->T_rate = val;
- }
- else if ( strcmp(tmp[0], "t_freq") == 0 )
- {
- val = atof(tmp[1]);
- control->T_freq = val;
- }
- else if ( strcmp(tmp[0], "pressure") == 0 )
- {
- if ( control->ensemble == iNPT )
+ if ( control->T_init < 0.1 )
+ {
+ control->T_init = 0.1;
+ }
+ }
+ else if ( strcmp(tmp[0], "temp_final") == 0 )
{
- control->P[0] = control->P[1] = control->P[2] = atof(tmp[1]);
+ val = atof(tmp[1]);
+ control->T_final = val;
+
+ if ( control->T_final < 0.1 )
+ {
+ control->T_final = 0.1;
+ }
}
- else if ( control->ensemble == sNPT )
+ else if ( strcmp(tmp[0], "t_mass") == 0 )
{
- control->P[0] = atof(tmp[1]);
- control->P[1] = atof(tmp[2]);
- control->P[2] = atof(tmp[3]);
+ val = atof(tmp[1]);
+ control->Tau_T = val * 1.e-3; // convert t_mass from fs to ps
}
- }
- else if ( strcmp(tmp[0], "p_mass") == 0 )
- {
- // convert p_mass from fs to ps
- if ( control->ensemble == iNPT )
+ else if ( strcmp(tmp[0], "t_mode") == 0 )
{
- control->Tau_P[0] = control->Tau_P[1] = control->Tau_P[2] =
- atof(tmp[1]) * 1.e-3;
+ ival = atoi(tmp[1]);
+ control->T_mode = ival;
}
- else if ( control->ensemble == sNPT )
+ else if ( strcmp(tmp[0], "t_rate") == 0 )
{
- control->Tau_P[0] = atof(tmp[1]) * 1.e-3;
- control->Tau_P[1] = atof(tmp[2]) * 1.e-3;
- control->Tau_P[2] = atof(tmp[3]) * 1.e-3;
+ val = atof(tmp[1]);
+ control->T_rate = val;
+ }
+ else if ( strcmp(tmp[0], "t_freq") == 0 )
+ {
+ val = atof(tmp[1]);
+ control->T_freq = val;
+ }
+ else if ( strcmp(tmp[0], "pressure") == 0 )
+ {
+ if ( control->ensemble == iNPT )
+ {
+ control->P[0] = control->P[1] = control->P[2] = atof(tmp[1]);
+ }
+ else if ( control->ensemble == sNPT )
+ {
+ control->P[0] = atof(tmp[1]);
+ control->P[1] = atof(tmp[2]);
+ control->P[2] = atof(tmp[3]);
+ }
+ }
+ else if ( strcmp(tmp[0], "p_mass") == 0 )
+ {
+ // convert p_mass from fs to ps
+ if ( control->ensemble == iNPT )
+ {
+ control->Tau_P[0] = control->Tau_P[1] = control->Tau_P[2] =
+ atof(tmp[1]) * 1.e-3;
+ }
+ else if ( control->ensemble == sNPT )
+ {
+ control->Tau_P[0] = atof(tmp[1]) * 1.e-3;
+ control->Tau_P[1] = atof(tmp[2]) * 1.e-3;
+ control->Tau_P[2] = atof(tmp[3]) * 1.e-3;
+ }
+ }
+ else if ( strcmp(tmp[0], "pt_mass") == 0 )
+ {
+ val = atof(tmp[1]);
+ control->Tau_PT[0] = control->Tau_PT[1] = control->Tau_PT[2] =
+ val * 1.e-3; // convert pt_mass from fs to ps
+ }
+ else if ( strcmp(tmp[0], "compress") == 0 )
+ {
+ val = atof(tmp[1]);
+ control->compressibility = val;
+ }
+ else if ( strcmp(tmp[0], "press_mode") == 0 )
+ {
+ ival = atoi(tmp[1]);
+ control->press_mode = ival;
+ }
+ else if ( strcmp(tmp[0], "geo_format") == 0 )
+ {
+ ival = atoi( tmp[1] );
+ control->geo_format = ival;
+ }
+ else if ( strcmp(tmp[0], "write_freq") == 0 )
+ {
+ ival = atoi(tmp[1]);
+ out_control->write_steps = ival;
+ }
+ else if ( strcmp(tmp[0], "traj_compress") == 0 )
+ {
+ ival = atoi(tmp[1]);
+ out_control->traj_compress = ival;
+ }
+ else if ( strcmp(tmp[0], "traj_method") == 0 )
+ {
+ ival = atoi(tmp[1]);
+ out_control->traj_method = ival;
+ }
+ else if ( strcmp(tmp[0], "traj_title") == 0 )
+ {
+ strcpy( out_control->traj_title, tmp[1] );
+ }
+ else if ( strcmp(tmp[0], "atom_info") == 0 )
+ {
+ ival = atoi(tmp[1]);
+ out_control->atom_info += ival * 4;
+ }
+ else if ( strcmp(tmp[0], "atom_velocities") == 0 )
+ {
+ ival = atoi(tmp[1]);
+ out_control->atom_info += ival * 2;
+ }
+ else if ( strcmp(tmp[0], "atom_forces") == 0 )
+ {
+ ival = atoi(tmp[1]);
+ out_control->atom_info += ival * 1;
+ }
+ else if ( strcmp(tmp[0], "bond_info") == 0 )
+ {
+ ival = atoi(tmp[1]);
+ out_control->bond_info = ival;
+ }
+ else if ( strcmp(tmp[0], "angle_info") == 0 )
+ {
+ ival = atoi(tmp[1]);
+ out_control->angle_info = ival;
+ }
+ else if ( strcmp(tmp[0], "molecular_analysis") == 0 )
+ {
+ ival = atoi(tmp[1]);
+ control->molecular_analysis = ival;
+ }
+ else if ( strcmp(tmp[0], "ignore") == 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 )
+ {
+ ival = atoi(tmp[1]);
+ control->dipole_anal = ival;
+ }
+ else if ( strcmp(tmp[0], "freq_dipole_anal") == 0 )
+ {
+ ival = atoi(tmp[1]);
+ control->freq_dipole_anal = ival;
+ }
+ else if ( strcmp(tmp[0], "diffusion_coef") == 0 )
+ {
+ ival = atoi(tmp[1]);
+ control->diffusion_coef = ival;
+ }
+ else if ( strcmp(tmp[0], "freq_diffusion_coef") == 0 )
+ {
+ ival = atoi(tmp[1]);
+ control->freq_diffusion_coef = ival;
+ }
+ else if ( strcmp(tmp[0], "restrict_type") == 0 )
+ {
+ ival = atoi(tmp[1]);
+ control->restrict_type = ival;
+ }
+ else
+ {
+ fprintf( stderr, "[ERROR] unknown control file parameter (%s)\n", tmp[0] );
+ MPI_Abort( MPI_COMM_WORLD, INVALID_INPUT );
}
- }
- else if ( strcmp(tmp[0], "pt_mass") == 0 )
- {
- val = atof(tmp[1]);
- control->Tau_PT[0] = control->Tau_PT[1] = control->Tau_PT[2] =
- val * 1.e-3; // convert pt_mass from fs to ps
- }
- else if ( strcmp(tmp[0], "compress") == 0 )
- {
- val = atof(tmp[1]);
- control->compressibility = val;
- }
- else if ( strcmp(tmp[0], "press_mode") == 0 )
- {
- ival = atoi(tmp[1]);
- control->press_mode = ival;
- }
- else if ( strcmp(tmp[0], "geo_format") == 0 )
- {
- ival = atoi( tmp[1] );
- control->geo_format = ival;
- }
- else if ( strcmp(tmp[0], "write_freq") == 0 )
- {
- ival = atoi(tmp[1]);
- out_control->write_steps = ival;
- }
- else if ( strcmp(tmp[0], "traj_compress") == 0 )
- {
- ival = atoi(tmp[1]);
- out_control->traj_compress = ival;
- }
- else if ( strcmp(tmp[0], "traj_method") == 0 )
- {
- ival = atoi(tmp[1]);
- out_control->traj_method = ival;
- }
- else if ( strcmp(tmp[0], "traj_title") == 0 )
- {
- strcpy( out_control->traj_title, tmp[1] );
- }
- else if ( strcmp(tmp[0], "atom_info") == 0 )
- {
- ival = atoi(tmp[1]);
- out_control->atom_info += ival * 4;
- }
- else if ( strcmp(tmp[0], "atom_velocities") == 0 )
- {
- ival = atoi(tmp[1]);
- out_control->atom_info += ival * 2;
- }
- else if ( strcmp(tmp[0], "atom_forces") == 0 )
- {
- ival = atoi(tmp[1]);
- out_control->atom_info += ival * 1;
- }
- else if ( strcmp(tmp[0], "bond_info") == 0 )
- {
- ival = atoi(tmp[1]);
- out_control->bond_info = ival;
- }
- else if ( strcmp(tmp[0], "angle_info") == 0 )
- {
- ival = atoi(tmp[1]);
- out_control->angle_info = ival;
- }
- else if ( strcmp(tmp[0], "molecular_analysis") == 0 )
- {
- ival = atoi(tmp[1]);
- control->molecular_analysis = ival;
- }
- else if ( strcmp(tmp[0], "ignore") == 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 )
- {
- ival = atoi(tmp[1]);
- control->dipole_anal = ival;
- }
- else if ( strcmp(tmp[0], "freq_dipole_anal") == 0 )
- {
- ival = atoi(tmp[1]);
- control->freq_dipole_anal = ival;
- }
- else if ( strcmp(tmp[0], "diffusion_coef") == 0 )
- {
- ival = atoi(tmp[1]);
- control->diffusion_coef = ival;
- }
- else if ( strcmp(tmp[0], "freq_diffusion_coef") == 0 )
- {
- ival = atoi(tmp[1]);
- control->freq_diffusion_coef = ival;
- }
- else if ( strcmp(tmp[0], "restrict_type") == 0 )
- {
- ival = atoi(tmp[1]);
- control->restrict_type = ival;
- }
- else
- {
- fprintf( stderr, "WARNING: unknown parameter %s\n", tmp[0] );
- MPI_Abort( MPI_COMM_WORLD, 15 );
}
}
+ if ( ferror( fp ) )
+ {
+ fprintf( stderr, "[ERROR] parsing control file failed (I/O error). TERMINATING...\n" );
+ MPI_Abort( MPI_COMM_WORLD, INVALID_INPUT );
+ }
+
/* determine target T */
if ( control->T_mode == 0 )
+ {
control->T = control->T_final;
- else control->T = control->T_init;
+ }
+ else
+ {
+ control->T = control->T_init;
+ }
/* free memory allocations at the top */
for ( i = 0; i < MAX_TOKENS; i++ )
- sfree( tmp[i], "tmp[i]" );
- sfree( tmp, "tmp" );
- sfree( s, "s" );
+ {
+ sfree( tmp[i], "Read_Control_File::tmp[i]" );
+ }
+ sfree( tmp, "Read_Control_File::tmp" );
+ sfree( s, "Read_Control_File::s" );
// fprintf( stderr,"%d %d %10.5f %d %10.5f %10.5f\n",
// control->ensemble, control->nsteps, control->dt,
diff --git a/PuReMD/src/forces.c b/PuReMD/src/forces.c
index d438a0eb..c406417b 100644
--- a/PuReMD/src/forces.c
+++ b/PuReMD/src/forces.c
@@ -926,7 +926,7 @@ void Compute_Forces( reax_system *system, control_params *control,
comm = mpi_data->world;
/********* init forces ************/
#if defined(PURE_REAX)
- if ( control->qeq_freq && (data->step - data->prev_steps) % control->qeq_freq == 0 )
+ if ( control->charge_freq && (data->step - data->prev_steps) % control->charge_freq == 0 )
qeq_flag = 1;
else qeq_flag = 0;
#elif defined(LAMMPS_REAX)
@@ -970,7 +970,7 @@ void Compute_Forces( reax_system *system, control_params *control,
#if defined(LOG_PERFORMANCE)
//MPI_Barrier( mpi_data->world );
if ( system->my_rank == MASTER_NODE )
- Update_Timing_Info( &t_start, &(data->timing.qEq) );
+ Update_Timing_Info( &t_start, &data->timing.cm );
#endif
#if defined(DEBUG_FOCUS)
fprintf(stderr, "p%d @ step%d: qeq completed\n", system->my_rank, data->step);
diff --git a/PuReMD/src/io_tools.c b/PuReMD/src/io_tools.c
index 560bb5a1..7f3c4c02 100644
--- a/PuReMD/src/io_tools.c
+++ b/PuReMD/src/io_tools.c
@@ -1274,8 +1274,8 @@ void Output_Results( reax_system *system, control_params *control,
data->timing.init_forces * denom,
data->timing.bonded * denom,
data->timing.nonb * denom,
- data->timing.qEq * denom,
- (int)((data->timing.s_matvecs + data->timing.t_matvecs)*denom) );
+ data->timing.cm * denom,
+ (int)(data->timing.cm_solver_iters * denom) );
Reset_Timing( &(data->timing) );
fflush( out_control->log );
diff --git a/PuReMD/src/qEq.c b/PuReMD/src/qEq.c
index c8a348c7..15cc0249 100644
--- a/PuReMD/src/qEq.c
+++ b/PuReMD/src/qEq.c
@@ -382,18 +382,18 @@ void QEq( reax_system *system, control_params *control, simulation_data *data,
#endif
//s_matvecs = dual_CG(system, workspace, workspace->H, workspace->b,
- // control->q_err, workspace->x, mpi_data, out_control->log);
+ // control->cm_solver_q_err, workspace->x, mpi_data, out_control->log);
//t_matvecs = 0;
for ( j = 0; j < system->n; ++j )
workspace->s[j] = workspace->x[j][0];
s_matvecs = CG(system, workspace, workspace->H, workspace->b_s,//newQEq sCG
- control->q_err, workspace->s, mpi_data, out_control->log );
+ control->cm_solver_q_err, workspace->s, mpi_data, out_control->log );
for ( j = 0; j < system->n; ++j )
workspace->x[j][0] = workspace->s[j];
//s_matvecs = PCG( system, workspace, workspace->H, workspace->b_s,
- // control->q_err, workspace->L, workspace->U, workspace->s,
+ // control->cm_solver_q_err, workspace->L, workspace->U, workspace->s,
// mpi_data, out_control->log );
#if defined(DEBUG)
fprintf( stderr, "p%d: first CG completed\n", system->my_rank );
@@ -402,12 +402,12 @@ void QEq( reax_system *system, control_params *control, simulation_data *data,
for ( j = 0; j < system->n; ++j )
workspace->t[j] = workspace->x[j][1];
t_matvecs = CG(system, workspace, workspace->H, workspace->b_t,//newQEq sCG
- control->q_err, workspace->t, mpi_data, out_control->log );
+ control->cm_solver_q_err, workspace->t, mpi_data, out_control->log );
for ( j = 0; j < system->n; ++j )
workspace->x[j][1] = workspace->t[j];
//t_matvecs = PCG( system, workspace, workspace->H, workspace->b_t,
- // control->q_err, workspace->L, workspace->U, workspace->t,
+ // control->cm_solver_q_err, workspace->L, workspace->U, workspace->t,
// mpi_data, out_control->log );
#if defined(DEBUG)
fprintf( stderr, "p%d: second CG completed\n", system->my_rank );
@@ -422,8 +422,7 @@ void QEq( reax_system *system, control_params *control, simulation_data *data,
#if defined(LOG_PERFORMANCE)
if ( system->my_rank == MASTER_NODE )
{
- data->timing.s_matvecs += s_matvecs;
- data->timing.t_matvecs += t_matvecs;
+ data->timing.cm_solver_iters += s_matvecs + t_matvecs;
}
#endif
}
diff --git a/PuReMD/src/reax_types.h b/PuReMD/src/reax_types.h
index 152c0c54..dfaa9c37 100644
--- a/PuReMD/src/reax_types.h
+++ b/PuReMD/src/reax_types.h
@@ -22,6 +22,11 @@
#ifndef __REAX_TYPES_H_
#define __REAX_TYPES_H_
+#if (defined(HAVE_CONFIG_H) && !defined(__CONFIG_H_))
+ #define __CONFIG_H_
+ #include "config.h"
+#endif
+
#include <ctype.h>
#include <math.h>
#include <mpi.h>
@@ -36,7 +41,6 @@
#define PURE_REAX
//#define LAMMPS_REAX
-
//#define DEBUG
//#define DEBUG_FOCUS
//#define TEST_ENERGY
@@ -47,12 +51,13 @@
//#define OLD_BOUNDARIES
//#define MIDPOINT_BOUNDARIES
-#define REAX_MAX_STR 1024
-#define REAX_MAX_NBRS 6
-#define REAX_MAX_3BODY_PARAM 5
-#define REAX_MAX_4BODY_PARAM 5
-#define REAX_MAX_ATOM_TYPES 25
-#define REAX_MAX_MOLECULE_SIZE 20
+#define ZERO (0.000000000000000e+00)
+#define REAX_MAX_STR (1024)
+#define REAX_MAX_NBRS (6)
+#define REAX_MAX_3BODY_PARAM (5)
+#define REAX_MAX_4BODY_PARAM (5)
+#define REAX_MAX_ATOM_TYPES (25)
+#define REAX_MAX_MOLECULE_SIZE (20)
/* NaN IEEE 754 representation for C99 in math.h
* Note: function choice must match REAL typedef below */
@@ -63,6 +68,68 @@
#define NAN_REAL(a) (0)
#endif
+
+/* method for computing atomic charges */
+enum charge_method
+{
+ QEQ_CM = 0,
+ EE_CM = 1,
+ ACKS2_CM = 2,
+};
+
+/* linear solver type used in charge method */
+enum solver
+{
+ GMRES_S = 0,
+ GMRES_H_S = 1,
+ CG_S = 2,
+ SDM_S = 3,
+ BiCGStab_S = 4,
+};
+
+/* preconditioner computation type for charge method linear solver */
+enum pre_comp
+{
+ NONE_PC = 0,
+ JACOBI_PC = 1,
+ ICHOLT_PC = 2,
+ ILUT__PC = 3,
+ ILUTP_PC = 4,
+ FG_ILUT_PC = 5,
+ SAI_PC = 6,
+};
+
+/* preconditioner application type for ICHOL/ILU preconditioners,
+ * used for charge method linear solver */
+enum pre_app
+{
+ TRI_SOLVE_PA = 0,
+ TRI_SOLVE_LEVEL_SCHED_PA = 1,
+ TRI_SOLVE_GC_PA = 2,
+ JACOBI_ITER_PA = 3,
+};
+
+/* interaction list (reax_list) storage format */
+enum reax_list_format
+{
+ /* store half of interactions, when i < j (atoms i and j) */
+ HALF_LIST = 0,
+ /* store all interactions */
+ FULL_LIST = 1,
+};
+
+/* sparse matrix (sparse_matrix) storage format */
+enum sparse_matrix_format
+{
+ /* store upper half of nonzeros in a symmetric matrix (a_{ij}, i >= j) */
+ SYM_HALF_MATRIX = 0,
+ /* store all nonzeros in a symmetric matrix */
+ SYM_FULL_MATRIX = 1,
+ /* store all nonzeros in a matrix */
+ FULL_MATRIX = 2,
+};
+
+
/********************** TYPE DEFINITIONS ********************/
typedef int ivec[3];
typedef double real;
@@ -148,15 +215,8 @@ typedef struct
MPI_Datatype bond_view;
MPI_Datatype angle_line;
MPI_Datatype angle_view;
-
- //MPI_Request send_req1[REAX_MAX_NBRS];
- //MPI_Request send_req2[REAX_MAX_NBRS];
- //MPI_Status send_stat1[REAX_MAX_NBRS];
- //MPI_Status send_stat2[REAX_MAX_NBRS];
- //MPI_Status recv_stat1[REAX_MAX_NBRS];
- //MPI_Status recv_stat2[REAX_MAX_NBRS];
-
mpi_out_data out_buffers[REAX_MAX_NBRS];
+
void *in1_buffer;
void *in2_buffer;
} mpi_datatypes;
@@ -493,64 +553,161 @@ typedef struct
/* system control parameters */
typedef struct
{
+ /* simulation name, as supplied via control file */
char sim_name[REAX_MAX_STR];
- int nprocs;
+ /* number of MPI processors, as supplied via control file */
+ int nprocs;
+ /* MPI processors per each simulation dimension (cartesian topology),
+ * as supplied via control file */
ivec procs_by_dim;
- /* ensemble values:
- 0 : NVE
- 1 : bNVT (Berendsen)
- 2 : nhNVT (Nose-Hoover)
- 3 : sNPT (Parrinello-Rehman-Nose-Hoover) semiisotropic
- 4 : iNPT (Parrinello-Rehman-Nose-Hoover) isotropic
- 5 : NPT (Parrinello-Rehman-Nose-Hoover) Anisotropic*/
- int ensemble;
- int nsteps;
+ /* ensemble type for simulation, values:
+ * 0 : NVE
+ * 1 : bNVT (Berendsen)
+ * 2 : nhNVT (Nose-Hoover)
+ * 3 : sNPT (Parrinello-Rehman-Nose-Hoover) semiisotropic
+ * 4 : iNPT (Parrinello-Rehman-Nose-Hoover) isotropic
+ * 5 : NPT (Parrinello-Rehman-Nose-Hoover) Anisotropic */
+ int ensemble;
+ /* num. of simulation time steps */
+ int nsteps;
+ /* length of time step, in femtoseconds */
real dt;
- int geo_format;
- int restart;
-
- int restrict_bonds;
- int remove_CoM_vel;
- int random_vel;
- int reposition_atoms;
-
- int reneighbor;
+ /* format of geometry input file */
+ int geo_format;
+ /* format of restart file */
+ int restart;
+
+ /**/
+ int restrict_bonds;
+ /* flag to control if center of mass velocity is removed */
+ int remove_CoM_vel;
+ /* flag to control if atomic initial velocity is randomly assigned */
+ int random_vel;
+ /* flag to control how atom repositioning is performed, values:
+ * 0: fit to periodic box
+ * 1: put center of mass to box center
+ * 2: put center of mass to box origin */
+ int reposition_atoms;
+
+ /* flag to control the frequency (in terms of simulation time stesp)
+ * at which atom reneighboring is performed */
+ int reneighbor;
+ /* far neighbor (Verlet list) interaction cutoff, in Angstroms */
real vlist_cut;
+ /* bond interaction cutoff, in Angstroms */
real bond_cut;
- real nonb_cut, nonb_low;
+ /* non-bonded interaction cutoff, in Angstroms */
+ real nonb_cut;
+ /* ???, as supplied by force field parameters, in Angstroms */
+ real nonb_low;
+ /* hydrogen bond interaction cutoff, in Angstroms */
real hbond_cut;
+ /* ghost region cutoff (user-supplied via control file), in Angstroms */
real user_ghost_cut;
+ /* bond graph cutoff, as supplied by control file, in Angstroms */
real bg_cut;
+ /* bond order cutoff, as supplied by force field parameters, in Angstroms */
real bo_cut;
+ /* three body interaction cutoff, as supplied by control file, in Angstroms */
real thb_cut;
+ /* flag to control if force computations are tablulated */
int tabulate;
- int qeq_freq;
- real q_err;
- int refactor;
- real droptol;
-
- real T_init, T_final, T;
+ /* method for computing atomic charges */
+ unsigned int charge_method;
+ /* frequency (in terms of simulation time steps) at which to
+ * re-compute atomic charge distribution */
+ int charge_freq;
+ /* iterative linear solver type */
+ unsigned int cm_solver_type;
+ /* system net charge */
+ real cm_q_net;
+ /* max. iterations for linear solver */
+ unsigned int cm_solver_max_iters;
+ /* max. iterations before restarting in specific solvers, e.g., GMRES(k) */
+ unsigned int cm_solver_restart;
+ /* error tolerance of solution produced by charge distribution
+ * sparse iterative linear solver */
+ real cm_solver_q_err;
+ /* ratio used in computing sparser charge matrix,
+ * between 0.0 and 1.0 */
+ real cm_domain_sparsity;
+ /* TRUE if enabled, FALSE otherwise */
+ unsigned int cm_domain_sparsify_enabled;
+ /* order of spline extrapolation used for computing initial guess
+ * to linear solver */
+ unsigned int cm_init_guess_extrap1;
+ /* order of spline extrapolation used for computing initial guess
+ * to linear solver */
+ unsigned int cm_init_guess_extrap2;
+ /* preconditioner type for linear solver */
+ unsigned int cm_solver_pre_comp_type;
+ /* frequency (in terms of simulation time steps) at which to recompute
+ * incomplete factorizations */
+ unsigned int cm_solver_pre_comp_refactor;
+ /* drop tolerance of incomplete factorization schemes (ILUT, ICHOLT, etc.)
+ * used for preconditioning the iterative linear solver used in charge distribution */
+ real cm_solver_pre_comp_droptol;
+ /* num. of sweeps for computing preconditioner factors
+ * in fine-grained iterative methods (FG-ICHOL, FG-ILU) */
+ unsigned int cm_solver_pre_comp_sweeps;
+ /* relative num. of non-zeros to charge matrix used to
+ * compute the sparse approximate inverse preconditioner,
+ * between 0.0 and 1.0 */
+ real cm_solver_pre_comp_sai_thres;
+ /* preconditioner application type */
+ unsigned int cm_solver_pre_app_type;
+ /* num. of iterations used to apply preconditioner via
+ * Jacobi relaxation scheme (truncated Neumann series) */
+ unsigned int cm_solver_pre_app_jacobi_iters;
+
+ /* initial temperature of simulation, in Kelvin */
+ real T_init;
+ /* final temperature of simulation, in Kelvin */
+ real T_final;
+ /* current temperature of simulation, in Kelvin */
+ real T;
+ /**/
real Tau_T;
+ /**/
int T_mode;
- real T_rate, T_freq;
+ /**/
+ real T_rate;
+ /**/
+ real T_freq;
+ /**/
int virial;
- rvec P, Tau_P, Tau_PT;
- int press_mode;
+ /**/
+ rvec P;
+ /**/
+ rvec Tau_P;
+ /**/
+ rvec Tau_PT;
+ /**/
+ int press_mode;
+ /**/
real compressibility;
- int molecular_analysis;
- int num_ignored;
+ /**/
+ int molecular_analysis;
+ /**/
+ int num_ignored;
+ /**/
int ignore[REAX_MAX_ATOM_TYPES];
- int dipole_anal;
- int freq_dipole_anal;
- int diffusion_coef;
- int freq_diffusion_coef;
- int restrict_type;
+ /**/
+ int dipole_anal;
+ /**/
+ int freq_dipole_anal;
+ /**/
+ int diffusion_coef;
+ /**/
+ int freq_diffusion_coef;
+ /**/
+ int restrict_type;
} control_params;
@@ -596,19 +753,46 @@ typedef struct
typedef struct
{
+ /* start time of event */
real start;
+ /* end time of event */
real end;
+ /* total elapsed time of event */
real elapsed;
-
+ /* total simulation time */
real total;
+ /* communication time */
real comm;
+ /* neighbor list generation time */
real nbrs;
+ /* force initialization time */
real init_forces;
+ /* bonded force calculation time */
real bonded;
+ /* non-bonded force calculation time */
real nonb;
- real qEq;
- int s_matvecs;
- int t_matvecs;
+ /* atomic charge distribution calculation time */
+ real cm;
+ /**/
+ real cm_sort_mat_rows;
+ /**/
+ real cm_solver_comm;
+ /**/
+ real cm_solver_allreduce;
+ /**/
+ real cm_solver_pre_comp;
+ /**/
+ real cm_solver_pre_app; // update CG()
+ /* num. of steps in iterative linear solver for charge distribution */
+ int cm_solver_iters;
+ /**/
+ real cm_solver_spmv; // update CG()
+ /**/
+ real cm_solver_vector_ops; // update CG()
+ /**/
+ real cm_solver_orthog;
+ /**/
+ real cm_solver_tri_solve;
} reax_timing;
@@ -763,6 +947,8 @@ typedef struct
typedef struct
{
+ /* matrix storage format */
+ int format;
int cap, n, m;
int *start, *end;
sparse_matrix_entry *entries;
@@ -894,6 +1080,9 @@ typedef struct
int type;
// list_type select;
+ /* list storage format (half or full) */
+ int format;
+
void *v;
three_body_interaction_data *three_body_list;
bond_data *bond_list;
diff --git a/PuReMD/src/reset_tools.c b/PuReMD/src/reset_tools.c
index eb7c6465..1acbc002 100644
--- a/PuReMD/src/reset_tools.c
+++ b/PuReMD/src/reset_tools.c
@@ -107,14 +107,13 @@ void Reset_Simulation_Data( simulation_data* data, int virial )
void Reset_Timing( reax_timing *rt )
{
rt->total = Get_Time();
- rt->comm = 0;
- rt->nbrs = 0;
- rt->init_forces = 0;
- rt->bonded = 0;
- rt->nonb = 0;
- rt->qEq = 0;
- rt->s_matvecs = 0;
- rt->t_matvecs = 0;
+ rt->comm = 0.0;
+ rt->nbrs = 0.0;
+ rt->init_forces = 0.0;
+ rt->bonded = 0.0;
+ rt->nonb = 0.0;
+ rt->cm = 0.0;
+ rt->cm_solver_iters = 0;
}
#ifdef TEST_FORCES
diff --git a/PuReMD/src/tool_box.c b/PuReMD/src/tool_box.c
index a0892870..d898f219 100644
--- a/PuReMD/src/tool_box.c
+++ b/PuReMD/src/tool_box.c
@@ -486,3 +486,65 @@ void sfree( void *ptr, char *name )
free( ptr );
ptr = NULL;
}
+
+
+/* Safe wrapper around libc fopen
+ *
+ * fname: name of file to be opened
+ * mode: mode in which to open file
+ * msg: message to be printed in case of error
+ * */
+FILE * sfopen( const char * fname, const char * mode, const char * msg )
+{
+ FILE * ptr;
+
+ if ( fname == NULL )
+ {
+ fprintf( stderr, "[ERROR] trying to open file: NULL file name (%s). Terminating...\n",
+ msg );
+ exit( INVALID_INPUT );
+ }
+ if ( mode == NULL )
+ {
+ fprintf( stderr, "[ERROR] trying to open file: NULL mode (%s). Terminating...\n",
+ msg );
+ exit( INVALID_INPUT );
+ }
+
+ ptr = fopen( fname, mode );
+
+ if ( ptr == NULL )
+ {
+ fprintf( stderr, "[ERROR] failed to open file %s with mode %s (%s)\n",
+ fname, mode, msg );
+ 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 (%s). Returning...\n", msg );
+ return;
+ }
+
+ ret = fclose( fp );
+
+ if ( ret != 0 )
+ {
+ fprintf( stderr, "[ERROR] error detected when closing file (%s). Terminating...\n", msg );
+ exit( INVALID_INPUT );
+ }
+}
diff --git a/PuReMD/src/tool_box.h b/PuReMD/src/tool_box.h
index bf25c64f..22cbc195 100644
--- a/PuReMD/src/tool_box.h
+++ b/PuReMD/src/tool_box.h
@@ -66,4 +66,8 @@ void *smalloc( long, char*, MPI_Comm );
void *scalloc( int, int, char*, MPI_Comm );
void sfree( void*, char* );
+FILE * sfopen( const char *, const char *, const char * );
+
+void sfclose( FILE *, const char * );
+
#endif
diff --git a/PuReMD/src/traj.c b/PuReMD/src/traj.c
index 40df8bbd..d03a8289 100644
--- a/PuReMD/src/traj.c
+++ b/PuReMD/src/traj.c
@@ -261,7 +261,7 @@ int Write_Header( reax_system *system, control_params *control,
control->thb_cut );
strncat( out_control->buffer, out_control->line, HEADER_LINE_LEN + 1 );
- sprintf( out_control->line, SCI_LINE, "QEq_tolerance:", control->q_err );
+ sprintf( out_control->line, SCI_LINE, "QEq_tolerance:", control->cm_solver_q_err );
strncat( out_control->buffer, out_control->line, HEADER_LINE_LEN + 1 );
/* temperature controls */
--
GitLab