/*----------------------------------------------------------------------
PuReMD - Purdue ReaxFF Molecular Dynamics Program
Copyright (2010) Purdue University
Hasan Metin Aktulga, haktulga@cs.purdue.edu
Joseph Fogarty, jcfogart@mail.usf.edu
Sagar Pandit, pandit@usf.edu
Ananth Y Grama, ayg@cs.purdue.edu
This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License as
published by the Free Software Foundation; either version 2 of
the License, or (at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
See the GNU General Public License for more details:
<http://www.gnu.org/licenses/>.
----------------------------------------------------------------------*/
#if (defined(HAVE_CONFIG_H) && !defined(__CONFIG_H_))
#define __CONFIG_H_
#include "../../common/include/config.h"
#endif
#if defined(PURE_REAX)
#include "control.h"
#include "tool_box.h"
#elif defined(LAMMPS_REAX)
#include "reax_control.h"
#include "reax_tool_box.h"
#endif
void Read_Control_File( const char * const control_file, control_params * const control,
output_controls * const out_control )
{
FILE *fp;
char *s, **tmp;
int c, i, ival;
real val;
fp = sfopen( control_file, "r", "Read_Control_File::fp" );
/* assign default values */
strncpy( control->sim_name, "default.sim", sizeof(control->sim_name) - 1 );
control->sim_name[sizeof(control->sim_name) - 1] = '\0';
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->gpus_per_node = 1;
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;
out_control->energy_update_freq = 0;
control->reneighbor = 1;
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->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 = 1;
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->polarization_energy_enabled = TRUE;
control->T_init = 0.0;
control->T_final = 300.;
control->Tau_T = 500.0;
control->T_mode = 0;
control->T_rate = 1.0;
control->T_freq = 1.0;
control->P[0] = 0.000101325;
control->P[1] = 0.000101325;
control->P[2] = 0.000101325;
control->Tau_P[0] = 500.0;
control->Tau_P[1] = 500.0;
control->Tau_P[2] = 500.0;
control->Tau_PT[0] = 500.0;
control->Tau_PT[1] = 500.0;
control->Tau_PT[2] = 500.0;
control->compressibility = 1.0;
control->press_mode = 0;
out_control->write_steps = 0;
out_control->traj_compress = 0;
out_control->traj_method = REG_TRAJ;
strncpy( out_control->traj_title, "default_title", sizeof(out_control->traj_title) - 1 );
out_control->traj_title[sizeof(out_control->traj_title) - 1] = '\0';
out_control->atom_info = 0;
out_control->bond_info = 0;
out_control->angle_info = 0;
control->molecular_analysis = 0;
control->dipole_anal = 0;
control->freq_dipole_anal = 0;
control->diffusion_coef = 0;
control->freq_diffusion_coef = 0;
control->restrict_type = 0;
/* memory allocations */
s = (char*) smalloc( sizeof(char) * MAX_LINE, "Read_Control_File::s" );
tmp = (char**) smalloc( sizeof(char*) * MAX_TOKENS, "Read_Control_File::tmp" );
for ( i = 0; i < MAX_TOKENS; i++ )
{
tmp[i] = (char*) smalloc( sizeof(char) * MAX_LINE, "Read_Control_File::tmp[i]" );
}
/* read control parameters file */
while ( fgets( s, MAX_LINE, fp ) )
{
c = Tokenize( s, &tmp, MAX_LINE );
if ( c > 0 )
{
if ( strncmp(tmp[0], "simulation_name", MAX_LINE) == 0 )
{
strncpy( control->sim_name, tmp[1], sizeof(control->sim_name) - 1 );
control->sim_name[sizeof(control->sim_name) - 1] = '\0';
}
else if ( strncmp(tmp[0], "ensemble_type", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
control->ensemble = ival;
}
else if ( strncmp(tmp[0], "nsteps", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
control->nsteps = ival;
}
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 ( strncmp(tmp[0], "gpus_per_node", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
control->gpus_per_node = ival;
}
else if ( strncmp(tmp[0], "proc_by_dim", MAX_LINE) == 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 );
}
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;
control->nprocs = control->procs_by_dim[0] * control->procs_by_dim[1] *
control->procs_by_dim[2];
}
else if ( strncmp(tmp[0], "periodic_boundaries", MAX_LINE) == 0 )
{
// skip since not supported in distributed memory code
;
}
// else if( strncmp(tmp[0], "restart", MAX_LINE) == 0 )
// {
// ival = atoi(tmp[1]);
// control->restart = ival;
// }
// else if( strncmp(tmp[0], "restart_from", MAX_LINE) == 0 )
// {
// strncpy( control->restart_from, tmp[1], sizeof(control->restart_from) - 1 );
// control->restart_from[sizeof(control->restart_from) - 1] = '\0';
// }
else if ( strncmp(tmp[0], "random_vel", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
control->random_vel = ival;
}
else if ( strncmp(tmp[0], "restart_format", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
out_control->restart_format = ival;
}
else if ( strncmp(tmp[0], "restart_freq", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
out_control->restart_freq = ival;
}
else if ( strncmp(tmp[0], "reposition_atoms", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
control->reposition_atoms = ival;
}
else if ( strncmp(tmp[0], "restrict_bonds", MAX_LINE) == 0 )
{
ival = atoi( tmp[1] );
control->restrict_bonds = ival;
}
else if ( strncmp(tmp[0], "remove_CoM_vel", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
control->remove_CoM_vel = ival;
}
else if ( strncmp(tmp[0], "debug_level", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
out_control->debug_level = ival;
}
else if ( strncmp(tmp[0], "energy_update_freq", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
out_control->energy_update_freq = ival;
}
else if ( strncmp(tmp[0], "reneighbor", MAX_LINE) == 0 )
{
ival = atoi( tmp[1] );
control->reneighbor = ival;
}
else if ( strncmp(tmp[0], "vlist_buffer", MAX_LINE) == 0 )
{
val = atof(tmp[1]);
control->vlist_cut = val + control->nonb_cut;
}
else if ( strncmp(tmp[0], "nbrhood_cutoff", MAX_LINE) == 0 )
{
val = atof(tmp[1]);
control->bond_cut = val;
}
else if ( strncmp(tmp[0], "bond_graph_cutoff", MAX_LINE) == 0 )
{
val = atof(tmp[1]);
control->bg_cut = val;
}
else if ( strncmp(tmp[0], "thb_cutoff", MAX_LINE) == 0 )
{
val = atof(tmp[1]);
control->thb_cut = val;
}
else if ( strncmp(tmp[0], "hbond_cutoff", MAX_LINE) == 0 )
{
val = atof( tmp[1] );
control->hbond_cut = val;
}
else if ( strncmp(tmp[0], "ghost_cutoff", MAX_LINE) == 0 )
{
val = atof(tmp[1]);
control->user_ghost_cut = val;
}
else if ( strncmp(tmp[0], "tabulate_long_range", MAX_LINE) == 0 )
{
ival = atoi( tmp[1] );
control->tabulate = ival;
}
else if ( strncmp(tmp[0], "charge_method", MAX_LINE) == 0 )
{
ival = atoi( tmp[1] );
control->charge_method = ival;
}
else if ( strncmp(tmp[0], "charge_freq", MAX_LINE) == 0 )
{
ival = atoi( tmp[1] );
control->charge_freq = ival;
}
else if ( strncmp(tmp[0], "cm_q_net", MAX_LINE) == 0 )
{
val = atof( tmp[1] );
control->cm_q_net = val;
}
else if ( strncmp(tmp[0], "cm_solver_type", MAX_LINE) == 0 )
{
ival = atoi( tmp[1] );
control->cm_solver_type = ival;
}
else if ( strncmp(tmp[0], "cm_solver_max_iters", MAX_LINE) == 0 )
{
ival = atoi( tmp[1] );
control->cm_solver_max_iters = ival;
}
else if ( strncmp(tmp[0], "cm_solver_restart", MAX_LINE) == 0 )
{
ival = atoi( tmp[1] );
control->cm_solver_restart = ival;
}
else if ( strncmp(tmp[0], "cm_solver_q_err", MAX_LINE) == 0 )
{
val = atof( tmp[1] );
control->cm_solver_q_err = val;
}
else if ( strncmp(tmp[0], "cm_domain_sparsity", MAX_LINE) == 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 ( strncmp(tmp[0], "cm_solver_pre_comp_type", MAX_LINE) == 0 )
{
ival = atoi( tmp[1] );
control->cm_solver_pre_comp_type = ival;
}
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 ( strncmp(tmp[0], "cm_solver_pre_comp_droptol", MAX_LINE) == 0 )
{
val = atof( tmp[1] );
control->cm_solver_pre_comp_droptol = val;
}
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 ( 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 ( strncmp(tmp[0], "cm_solver_pre_app_type", MAX_LINE) == 0 )
{
ival = atoi( tmp[1] );
control->cm_solver_pre_app_type = ival;
}
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 ( strncmp(tmp[0], "include_polarization_energy", MAX_LINE) == 0 )
{
ival = atoi( tmp[1] );
control->polarization_energy_enabled = ival;
}
else if ( strncmp(tmp[0], "temp_init", MAX_LINE) == 0 )
{
val = atof(tmp[1]);
control->T_init = val;
if ( control->T_init < 0.001 )
{
control->T_init = 0.001;
}
}
else if ( strncmp(tmp[0], "temp_final", MAX_LINE) == 0 )
{
val = atof(tmp[1]);
control->T_final = val;
if ( control->T_final < 0.1 )
{
control->T_final = 0.1;
}
}
else if ( strncmp(tmp[0], "t_mass", MAX_LINE) == 0 )
{
val = atof(tmp[1]);
/* convert from fs to s */
control->Tau_T = val * 1.0e-15;
}
else if ( strncmp(tmp[0], "t_mode", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
control->T_mode = ival;
}
else if ( strncmp(tmp[0], "t_rate", MAX_LINE) == 0 )
{
val = atof(tmp[1]);
control->T_rate = val;
}
else if ( strncmp(tmp[0], "t_freq", MAX_LINE) == 0 )
{
val = atof(tmp[1]);
control->T_freq = val;
}
else if ( strncmp(tmp[0], "pressure", MAX_LINE) == 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 ( strncmp(tmp[0], "p_mass", MAX_LINE) == 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 ( strncmp(tmp[0], "pt_mass", MAX_LINE) == 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 ( strncmp(tmp[0], "compress", MAX_LINE) == 0 )
{
val = atof(tmp[1]);
control->compressibility = val;
}
else if ( strncmp(tmp[0], "press_mode", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
control->press_mode = ival;
}
else if ( strncmp(tmp[0], "geo_format", MAX_LINE) == 0 )
{
ival = atoi( tmp[1] );
control->geo_format = ival;
}
else if ( strncmp(tmp[0], "write_freq", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
out_control->write_steps = ival;
}
else if ( strncmp(tmp[0], "traj_compress", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
out_control->traj_compress = ival;
}
else if ( strncmp(tmp[0], "traj_format", MAX_LINE) == 0 )
{
// skip since not applicable to distributed memory code
;
}
else if ( strncmp(tmp[0], "traj_method", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
out_control->traj_method = ival;
}
else if ( strncmp(tmp[0], "traj_title", MAX_LINE) == 0 )
{
strncpy( out_control->traj_title, tmp[1], sizeof(out_control->traj_title) - 1 );
out_control->traj_title[sizeof(out_control->traj_title) - 1] = '\0';
}
else if ( strncmp(tmp[0], "atom_info", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
out_control->atom_info += ival * 4;
}
else if ( strncmp(tmp[0], "atom_velocities", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
out_control->atom_info += ival * 2;
}
else if ( strncmp(tmp[0], "atom_forces", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
out_control->atom_info += ival * 1;
}
else if ( strncmp(tmp[0], "bond_info", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
out_control->bond_info = ival;
}
else if ( strncmp(tmp[0], "angle_info", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
out_control->angle_info = ival;
}
else if ( strncmp(tmp[0], "test_forces", MAX_LINE) == 0 )
{
// skip since not supported in distributed memory code
;
}
else if ( strncmp(tmp[0], "molecular_analysis", MAX_LINE) == 0
|| strncmp(tmp[0], "molec_anal", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
control->molecular_analysis = ival;
}
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 ( strncmp(tmp[0], "freq_molec_anal", MAX_LINE) == 0 )
{
// skip since not supported in distributed memory code
;
}
else if ( strncmp(tmp[0], "dipole_anal", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
control->dipole_anal = ival;
}
else if ( strncmp(tmp[0], "freq_dipole_anal", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
control->freq_dipole_anal = ival;
}
else if ( strncmp(tmp[0], "diffusion_coef", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
control->diffusion_coef = ival;
}
else if ( strncmp(tmp[0], "freq_diffusion_coef", MAX_LINE) == 0 )
{
ival = atoi(tmp[1]);
control->freq_diffusion_coef = ival;
}
else if ( strncmp(tmp[0], "restrict_type", MAX_LINE) == 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 );
}
}
}
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;
}
/* free memory allocations at the top */
for ( i = 0; i < MAX_TOKENS; i++ )
{
sfree( tmp[i], "Read_Control_File::tmp[i]" );
}
sfree( tmp, "Read_Control_File::tmp" );
sfree( s, "Read_Control_File::s" );
sfclose( fp, "Read_Control_Field::fp" );
}