From b8d76adabce7f436edf6f30625698b3c5da6076d Mon Sep 17 00:00:00 2001
From: "Kurt A. O'Hearn" <ohearnku@msu.edu>
Date: Mon, 4 Jan 2021 14:11:33 -0500
Subject: [PATCH] sPuReMD: add QMMM library interface for external codes (e.g.,
Amber). Rework file I/O parsing to jive with multiple consecutive simulation
code changes. Improve file parser error checking (numeric types). Other code
clean-up.
---
sPuReMD/src/charges.c | 6 +-
sPuReMD/src/control.c | 843 +++++++++++++++++++--------------------
sPuReMD/src/control.h | 5 +-
sPuReMD/src/ffield.c | 213 +++++-----
sPuReMD/src/ffield.h | 3 +-
sPuReMD/src/forces.c | 2 +-
sPuReMD/src/geo_tools.c | 36 +-
sPuReMD/src/init_md.c | 30 +-
sPuReMD/src/io_tools.c | 26 +-
sPuReMD/src/lin_alg.c | 12 +-
sPuReMD/src/reax_types.h | 3 +-
sPuReMD/src/spuremd.c | 510 +++++++++++++++++++++--
sPuReMD/src/spuremd.h | 32 +-
sPuReMD/src/tool_box.c | 138 +++++++
sPuReMD/src/tool_box.h | 36 +-
sPuReMD/src/vector.h | 12 +-
tools/driver_qmmm.py | 524 ++++++++++++++++++++++++
17 files changed, 1777 insertions(+), 654 deletions(-)
create mode 100644 tools/driver_qmmm.py
diff --git a/sPuReMD/src/charges.c b/sPuReMD/src/charges.c
index a392fbfc..4695942e 100644
--- a/sPuReMD/src/charges.c
+++ b/sPuReMD/src/charges.c
@@ -588,7 +588,7 @@ static void Compute_Preconditioner_QEq( const reax_system * const system,
//if ( control->cm_solver_pre_comp_refactor == -1 )
//{
// data->timing.cm_last_pre_comp = data->timing.cm_solver_pre_comp;
- // data->timing.cm_total_loss = ZERO;
+ // data->timing.cm_total_loss = 0.0;
//}
#if defined(DEBUG)
@@ -710,7 +710,7 @@ static void Compute_Preconditioner_EE( const reax_system * const system,
//if ( control->cm_solver_pre_comp_refactor == -1 )
//{
// data->timing.cm_last_pre_comp = data->timing.cm_solver_pre_comp;
- // data->timing.cm_total_loss = ZERO;
+ // data->timing.cm_total_loss = 0.0;
//}
if ( control->cm_solver_pre_app_type == TRI_SOLVE_GC_PA )
@@ -851,7 +851,7 @@ static void Compute_Preconditioner_ACKS2( const reax_system * const system,
//if ( control->cm_solver_pre_comp_refactor == -1 )
//{
// data->timing.cm_last_pre_comp = data->timing.cm_solver_pre_comp;
- // data->timing.cm_total_loss = ZERO;
+ // data->timing.cm_total_loss = 0.0;
//}
if ( control->cm_solver_pre_app_type == TRI_SOLVE_GC_PA )
diff --git a/sPuReMD/src/control.c b/sPuReMD/src/control.c
index 28978fa2..8019737a 100644
--- a/sPuReMD/src/control.c
+++ b/sPuReMD/src/control.c
@@ -29,13 +29,416 @@
#include "tool_box.h"
-void Read_Control_File( FILE* fp, reax_system *system, control_params* control,
- output_controls *out_control )
+int Set_Control_Parameter( const char * const keyword,
+ const char ** const values, control_params * const control,
+ output_controls * const out_control )
{
- char *s, **tmp;
- int c, i, ival;
+ int i, ret;
real val;
+ ret = SUCCESS;
+
+ if ( strncmp(keyword, "simulation_name", MAX_LINE) == 0 )
+ {
+ strncpy( control->sim_name, values[0], sizeof(control->sim_name) - 1 );
+ control->sim_name[sizeof(control->sim_name) - 1] = '\0';
+ }
+ //else if( strncmp(keyword, "restart", MAX_LINE) == 0 ) {
+ // control->restart = sstrtol( values[0], __FILE__, __LINE__ );
+ //}
+ else if ( strncmp(keyword, "restart_format", MAX_LINE) == 0 )
+ {
+ out_control->restart_format = sstrtol( values[0], __FILE__, __LINE__ );
+ }
+ else if ( strncmp(keyword, "restart_freq", MAX_LINE) == 0 )
+ {
+ out_control->restart_freq = sstrtol( values[0], __FILE__, __LINE__ );
+ }
+ else if ( strncmp(keyword, "random_vel", MAX_LINE) == 0 )
+ {
+ control->random_vel = sstrtol( values[0], __FILE__, __LINE__ );
+ }
+ else if ( strncmp(keyword, "reposition_atoms", MAX_LINE) == 0 )
+ {
+ control->reposition_atoms = sstrtol( values[0], __FILE__, __LINE__ );
+ }
+ else if ( strncmp(keyword, "ensemble_type", MAX_LINE) == 0 )
+ {
+ control->ensemble = sstrtol( values[0], __FILE__, __LINE__ );
+ }
+ else if ( strncmp(keyword, "nsteps", MAX_LINE) == 0 )
+ {
+ control->nsteps = sstrtol( values[0], __FILE__, __LINE__ );
+ }
+ else if ( strncmp(keyword, "dt", MAX_LINE) == 0 )
+ {
+ val = sstrtod( values[0], __FILE__, __LINE__ );
+ /* convert from fs to ps */
+ control->dt = val * 1.0e-3;
+ }
+ else if ( strncmp(keyword, "gpus_per_node", MAX_LINE) == 0 )
+ {
+ // skip since not applicable to shared memory code
+ ;
+ }
+ else if ( strncmp(keyword, "proc_by_dim", MAX_LINE) == 0 )
+ {
+ // skip since not applicable to shared memory code
+ ;
+ }
+ else if ( strncmp(keyword, "periodic_boundaries", MAX_LINE) == 0 )
+ {
+ control->periodic_boundaries = sstrtol( values[0], __FILE__, __LINE__ );
+ }
+ else if ( strncmp(keyword, "geo_format", MAX_LINE) == 0 )
+ {
+ control->geo_format = sstrtol( values[0], __FILE__, __LINE__ );
+ }
+ else if ( strncmp(keyword, "restrict_bonds", MAX_LINE) == 0 )
+ {
+ control->restrict_bonds = sstrtol( values[0], __FILE__, __LINE__ );
+ }
+ else if ( strncmp(keyword, "tabulate_long_range", MAX_LINE) == 0 )
+ {
+ control->tabulate = sstrtol( values[0], __FILE__, __LINE__ );
+ }
+ else if ( strncmp(keyword, "reneighbor", MAX_LINE) == 0 )
+ {
+ control->reneighbor = sstrtol( values[0], __FILE__, __LINE__ );
+ }
+ else if ( strncmp(keyword, "vlist_buffer", MAX_LINE) == 0 )
+ {
+ val = sstrtod( values[0], __FILE__, __LINE__ );
+ control->vlist_cut = control->nonb_cut + val;
+ }
+ else if ( strncmp(keyword, "nbrhood_cutoff", MAX_LINE) == 0 )
+ {
+ val = sstrtod( values[0], __FILE__, __LINE__ );
+ control->bond_cut = val;
+ }
+ else if ( strncmp(keyword, "thb_cutoff", MAX_LINE) == 0 )
+ {
+ val = sstrtod( values[0], __FILE__, __LINE__ );
+ control->thb_cut = val;
+ }
+ else if ( strncmp(keyword, "hbond_cutoff", MAX_LINE) == 0 )
+ {
+ val = sstrtod( values[0], __FILE__, __LINE__ );
+ control->hbond_cut = val;
+ }
+ else if ( strncmp(keyword, "charge_method", MAX_LINE) == 0 )
+ {
+ control->charge_method = sstrtol( values[0], __FILE__, __LINE__ );
+ }
+ else if ( strncmp(keyword, "cm_q_net", MAX_LINE) == 0 )
+ {
+ val = sstrtod( values[0], __FILE__, __LINE__ );
+ control->cm_q_net = val;
+ }
+ else if ( strncmp(keyword, "cm_solver_type", MAX_LINE) == 0 )
+ {
+ control->cm_solver_type = sstrtol( values[0], __FILE__, __LINE__ );
+ }
+ else if ( strncmp(keyword, "cm_solver_max_iters", MAX_LINE) == 0 )
+ {
+ control->cm_solver_max_iters = sstrtol( values[0], __FILE__, __LINE__ );
+ }
+ else if ( strncmp(keyword, "cm_solver_restart", MAX_LINE) == 0 )
+ {
+ control->cm_solver_restart = sstrtol( values[0], __FILE__, __LINE__ );
+ }
+ else if ( strncmp(keyword, "cm_solver_q_err", MAX_LINE) == 0 )
+ {
+ val = sstrtod( values[0], __FILE__, __LINE__ );
+ control->cm_solver_q_err = val;
+ }
+ else if ( strncmp(keyword, "cm_domain_sparsity", MAX_LINE) == 0 )
+ {
+ val = sstrtod( values[0], __FILE__, __LINE__ );
+ control->cm_domain_sparsity = val;
+ if ( val < 1.0 )
+ {
+ control->cm_domain_sparsify_enabled = TRUE;
+ }
+ }
+ else if ( strncmp(keyword, "cm_init_guess_type", MAX_LINE) == 0 )
+ {
+ control->cm_init_guess_type = sstrtol( values[0], __FILE__, __LINE__ );
+ }
+ else if ( strncmp(keyword, "cm_init_guess_extrap1", MAX_LINE) == 0 )
+ {
+ control->cm_init_guess_extrap1 = sstrtol( values[0], __FILE__, __LINE__ );
+ }
+ else if ( strncmp(keyword, "cm_init_guess_extrap2", MAX_LINE) == 0 )
+ {
+ control->cm_init_guess_extrap2 = sstrtol( values[0], __FILE__, __LINE__ );
+ }
+ else if ( strncmp(keyword, "cm_init_guess_gd_model", MAX_LINE) == 0 )
+ {
+ strncpy( control->cm_init_guess_gd_model, values[0], sizeof(control->cm_init_guess_gd_model) - 1 );
+ control->cm_init_guess_gd_model[sizeof(control->cm_init_guess_gd_model) - 1] = '\0';
+ }
+ else if ( strncmp(keyword, "cm_init_guess_win_size", MAX_LINE) == 0 )
+ {
+ control->cm_init_guess_win_size = sstrtol( values[0], __FILE__, __LINE__ );
+ }
+ else if ( strncmp(keyword, "cm_solver_pre_comp_type", MAX_LINE) == 0 )
+ {
+ control->cm_solver_pre_comp_type = sstrtol( values[0], __FILE__, __LINE__ );
+ }
+ else if ( strncmp(keyword, "cm_solver_pre_comp_refactor", MAX_LINE) == 0 )
+ {
+ control->cm_solver_pre_comp_refactor = sstrtol( values[0], __FILE__, __LINE__ );
+ }
+ else if ( strncmp(keyword, "cm_solver_pre_comp_droptol", MAX_LINE) == 0 )
+ {
+ val = sstrtod( values[0], __FILE__, __LINE__ );
+ control->cm_solver_pre_comp_droptol = val;
+ }
+ else if ( strncmp(keyword, "cm_solver_pre_comp_sweeps", MAX_LINE) == 0 )
+ {
+ control->cm_solver_pre_comp_sweeps = sstrtol( values[0], __FILE__, __LINE__ );
+ }
+ else if ( strncmp(keyword, "cm_solver_pre_comp_sai_thres", MAX_LINE) == 0 )
+ {
+ val = sstrtod( values[0], __FILE__, __LINE__ );
+ control->cm_solver_pre_comp_sai_thres = val;
+ }
+ else if ( strncmp(keyword, "cm_solver_pre_app_type", MAX_LINE) == 0 )
+ {
+ control->cm_solver_pre_app_type = sstrtol( values[0], __FILE__, __LINE__ );
+ }
+ else if ( strncmp(keyword, "cm_solver_pre_app_jacobi_iters", MAX_LINE) == 0 )
+ {
+ control->cm_solver_pre_app_jacobi_iters = sstrtol( values[0], __FILE__, __LINE__ );
+ }
+ else if ( strncmp(keyword, "temp_init", MAX_LINE) == 0 )
+ {
+ val = sstrtod( values[0], __FILE__, __LINE__ );
+ control->T_init = val;
+
+ if ( control->T_init < 0.001 )
+ {
+ control->T_init = 0.001;
+ }
+ }
+ else if ( strncmp(keyword, "temp_final", MAX_LINE) == 0 )
+ {
+ val = sstrtod( values[0], __FILE__, __LINE__ );
+ control->T_final = val;
+
+ if ( control->T_final < 0.1 )
+ {
+ control->T_final = 0.1;
+ }
+ }
+ else if ( strncmp(keyword, "t_mass", MAX_LINE) == 0 )
+ {
+ val = sstrtod( values[0], __FILE__, __LINE__ );
+ /* convert from fs to s */
+ control->Tau_T = val * 1.0e-15;
+ }
+ else if ( strncmp(keyword, "t_mode", MAX_LINE) == 0 )
+ {
+ control->T_mode = sstrtol( values[0], __FILE__, __LINE__ );
+ }
+ else if ( strncmp(keyword, "t_rate", MAX_LINE) == 0 )
+ {
+ val = sstrtod( values[0], __FILE__, __LINE__ );
+ control->T_rate = val;
+ }
+ else if ( strncmp(keyword, "t_freq", MAX_LINE) == 0 )
+ {
+ val = sstrtod( values[0], __FILE__, __LINE__ );
+ control->T_freq = val;
+ }
+ else if ( strncmp(keyword, "pressure", MAX_LINE) == 0 )
+ {
+ if ( control->ensemble == iNPT )
+ {
+ val = sstrtod( values[0], __FILE__, __LINE__ );
+ control->P[0] = val;
+ control->P[1] = val;
+ control->P[2] = val;
+ }
+ else if ( control->ensemble == sNPT )
+ {
+ val = sstrtod( values[0], __FILE__, __LINE__ );
+ control->P[0] = val;
+
+ val = sstrtod( values[1], __FILE__, __LINE__ );
+ control->P[1] = val;
+
+ val = sstrtod( values[2], __FILE__, __LINE__ );
+ control->P[2] = val;
+ }
+ }
+ else if ( strncmp(keyword, "p_mass", MAX_LINE) == 0 )
+ {
+ if ( control->ensemble == iNPT )
+ {
+ val = sstrtod( values[0], __FILE__, __LINE__ );
+ control->Tau_P[0] = val * 1.0e-3; // convert p_mass from fs to ps
+ }
+ else if ( control->ensemble == sNPT )
+ {
+ val = sstrtod( values[0], __FILE__, __LINE__ );
+ control->Tau_P[0] = val * 1.0e-3; // convert p_mass from fs to ps
+
+ val = sstrtod( values[1], __FILE__, __LINE__ );
+ control->Tau_P[1] = val * 1.0e-3; // convert p_mass from fs to ps
+
+ val = sstrtod( values[2], __FILE__, __LINE__ );
+ control->Tau_P[2] = val * 1.0e-3; // convert p_mass from fs to ps
+ }
+ }
+ else if ( strncmp(keyword, "pt_mass", MAX_LINE) == 0 )
+ {
+ val = sstrtod( values[0], __FILE__, __LINE__ );
+ control->Tau_PT = val * 1.0e-3; // convert pt_mass from fs to ps
+ }
+ else if ( strncmp(keyword, "compress", MAX_LINE) == 0 )
+ {
+ val = sstrtod( values[0], __FILE__, __LINE__ );
+ control->compressibility = val;
+ }
+ else if ( strncmp(keyword, "press_mode", MAX_LINE) == 0 )
+ {
+ control->press_mode = sstrtol( values[0], __FILE__, __LINE__ );
+ }
+ else if ( strncmp(keyword, "compute_pressure", MAX_LINE) == 0 )
+ {
+ control->compute_pressure = sstrtol( values[0], __FILE__, __LINE__ );
+ }
+ else if ( strncmp(keyword, "remove_CoM_vel", MAX_LINE) == 0 )
+ {
+ control->remove_CoM_vel = sstrtol( values[0], __FILE__, __LINE__ );
+ }
+ else if ( strncmp(keyword, "energy_update_freq", MAX_LINE) == 0 )
+ {
+ out_control->log_update_freq = sstrtol( values[0], __FILE__, __LINE__ );
+ }
+ else if ( strncmp(keyword, "write_freq", MAX_LINE) == 0 )
+ {
+ out_control->write_steps = sstrtol( values[0], __FILE__, __LINE__ );
+ }
+ else if ( strncmp(keyword, "traj_compress", MAX_LINE) == 0 )
+ {
+ out_control->traj_compress = sstrtol( values[0], __FILE__, __LINE__ );
+
+ if ( out_control->traj_compress )
+ {
+ out_control->write = (int (*)(FILE *, const char *, ...)) &gzprintf;
+ }
+ else
+ {
+ out_control->write = &fprintf;
+ }
+ }
+ else if ( strncmp(keyword, "traj_format", MAX_LINE) == 0 )
+ {
+ out_control->traj_format = sstrtol( values[0], __FILE__, __LINE__ );
+
+ if ( out_control->traj_format == 0 )
+ {
+ out_control->write_header = &Write_Custom_Header;
+ out_control->append_traj_frame = &Append_Custom_Frame;
+ }
+ else if ( out_control->traj_format == 1 )
+ {
+ out_control->write_header = &Write_xyz_Header;
+ out_control->append_traj_frame = &Append_xyz_Frame;
+ }
+ }
+ else if ( strncmp(keyword, "traj_title", MAX_LINE) == 0 )
+ {
+ strncpy( out_control->traj_title, values[0], sizeof(out_control->traj_title) - 1 );
+ out_control->traj_title[sizeof(out_control->traj_title) - 1] = '\0';
+ }
+ else if ( strncmp(keyword, "atom_info", MAX_LINE) == 0 )
+ {
+ out_control->atom_format +=
+ sstrtol( values[0], __FILE__, __LINE__ ) * 4;
+ }
+ else if ( strncmp(keyword, "atom_velocities", MAX_LINE) == 0 )
+ {
+ out_control->atom_format +=
+ sstrtol( values[0], __FILE__, __LINE__ ) * 2;
+ }
+ else if ( strncmp(keyword, "atom_forces", MAX_LINE) == 0 )
+ {
+ out_control->atom_format +=
+ sstrtol( values[0], __FILE__, __LINE__ );
+ }
+ else if ( strncmp(keyword, "bond_info", MAX_LINE) == 0 )
+ {
+ out_control->bond_info = sstrtol( values[0], __FILE__, __LINE__ );
+ }
+ else if ( strncmp(keyword, "angle_info", MAX_LINE) == 0 )
+ {
+ out_control->angle_info = sstrtol( values[0], __FILE__, __LINE__ );
+ }
+ else if ( strncmp(keyword, "test_forces", MAX_LINE) == 0 )
+ {
+ ;
+ }
+ else if ( strncmp(keyword, "molec_anal", MAX_LINE) == 0 )
+ {
+ control->molec_anal = sstrtol( values[0], __FILE__, __LINE__ );
+ }
+ else if ( strncmp(keyword, "freq_molec_anal", MAX_LINE) == 0 )
+ {
+ control->freq_molec_anal = sstrtol( values[0], __FILE__, __LINE__ );
+ }
+ else if ( strncmp(keyword, "bond_graph_cutoff", MAX_LINE) == 0 )
+ {
+ val = sstrtod( values[0], __FILE__, __LINE__ );
+ control->bg_cut = val;
+ }
+ else if ( strncmp(keyword, "ignore", MAX_LINE) == 0 )
+ {
+ control->num_ignored = sstrtol( values[0], __FILE__, __LINE__ );
+ for ( i = 0; i < control->num_ignored; ++i )
+ control->ignore[ sstrtol( values[i + 1], __FILE__, __LINE__ ) ] = 1;
+ }
+ else if ( strncmp(keyword, "dipole_anal", MAX_LINE) == 0 )
+ {
+ control->dipole_anal = sstrtol( values[0], __FILE__, __LINE__ );
+ }
+ else if ( strncmp(keyword, "freq_dipole_anal", MAX_LINE) == 0 )
+ {
+ control->freq_dipole_anal = sstrtol( values[0], __FILE__, __LINE__ );
+ }
+ else if ( strncmp(keyword, "diffusion_coef", MAX_LINE) == 0 )
+ {
+ control->diffusion_coef = sstrtol( values[0], __FILE__, __LINE__ );
+ }
+ else if ( strncmp(keyword, "freq_diffusion_coef", MAX_LINE) == 0 )
+ {
+ control->freq_diffusion_coef = sstrtol( values[0], __FILE__, __LINE__ );
+ }
+ else if ( strncmp(keyword, "restrict_type", MAX_LINE) == 0 )
+ {
+ control->restrict_type = sstrtol( values[0], __FILE__, __LINE__ );
+ }
+ else
+ {
+ ret = FAILURE;
+ }
+
+ return ret;
+}
+
+
+void Read_Control_File( const char * const control_file, reax_system * const system,
+ control_params * const control, output_controls * const out_control )
+{
+ char *s, **tmp;
+ int c, i, ret;
+ FILE *fp;
+
+ fp = sfopen( control_file, "r" );
+
assert( fp != NULL );
/* assign default values */
@@ -166,432 +569,10 @@ void Read_Control_File( FILE* fp, reax_system *system, control_params* control,
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], "restart", MAX_LINE) == 0 ) {
- // ival = atoi(tmp[1]);
- // control->restart = 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], "random_vel", MAX_LINE) == 0 )
- {
- ival = atoi(tmp[1]);
- control->random_vel = ival;
- }
- else if ( strncmp(tmp[0], "reposition_atoms", MAX_LINE) == 0 )
- {
- ival = atoi(tmp[1]);
- control->reposition_atoms = ival;
- }
- 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]);
- /* convert from fs to ps */
- control->dt = val * 1.0e-3;
- }
- else if ( strncmp(tmp[0], "gpus_per_node", MAX_LINE) == 0 )
- {
- // skip since not applicable to shared memory code
- ;
- }
- else if ( strncmp(tmp[0], "proc_by_dim", MAX_LINE) == 0 )
- {
- // skip since not applicable to shared memory code
- ;
- }
- else if ( strncmp(tmp[0], "periodic_boundaries", MAX_LINE) == 0 )
- {
- ival = atoi( tmp[1] );
- control->periodic_boundaries = ival;
- }
- else if ( strncmp(tmp[0], "geo_format", MAX_LINE) == 0 )
- {
- ival = atoi( tmp[1] );
- control->geo_format = ival;
- }
- else if ( strncmp(tmp[0], "restrict_bonds", MAX_LINE) == 0 )
- {
- ival = atoi( tmp[1] );
- control->restrict_bonds = ival;
- }
- else if ( strncmp(tmp[0], "tabulate_long_range", MAX_LINE) == 0 )
- {
- ival = atoi( tmp[1] );
- control->tabulate = 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 = control->nonb_cut + val;
- }
- else if ( strncmp(tmp[0], "nbrhood_cutoff", MAX_LINE) == 0 )
- {
- val = atof(tmp[1]);
- control->bond_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], "charge_method", MAX_LINE) == 0 )
- {
- ival = atoi( tmp[1] );
- control->charge_method = 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_type", MAX_LINE) == 0 )
- {
- ival = atoi( tmp[1] );
- control->cm_init_guess_type = ival;
- }
- 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_init_guess_gd_model", MAX_LINE) == 0 )
- {
- strncpy( control->cm_init_guess_gd_model, tmp[1], sizeof(control->cm_init_guess_gd_model) - 1 );
- control->cm_init_guess_gd_model[sizeof(control->cm_init_guess_gd_model) - 1] = '\0';
- }
- else if ( strncmp(tmp[0], "cm_init_guess_win_size", MAX_LINE) == 0 )
- {
- ival = atoi( tmp[1] );
- control->cm_init_guess_win_size = 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], "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;
+ ret = Set_Control_Parameter( tmp[0],
+ (const char ** const) &tmp[1], control, out_control );
- 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 )
- {
- val = atof(tmp[1]);
- control->P[0] = val;
- control->P[1] = val;
- control->P[2] = val;
- }
- else if ( control->ensemble == sNPT )
- {
- val = atof(tmp[1]);
- control->P[0] = val;
-
- val = atof(tmp[2]);
- control->P[1] = val;
-
- val = atof(tmp[3]);
- control->P[2] = val;
- }
- }
- else if ( strncmp(tmp[0], "p_mass", MAX_LINE) == 0 )
- {
- if ( control->ensemble == iNPT )
- {
- val = atof(tmp[1]);
- control->Tau_P[0] = val * 1.0e-3; // convert p_mass from fs to ps
- }
- else if ( control->ensemble == sNPT )
- {
- val = atof(tmp[1]);
- control->Tau_P[0] = val * 1.0e-3; // convert p_mass from fs to ps
-
- val = atof(tmp[2]);
- control->Tau_P[1] = val * 1.0e-3; // convert p_mass from fs to ps
-
- val = atof(tmp[3]);
- control->Tau_P[2] = val * 1.0e-3; // convert p_mass from fs to ps
- }
- }
- else if ( strncmp(tmp[0], "pt_mass", MAX_LINE) == 0 )
- {
- val = atof(tmp[1]);
- control->Tau_PT = val * 1.0e-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 )
- {
- val = atoi(tmp[1]);
- control->press_mode = val;
- }
- else if ( strncmp(tmp[0], "compute_pressure", MAX_LINE) == 0 )
- {
- val = atoi(tmp[1]);
- control->compute_pressure = val;
- }
- else if ( strncmp(tmp[0], "remove_CoM_vel", MAX_LINE) == 0 )
- {
- val = atoi(tmp[1]);
- control->remove_CoM_vel = val;
- }
- else if ( strncmp(tmp[0], "energy_update_freq", MAX_LINE) == 0 )
- {
- ival = atoi(tmp[1]);
- out_control->log_update_freq = 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;
-
- if ( out_control->traj_compress )
- {
- out_control->write = (int (*)(FILE *, const char *, ...)) &gzprintf;
- }
- else
- {
- out_control->write = &fprintf;
- }
- }
- else if ( strncmp(tmp[0], "traj_format", MAX_LINE) == 0 )
- {
- ival = atoi(tmp[1]);
- out_control->traj_format = ival;
-
- if ( out_control->traj_format == 0 )
- {
- out_control->write_header = &Write_Custom_Header;
- out_control->append_traj_frame = &Append_Custom_Frame;
- }
- else if ( out_control->traj_format == 1 )
- {
- out_control->write_header = &Write_xyz_Header;
- out_control->append_traj_frame = &Append_xyz_Frame;
- }
- }
- 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_format += ival * 4;
- }
- else if ( strncmp(tmp[0], "atom_velocities", MAX_LINE) == 0 )
- {
- ival = atoi(tmp[1]);
- out_control->atom_format += ival * 2;
- }
- else if ( strncmp(tmp[0], "atom_forces", MAX_LINE) == 0 )
- {
- ival = atoi(tmp[1]);
- out_control->atom_format += 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 )
- {
- ival = atoi(tmp[1]);
- }
- else if ( strncmp(tmp[0], "molec_anal", MAX_LINE) == 0 )
- {
- ival = atoi(tmp[1]);
- control->molec_anal = ival;
- }
- else if ( strncmp(tmp[0], "freq_molec_anal", MAX_LINE) == 0 )
- {
- ival = atoi(tmp[1]);
- control->freq_molec_anal = ival;
- }
- else if ( strncmp(tmp[0], "bond_graph_cutoff", MAX_LINE) == 0 )
- {
- val = atof(tmp[1]);
- control->bg_cut = val;
- }
- 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], "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
+ if ( ret == FAILURE )
{
fprintf( stderr, "WARNING: unknown parameter %s\n", tmp[0] );
exit( UNKNOWN_OPTION );
@@ -651,6 +632,8 @@ void Read_Control_File( FILE* fp, reax_system *system, control_params* control,
control->ensemble, control->nsteps, control->dt, control->tabulate,
control->T, control->P[0], control->P[1], control->P[2] );
- fprintf(stderr, "control file read\n" );
+ fprintf( stderr, "control file read\n" );
#endif
+
+ sfclose( fp, "Read_Control_File::fp" );
}
diff --git a/sPuReMD/src/control.h b/sPuReMD/src/control.h
index 14fe5493..13ae7da0 100644
--- a/sPuReMD/src/control.h
+++ b/sPuReMD/src/control.h
@@ -24,7 +24,10 @@
#include "reax_types.h"
+int Set_Control_Parameter( const char * const, const char ** const,
+ control_params * const, output_controls * const );
-void Read_Control_File( FILE*, reax_system*, control_params*, output_controls* );
+void Read_Control_File( const char * const, reax_system * const, control_params * const,
+ output_controls * const );
#endif
diff --git a/sPuReMD/src/ffield.c b/sPuReMD/src/ffield.c
index 84028b33..569b6844 100644
--- a/sPuReMD/src/ffield.c
+++ b/sPuReMD/src/ffield.c
@@ -26,14 +26,17 @@
#include "tool_box.h"
-void Read_Force_Field( FILE *fp, reax_system * const system,
- reax_interaction * const reax )
+void Read_Force_Field( const char * const ffield_file,
+ reax_system * const system, reax_interaction * const reax )
{
char *s;
char **tmp;
char ****tor_flag;
int i, j, k, l, m, n, o, p, cnt;
real val;
+ FILE *fp;
+
+ fp = sfopen( ffield_file, "r" );
assert( fp != NULL );
@@ -54,7 +57,7 @@ void Read_Force_Field( FILE *fp, reax_system * const system,
Tokenize( s, &tmp, MAX_TOKEN_LEN );
/* reading the number of global parameters */
- n = atoi(tmp[0]);
+ n = sstrtol( tmp[0], __FILE__, __LINE__ );
if ( n < 1 )
{
fprintf( stderr, "[WARNING] number of globals in ffield file is 0!\n" );
@@ -83,7 +86,7 @@ void Read_Force_Field( FILE *fp, reax_system * const system,
fgets( s, MAX_LINE, fp );
Tokenize( s, &tmp, MAX_TOKEN_LEN );
- val = (real) atof(tmp[0]);
+ val = (real) sstrtod( tmp[0], __FILE__, __LINE__ );
reax->gp.l[i] = val;
}
@@ -91,7 +94,7 @@ void Read_Force_Field( FILE *fp, reax_system * const system,
/* next line is number of atom types and some comments */
fgets( s, MAX_LINE, fp );
Tokenize( s, &tmp, MAX_TOKEN_LEN );
- n = atoi(tmp[0]);
+ n = sstrtol( tmp[0], __FILE__, __LINE__ );
/* 3 lines of comments */
fgets( s, MAX_LINE, fp );
@@ -260,21 +263,21 @@ void Read_Force_Field( FILE *fp, reax_system * const system,
reax->sbp[i].name[j] = toupper( reax->sbp[i].name[j] );
}
- val = atof(tmp[1]);
+ val = sstrtod( tmp[1], __FILE__, __LINE__ );
reax->sbp[i].r_s = val;
- val = atof(tmp[2]);
+ val = sstrtod( tmp[2], __FILE__, __LINE__ );
reax->sbp[i].valency = val;
- val = atof(tmp[3]);
+ val = sstrtod( tmp[3], __FILE__, __LINE__ );
reax->sbp[i].mass = val;
- val = atof(tmp[4]);
+ val = sstrtod( tmp[4], __FILE__, __LINE__ );
reax->sbp[i].r_vdw = val;
- val = atof(tmp[5]);
+ val = sstrtod( tmp[5], __FILE__, __LINE__ );
reax->sbp[i].epsilon = val;
- val = atof(tmp[6]);
+ val = sstrtod( tmp[6], __FILE__, __LINE__ );
reax->sbp[i].gamma = val;
- val = atof(tmp[7]);
+ val = sstrtod( tmp[7], __FILE__, __LINE__ );
reax->sbp[i].r_pi = val;
- val = atof(tmp[8]);
+ val = sstrtod( tmp[8], __FILE__, __LINE__ );
reax->sbp[i].valency_e = val;
reax->sbp[i].nlp_opt = 0.5 * (reax->sbp[i].valency_e - reax->sbp[i].valency);
@@ -282,23 +285,23 @@ void Read_Force_Field( FILE *fp, reax_system * const system,
fgets( s, MAX_LINE, fp );
Tokenize( s, &tmp, MAX_TOKEN_LEN );
- val = atof(tmp[0]);
+ val = sstrtod( tmp[0], __FILE__, __LINE__ );
reax->sbp[i].alpha = val;
- val = atof(tmp[1]);
+ val = sstrtod( tmp[1], __FILE__, __LINE__ );
reax->sbp[i].gamma_w = val;
- val = atof(tmp[2]);
+ val = sstrtod( tmp[2], __FILE__, __LINE__ );
reax->sbp[i].valency_boc = val;
- val = atof(tmp[3]);
+ val = sstrtod( tmp[3], __FILE__, __LINE__ );
reax->sbp[i].p_ovun5 = val;
- val = atof(tmp[4]);
- val = atof(tmp[5]);
+ val = sstrtod( tmp[4], __FILE__, __LINE__ );
+ val = sstrtod( tmp[5], __FILE__, __LINE__ );
reax->sbp[i].chi = val;
- val = atof(tmp[6]);
+ val = sstrtod( tmp[6], __FILE__, __LINE__ );
reax->sbp[i].eta = 2.0 * val;
/* this is the parameter that is used to determine
* which type of atoms participate in h-bonds.
* 1 is for H - 2 for O, N, S - 0 for all others.*/
- val = atof(tmp[7]);
+ val = sstrtod( tmp[7], __FILE__, __LINE__ );
reax->sbp[i].p_hbond = (int)(val + 0.1);
//0.1 is to avoid from truncating down!
@@ -306,39 +309,39 @@ void Read_Force_Field( FILE *fp, reax_system * const system,
fgets( s, MAX_LINE, fp );
Tokenize( s, &tmp, MAX_TOKEN_LEN );
- val = atof(tmp[0]);
+ val = sstrtod( tmp[0], __FILE__, __LINE__ );
reax->sbp[i].r_pi_pi = val;
- val = atof(tmp[1]);
+ val = sstrtod( tmp[1], __FILE__, __LINE__ );
reax->sbp[i].p_lp2 = val;
- val = atof(tmp[2]);
- val = atof(tmp[3]);
+ val = sstrtod( tmp[2], __FILE__, __LINE__ );
+ val = sstrtod( tmp[3], __FILE__, __LINE__ );
reax->sbp[i].b_o_131 = val;
- val = atof(tmp[4]);
+ val = sstrtod( tmp[4], __FILE__, __LINE__ );
reax->sbp[i].b_o_132 = val;
- val = atof(tmp[5]);
+ val = sstrtod( tmp[5], __FILE__, __LINE__ );
reax->sbp[i].b_o_133 = val;
- val = atof(tmp[6]);
+ val = sstrtod( tmp[6], __FILE__, __LINE__ );
reax->sbp[i].b_s_acks2 = val;
- val = atof(tmp[7]);
+ val = sstrtod( tmp[7], __FILE__, __LINE__ );
/* line 4 */
fgets( s, MAX_LINE, fp );
Tokenize( s, &tmp, MAX_TOKEN_LEN );
- val = atof(tmp[0]);
+ val = sstrtod( tmp[0], __FILE__, __LINE__ );
reax->sbp[i].p_ovun2 = val;
- val = atof(tmp[1]);
+ val = sstrtod( tmp[1], __FILE__, __LINE__ );
reax->sbp[i].p_val3 = val;
- val = atof(tmp[2]);
- val = atof(tmp[3]);
+ val = sstrtod( tmp[2], __FILE__, __LINE__ );
+ val = sstrtod( tmp[3], __FILE__, __LINE__ );
reax->sbp[i].valency_val = val;
- val = atof(tmp[4]);
+ val = sstrtod( tmp[4], __FILE__, __LINE__ );
reax->sbp[i].p_val5 = val;
- val = atof(tmp[5]);
+ val = sstrtod( tmp[5], __FILE__, __LINE__ );
reax->sbp[i].rcore2 = val;
- val = atof(tmp[6]);
+ val = sstrtod( tmp[6], __FILE__, __LINE__ );
reax->sbp[i].ecore2 = val;
- val = atof(tmp[7]);
+ val = sstrtod( tmp[7], __FILE__, __LINE__ );
reax->sbp[i].acore2 = val;
/* inner-wall parameters present */
@@ -442,7 +445,7 @@ void Read_Force_Field( FILE *fp, reax_system * const system,
/* next line is number of two body combination and some comments */
fgets( s, MAX_LINE, fp );
Tokenize( s, &tmp, MAX_TOKEN_LEN );
- l = atoi(tmp[0]);
+ l = sstrtol( tmp[0], __FILE__, __LINE__ );
/* a line of comments */
fgets(s, MAX_LINE, fp);
@@ -453,35 +456,35 @@ void Read_Force_Field( FILE *fp, reax_system * const system,
fgets( s, MAX_LINE, fp );
Tokenize( s, &tmp, MAX_TOKEN_LEN );
- j = atoi(tmp[0]) - 1;
- k = atoi(tmp[1]) - 1;
+ j = sstrtol( tmp[0], __FILE__, __LINE__ ) - 1;
+ k = sstrtol( tmp[1], __FILE__, __LINE__ ) - 1;
if ( j < reax->num_atom_types && k < reax->num_atom_types )
{
- val = atof(tmp[2]);
+ val = sstrtod( tmp[2], __FILE__, __LINE__ );
reax->tbp[j][k].De_s = val;
reax->tbp[k][j].De_s = val;
- val = atof(tmp[3]);
+ val = sstrtod( tmp[3], __FILE__, __LINE__ );
reax->tbp[j][k].De_p = val;
reax->tbp[k][j].De_p = val;
- val = atof(tmp[4]);
+ val = sstrtod( tmp[4], __FILE__, __LINE__ );
reax->tbp[j][k].De_pp = val;
reax->tbp[k][j].De_pp = val;
- val = atof(tmp[5]);
+ val = sstrtod( tmp[5], __FILE__, __LINE__ );
reax->tbp[j][k].p_be1 = val;
reax->tbp[k][j].p_be1 = val;
- val = atof(tmp[6]);
+ val = sstrtod( tmp[6], __FILE__, __LINE__ );
reax->tbp[j][k].p_bo5 = val;
reax->tbp[k][j].p_bo5 = val;
- val = atof(tmp[7]);
+ val = sstrtod( tmp[7], __FILE__, __LINE__ );
reax->tbp[j][k].v13cor = val;
reax->tbp[k][j].v13cor = val;
- val = atof(tmp[8]);
+ val = sstrtod( tmp[8], __FILE__, __LINE__ );
reax->tbp[j][k].p_bo6 = val;
reax->tbp[k][j].p_bo6 = val;
- val = atof(tmp[9]);
+ val = sstrtod( tmp[9], __FILE__, __LINE__ );
reax->tbp[j][k].p_ovun1 = val;
reax->tbp[k][j].p_ovun1 = val;
@@ -489,28 +492,28 @@ void Read_Force_Field( FILE *fp, reax_system * const system,
fgets( s, MAX_LINE, fp );
Tokenize( s, &tmp, MAX_TOKEN_LEN );
- val = atof(tmp[0]);
+ val = sstrtod( tmp[0], __FILE__, __LINE__ );
reax->tbp[j][k].p_be2 = val;
reax->tbp[k][j].p_be2 = val;
- val = atof(tmp[1]);
+ val = sstrtod( tmp[1], __FILE__, __LINE__ );
reax->tbp[j][k].p_bo3 = val;
reax->tbp[k][j].p_bo3 = val;
- val = atof(tmp[2]);
+ val = sstrtod( tmp[2], __FILE__, __LINE__ );
reax->tbp[j][k].p_bo4 = val;
reax->tbp[k][j].p_bo4 = val;
- val = atof(tmp[3]);
+ val = sstrtod( tmp[3], __FILE__, __LINE__ );
- val = atof(tmp[4]);
+ val = sstrtod( tmp[4], __FILE__, __LINE__ );
reax->tbp[j][k].p_bo1 = val;
reax->tbp[k][j].p_bo1 = val;
- val = atof(tmp[5]);
+ val = sstrtod( tmp[5], __FILE__, __LINE__ );
reax->tbp[j][k].p_bo2 = val;
reax->tbp[k][j].p_bo2 = val;
- val = atof(tmp[6]);
+ val = sstrtod( tmp[6], __FILE__, __LINE__ );
reax->tbp[j][k].ovc = val;
reax->tbp[k][j].ovc = val;
- val = atof(tmp[7]);
+ val = sstrtod( tmp[7], __FILE__, __LINE__ );
}
}
@@ -568,54 +571,54 @@ void Read_Force_Field( FILE *fp, reax_system * const system,
* combination rules defined above */
fgets( s, MAX_LINE, fp );
Tokenize( s, &tmp, MAX_TOKEN_LEN );
- l = atoi(tmp[0]);
+ l = sstrtol( tmp[0], __FILE__, __LINE__ );
for ( i = 0; i < l; i++ )
{
fgets( s, MAX_LINE, fp );
Tokenize( s, &tmp, MAX_TOKEN_LEN );
- j = atoi(tmp[0]) - 1;
- k = atoi(tmp[1]) - 1;
+ j = sstrtol( tmp[0], __FILE__, __LINE__ ) - 1;
+ k = sstrtol( tmp[1], __FILE__, __LINE__ ) - 1;
if ( j < reax->num_atom_types && k < reax->num_atom_types )
{
- val = atof(tmp[2]);
+ val = sstrtod( tmp[2], __FILE__, __LINE__ );
if ( val > 0.0 )
{
reax->tbp[j][k].D = val;
reax->tbp[k][j].D = val;
}
- val = atof(tmp[3]);
+ val = sstrtod( tmp[3], __FILE__, __LINE__ );
if ( val > 0.0 )
{
reax->tbp[j][k].r_vdW = 2.0 * val;
reax->tbp[k][j].r_vdW = 2.0 * val;
}
- val = atof(tmp[4]);
+ val = sstrtod( tmp[4], __FILE__, __LINE__ );
if ( val > 0.0 )
{
reax->tbp[j][k].alpha = val;
reax->tbp[k][j].alpha = val;
}
- val = atof(tmp[5]);
+ val = sstrtod( tmp[5], __FILE__, __LINE__ );
if ( val > 0.0 )
{
reax->tbp[j][k].r_s = val;
reax->tbp[k][j].r_s = val;
}
- val = atof(tmp[6]);
+ val = sstrtod( tmp[6], __FILE__, __LINE__ );
if ( val > 0.0 )
{
reax->tbp[j][k].r_p = val;
reax->tbp[k][j].r_p = val;
}
- val = atof(tmp[7]);
+ val = sstrtod( tmp[7], __FILE__, __LINE__ );
if ( val > 0.0 )
{
reax->tbp[j][k].r_pp = val;
@@ -641,16 +644,16 @@ void Read_Force_Field( FILE *fp, reax_system * const system,
/* next line is number of 3-body params and some comments */
fgets( s, MAX_LINE, fp );
Tokenize( s, &tmp, MAX_TOKEN_LEN );
- l = atoi( tmp[0] );
+ l = sstrtol( tmp[0], __FILE__, __LINE__ );
for ( i = 0; i < l; i++ )
{
fgets( s, MAX_LINE, fp );
Tokenize( s, &tmp, MAX_TOKEN_LEN );
- j = atoi(tmp[0]) - 1;
- k = atoi(tmp[1]) - 1;
- m = atoi(tmp[2]) - 1;
+ j = sstrtol( tmp[0], __FILE__, __LINE__ ) - 1;
+ k = sstrtol( tmp[1], __FILE__, __LINE__ ) - 1;
+ m = sstrtol( tmp[2], __FILE__, __LINE__ ) - 1;
if ( j < reax->num_atom_types
&& k < reax->num_atom_types
@@ -660,31 +663,31 @@ void Read_Force_Field( FILE *fp, reax_system * const system,
reax->thbp[j][k][m].cnt++;
reax->thbp[m][k][j].cnt++;
- val = atof(tmp[3]);
+ val = sstrtod( tmp[3], __FILE__, __LINE__ );
reax->thbp[j][k][m].prm[cnt].theta_00 = val;
reax->thbp[m][k][j].prm[cnt].theta_00 = val;
- val = atof(tmp[4]);
+ val = sstrtod( tmp[4], __FILE__, __LINE__ );
reax->thbp[j][k][m].prm[cnt].p_val1 = val;
reax->thbp[m][k][j].prm[cnt].p_val1 = val;
- val = atof(tmp[5]);
+ val = sstrtod( tmp[5], __FILE__, __LINE__ );
reax->thbp[j][k][m].prm[cnt].p_val2 = val;
reax->thbp[m][k][j].prm[cnt].p_val2 = val;
- val = atof(tmp[6]);
+ val = sstrtod( tmp[6], __FILE__, __LINE__ );
reax->thbp[j][k][m].prm[cnt].p_coa1 = val;
reax->thbp[m][k][j].prm[cnt].p_coa1 = val;
- val = atof(tmp[7]);
+ val = sstrtod( tmp[7], __FILE__, __LINE__ );
reax->thbp[j][k][m].prm[cnt].p_val7 = val;
reax->thbp[m][k][j].prm[cnt].p_val7 = val;
- val = atof(tmp[8]);
+ val = sstrtod( tmp[8], __FILE__, __LINE__ );
reax->thbp[j][k][m].prm[cnt].p_pen1 = val;
reax->thbp[m][k][j].prm[cnt].p_pen1 = val;
- val = atof(tmp[9]);
+ val = sstrtod( tmp[9], __FILE__, __LINE__ );
reax->thbp[j][k][m].prm[cnt].p_val4 = val;
reax->thbp[m][k][j].prm[cnt].p_val4 = val;
}
@@ -717,17 +720,17 @@ void Read_Force_Field( FILE *fp, reax_system * const system,
/* next line is number of 4-body params and some comments */
fgets( s, MAX_LINE, fp );
Tokenize( s, &tmp, MAX_TOKEN_LEN );
- l = atoi( tmp[0] );
+ l = sstrtol( tmp[0], __FILE__, __LINE__ );
for ( i = 0; i < l; i++ )
{
fgets( s, MAX_LINE, fp );
Tokenize( s, &tmp, MAX_TOKEN_LEN );
- j = atoi(tmp[0]) - 1;
- k = atoi(tmp[1]) - 1;
- m = atoi(tmp[2]) - 1;
- n = atoi(tmp[3]) - 1;
+ j = sstrtol( tmp[0], __FILE__, __LINE__ ) - 1;
+ k = sstrtol( tmp[1], __FILE__, __LINE__ ) - 1;
+ m = sstrtol( tmp[2], __FILE__, __LINE__ ) - 1;
+ n = sstrtol( tmp[3], __FILE__, __LINE__ ) - 1;
/* this means the entry is not in compact form */
if ( j >= 0 && n >= 0 )
@@ -749,23 +752,23 @@ void Read_Force_Field( FILE *fp, reax_system * const system,
// reax->fbp[j][k][m][n].cnt++;
// reax->fbp[n][m][k][j].cnt++;
- val = atof(tmp[4]);
+ val = sstrtod( tmp[4], __FILE__, __LINE__ );
reax->fbp[j][k][m][n].prm[0].V1 = val;
reax->fbp[n][m][k][j].prm[0].V1 = val;
- val = atof(tmp[5]);
+ val = sstrtod( tmp[5], __FILE__, __LINE__ );
reax->fbp[j][k][m][n].prm[0].V2 = val;
reax->fbp[n][m][k][j].prm[0].V2 = val;
- val = atof(tmp[6]);
+ val = sstrtod( tmp[6], __FILE__, __LINE__ );
reax->fbp[j][k][m][n].prm[0].V3 = val;
reax->fbp[n][m][k][j].prm[0].V3 = val;
- val = atof(tmp[7]);
+ val = sstrtod( tmp[7], __FILE__, __LINE__ );
reax->fbp[j][k][m][n].prm[0].p_tor1 = val;
reax->fbp[n][m][k][j].prm[0].p_tor1 = val;
- val = atof(tmp[8]);
+ val = sstrtod( tmp[8], __FILE__, __LINE__ );
reax->fbp[j][k][m][n].prm[0].p_cot1 = val;
reax->fbp[n][m][k][j].prm[0].p_cot1 = val;
}
@@ -788,20 +791,20 @@ void Read_Force_Field( FILE *fp, reax_system * const system,
if ( tor_flag[p][k][m][o] == 0 )
{
- reax->fbp[p][k][m][o].prm[0].V1 = atof(tmp[4]);
- reax->fbp[p][k][m][o].prm[0].V2 = atof(tmp[5]);
- reax->fbp[p][k][m][o].prm[0].V3 = atof(tmp[6]);
- reax->fbp[p][k][m][o].prm[0].p_tor1 = atof(tmp[7]);
- reax->fbp[p][k][m][o].prm[0].p_cot1 = atof(tmp[8]);
+ reax->fbp[p][k][m][o].prm[0].V1 = sstrtod( tmp[4], __FILE__, __LINE__ );
+ reax->fbp[p][k][m][o].prm[0].V2 = sstrtod( tmp[5], __FILE__, __LINE__ );
+ reax->fbp[p][k][m][o].prm[0].V3 = sstrtod( tmp[6], __FILE__, __LINE__ );
+ reax->fbp[p][k][m][o].prm[0].p_tor1 = sstrtod( tmp[7], __FILE__, __LINE__ );
+ reax->fbp[p][k][m][o].prm[0].p_cot1 = sstrtod( tmp[8], __FILE__, __LINE__ );
}
if ( tor_flag[o][m][k][p] == 0 )
{
- reax->fbp[o][m][k][p].prm[0].V1 = atof(tmp[4]);
- reax->fbp[o][m][k][p].prm[0].V2 = atof(tmp[5]);
- reax->fbp[o][m][k][p].prm[0].V3 = atof(tmp[6]);
- reax->fbp[o][m][k][p].prm[0].p_tor1 = atof(tmp[7]);
- reax->fbp[o][m][k][p].prm[0].p_cot1 = atof(tmp[8]);
+ reax->fbp[o][m][k][p].prm[0].V1 = sstrtod( tmp[4], __FILE__, __LINE__ );
+ reax->fbp[o][m][k][p].prm[0].V2 = sstrtod( tmp[5], __FILE__, __LINE__ );
+ reax->fbp[o][m][k][p].prm[0].V3 = sstrtod( tmp[6], __FILE__, __LINE__ );
+ reax->fbp[o][m][k][p].prm[0].p_tor1 = sstrtod( tmp[7], __FILE__, __LINE__ );
+ reax->fbp[o][m][k][p].prm[0].p_cot1 = sstrtod( tmp[8], __FILE__, __LINE__ );
}
}
}
@@ -812,29 +815,29 @@ void Read_Force_Field( FILE *fp, reax_system * const system,
/* next line is number of hydrogen bond params and some comments */
fgets( s, MAX_LINE, fp );
Tokenize( s, &tmp, MAX_TOKEN_LEN );
- l = atoi( tmp[0] );
+ l = sstrtol( tmp[0], __FILE__, __LINE__ );
for ( i = 0; i < l; i++ )
{
fgets( s, MAX_LINE, fp );
Tokenize( s, &tmp, MAX_TOKEN_LEN );
- j = atoi(tmp[0]) - 1;
- k = atoi(tmp[1]) - 1;
- m = atoi(tmp[2]) - 1;
+ j = sstrtol( tmp[0], __FILE__, __LINE__ ) - 1;
+ k = sstrtol( tmp[1], __FILE__, __LINE__ ) - 1;
+ m = sstrtol( tmp[2], __FILE__, __LINE__ ) - 1;
if ( j < reax->num_atom_types && m < reax->num_atom_types )
{
- val = atof(tmp[3]);
+ val = sstrtod( tmp[3], __FILE__, __LINE__ );
reax->hbp[j][k][m].r0_hb = val;
- val = atof(tmp[4]);
+ val = sstrtod( tmp[4], __FILE__, __LINE__ );
reax->hbp[j][k][m].p_hb1 = val;
- val = atof(tmp[5]);
+ val = sstrtod( tmp[5], __FILE__, __LINE__ );
reax->hbp[j][k][m].p_hb2 = val;
- val = atof(tmp[6]);
+ val = sstrtod( tmp[6], __FILE__, __LINE__ );
reax->hbp[j][k][m].p_hb3 = val;
}
@@ -866,4 +869,6 @@ void Read_Force_Field( FILE *fp, reax_system * const system,
sfree( tor_flag, "Read_Force_Field::tor_flag" );
}
+
+ sfclose( fp, "Read_Force_Field::fp" );
}
diff --git a/sPuReMD/src/ffield.h b/sPuReMD/src/ffield.h
index 692e772c..88bb91a9 100644
--- a/sPuReMD/src/ffield.h
+++ b/sPuReMD/src/ffield.h
@@ -25,7 +25,8 @@
#include "reax_types.h"
-void Read_Force_Field( FILE *, reax_system * const, reax_interaction * const );
+void Read_Force_Field( const char * const, reax_system * const,
+ reax_interaction * const );
#endif
diff --git a/sPuReMD/src/forces.c b/sPuReMD/src/forces.c
index cf3efeca..3b6afe5b 100644
--- a/sPuReMD/src/forces.c
+++ b/sPuReMD/src/forces.c
@@ -168,7 +168,7 @@ static void Compute_NonBonded_Forces( reax_system *system, control_params *contr
+ data->timing.cm_solver_vector_ops
+ data->timing.cm_solver_orthog
+ data->timing.cm_solver_tri_solve;
- data->timing.cm_total_loss = ZERO;
+ data->timing.cm_total_loss = 0.0;
}
else
{
diff --git a/sPuReMD/src/geo_tools.c b/sPuReMD/src/geo_tools.c
index 068f1408..88f8799c 100644
--- a/sPuReMD/src/geo_tools.c
+++ b/sPuReMD/src/geo_tools.c
@@ -351,12 +351,9 @@ void Read_PDB( const char * const pdb_file, reax_system* system, control_params
char occupancy[7], temp_factor[7];
char seg_id[5], element[3], charge[3];
char alt_loc, chain_id, icode;
- char *endptr;
rvec x;
reax_atom *atom;
- endptr = NULL;
-
pdb = sfopen( pdb_file, "r" );
Allocate_Tokenizer_Space( &s, MAX_LINE, &s1, MAX_LINE,
@@ -459,8 +456,9 @@ void Read_PDB( const char * const pdb_file, reax_system* system, control_params
}
/* if the point is inside my_box, add it to my lists */
- Make_Point( strtod( s_x, &endptr ), strtod( s_y, &endptr ),
- strtod( s_z, &endptr ), &x );
+ Make_Point( sstrtod( s_x, __FILE__, __LINE__ ),
+ sstrtod( s_y, __FILE__, __LINE__ ),
+ sstrtod( s_z, __FILE__, __LINE__ ), &x );
Fit_to_Periodic_Box( &system->box, x );
@@ -468,7 +466,7 @@ void Read_PDB( const char * const pdb_file, reax_system* system, control_params
{
/* store orig_id, type, name and coord info of the new atom */
atom = &system->atoms[top];
- pdb_serial = (int) strtod( serial, &endptr );
+ pdb_serial = (int) sstrtod( serial, __FILE__, __LINE__ );
workspace->orig_id[top] = pdb_serial;
Trim_Spaces( element, sizeof(element) );
@@ -508,13 +506,13 @@ void Read_PDB( const char * const pdb_file, reax_system* system, control_params
// "CONECT line exceeds max num restrictions allowed.\n" );
/* read bond restrictions */
- // if( is_Valid_Serial( workspace, pdb_serial = atoi(tmp[1]) ) )
+ // if( is_Valid_Serial( workspace, pdb_serial = sstrtol( tmp[1]), __FILE__, __LINE__ ) )
// ratom = workspace->map_serials[ pdb_serial ];
// workspace->restricted[ ratom ] = c1 - 2;
// for( i = 2; i < c1; ++i )
// {
- // if( is_Valid_Serial(workspace, pdb_serial = atoi(tmp[i])) )
+ // if( is_Valid_Serial(workspace, pdb_serial = sstrtol( tmp[i], __FILE__, __LINE__ )) )
// workspace->restricted_list[ ratom ][ i-2 ] =
// workspace->map_serials[ pdb_serial ];
// }
@@ -625,10 +623,8 @@ void Read_BGF( const char * const bgf_file, reax_system* system, control_params
char element[6], charge[9];
char chain_id;
char s_a[12], s_b[12], s_c[12], s_alpha[12], s_beta[12], s_gamma[12];
- char *endptr;
int i, n, atom_cnt, token_cnt, bgf_serial, ratom, crystx_found;
- endptr = NULL;
ratom = 0;
crystx_found = FALSE;
@@ -727,15 +723,15 @@ void Read_BGF( const char * const bgf_file, reax_system* system, control_params
// occupancy, temp_factor, charge );
/* add to mapping */
- bgf_serial = strtod( serial, &endptr );
+ bgf_serial = sstrtod( serial, __FILE__, __LINE__ );
Check_Input_Range( bgf_serial, 0, MAX_ATOM_ID, "Invalid bgf serial" );
workspace->map_serials[ bgf_serial ] = atom_cnt;
workspace->orig_id[ atom_cnt ] = bgf_serial;
/* copy atomic positions */
- system->atoms[atom_cnt].x[0] = strtod( s_x, &endptr );
- system->atoms[atom_cnt].x[1] = strtod( s_y, &endptr );
- system->atoms[atom_cnt].x[2] = strtod( s_z, &endptr );
+ system->atoms[atom_cnt].x[0] = sstrtod( s_x, __FILE__, __LINE__ );
+ system->atoms[atom_cnt].x[1] = sstrtod( s_y, __FILE__, __LINE__ );
+ system->atoms[atom_cnt].x[2] = sstrtod( s_z, __FILE__, __LINE__ );
/* atom name and type */
strncpy( system->atoms[atom_cnt].name, atom_name,
@@ -763,8 +759,12 @@ void Read_BGF( const char * const bgf_file, reax_system* system, control_params
s_a, s_b, s_c, s_alpha, s_beta, s_gamma );
/* compute full volume tensor from the angles */
- Setup_Box( atof(s_a), atof(s_b), atof(s_c),
- atof(s_alpha), atof(s_beta), atof(s_gamma),
+ Setup_Box( sstrtod( s_a, __FILE__, __LINE__ ),
+ sstrtod( s_b, __FILE__, __LINE__ ),
+ sstrtod( s_c, __FILE__, __LINE__ ),
+ sstrtod( s_alpha, __FILE__, __LINE__ ),
+ sstrtod( s_beta, __FILE__, __LINE__ ),
+ sstrtod( s_gamma, __FILE__, __LINE__ ),
&system->box );
crystx_found = TRUE;
@@ -778,7 +778,7 @@ void Read_BGF( const char * const bgf_file, reax_system* system, control_params
"CONECT line exceeds max restrictions allowed.\n" );
/* read bond restrictions */
- bgf_serial = atoi(tokens[1]);
+ bgf_serial = sstrtol( tokens[1], __FILE__, __LINE__ );
if ( is_Valid_Serial( workspace, bgf_serial ) )
{
ratom = workspace->map_serials[ bgf_serial ];
@@ -787,7 +787,7 @@ void Read_BGF( const char * const bgf_file, reax_system* system, control_params
workspace->restricted[ ratom ] = token_cnt - 2;
for ( i = 2; i < token_cnt; ++i )
{
- if ( is_Valid_Serial( workspace, bgf_serial = atoi(tokens[i]) ) )
+ if ( is_Valid_Serial( workspace, bgf_serial = sstrtol( tokens[i], __FILE__, __LINE__ )) )
{
workspace->restricted_list[ratom][i - 2] =
workspace->map_serials[bgf_serial];
diff --git a/sPuReMD/src/init_md.c b/sPuReMD/src/init_md.c
index caf475ee..0b13e1bc 100644
--- a/sPuReMD/src/init_md.c
+++ b/sPuReMD/src/init_md.c
@@ -259,22 +259,22 @@ static void Init_Simulation_Data( reax_system *system, control_params *control,
/* init timing info */
data->timing.start = Get_Time( );
data->timing.total = data->timing.start;
- data->timing.nbrs = 0;
- data->timing.init_forces = 0;
- data->timing.bonded = 0;
- data->timing.nonb = 0;
- data->timing.cm = ZERO;
- data->timing.cm_sort_mat_rows = ZERO;
- data->timing.cm_solver_pre_comp = ZERO;
- data->timing.cm_solver_pre_app = ZERO;
+ data->timing.nbrs = 0.0;
+ data->timing.init_forces = 0.0;
+ data->timing.bonded = 0.0;
+ data->timing.nonb = 0.0;
+ data->timing.cm = 0.0;
+ data->timing.cm_sort_mat_rows = 0.0;
+ data->timing.cm_solver_pre_comp = 0.0;
+ data->timing.cm_solver_pre_app = 0.0;
data->timing.cm_solver_iters = 0;
- data->timing.cm_solver_spmv = ZERO;
- data->timing.cm_solver_vector_ops = ZERO;
- data->timing.cm_solver_orthog = ZERO;
- data->timing.cm_solver_tri_solve = ZERO;
- data->timing.cm_last_pre_comp = ZERO;
- data->timing.cm_total_loss = ZERO;
- data->timing.cm_optimum = ZERO;
+ data->timing.cm_solver_spmv = 0.0;
+ data->timing.cm_solver_vector_ops = 0.0;
+ data->timing.cm_solver_orthog = 0.0;
+ data->timing.cm_solver_tri_solve = 0.0;
+ data->timing.cm_last_pre_comp = 0.0;
+ data->timing.cm_total_loss = 0.0;
+ data->timing.cm_optimum = 0.0;
}
diff --git a/sPuReMD/src/io_tools.c b/sPuReMD/src/io_tools.c
index ad19a383..819d7678 100644
--- a/sPuReMD/src/io_tools.c
+++ b/sPuReMD/src/io_tools.c
@@ -642,7 +642,7 @@ void Output_Results( reax_system *system, control_params *control,
data->timing.nonb * f_update,
data->timing.cm * f_update,
data->timing.cm_sort_mat_rows * f_update,
- (double)data->timing.cm_solver_iters * f_update,
+ (double) data->timing.cm_solver_iters * f_update,
data->timing.cm_solver_pre_comp * f_update,
data->timing.cm_solver_pre_app * f_update,
data->timing.cm_solver_spmv * f_update,
@@ -651,19 +651,19 @@ void Output_Results( reax_system *system, control_params *control,
data->timing.cm_solver_tri_solve * f_update );
data->timing.total = Get_Time( );
- data->timing.nbrs = 0;
- data->timing.init_forces = 0;
- data->timing.bonded = 0;
- data->timing.nonb = 0;
- data->timing.cm = ZERO;
- data->timing.cm_sort_mat_rows = ZERO;
- data->timing.cm_solver_pre_comp = ZERO;
- data->timing.cm_solver_pre_app = ZERO;
+ data->timing.nbrs = 0.0;
+ data->timing.init_forces = 0.0;
+ data->timing.bonded = 0.0;
+ data->timing.nonb = 0.0;
+ data->timing.cm = 0.0;
+ data->timing.cm_sort_mat_rows = 0.0;
+ data->timing.cm_solver_pre_comp = 0.0;
+ data->timing.cm_solver_pre_app = 0.0;
data->timing.cm_solver_iters = 0;
- data->timing.cm_solver_spmv = ZERO;
- data->timing.cm_solver_vector_ops = ZERO;
- data->timing.cm_solver_orthog = ZERO;
- data->timing.cm_solver_tri_solve = ZERO;
+ data->timing.cm_solver_spmv = 0.0;
+ data->timing.cm_solver_vector_ops = 0.0;
+ data->timing.cm_solver_orthog = 0.0;
+ data->timing.cm_solver_tri_solve = 0.0;
fflush( out_control->out );
fflush( out_control->pot );
diff --git a/sPuReMD/src/lin_alg.c b/sPuReMD/src/lin_alg.c
index 7fe4ec4b..53f736e3 100644
--- a/sPuReMD/src/lin_alg.c
+++ b/sPuReMD/src/lin_alg.c
@@ -1153,7 +1153,7 @@ real FG_ICHOLT( const sparse_matrix * const A, const real * droptol,
#endif
for ( i = 0; i < A->n; ++i )
{
- if ( A->val[A->start[i + 1] - 1] < ZERO )
+ if ( A->val[A->start[i + 1] - 1] < 0.0 )
{
gamma[i] = -1.0;
}
@@ -1224,7 +1224,7 @@ real FG_ICHOLT( const sparse_matrix * const A, const real * droptol,
y = U_T_temp.start[U_T_temp.j[pj]];
ei_y = U_T_temp.start[U_T_temp.j[pj] + 1];
- sum = ZERO;
+ sum = 0.0;
/* sparse vector-sparse vector inner product for nonzero (i, j):
* dot( U^T(i,1:j-1), U^T(j,1:j-1) ) */
@@ -1252,7 +1252,7 @@ real FG_ICHOLT( const sparse_matrix * const A, const real * droptol,
{
#if defined(DEBUG_FOCUS)
/* sanity check */
- if ( sum < ZERO )
+ if ( sum < 0.0 )
{
fprintf( stderr, "[ERROR] Numeric breakdown in FG_ICHOLT. Terminating.\n");
fprintf( stderr, " [INFO] DAD(%5d,%5d) = %10.3f\n",
@@ -1362,7 +1362,7 @@ real FG_ILUT( const sparse_matrix * const A, const real * droptol,
#endif
for ( i = 0; i < A->n; ++i )
{
- if ( A->val[A->start[i + 1] - 1] < ZERO )
+ if ( A->val[A->start[i + 1] - 1] < 0.0 )
{
gamma[i] = -1.0;
}
@@ -1452,7 +1452,7 @@ real FG_ILUT( const sparse_matrix * const A, const real * droptol,
y = L_temp.start[L_temp.j[pj]];
ei_y = L_temp.start[L_temp.j[pj] + 1];
- sum = ZERO;
+ sum = 0.0;
/* sparse vector-sparse vector inner products for nonzero (i, j):
* dot( L(i,1:j-1), U^T(j,1:j-1) )
@@ -1505,7 +1505,7 @@ real FG_ILUT( const sparse_matrix * const A, const real * droptol,
y = L_temp.start[L_temp.j[pj]];
ei_y = L_temp.start[L_temp.j[pj] + 1];
- sum = ZERO;
+ sum = 0.0;
/* sparse vector-sparse vector inner products for nonzero (i, j):
* dot( L(j,1:j-1), U^T(i,1:j-1) )
diff --git a/sPuReMD/src/reax_types.h b/sPuReMD/src/reax_types.h
index fa5401d1..097c99c9 100644
--- a/sPuReMD/src/reax_types.h
+++ b/sPuReMD/src/reax_types.h
@@ -195,7 +195,6 @@
/* min. temperature scaler for atomic positions and velocities in NPT ensembles */
#define MIN_dT (0.0)
-#define ZERO (0.000000000000000e+00)
#define ALMOST_ZERO (1.0e-10)
#define NEG_INF (-1.0e10)
#define NO_BOND (1.0e-3)
@@ -1073,7 +1072,7 @@ struct control_params
struct thermostat
{
- /* temperature scaler of the system at the current simulation step */
+ /* temperature scalar of the system at the current simulation step */
real T;
/* temperature tensor of the system at the current simulation step */
rtensor Temp;
diff --git a/sPuReMD/src/spuremd.c b/sPuReMD/src/spuremd.c
index 41fa92f7..694dfc78 100644
--- a/sPuReMD/src/spuremd.c
+++ b/sPuReMD/src/spuremd.c
@@ -21,10 +21,9 @@
#include "spuremd.h"
+#include "allocate.h"
#include "analyze.h"
-#if defined(DEBUG_FOCUS)
- #include "box.h"
-#endif
+#include "box.h"
#include "control.h"
#include "ffield.h"
#include "forces.h"
@@ -101,49 +100,196 @@ static void Read_Input_Files( const char * const geo_file,
simulation_data * const data, static_storage * const workspace,
output_controls * const out_control )
{
- FILE *ffield, *ctrl;
-
- ffield = sfopen( ffield_file, "r" );
- ctrl = sfopen( control_file, "r" );
-
- Read_Force_Field( ffield, system, &system->reax_param );
-
- Read_Control_File( ctrl, system, control, out_control );
-
- if ( control->geo_format == CUSTOM )
+ if ( ffield_file != NULL )
{
- Read_Geo( geo_file, system, control, data, workspace );
+ Read_Force_Field( ffield_file, system, &system->reax_param );
}
- else if ( control->geo_format == PDB )
+
+ if ( control_file != NULL )
{
- Read_PDB( geo_file, system, control, data, workspace );
+ Read_Control_File( control_file, system, control, out_control );
}
- else if ( control->geo_format == BGF )
+
+ if ( geo_file != NULL )
{
- Read_BGF( geo_file, system, control, data, workspace );
+ if ( control->geo_format == CUSTOM )
+ {
+ Read_Geo( geo_file, system, control, data, workspace );
+ }
+ else if ( control->geo_format == PDB )
+ {
+ Read_PDB( geo_file, system, control, data, workspace );
+ }
+ else if ( control->geo_format == BGF )
+ {
+ Read_BGF( geo_file, system, control, data, workspace );
+ }
+ else if ( control->geo_format == ASCII_RESTART )
+ {
+ Read_ASCII_Restart( geo_file, system, control, data, workspace );
+ control->restart = TRUE;
+ }
+ else if ( control->geo_format == BINARY_RESTART )
+ {
+ Read_Binary_Restart( geo_file, system, control, data, workspace );
+ control->restart = TRUE;
+ }
+ else
+ {
+ fprintf( stderr, "[ERROR] unknown geo file format. terminating!\n" );
+ exit( INVALID_GEO );
+ }
}
- else if ( control->geo_format == ASCII_RESTART )
+
+#if defined(DEBUG_FOCUS)
+ Print_Box( &system->box, stderr );
+#endif
+}
+
+
+/* Allocate top-level data structures and parse input files
+ * for the first simulation
+ *
+ * qm_num_atoms: num. atoms in the QM region
+ * qm_types: element types for QM atoms
+ * qm_pos_x: x-coordinate of QM atom positions, in Angstroms
+ * qm_pos_y: y-coordinate of QM atom positions, in Angstroms
+ * qm_pos_z: z-coordinate of QM atom positions, in Angstroms
+ * mm_num_atoms: num. atoms in the MM region
+ * mm_types: element types for MM atoms
+ * mm_pos_x: x-coordinate of MM atom positions, in Angstroms
+ * mm_pos_y: y-coordinate of MM atom positions, in Angstroms
+ * mm_pos_z: z-coordinate of MM atom positions, in Angstroms
+ * mm_q: charge of MM atom, in Coulombs
+ * sim_box: simulation box information, where the entries are
+ * - box length per dimension (3 entries)
+ * - angles per dimension (3 entries)
+ * ffield_file: file containing force field parameters
+ * control_file: file containing simulation parameters
+ */
+void * setup_qmmm_( int qm_num_atoms, const int * const qm_types,
+ const double * const qm_pos_x, const double * const qm_pos_y,
+ const double * const qm_pos_z,
+ int mm_num_atoms, const int * const mm_types,
+ const double * const mm_pos_x, const double * const mm_pos_y,
+ const double * const mm_pos_z, const double * const mm_q,
+ const double * const sim_box,
+ const char * const ffield_file,
+ const char * const control_file )
+{
+ int i;
+// char atom_name[9];
+ rvec x;
+ spuremd_handle *spmd_handle;
+
+ /* top-level allocation */
+ spmd_handle = (spuremd_handle*) smalloc( sizeof(spuremd_handle),
+ "setup::spmd_handle" );
+
+ /* second-level allocations */
+ spmd_handle->system = smalloc( sizeof(reax_system),
+ "Setup::spmd_handle->system" );
+ spmd_handle->system->prealloc_allocated = FALSE;
+ spmd_handle->system->ffield_params_allocated = FALSE;
+ spmd_handle->system->g.allocated = FALSE;
+
+ spmd_handle->control = smalloc( sizeof(control_params),
+ "Setup::spmd_handle->control" );
+
+ spmd_handle->data = smalloc( sizeof(simulation_data),
+ "Setup::spmd_handle->data" );
+
+ spmd_handle->workspace = smalloc( sizeof(static_storage),
+ "Setup::spmd_handle->workspace" );
+ spmd_handle->workspace->H.allocated = FALSE;
+ spmd_handle->workspace->H_full.allocated = FALSE;
+ spmd_handle->workspace->H_sp.allocated = FALSE;
+ spmd_handle->workspace->H_p.allocated = FALSE;
+ spmd_handle->workspace->H_spar_patt.allocated = FALSE;
+ spmd_handle->workspace->H_spar_patt_full.allocated = FALSE;
+ spmd_handle->workspace->H_app_inv.allocated = FALSE;
+ spmd_handle->workspace->L.allocated = FALSE;
+ spmd_handle->workspace->U.allocated = FALSE;
+
+ spmd_handle->lists = smalloc( sizeof(reax_list *) * LIST_N,
+ "Setup::spmd_handle->lists" );
+ for ( i = 0; i < LIST_N; ++i )
{
- Read_ASCII_Restart( geo_file, system, control, data, workspace );
- control->restart = TRUE;
+ spmd_handle->lists[i] = smalloc( sizeof(reax_list),
+ "Setup::spmd_handle->lists[i]" );
+ spmd_handle->lists[i]->allocated = FALSE;
}
- else if ( control->geo_format == BINARY_RESTART )
+ spmd_handle->out_control = smalloc( sizeof(output_controls),
+ "Setup::spmd_handle->out_control" );
+
+ spmd_handle->output_enabled = FALSE;
+ spmd_handle->realloc = TRUE;
+ spmd_handle->callback = NULL;
+ spmd_handle->data->sim_id = 0;
+
+ spmd_handle->system->N = qm_num_atoms + mm_num_atoms;
+
+ PreAllocate_Space( spmd_handle->system, spmd_handle->control,
+ spmd_handle->workspace, spmd_handle->system->N );
+
+ Setup_Box( sim_box[0], sim_box[1], sim_box[2],
+ sim_box[3], sim_box[4], sim_box[5],
+ &spmd_handle->system->box );
+
+ for ( i = 0; i < qm_num_atoms; ++i )
{
- Read_Binary_Restart( geo_file, system, control, data, workspace );
- control->restart = TRUE;
+ x[0] = qm_pos_x[i];
+ x[1] = qm_pos_y[i];
+ x[2] = qm_pos_z[i];
+
+ Fit_to_Periodic_Box( &spmd_handle->system->box, x );
+
+ spmd_handle->workspace->orig_id[i] = i + 1;
+// spmd_handle->system->atoms[i].type = Get_Atom_Type( &system->reax_param,
+// element, sizeof(element) );
+ spmd_handle->system->atoms[i].type = qm_types[i];
+// strncpy( spmd_handle->system->atoms[i].name, atom_name,
+// sizeof(spmd_handle->system->atoms[i].name) - 1 );
+// spmd_handle->system->atoms[i].name[sizeof(spmd_handle->system->atoms[i].name) - 1] = '\0';
+ rvec_Copy( spmd_handle->system->atoms[i].x, x );
+ rvec_MakeZero( spmd_handle->system->atoms[i].v );
+ rvec_MakeZero( spmd_handle->system->atoms[i].f );
+ spmd_handle->system->atoms[i].q = 0.0;
+
+// mask[i] = 1;
}
- else
+
+ for ( i = qm_num_atoms; i < qm_num_atoms + mm_num_atoms; ++i )
{
- fprintf( stderr, "[ERROR] unknown geo file format. terminating!\n" );
- exit( INVALID_GEO );
+ x[0] = mm_pos_x[i - qm_num_atoms];
+ x[1] = mm_pos_y[i - qm_num_atoms];
+ x[2] = mm_pos_z[i - qm_num_atoms];
+
+ Fit_to_Periodic_Box( &spmd_handle->system->box, x );
+
+ spmd_handle->workspace->orig_id[i] = i + 1;
+// spmd_handle->system->atoms[i].type = Get_Atom_Type( &system->reax_param,
+// element, sizeof(element) );
+ spmd_handle->system->atoms[i].type = mm_types[i - qm_num_atoms];
+// strncpy( spmd_handle->system->atoms[i].name, atom_name,
+// sizeof(spmd_handle->system->atoms[i].name) - 1 );
+// spmd_handle->system->atoms[i].name[sizeof(spmd_handle->system->atoms[i].name) - 1] = '\0';
+ rvec_Copy( spmd_handle->system->atoms[i].x, x );
+ rvec_MakeZero( spmd_handle->system->atoms[i].v );
+ rvec_MakeZero( spmd_handle->system->atoms[i].f );
+ spmd_handle->system->atoms[i].q = mm_q[i - qm_num_atoms];
+
+// mask[i] = 0;
}
- sfclose( ffield, "Read_Input_Files::ffield" );
- sfclose( ctrl, "Read_Input_Files::ctrl" );
+ Read_Input_Files( NULL, ffield_file, control_file,
+ spmd_handle->system, spmd_handle->control,
+ spmd_handle->data, spmd_handle->workspace,
+ spmd_handle->out_control );
-#if defined(DEBUG_FOCUS)
- Print_Box( &system->box, stderr );
-#endif
+ spmd_handle->system->N_max = (int) CEIL( SAFE_ZONE * spmd_handle->system->N );
+
+ return (void *) spmd_handle;
}
@@ -154,7 +300,7 @@ static void Read_Input_Files( const char * const geo_file,
* ffield_file: file containing force field parameters
* control_file: file containing simulation parameters
*/
-void* setup( const char * const geo_file, const char * const ffield_file,
+void * setup( const char * const geo_file, const char * const ffield_file,
const char * const control_file )
{
int i;
@@ -212,7 +358,7 @@ void* setup( const char * const geo_file, const char * const ffield_file,
spmd_handle->system->N_max = (int) CEIL( SAFE_ZONE * spmd_handle->system->N );
- return (void*) spmd_handle;
+ return (void *) spmd_handle;
}
@@ -220,6 +366,8 @@ void* setup( const char * const geo_file, const char * const ffield_file,
*
* handle: pointer to wrapper struct with top-level data structures
* callback: function pointer to attach for callback
+ *
+ * returns: SPUREMD_SUCCESS upon success, SPUREMD_FAILURE otherwise
*/
int setup_callback( const void * const handle, const callback_function callback )
{
@@ -243,6 +391,8 @@ int setup_callback( const void * const handle, const callback_function callback
/* Run the simulation according to the prescribed parameters
*
* handle: pointer to wrapper struct with top-level data structures
+ *
+ * returns: SPUREMD_SUCCESS upon success, SPUREMD_FAILURE otherwise
*/
int simulate( const void * const handle )
{
@@ -391,6 +541,8 @@ int simulate( const void * const handle )
/* Deallocate all data structures post-simulation
*
* handle: pointer to wrapper struct with top-level data structures
+ *
+ * returns: SPUREMD_SUCCESS upon success, SPUREMD_FAILURE otherwise
*/
int cleanup( const void * const handle )
{
@@ -427,6 +579,139 @@ int cleanup( const void * const handle )
}
+/* Reset for the next simulation by parsing input files and triggering
+ * reallocation if more space is needed
+ *
+ * handle: pointer to wrapper struct with top-level data structures
+ * qm_num_atoms: num. atoms in the QM region
+ * qm_types: element types for QM atoms
+ * qm_pos_x: x-coordinate of QM atom positions, in Angstroms
+ * qm_pos_y: y-coordinate of QM atom positions, in Angstroms
+ * qm_pos_z: z-coordinate of QM atom positions, in Angstroms
+ * mm_num_atoms: num. atoms in the MM region
+ * mm_types: element types for MM atoms
+ * mm_pos_x: x-coordinate of MM atom positions, in Angstroms
+ * mm_pos_y: y-coordinate of MM atom positions, in Angstroms
+ * mm_pos_z: z-coordinate of MM atom positions, in Angstroms
+ * mm_q: charge of MM atom, in Coulombs
+ * sim_box: simulation box information, where the entries are
+ * - box length per dimension (3 entries)
+ * - angles per dimension (3 entries)
+ * ffield_file: file containing force field parameters
+ * control_file: file containing simulation parameters
+ *
+ * returns: SPUREMD_SUCCESS upon success, SPUREMD_FAILURE otherwise
+ */
+int reset_qmmm_( const void * const handle,
+ int qm_num_atoms, const int * const qm_types,
+ const double * const qm_pos_x, const double * const qm_pos_y,
+ const double * const qm_pos_z,
+ int mm_num_atoms, const int * const mm_types,
+ const double * const mm_pos_x, const double * const mm_pos_y,
+ const double * const mm_pos_z, const double * const mm_q,
+ const double * const sim_box,
+ const char * const ffield_file, const char * const control_file )
+{
+ int i, ret;
+ rvec x;
+ spuremd_handle *spmd_handle;
+
+ ret = SPUREMD_FAILURE;
+
+ if ( handle != NULL )
+ {
+ spmd_handle = (spuremd_handle*) handle;
+
+ /* close files used in previous simulation */
+ if ( spmd_handle->output_enabled == TRUE )
+ {
+ Finalize_Out_Controls( spmd_handle->system, spmd_handle->control,
+ spmd_handle->workspace, spmd_handle->out_control );
+ }
+
+ spmd_handle->realloc = FALSE;
+ spmd_handle->data->sim_id++;
+
+ spmd_handle->system->N = qm_num_atoms + mm_num_atoms;
+
+ PreAllocate_Space( spmd_handle->system, spmd_handle->control,
+ spmd_handle->workspace, spmd_handle->system->N );
+
+ Setup_Box( sim_box[0], sim_box[1], sim_box[2],
+ sim_box[3], sim_box[4], sim_box[5],
+ &spmd_handle->system->box );
+
+ for ( i = 0; i < qm_num_atoms; ++i )
+ {
+ x[0] = qm_pos_x[i];
+ x[1] = qm_pos_y[i];
+ x[2] = qm_pos_z[i];
+
+ Fit_to_Periodic_Box( &spmd_handle->system->box, x );
+
+ spmd_handle->workspace->orig_id[i] = i + 1;
+// spmd_handle->system->atoms[i].type = Get_Atom_Type( &system->reax_param,
+// element, sizeof(element) );
+ spmd_handle->system->atoms[i].type = qm_types[i];
+// strncpy( spmd_handle->system->atoms[i].name, atom_name,
+// sizeof(spmd_handle->system->atoms[i].name) - 1 );
+// spmd_handle->system->atoms[i].name[sizeof(spmd_handle->system->atoms[i].name) - 1] = '\0';
+ rvec_Copy( spmd_handle->system->atoms[i].x, x );
+ rvec_MakeZero( spmd_handle->system->atoms[i].v );
+ rvec_MakeZero( spmd_handle->system->atoms[i].f );
+ spmd_handle->system->atoms[i].q = 0.0;
+
+// mask[i] = 1;
+ }
+
+ for ( i = qm_num_atoms; i < qm_num_atoms + mm_num_atoms; ++i )
+ {
+ x[0] = mm_pos_x[i - qm_num_atoms];
+ x[1] = mm_pos_y[i - qm_num_atoms];
+ x[2] = mm_pos_z[i - qm_num_atoms];
+
+ Fit_to_Periodic_Box( &spmd_handle->system->box, x );
+
+ spmd_handle->workspace->orig_id[i] = i + 1;
+// spmd_handle->system->atoms[i].type = Get_Atom_Type( &system->reax_param,
+// element, sizeof(element) );
+ spmd_handle->system->atoms[i].type = mm_types[i - qm_num_atoms];
+// strncpy( spmd_handle->system->atoms[i].name, atom_name,
+// sizeof(spmd_handle->system->atoms[i].name) - 1 );
+// spmd_handle->system->atoms[i].name[sizeof(spmd_handle->system->atoms[i].name) - 1] = '\0';
+ rvec_Copy( spmd_handle->system->atoms[i].x, x );
+ rvec_MakeZero( spmd_handle->system->atoms[i].v );
+ rvec_MakeZero( spmd_handle->system->atoms[i].f );
+ spmd_handle->system->atoms[i].q = mm_q[i - qm_num_atoms];
+
+// mask[i] = 0;
+ }
+
+ Read_Input_Files( NULL, ffield_file, control_file,
+ spmd_handle->system, spmd_handle->control,
+ spmd_handle->data, spmd_handle->workspace,
+ spmd_handle->out_control );
+
+ if ( spmd_handle->system->N > spmd_handle->system->N_max )
+ {
+ /* deallocate everything which needs more space
+ * (i.e., structures whose space is a function of the number of atoms),
+ * except for data structures allocated while parsing input files */
+ Finalize( spmd_handle->system, spmd_handle->control, spmd_handle->data,
+ spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control,
+ spmd_handle->output_enabled, TRUE );
+
+ spmd_handle->system->N_max = (int) CEIL( SAFE_ZONE * spmd_handle->system->N );
+ spmd_handle->realloc = TRUE;
+ }
+
+ ret = SPUREMD_SUCCESS;
+ }
+
+ return ret;
+}
+
+
/* Reset for the next simulation by parsing input files and triggering
* reallocation if more space is needed
*
@@ -434,6 +719,8 @@ int cleanup( const void * const handle )
* geo_file: file containing geometry info of the structure to simulate
* ffield_file: file containing force field parameters
* control_file: file containing simulation parameters
+ *
+ * returns: SPUREMD_SUCCESS upon success, SPUREMD_FAILURE otherwise
*/
int reset( const void * const handle, const char * const geo_file,
const char * const ffield_file, const char * const control_file )
@@ -485,9 +772,11 @@ int reset( const void * const handle, const char * const geo_file,
/* Getter for atom positions
*
* handle: pointer to wrapper struct with top-level data structures
- * pos_x: x-coordinate of atom positions (allocated by caller)
- * pos_y: y-coordinate of atom positions (allocated by caller)
- * pos_z: z-coordinate of atom positions (allocated by caller)
+ * pos_x: x-coordinate of atom positions, in Angstroms (allocated by caller)
+ * pos_y: y-coordinate of atom positions, in Angstroms (allocated by caller)
+ * pos_z: z-coordinate of atom positions, in Angstroms (allocated by caller)
+ *
+ * returns: SPUREMD_SUCCESS upon success, SPUREMD_FAILURE otherwise
*/
int get_atom_positions( const void * const handle, double * const pos_x,
double * const pos_y, double * const pos_z )
@@ -515,10 +804,82 @@ int get_atom_positions( const void * const handle, double * const pos_x,
}
+/* Getter for atom velocities
+ *
+ * handle: pointer to wrapper struct with top-level data structures
+ * vel_x: x-coordinate of atom velocities, in Angstroms / ps (allocated by caller)
+ * vel_y: y-coordinate of atom velocities, in Angstroms / ps (allocated by caller)
+ * vel_z: z-coordinate of atom velocities, in Angstroms / ps (allocated by caller)
+ *
+ * returns: SPUREMD_SUCCESS upon success, SPUREMD_FAILURE otherwise
+ */
+int get_atom_velocities( const void * const handle, double * const vel_x,
+ double * const vel_y, double * const vel_z )
+{
+ int i, ret;
+ spuremd_handle *spmd_handle;
+
+ ret = SPUREMD_FAILURE;
+
+ if ( handle != NULL )
+ {
+ spmd_handle = (spuremd_handle*) handle;
+
+ for ( i = 0; i < spmd_handle->system->N; ++i )
+ {
+ vel_x[i] = spmd_handle->system->atoms[i].v[0];
+ vel_y[i] = spmd_handle->system->atoms[i].v[1];
+ vel_z[i] = spmd_handle->system->atoms[i].v[2];
+ }
+
+ ret = SPUREMD_SUCCESS;
+ }
+
+ return ret;
+}
+
+
+/* Getter for atom forces
+ *
+ * handle: pointer to wrapper struct with top-level data structures
+ * f_x: x-coordinate of atom forces, in Angstroms * Daltons / ps^2 (allocated by caller)
+ * f_y: y-coordinate of atom forces, in Angstroms * Daltons / ps^2 (allocated by caller)
+ * f_z: z-coordinate of atom forces, in Angstroms * Daltons / ps^2 (allocated by caller)
+ *
+ * returns: SPUREMD_SUCCESS upon success, SPUREMD_FAILURE otherwise
+ */
+int get_atom_forces( const void * const handle, double * const f_x,
+ double * const f_y, double * const f_z )
+{
+ int i, ret;
+ spuremd_handle *spmd_handle;
+
+ ret = SPUREMD_FAILURE;
+
+ if ( handle != NULL )
+ {
+ spmd_handle = (spuremd_handle*) handle;
+
+ for ( i = 0; i < spmd_handle->system->N; ++i )
+ {
+ f_x[i] = spmd_handle->system->atoms[i].f[0];
+ f_y[i] = spmd_handle->system->atoms[i].f[1];
+ f_z[i] = spmd_handle->system->atoms[i].f[2];
+ }
+
+ ret = SPUREMD_SUCCESS;
+ }
+
+ return ret;
+}
+
+
/* Getter for atom charges
*
* handle: pointer to wrapper struct with top-level data structures
- * q: atom charges (allocated by caller)
+ * q: atom charges, in Coulombs (allocated by caller)
+ *
+ * returns: SPUREMD_SUCCESS upon success, SPUREMD_FAILURE otherwise
*/
int get_atom_charges( const void * const handle, double * const q )
{
@@ -543,10 +904,52 @@ int get_atom_charges( const void * const handle, double * const q )
}
+/* Getter for system energies
+ *
+ * handle: pointer to wrapper struct with top-level data structures
+ * e_pot: system potential energy, in kcal / mol (reference from caller)
+ * e_kin: system kinetic energy, in kcal / mol (reference from caller)
+ * e_tot: system total energy, in kcal / mol (reference from caller)
+ * t_scalar: temperature scalar, in K (reference from caller)
+ * vol: volume of the simulation box, in Angstroms^3 (reference from caller)
+ * pres: average pressure, in K (reference from caller)
+ *
+ * returns: SPUREMD_SUCCESS upon success, SPUREMD_FAILURE otherwise
+ */
+int get_system_info( const void * const handle, double * const e_pot,
+ double * const e_kin, double * const e_tot, double * const temp,
+ double * const vol, double * const pres )
+{
+ int ret;
+ spuremd_handle *spmd_handle;
+
+ ret = SPUREMD_FAILURE;
+
+ if ( handle != NULL )
+ {
+ spmd_handle = (spuremd_handle*) handle;
+
+ *e_pot = spmd_handle->data->E_Pot;
+ *e_kin = spmd_handle->data->E_Kin;
+ *e_tot = spmd_handle->data->E_Tot;
+ *temp = spmd_handle->data->therm.T;
+ *vol = spmd_handle->system->box.volume;
+ *pres = (spmd_handle->control->P[0] + spmd_handle->control->P[1]
+ + spmd_handle->control->P[2]) / 3.0;
+
+ ret = SPUREMD_SUCCESS;
+ }
+
+ return ret;
+}
+
+
/* Setter for writing output to files
*
* handle: pointer to wrapper struct with top-level data structures
* enabled: TRUE enables writing output to files, FALSE otherwise
+ *
+ * returns: SPUREMD_SUCCESS upon success, SPUREMD_FAILURE otherwise
*/
int set_output_enabled( const void * const handle, const int enabled )
{
@@ -564,3 +967,34 @@ int set_output_enabled( const void * const handle, const int enabled )
return ret;
}
+
+
+/* Setter for simulation parameter values as defined in the input control file
+ *
+ * handle: pointer to wrapper struct with top-level data structures
+ * control_keyword: keyword from the control file to set the value for
+ * control_value: value to set
+ *
+ * returns: SPUREMD_SUCCESS upon success, SPUREMD_FAILURE otherwise
+ */
+int set_control_parameter( const void * const handle, const char * const keyword,
+ const char ** const values )
+{
+ int ret, ret_;
+ spuremd_handle *spmd_handle;
+
+ ret = SPUREMD_FAILURE;
+
+ if ( handle != NULL )
+ {
+ spmd_handle = (spuremd_handle*) handle;
+ ret_ = Set_Control_Parameter( keyword, values, spmd_handle->control,
+ spmd_handle->out_control );
+ if ( ret_ == SUCCESS )
+ {
+ ret = SPUREMD_SUCCESS;
+ }
+ }
+
+ return ret;
+}
diff --git a/sPuReMD/src/spuremd.h b/sPuReMD/src/spuremd.h
index d84c1875..3ee91704 100644
--- a/sPuReMD/src/spuremd.h
+++ b/sPuReMD/src/spuremd.h
@@ -33,7 +33,15 @@
extern "C" {
#endif
-void* setup( const char * const, const char * const,
+void * setup_qmmm_( int, const int * const,
+ const double * const, const double * const,
+ const double * const, int, const int * const,
+ const double * const, const double * const,
+ const double * const, const double * const,
+ const double * const,
+ const char * const, const char * const );
+
+void * setup( const char * const, const char * const,
const char * const );
int setup_callback( const void * const, const callback_function );
@@ -42,16 +50,38 @@ int simulate( const void * const );
int cleanup( const void * const );
+int reset_qmmm_( const void * const, int,
+ const int * const,
+ const double * const, const double * const,
+ const double * const, int, const int * const,
+ const double * const, const double * const,
+ const double * const, const double * const,
+ const double * const,
+ const char * const, const char * const );
+
int reset( const void * const, const char * const,
const char * const, const char * const );
int get_atom_positions( const void * const, double * const,
double * const, double * const );
+int get_atom_velocities( const void * const, double * const,
+ double * const, double * const );
+
+int get_atom_forces( const void * const, double * const,
+ double * const, double * const );
+
int get_atom_charges( const void * const, double * const );
+int get_system_info( const void * const, double * const,
+ double * const, double * const, double * const,
+ double * const, double * const );
+
int set_output_enabled( const void * const, const int );
+int set_control_parameter( const void * const, const char * const,
+ const char ** const );
+
#if defined(__cplusplus)
}
#endif
diff --git a/sPuReMD/src/tool_box.c b/sPuReMD/src/tool_box.c
index b9ed9cbe..1766879c 100644
--- a/sPuReMD/src/tool_box.c
+++ b/sPuReMD/src/tool_box.c
@@ -21,9 +21,15 @@
#include "tool_box.h"
+#include <stdlib.h>
#include <ctype.h>
+#include <limits.h>
+#include <errno.h>
#include <time.h>
+/* base 10 for result of string-to-integer conversion */
+#define INTBASE (10)
+
/************** taken from box.c **************/
/* Applies transformation on atomic position between Cartesian and
@@ -553,3 +559,135 @@ void sfclose( FILE * fp, const char * msg )
exit( INVALID_INPUT );
}
}
+
+
+/* Safe wrapper around strtol
+ *
+ * str: string to be converted
+ * filename: source filename of caller
+ * line: source line of caller
+ *
+ * returns: result of conversion (integer)
+ * */
+int sstrtol( const char * const str,
+ const char * const filename, int line )
+{
+ long ret;
+ char *endptr;
+
+ if ( str[0] == '\0' )
+ {
+ fprintf( stderr, "[ERROR] sstrtol: NULL string\n" );
+ /* strlen safe here only if filename is NULL-terminated
+ * before calling sconvert_string_to_int */
+ fprintf( stderr, " [INFO] At line %d in file %.*s\n",
+ line, (int) strlen(filename), filename );
+ exit( INVALID_INPUT );
+ }
+
+ errno = 0;
+ ret = strtol( str, &endptr, INTBASE );
+
+ if ( (errno == ERANGE && (ret == LONG_MAX || ret == LONG_MIN) )
+ || (errno != 0 && ret == 0) )
+ {
+ fprintf( stderr, "[ERROR] strtol: invalid string\n" );
+ /* strlen safe here only if filename is NULL-terminated
+ * before calling sconvert_string_to_int */
+ fprintf( stderr, " [INFO] At line %d in file %.*s\n",
+ line, (int) strlen(filename), filename );
+ fprintf( stderr, " [INFO] str: %.*s\n",
+ (int) strlen(str), str );
+ exit( INVALID_INPUT );
+ }
+ else if ( endptr == str )
+ {
+ fprintf( stderr, "[ERROR] strtol: no digits found\n" );
+ /* strlen safe here only if filename is NULL-terminated
+ * before calling sconvert_string_to_int */
+ fprintf( stderr, " [INFO] At line %d in file %.*s\n",
+ line, (int) strlen(filename), filename );
+ fprintf( stderr, " [INFO] str: %.*s\n",
+ (int) strlen(str), str );
+ exit( INVALID_INPUT );
+ }
+ else if ( *endptr != '\0' )
+ {
+ fprintf( stderr, "[ERROR] strtol: non-numeric trailing characters\n" );
+ /* strlen safe here only if filename is NULL-terminated
+ * before calling sconvert_string_to_int */
+ fprintf( stderr, " [INFO] At line %d in file %.*s\n",
+ line, (int) strlen(filename), filename );
+ fprintf( stderr, " [INFO] str: %.*s\n",
+ (int) strlen(str), str );
+ exit( INVALID_INPUT );
+ }
+
+ return (int) ret;
+}
+
+
+/* Safe wrapper around strtod
+ *
+ * str: string to be converted
+ * filename: source filename of caller
+ * line: source line of caller
+ *
+ * returns: result of conversion (double)
+ * */
+double sstrtod( const char * const str,
+ const char * const filename, int line )
+{
+ double ret;
+ char *endptr;
+
+ if ( str[0] == '\0' )
+ {
+ fprintf( stderr, "[ERROR] sstrtod: NULL string\n" );
+ /* strlen safe here only if filename is NULL-terminated
+ * before calling sconvert_string_to_int */
+ fprintf( stderr, " [INFO] At line %d in file %.*s\n",
+ line, (int) strlen(filename), filename );
+ exit( INVALID_INPUT );
+ }
+
+ errno = 0;
+ ret = strtod( str, &endptr );
+
+ if ( (errno == ERANGE && (ret == LONG_MAX || ret == LONG_MIN) )
+ || (errno != 0 && ret == 0) )
+ {
+ fprintf( stderr, "[ERROR] strtod: invalid string\n" );
+ /* strlen safe here only if filename is NULL-terminated
+ * before calling sconvert_string_to_int */
+ fprintf( stderr, " [INFO] At line %d in file %.*s\n",
+ line, (int) strlen(filename), filename );
+ fprintf( stderr, " [INFO] str: %.*s\n",
+ (int) strlen(str), str );
+ exit( INVALID_INPUT );
+ }
+ else if ( endptr == str )
+ {
+ fprintf( stderr, "[ERROR] strtod: no digits found\n" );
+ /* strlen safe here only if filename is NULL-terminated
+ * before calling sconvert_string_to_int */
+ fprintf( stderr, " [INFO] At line %d in file %.*s\n",
+ line, (int) strlen(filename), filename );
+ fprintf( stderr, " [INFO] str: %.*s\n",
+ (int) strlen(str), str );
+ exit( INVALID_INPUT );
+ }
+ else if ( *endptr != '\0' )
+ {
+ fprintf( stderr, "[ERROR] strtod: non-numeric trailing characters\n" );
+ /* strlen safe here only if filename is NULL-terminated
+ * before calling sconvert_string_to_int */
+ fprintf( stderr, " [INFO] At line %d in file %.*s\n",
+ line, (int) strlen(filename), filename );
+ fprintf( stderr, " [INFO] str: %.*s\n",
+ (int) strlen(str), str );
+ exit( INVALID_INPUT );
+ }
+
+ return ret;
+}
diff --git a/sPuReMD/src/tool_box.h b/sPuReMD/src/tool_box.h
index a76bd32d..2c232b1b 100644
--- a/sPuReMD/src/tool_box.h
+++ b/sPuReMD/src/tool_box.h
@@ -26,20 +26,20 @@
/* from box.h */
-void Transform( rvec, simulation_box*, int, rvec );
+void Transform( rvec, simulation_box *, int, rvec );
-void Transform_to_UnitBox( rvec, simulation_box*, int, rvec );
+void Transform_to_UnitBox( rvec, simulation_box *, int, rvec );
-void Fit_to_Periodic_Box( simulation_box*, rvec );
+void Fit_to_Periodic_Box( simulation_box *, rvec );
-int is_Inside_Box( simulation_box*, rvec );
+int is_Inside_Box( simulation_box *, rvec );
/* from geo_tools.h */
-void Make_Point( real, real, real, rvec* );
+void Make_Point( real, real, real, rvec * );
-int is_Valid_Serial( static_storage*, int );
+int is_Valid_Serial( static_storage *, int );
-int Check_Input_Range( int, int, int, char* );
+int Check_Input_Range( int, int, int, char * );
void Trim_Spaces( char * const, const size_t );
@@ -49,26 +49,26 @@ real Get_Time( );
real Get_Timing_Info( real );
/* from io_tools.h */
-int Get_Atom_Type( reax_interaction*, char*, size_t );
+int Get_Atom_Type( reax_interaction *, char *, size_t );
-char *Get_Element( reax_system*, int );
+char * Get_Element( reax_system *, int );
-char *Get_Atom_Name( reax_system*, int );
+char * Get_Atom_Name( reax_system *, int );
-void Allocate_Tokenizer_Space( char**, size_t, char**, size_t, char***,
+void Allocate_Tokenizer_Space( char **, size_t, char **, size_t, char ***,
size_t, size_t );
void Deallocate_Tokenizer_Space( char **, char **, char ***,
size_t );
-int Tokenize( char*, char***, size_t );
+int Tokenize( char *, char ***, size_t );
/* from lammps */
-void *smalloc( size_t, const char * );
+void * smalloc( size_t, const char * );
-void* srealloc( void *, size_t, const char * );
+void * srealloc( void *, size_t, const char * );
-void* scalloc( size_t, size_t, const char * );
+void * scalloc( size_t, size_t, const char * );
void sfree( void *, const char * );
@@ -76,4 +76,10 @@ FILE * sfopen( const char *, const char * );
void sfclose( FILE *, const char * );
+int sstrtol( const char * const,
+ const char * const, int );
+
+double sstrtod( const char * const,
+ const char * const, int );
+
#endif
diff --git a/sPuReMD/src/vector.h b/sPuReMD/src/vector.h
index f713804a..2f3dd52e 100644
--- a/sPuReMD/src/vector.h
+++ b/sPuReMD/src/vector.h
@@ -73,7 +73,7 @@ static inline void Vector_MakeZero( real * const v, const unsigned int k )
#endif
for ( i = 0; i < k; ++i )
{
- v[i] = ZERO;
+ v[i] = 0.0;
}
}
@@ -146,7 +146,7 @@ static inline real Dot( const real * const v1, const real * const v2,
#pragma omp single
#endif
{
- ret2_omp = ZERO;
+ ret2_omp = 0.0;
}
#if defined(_OPENMP)
@@ -169,7 +169,7 @@ static inline real Norm( const real * const v1, const unsigned int k )
#pragma omp single
#endif
{
- ret2_omp = ZERO;
+ ret2_omp = 0.0;
}
#if defined(_OPENMP)
@@ -478,7 +478,7 @@ static inline void rvec_MakeZero( rvec v )
#endif
for ( i = 0; i < 3; ++i )
{
- v[i] = ZERO;
+ v[i] = 0.0;
}
}
@@ -662,7 +662,7 @@ static inline void rtensor_Identity( rtensor t )
{
for ( j = 0; j < 3; ++j )
{
- t[i][j] = (i == j ? 1.0 : ZERO);
+ t[i][j] = (i == j ? 1.0 : 0.0);
}
}
}
@@ -676,7 +676,7 @@ static inline void rtensor_MakeZero( rtensor t )
{
for ( j = 0; j < 3; ++j )
{
- t[i][j] = ZERO;
+ t[i][j] = 0.0;
}
}
}
diff --git a/tools/driver_qmmm.py b/tools/driver_qmmm.py
new file mode 100644
index 00000000..ec315724
--- /dev/null
+++ b/tools/driver_qmmm.py
@@ -0,0 +1,524 @@
+#!/bin/python3
+
+from ctypes import c_int, c_double, c_char, c_char_p, c_void_p, \
+ Structure, Union, POINTER, CFUNCTYPE, cdll
+import sqlite3 as sq3
+from os import path
+
+
+class BondOrderData(Structure):
+ _fields_ = [
+ ("BO", c_double),
+ ("BO_s", c_double),
+ ("BO_pi", c_double),
+ ("BO_pi2", c_double),
+ ("Cdbo", c_double),
+ ("Cdbopi", c_double),
+ ("Cdbopi2", c_double),
+ ("C1dbo", c_double),
+ ("C2dbo", c_double),
+ ("C3dbo", c_double),
+ ("C1dbopi", c_double),
+ ("C2dbopi", c_double),
+ ("C3dbopi", c_double),
+ ("C4dbopi", c_double),
+ ("C1dbopi2", c_double),
+ ("C2dbopi2", c_double),
+ ("C3dbopi2", c_double),
+ ("C4dbopi2", c_double),
+ ("dBOp", c_double * 3),
+ ("dln_BOp_s", c_double * 3),
+ ("dln_BOp_pi", c_double * 3),
+ ("dln_BOp_pi2", c_double * 3),
+ ]
+
+
+class ThreeBodyData(Structure):
+ _fields_ = [
+ ("thb", c_int),
+ ("pthb", c_int),
+ ("theta", c_double),
+ ("cos_theta", c_double),
+ ("dcos_di", c_double * 3),
+ ("dcos_dj", c_double * 3),
+ ("dcos_dk", c_double * 3),
+ ]
+
+
+class BondData(Structure):
+ _fields_ = [
+ ("nbr", c_int),
+ ("sym_index", c_int),
+ ("dbond_index", c_int),
+ ("rel_box", c_int * 3),
+ ("d", c_double),
+ ("dvec", c_double * 3),
+ ("bo_data", BondOrderData),
+ ]
+
+
+class DBondData(Structure):
+ _fields_ = [
+ ("wrt", c_int),
+ ("dBO", c_double * 3),
+ ("dBOpi", c_double * 3),
+ ("dBOpi2", c_double * 3),
+ ]
+
+
+class DDeltaData(Structure):
+ _fields_ = [
+ ("wrt", c_int),
+ ("dVal", c_double * 3),
+ ]
+
+
+class FarNbrData(Structure):
+ _fields_ = [
+ ("nbr", c_int),
+ ("rel_box", c_int * 3),
+ ("d", c_double),
+ ("dvec", c_double * 3),
+ ]
+
+
+class NearNbrData(Structure):
+ _fields_ = [
+ ("nbr", c_int),
+ ("rel_box", c_int * 3),
+ ("d", c_double),
+ ("dvec", c_double * 3),
+ ]
+
+
+class HBondData(Structure):
+ _fields_ = [
+ ("nbr", c_int),
+ ("scl", c_int),
+ ("ptr", POINTER(FarNbrData)),
+ ]
+
+
+class Thermostat(Structure):
+ _fields_ = [
+ ("T", c_double),
+ ("xi", c_double),
+ ("v_xi", c_double),
+ ("v_xi_old", c_double),
+ ("G_xi", c_double),
+ ]
+
+
+class IsotropicBarostat(Structure):
+ _fields_ = [
+ ("P", c_double),
+ ("eps", c_double),
+ ("v_eps", c_double),
+ ("v_eps_old", c_double),
+ ("a_eps", c_double),
+ ]
+
+
+class FlexibleBarostat(Structure):
+ _fields_ = [
+ ("P", c_double * 9),
+ ("P_scalar", c_double),
+ ("eps", c_double),
+ ("v_eps", c_double),
+ ("v_eps_old", c_double),
+ ("a_eps", c_double),
+ ("h0", c_double * 9),
+ ("v_g0", c_double * 9),
+ ("v_g0_old", c_double * 9),
+ ("a_g0", c_double * 9),
+ ]
+
+
+class ReaxTiming(Structure):
+ _fields_ = [
+ ("start", c_double),
+ ("end", c_double),
+ ("elapsed", c_double),
+ ("total", c_double),
+ ("nbrs", c_double),
+ ("init_forces", c_double),
+ ("bonded", c_double),
+ ("nonb", c_double),
+ ("cm", c_double),
+ ("cm_sort_mat_rows", c_double),
+ ("cm_solver_pre_comp", c_double),
+ ("cm_solver_pre_app", c_double),
+ ("cm_solver_iters", c_int),
+ ("cm_solver_spmv", c_double),
+ ("cm_solver_vector_ops", c_double),
+ ("cm_solver_orthog", c_double),
+ ("cm_solver_tri_solve", c_double),
+ ("cm_last_pre_comp", c_double),
+ ("cm_total_loss", c_double),
+ ("cm_optimum", c_double),
+ ("num_retries", c_int),
+ ]
+
+
+class SimulationData(Structure):
+ _fields_ = [
+ ("sim_id", c_int),
+ ("step", c_int),
+ ("prev_step", c_int),
+ ("time", c_double),
+ ("M", c_double),
+ ("inv_M", c_double),
+ ("xcm", c_double * 3),
+ ("vcm", c_double * 3),
+ ("fcm", c_double * 3),
+ ("amcm", c_double * 3),
+ ("avcm", c_double * 3),
+ ("etran_cm", c_double),
+ ("erot_cm", c_double),
+ ("kinetic", c_double * 9),
+ ("virial", c_double * 9),
+ ("E_Tot", c_double),
+ ("E_Kin", c_double),
+ ("E_Pot", c_double),
+ ("E_BE", c_double),
+ ("E_Ov", c_double),
+ ("E_Un", c_double),
+ ("E_Lp", c_double),
+ ("E_Ang", c_double),
+ ("E_Pen", c_double),
+ ("E_Coa", c_double),
+ ("E_HB", c_double),
+ ("E_Tor", c_double),
+ ("E_Con", c_double),
+ ("E_vdW", c_double),
+ ("E_Ele", c_double),
+ ("E_Pol", c_double),
+ ("N_f", c_double),
+ ("t_scale", c_double * 3),
+ ("p_scale", c_double * 9),
+ ("therm", Thermostat),
+ ("iso_bar", IsotropicBarostat),
+ ("flex_bar", FlexibleBarostat),
+ ("inv_W", c_double),
+ ("int_press", c_double * 3),
+ ("ext_press", c_double * 3),
+ ("kin_press", c_double),
+ ("tot_press", c_double * 3),
+ ("timing", ReaxTiming),
+ ]
+
+
+class ReaxAtom(Structure):
+ _fields_ = [
+ ("type", c_int),
+ ("rel_map", c_int * 3),
+ ("name", c_char * 9),
+ ("x", c_double * 3),
+ ("v", c_double * 3),
+ ("f", c_double * 3),
+ ("q", c_double),
+ ]
+
+
+def create_db(name='spuremd.db'):
+ conn = sq3.connect(name)
+
+ conn.executescript("""
+ CREATE TABLE simulation(
+ id integer,
+ date text,
+ name text,
+ ensemble_type integer,
+ steps integer,
+ time_step integer,
+ restart_format integer,
+ random_velocity integer,
+ reposition_atoms integer,
+ peroidic_boundary integer,
+ geo_format integer,
+ restrict_bonds integer,
+ tabulate_long_range integer,
+ reneighbor integer,
+ vlist_cutoff real,
+ neighbor_cutoff real,
+ three_body_cutoff real,
+ hydrogen_bond_cutoff real,
+ bond_graph_cutoff real,
+ charge_method integer,
+ cm_q_net real,
+ cm_solver_type integer,
+ cm_solver_max_iters integer,
+ cm_solver_restart integer,
+ cm_solver_q_err real,
+ cm_domain_sparsity real,
+ cm_solver_pre_comp_type integer,
+ cm_solver_pre_comp_refactor integer,
+ cm_solver_pre_comp_droptol real,
+ cm_solver_pre_comp_sweeps integer,
+ cm_solver_pre_comp_sai_thres real,
+ cm_solver_pre_app_type integer,
+ cm_solver_pre_app_jacobi_iters integer,
+ temp_init real,
+ temp_final real,
+ temp_mass real,
+ temp_mode integer,
+ temp_rate real,
+ temp_freq integer,
+ pressure real,
+ pressure_mass real,
+ compress integer,
+ pressure_mode integer,
+ remove_center_of_mass integer,
+ debug_level integer,
+ write_freq integer,
+ traj_compress integer,
+ traj_format integer,
+ traj_title text,
+ atom_info integer,
+ atom_velocities integer,
+ atom_forces integer,
+ bond_info integer,
+ angle_info integer,
+ test_forces integer,
+ molecule_analysis integer,
+ freq_molecule_analysis integer,
+ ignore integer,
+ dipole_analysis integer,
+ freq_dipole_analysis integer,
+ diffusion_coefficient integer,
+ freq_diffusion_coefficient integer,
+ restrict_type integer,
+ PRIMARY KEY (id)
+ );
+
+ CREATE TABLE system_properties(
+ id integer,
+ step integer,
+ total_energy real,
+ potential_energy real,
+ kinetic_energy real,
+ temperature real,
+ volume real,
+ pressure real,
+ PRIMARY KEY (id, step)
+ );
+
+ CREATE TABLE potential(
+ id integer,
+ step integer,
+ bond_energy real,
+ atom_energy real,
+ lone_pair_energy real,
+ angle_energy real,
+ coa_energy real,
+ hydrogen_bond_energy real,
+ torsion_energy real,
+ conjugation_energy real,
+ van_der_waals_energy real,
+ coulombic_energy real,
+ polarization_energy real,
+ PRIMARY KEY (id, step)
+ );
+
+ CREATE TABLE trajectory(
+ id integer,
+ step integer,
+ atom_id integer,
+ position_x real,
+ position_y real,
+ position_z real,
+ charge real,
+ PRIMARY KEY (id, step, atom_id)
+ );
+
+ CREATE TABLE performance(
+ id integer,
+ step integer,
+ time_total real,
+ time_nbrs real,
+ time_init real,
+ time_bonded real,
+ time_nonbonded real,
+ time_cm real,
+ time_cm_sort real,
+ cm_solver_iters integer,
+ time_cm_pre_comp real,
+ time_cm_pre_app real,
+ time_cm_solver_spmv real,
+ time_cm_solver_vec_ops real,
+ time_cm_solver_orthog real,
+ time_cm_solver_tri_solve real,
+ PRIMARY KEY (id, step)
+ );
+ """)
+
+ conn.close()
+
+
+if __name__ == '__main__':
+ lib = cdll.LoadLibrary("libspuremd.so.1")
+
+ setup_qmmm = lib.setup_qmmm_
+ setup_qmmm.argtypes = [c_int, POINTER(c_int),
+ POINTER(c_double), POINTER(c_double), POINTER(c_double),
+ c_int, POINTER(c_int), POINTER(c_double), POINTER(c_double),
+ POINTER(c_double), POINTER(c_double), POINTER(c_double),
+ c_char_p, c_char_p]
+ setup_qmmm.restype = c_void_p
+
+ simulate = lib.simulate
+ simulate.argtypes = [c_void_p]
+ simulate.restype = c_int
+
+ cleanup = lib.cleanup
+ cleanup.argtypes = [c_void_p]
+ cleanup.restype = c_int
+
+ reset_qmmm = lib.reset_qmmm_
+ reset_qmmm.argtypes = [c_void_p, c_int, POINTER(c_int),
+ POINTER(c_double), POINTER(c_double), POINTER(c_double),
+ c_int, POINTER(c_int), POINTER(c_double), POINTER(c_double),
+ POINTER(c_double), POINTER(c_double), POINTER(c_double),
+ c_char_p, c_char_p]
+ reset_qmmm.restype = c_int
+
+ CALLBACKFUNC = CFUNCTYPE(None, c_int, POINTER(ReaxAtom),
+ POINTER(SimulationData))
+
+ setup_callback = lib.setup_callback
+ setup_callback.argtypes = [c_void_p, CALLBACKFUNC]
+ setup_callback.restype = c_int
+
+ set_control_parameter = lib.set_control_parameter
+ set_control_parameter.argtypes = [c_void_p, c_char_p, POINTER(c_char_p)]
+ set_control_parameter.restype = c_int
+
+ get_atom_positions = lib.get_atom_positions
+ get_atom_positions.argtypes = [c_void_p, POINTER(c_double),
+ POINTER(c_double), POINTER(c_double)]
+ get_atom_positions.restype = c_int
+
+ get_atom_forces = lib.get_atom_forces
+ get_atom_forces.argtypes = [c_void_p, POINTER(c_double),
+ POINTER(c_double), POINTER(c_double)]
+ get_atom_forces.restype = c_int
+
+ get_atom_charges = lib.get_atom_charges
+ get_atom_charges.argtypes = [c_void_p, POINTER(c_double)]
+ get_atom_charges.restype = c_int
+
+ def get_simulation_step_results(num_atoms, atoms, data):
+ print("{0:24.15f} {1:24.15f} {2:24.15f}".format(
+ data[0].E_Tot, data[0].E_Kin, data[0].E_Pot))
+
+ # bulk water
+ sim_box = (c_double * 6)(40.299, 40.299, 40.299, 90.0, 90.0, 90.0)
+ num_qm_atoms = 10
+ num_mm_atoms = 10
+ num_atoms = num_qm_atoms + num_mm_atoms
+ qm_types = (c_int * num_qm_atoms)(2, 1, 1, 2, 1, 1, 2, 1, 1, 2)
+ qm_p_x = (c_double * num_qm_atoms)(5.690, 4.760, 5.800, 15.551, 14.981, 14.961, 17.431, 17.761, 17.941, 11.351)
+ qm_p_y = (c_double * num_qm_atoms)(12.751, 12.681, 13.641, 15.111, 14.951, 15.211, 6.180, 7.120, 5.640, 7.030)
+ qm_p_z = (c_double * num_qm_atoms)(11.651, 11.281, 12.091, 7.030, 7.840, 6.230, 8.560, 8.560, 9.220, 7.170)
+ mm_types = (c_int * num_qm_atoms)(1, 1, 2, 1, 1, 2, 1, 1, 2, 1)
+ mm_p_x = (c_double * num_mm_atoms)(11.921, 10.751, 17.551, 17.431, 17.251, 7.680, 6.900, 8.020, 8.500, 8.460)
+ mm_p_y = (c_double * num_mm_atoms)(7.810, 7.290, 6.070, 5.940, 5.260, 11.441, 11.611, 12.311, 7.980, 8.740)
+ mm_p_z = (c_double * num_mm_atoms)(6.920, 7.930, 2.310, 1.320, 2.800, 10.231, 10.831, 9.871, 18.231, 18.881)
+ mm_q = (c_double * num_mm_atoms)(-2.0, 1.0, 1.0, -2.0, 1.0, 1.0, -2.0, 1.0, 1.0, -2.0)
+
+ handle = setup_qmmm(c_int(num_qm_atoms), qm_types, qm_p_x, qm_p_y, qm_p_z,
+ c_int(num_mm_atoms), mm_types, mm_p_x, mm_p_y, mm_p_z, mm_q, sim_box,
+ b"data/benchmarks/water/ffield.water",
+ b"environ/control_water")
+
+ ret = setup_callback(handle, CALLBACKFUNC(get_simulation_step_results))
+
+ if ret != 0:
+ print("[ERROR] setup_callback returned {0}".format(ret))
+
+ keyword = b"nsteps"
+ values = (c_char_p)(b"10")
+ ret = set_control_parameter(handle, keyword, values)
+
+ if ret != 0:
+ print("[ERROR] set_control_parameter returned {0}".format(ret))
+
+ print("{0:24}|{1:24}|{2:24}".format("Total Energy", "Kinetic Energy", "Potential Energy"))
+
+ ret = simulate(handle)
+
+ if ret != 0:
+ print("[ERROR] simulate returned {0}".format(ret))
+
+ p_x = (c_double * num_atoms)()
+ p_y = (c_double * num_atoms)()
+ p_z = (c_double * num_atoms)()
+ ret = get_atom_positions(handle, p_x, p_y, p_z)
+
+ if ret != 0:
+ print("[ERROR] get_atom_positions returned {0}".format(ret))
+
+ f_x = (c_double * num_atoms)()
+ f_y = (c_double * num_atoms)()
+ f_z = (c_double * num_atoms)()
+ ret = get_atom_forces(handle, f_x, f_y, f_z)
+
+ if ret != 0:
+ print("[ERROR] get_atom_forces returned {0}".format(ret))
+
+ q = (c_double * num_atoms)()
+ ret = get_atom_charges(handle, q)
+
+ if ret != 0:
+ print("[ERROR] get_atom_charges returned {0}".format(ret))
+
+ # silica
+ sim_box = (c_double * 6)(36.477, 50.174, 52.110, 90.0, 90.0, 90.0)
+ num_qm_atoms = 15
+ num_mm_atoms = 15
+ num_atoms = num_qm_atoms + num_mm_atoms
+ qm_types = (c_int * num_qm_atoms)(2, 2, 2, 2, 10, 2, 10, 2, 2, 10, 2, 10, 2, 2)
+ qm_p_x = (c_double * num_qm_atoms)(56.987, 32.795, 26.543, 27.616, 26.560, 54.035, 54.425, 29.979, 38.008, 48.769, 57.113, 26.458, 52.299, 55.789, 45.752)
+ qm_p_y = (c_double * num_qm_atoms)(39.868, 24.104, 26.261, 27.534, 39.146, 39.112, 37.117, 43.558, 47.170, 42.454, 35.565, 35.477, 39.113, 41.444, 44.237)
+ qm_p_z = (c_double * num_qm_atoms)(41.795, 25.968, 36.254, 24.459, 32.281, 23.745, 32.278, 21.696, 27.275, 20.461, 31.366, 23.522, 26.519, 30.466, 10.521)
+ mm_types = (c_int * num_qm_atoms)(10, 2, 10, 2, 2, 2, 10, 2, 10, 2, 2, 10, 2, 2, 10)
+ mm_p_x = (c_double * num_mm_atoms)(49.617, 26.736, 58.166, 36.655, 55.773, 53.905, 52.254, 32.196, 53.817, 56.164, 28.678, 58.824, 24.807, 51.299, 28.950)
+ mm_p_y = (c_double * num_mm_atoms)(42.379, 37.815, 5.488, 48.615, 36.249, 1.450, 35.399, 6.595, 8.586, 1.585, 7.208, 36.067, 6.365, 25.234, 11.190)
+ mm_p_z = (c_double * num_mm_atoms)(21.790, 33.119, 0.945, 25.664, 31.954, 22.720, 33.665, 23.147, 31.684, 23.550, 30.801, 40.036, 24.849, 29.146, 24.542)
+ mm_q = (c_double * num_mm_atoms)(-1.0, -2.0, -1.0, -2.0, -2.0, -2.0, -1.0, -2.0, -1.0, -2.0, -2.0, -1.0, -2.0, -2.0, -1.0)
+
+ ret = reset_qmmm(handle, c_int(num_qm_atoms), qm_types, qm_p_x, qm_p_y, qm_p_z,
+ c_int(num_mm_atoms), mm_types, mm_p_x, mm_p_y, mm_p_z, mm_q, sim_box,
+ b"data/benchmarks/silica/ffield-bio",
+ b"environ/control_silica")
+
+ print("\n{0:24}|{1:24}|{2:24}".format("Total Energy", "Kinetic Energy", "Potential Energy"))
+
+ ret = simulate(handle)
+
+ if ret != 0:
+ print("[ERROR] simulate returned {0}".format(ret))
+
+ p_x = (c_double * num_atoms)()
+ p_y = (c_double * num_atoms)()
+ p_z = (c_double * num_atoms)()
+ ret = get_atom_positions(handle, p_x, p_y, p_z)
+
+ if ret != 0:
+ print("[ERROR] get_atom_positions returned {0}".format(ret))
+
+ f_x = (c_double * num_atoms)()
+ f_y = (c_double * num_atoms)()
+ f_z = (c_double * num_atoms)()
+ ret = get_atom_forces(handle, f_x, f_y, f_z)
+
+ if ret != 0:
+ print("[ERROR] get_atom_forces returned {0}".format(ret))
+
+ q = (c_double * num_atoms)()
+ ret = get_atom_charges(handle, q)
+
+ if ret != 0:
+ print("[ERROR] get_atom_charges returned {0}".format(ret))
+
+ cleanup(handle)
--
GitLab