diff --git a/PuReMD/src/control.c b/PuReMD/src/control.c
index f0f6ea0cb9780311237da174b07a743f6b29ca9b..e0e0e86d2999e7bdb1350b0c59921584b2655533 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,51 +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 = 0;
+ 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->charge_method = QEQ_CM;
+ control->charge_freq = 1;
+ control->cm_q_net = 0.0;
+ control->cm_solver_type = GMRES_S;
control->cm_solver_max_iters = 100;
- control->q_err = 1e-6;
- control->sai_thres = 0.1;
- control->refactor = 100;
- control->droptol = 1e-2;;
+ control->cm_solver_restart = 50;
+ 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.;
@@ -110,355 +118,435 @@ 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 ( 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;
- control->refactor = 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], "cm_solver_max_iters") == 0 )
- {
- ival = atoi( tmp[1] );
- control->cm_solver_max_iters = ival;
- }
- else if ( strcmp(tmp[0], "q_err") == 0 )
- {
- val = atof( tmp[1] );
- control->q_err = val;
- }
- else if ( strcmp(tmp[0], "cm_solver_pre_comp_type") == 0 )
- {
- val = atoi( tmp[1] );
- control->cm_solver_pre_comp_type = val;
- }
- else if ( strcmp(tmp[0], "sai_thres") == 0 )
- {
- val = atof( tmp[1] );
- control->sai_thres = 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 bbbf110adf7a3bd5be7d24c236f968612a8aa3bd..49e747a4c92042a1b2eba168a9c05417450b9332 100644
--- a/PuReMD/src/forces.c
+++ b/PuReMD/src/forces.c
@@ -978,7 +978,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)
diff --git a/PuReMD/src/linear_solvers.c b/PuReMD/src/linear_solvers.c
index bcb0432798fc5d4a3aff91c0848d025985eea50a..34b3f8e9b38157abd5797d9b33beace90a6df492 100644
--- a/PuReMD/src/linear_solvers.c
+++ b/PuReMD/src/linear_solvers.c
@@ -1227,7 +1227,7 @@ int CG( reax_system *system, control_params *control, simulation_data *data,
t_pa += Get_Timing_Info( t_start );
}
- else if ( control->cm_solver_pre_comp_type == DIAG_PC)
+ else if ( control->cm_solver_pre_comp_type == JACOBI_PC)
{
t_start = Get_Time( );
for ( j = 0; j < system->n; ++j )
@@ -1277,7 +1277,7 @@ int CG( reax_system *system, control_params *control, simulation_data *data,
t_pa += Get_Timing_Info( t_start );
}
- else if ( control->cm_solver_pre_comp_type == DIAG_PC)
+ else if ( control->cm_solver_pre_comp_type == JACOBI_PC)
{
t_start = Get_Time( );
for ( j = 0; j < system->n; ++j )
diff --git a/PuReMD/src/qEq.c b/PuReMD/src/qEq.c
index 969171981070ac646122b9e42844159c7f70cc02..3ec05786996bae3f0a25e0301c344c5111670d46 100644
--- a/PuReMD/src/qEq.c
+++ b/PuReMD/src/qEq.c
@@ -237,12 +237,12 @@ void Init_MatVec( reax_system *system, simulation_data *data,
int i; //, fillin;
reax_atom *atom;
- /*if( (data->step - data->prev_steps) % control->refactor == 0 ||
+ /*if( (data->step - data->prev_steps) % control->cm_solver_pre_comp_refactor == 0 ||
workspace->L == NULL ) {
//Print_Linear_System( system, control, workspace, data->step );
Sort_Matrix_Rows( workspace->H );
fprintf( stderr, "H matrix sorted\n" );
- Calculate_Droptol( workspace->H, workspace->droptol, control->droptol );
+ Calculate_Droptol( workspace->H, workspace->droptol, control->cm_solver_pre_comp_droptol );
fprintf( stderr, "drop tolerances calculated\n" );
if( workspace->L == NULL ) {
fillin = Estimate_LU_Fill( workspace->H, workspace->droptol );
@@ -376,7 +376,7 @@ static void Setup_Preconditioner_QEq( reax_system *system, control_params *contr
t_sort = Get_Timing_Info( time );
t_pc = setup_sparse_approx_inverse( system, data, workspace, mpi_data, workspace->H, &workspace->H_spar_patt,
- control->nprocs, control->sai_thres );
+ control->nprocs, control->cm_solver_pre_comp_sai_thres );
MPI_Reduce(&t_sort, &total_sort, 1, MPI_DOUBLE, MPI_SUM, MASTER_NODE, mpi_data->world);
@@ -424,7 +424,7 @@ void QEq( reax_system *system, control_params *control, simulation_data *data,
if( control->cm_solver_pre_comp_type == SAI_PC )
{
- if( control->refactor > 0 && ((data->step - data->prev_steps) % control->refactor == 0))
+ if( control->cm_solver_pre_comp_refactor > 0 && ((data->step - data->prev_steps) % control->cm_solver_pre_comp_refactor == 0))
{
Setup_Preconditioner_QEq( system, control, data, workspace, mpi_data );
@@ -436,7 +436,7 @@ void QEq( reax_system *system, control_params *control, simulation_data *data,
for ( j = 0; j < system->n; ++j )
workspace->s[j] = workspace->x[j][0];
iters = CG(system, control, data, workspace, workspace->H, workspace->b_s,
- control->q_err, workspace->s, mpi_data, out_control->log , control->nprocs );
+ control->cm_solver_q_err, workspace->s, mpi_data, out_control->log , control->nprocs );
for ( j = 0; j < system->n; ++j )
workspace->x[j][0] = workspace->s[j];
@@ -447,7 +447,7 @@ 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];
iters += CG(system, control, data, workspace, workspace->H, workspace->b_t,//newQEq sCG
- control->q_err, workspace->t, mpi_data, out_control->log, control->nprocs );
+ control->cm_solver_q_err, workspace->t, mpi_data, out_control->log, control->nprocs );
for ( j = 0; j < system->n; ++j )
workspace->x[j][1] = workspace->t[j];
diff --git a/PuReMD/src/reax_types.h b/PuReMD/src/reax_types.h
index 4335d66cb3e14ca6c0e4d6f2c566f8bb445795bf..f9933ffe9b3181c2555dd6708af2614a1608e506 100644
--- a/PuReMD/src/reax_types.h
+++ b/PuReMD/src/reax_types.h
@@ -70,18 +70,46 @@
#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,
- DIAG_PC = 1,
+ JACOBI_PC = 1,
ICHOLT_PC = 2,
- ILU_PAR_PC = 3,
- ILUT_PAR_PC = 4,
- ILU_SUPERLU_MT_PC = 5,
+ 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,
+};
+
/********************** TYPE DEFINITIONS ********************/
typedef int ivec[3];
@@ -513,67 +541,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;
+ /* format of geometry input file */
+ int geo_format;
+ /* format of restart file */
+ int restart;
- int reneighbor;
+ /**/
+ 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;
- int cm_solver_max_iters;
- real q_err;
- int cm_solver_pre_comp_type;
- real sai_thres;
- 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;
diff --git a/PuReMD/src/tool_box.c b/PuReMD/src/tool_box.c
index a089287001279e198fccbb9bbfd0aa660863c2d3..d898f2190d72de1424eea2730ed506e9e04ec610 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 bf25c64fc658e46d4762da484bd1fa7cbfac605c..22cbc1954e03cc7892a7c341e392dbb3603f51d7 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 40df8bbd1860eff00141ed97013e596430daf69c..d03a8289f22639477db5c9c7b3af919d37adba9e 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 */