tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

remove support for writing "native" trajectory files from USER-REAXC

Browse Source
This commit is contained in:
Axel Kohlmeyer
2021-04-18 05:16:41 -04:00
parent ab8d78c8f4
commit 6c88baceb7
17 changed files with 61 additions and 845 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
src/.gitignore vendored
View File

@ -1199,11 +1199,9 @@
/reaxc_bonds.cpp
/reaxc_control.cpp
/reaxc_ffield.cpp
/reaxc_ffield.h
/reaxc_forces.cpp
/reaxc_hydrogen_bonds.cpp
/reaxc_init_md.cpp
/reaxc_io_tools.cpp
/reaxc_list.cpp
/reaxc_lookup.cpp
/reaxc_multi_body.cpp
@ -1211,11 +1209,11 @@
/reaxc_reset_tools.cpp
/reaxc_tool_box.cpp
/reaxc_torsion_angles.cpp
/reaxc_traj.cpp
/reaxc_valence_angles.cpp
/reaxff_api.h
/reaxff_defs.h
/reaxff_inline.h
/reaxff_omp.h
/reaxff_types.h
/remap.cpp
/remap.h

2
src/Purge.list
View File

@ -62,6 +62,7 @@ reaxc_forces.h
reaxc_hydrogen_bonds.h
reaxc_init_md.h
reaxc_io_tools.h
reaxc_io_tools.cpp
reaxc_list.h
reaxc_lookup.h
reaxc_multi_body.h
@ -72,6 +73,7 @@ reaxc_system_props.cpp
reaxc_tool_box.h
reaxc_torsion_angles.h
reaxc_traj.h
reaxc_traj.cpp
reaxc_types.h
reaxc_valence_angles.h
reaxc_vector.h

8
src/USER-OMP/pair_reaxc_omp.cpp
View File

@ -211,7 +211,7 @@ void PairReaxCOMP::compute(int eflag, int vflag)
startTimeBase = MPI_Wtime();
#endif
Compute_ForcesOMP(api->system,api->control,api->data,api->workspace,&api->lists,api->out_control);
Compute_ForcesOMP(api->system,api->control,api->data,api->workspace,&api->lists);
read_reax_forces(vflag);
#ifdef OMP_TIMING
@ -270,8 +270,6 @@ void PairReaxCOMP::compute(int eflag, int vflag)
api->data->step = update->ntimestep;
Output_Results(api->system, api->control, api->data, &api->lists, api->out_control, world);
// populate tmpid and tmpbo arrays for fix reax/c/species
if (fixspecies_flag) {
@ -393,8 +391,8 @@ void PairReaxCOMP::setup()
write_reax_lists();
InitializeOMP(api->system, api->control, api->data, api->workspace,
&api->lists, api->out_control, world);
InitializeOMP(api->system, api->control, api->data,
api->workspace, &api->lists, world);
for (int k = 0; k < api->system->N; ++k) {
num_bonds[k] = api->system->my_atoms[k].num_bonds;

3
src/USER-OMP/reaxc_bond_orders_omp.cpp
View File

@ -366,8 +366,7 @@ namespace ReaxFF {
/* ---------------------------------------------------------------------- */
void BOOMP( reax_system *system, control_params * /* control */, simulation_data * /* data */,
storage *workspace, reax_list **lists, output_controls * /* out_control */)
void BOOMP( reax_system *system, storage *workspace, reax_list **lists)
{
double p_lp1 = system->reax_param.gp.l[15];
double p_boc1 = system->reax_param.gp.l[0];

9
src/USER-OMP/reaxc_forces_omp.cpp
View File

@ -43,10 +43,10 @@ namespace ReaxFF {
static void Compute_Bonded_ForcesOMP(reax_system *system, control_params *control,
simulation_data *data, storage *workspace,
reax_list **lists, output_controls *out_control)
reax_list **lists)
{
BOOMP(system, control, data, workspace, lists, out_control);
BOOMP(system, workspace, lists);
BondsOMP(system, data, workspace, lists);
Atom_EnergyOMP(system, data, workspace, lists);
Valence_AnglesOMP(system, control, data, workspace, lists);
@ -476,14 +476,13 @@ namespace ReaxFF {
void Compute_ForcesOMP(reax_system *system, control_params *control,
simulation_data *data, storage *workspace,
reax_list **lists, output_controls *out_control)
reax_list **lists)
{
// Init Forces
Init_Forces_noQEq_OMP(system, control, data, workspace, lists);
// Bonded Interactions
Compute_Bonded_ForcesOMP(system, control, data, workspace,
lists, out_control);
Compute_Bonded_ForcesOMP(system, control, data, workspace, lists);
// Nonbonded Interactions
Compute_NonBonded_ForcesOMP(system, control, data, workspace, lists);

5
src/USER-OMP/reaxc_init_md_omp.cpp
View File

@ -94,16 +94,13 @@ namespace ReaxFF {
// The only difference with the MPI-only function is calls to Init_ListsOMP and Init_Force_FunctionsOMP().
void InitializeOMP(reax_system *system, control_params *control,
simulation_data *data, storage *workspace,
reax_list **lists, output_controls *out_control,
MPI_Comm world)
reax_list **lists, MPI_Comm world)
{
Init_System(system,control);
Init_Simulation_Data(data);
Init_Workspace(system,control,workspace);
Init_ListsOMP(system,control,lists);
Init_Output_Files(system,control,out_control,world);
if (control->tabulate)
Init_Lookup_Tables(system,control,workspace,world);
}

11
src/USER-OMP/reaxff_omp.h
View File

@ -40,8 +40,7 @@ namespace ReaxFF
two_body_parameters *, int, double, double, double,
double, double, double, double);
extern void BOOMP(reax_system *, control_params *, simulation_data *,
storage *, reax_list **, output_controls *);
extern void BOOMP(reax_system *, storage *, reax_list **);
// bonds OpenMP
@ -51,8 +50,7 @@ namespace ReaxFF
// forces OpenMP
extern void Compute_ForcesOMP(reax_system *, control_params *,
simulation_data *, storage *,
reax_list **, output_controls *);
simulation_data *, storage *, reax_list **);
// hydrogen bonds
@ -61,9 +59,8 @@ namespace ReaxFF
// init md OpenMP
extern void InitializeOMP(reax_system *, control_params *,
simulation_data *, storage *, reax_list **,
output_controls *, MPI_Comm);
extern void InitializeOMP(reax_system *, control_params *, simulation_data *,
storage *, reax_list **,MPI_Comm);
// multi body

25
src/USER-REAXC/pair_reaxc.cpp
View File

@ -77,8 +77,6 @@ PairReaxC::PairReaxC(LAMMPS *lmp) : Pair(lmp)
memset(api->system,0,sizeof(reax_system));
api->control = new control_params;
memset(api->control,0,sizeof(control_params));
api->out_control = new output_controls;
memset(api->out_control,0,sizeof(output_controls));
api->data = new simulation_data;
api->workspace = new storage;
memory->create(api->lists, LIST_N,"reaxff:lists");
@ -124,7 +122,6 @@ PairReaxC::~PairReaxC()
delete[] fix_id;
if (setup_flag) {
Close_Output_Files(api->system,api->out_control);
// deallocate reax data-structures
@ -141,7 +138,6 @@ PairReaxC::~PairReaxC()
delete api->system;
delete api->control;
delete api->out_control;
delete api->data;
delete api->workspace;
memory->destroy(api->lists);
@ -192,8 +188,6 @@ void PairReaxC::settings(int narg, char **arg)
// read name of control file or use default controls
if (strcmp(arg[0],"NULL") == 0) {
strcpy(api->control->sim_name, "simulate");
api->out_control->energy_update_freq = 0;
api->control->tabulate = 0;
api->control->bond_cut = 5.;
@ -204,15 +198,9 @@ void PairReaxC::settings(int narg, char **arg)
api->control->nthreads = 1;
api->out_control->write_steps = 0;
strcpy(api->out_control->traj_title, "default_title");
api->out_control->atom_info = 0;
api->out_control->bond_info = 0;
api->out_control->angle_info = 0;
} else Read_Control_File(arg[0], api->control, api->out_control);
} else Read_Control_File(arg[0], api->control);
}
MPI_Bcast(api->control,sizeof(control_params),MPI_CHAR,0,world);
MPI_Bcast(api->out_control,sizeof(output_controls),MPI_CHAR,0,world);
// must reset these to local values after broadcast
api->control->me = comm->me;
@ -431,8 +419,8 @@ void PairReaxC::setup()
write_reax_lists();
api->system->wsize = comm->nprocs;
Initialize(api->system, api->control, api->data, api->workspace,
&api->lists, api->out_control, world);
Initialize(api->system, api->control, api->data,
api->workspace, &api->lists, world);
for (int k = 0; k < api->system->N; ++k) {
num_bonds[k] = api->system->my_atoms[k].num_bonds;
num_hbonds[k] = api->system->my_atoms[k].num_hbonds;
@ -500,7 +488,7 @@ void PairReaxC::compute(int eflag, int vflag)
// forces
Compute_Forces(api->system,api->control,api->data,api->workspace,&api->lists,api->out_control);
Compute_Forces(api->system,api->control,api->data,api->workspace,&api->lists);
read_reax_forces(vflag);
for (int k = 0; k < api->system->N; ++k) {
@ -526,9 +514,6 @@ void PairReaxC::compute(int eflag, int vflag)
ecoul += api->data->my_en.e_ele;
ecoul += api->data->my_en.e_pol;
// eng_vdwl += evdwl;
// eng_coul += ecoul;
// Store the different parts of the energy
// in a list for output by compute pair command
@ -554,8 +539,6 @@ void PairReaxC::compute(int eflag, int vflag)
api->data->step = update->ntimestep;
Output_Results(api->system, api->control, api->data, &api->lists, api->out_control, world);
// populate tmpid and tmpbo arrays for fix reax/c/species
int i, j;

5
src/USER-REAXC/reaxc_bond_orders.cpp
View File

@ -333,13 +333,10 @@ namespace ReaxFF {
return 1;
}
return 0;
}
void BO(reax_system *system, control_params * /*control*/, simulation_data * /*data*/,
storage *workspace, reax_list **lists, output_controls * /*out_control*/ )
void BO(reax_system *system, storage *workspace, reax_list **lists)
{
int i, j, pj, type_i, type_j;
int start_i, end_i, sym_index;

47
src/USER-REAXC/reaxc_control.cpp
View File

@ -41,7 +41,7 @@ using LAMMPS_NS::utils::logmesg;
using LAMMPS_NS::ValueTokenizer;
namespace ReaxFF {
static std::unordered_set<std::string> ignored_keywords = {
static std::unordered_set<std::string> inactive_keywords = {
"ensemble_type", "nsteps", "dt", "proc_by_dim", "random_vel",
"restart_format", "restart_freq", "reposition_atoms",
"restrict_bonds", "remove_CoM_vel", "debug_level", "reneighbor",
@ -50,8 +50,9 @@ namespace ReaxFF {
"t_freq", "pressure", "p_mass", "pt_mass", "compress", "press_mode",
"geo_format", "traj_compress", "traj_method", "molecular_analysis",
"ignore", "dipole_anal", "freq_dipole_anal", "diffusion_coef",
"freq_diffusion_coef", "restrict_type"
};
"freq_diffusion_coef", "restrict_type", "traj_title", "simulation_name",
"energy_update_freq", "atom_info", "atom_velocities", "atom_forces",
"bond_info", "angle_info" };
class parser_error : public std::exception {
std::string message;
@ -62,14 +63,11 @@ namespace ReaxFF {
// NOTE: this function is run on MPI rank 0 only
void Read_Control_File(const char *control_file, control_params *control,
output_controls *out_control)
void Read_Control_File(const char *control_file, control_params *control)
{
auto error = control->error_ptr;
auto lmp = control->lmp_ptr;
/* assign default values */
strcpy(control->sim_name, "simulate");
control->nthreads = 1;
control->tabulate = 0;
control->virial = 0;
@ -79,13 +77,6 @@ namespace ReaxFF {
control->thb_cutsq = 0.00001;
control->hbond_cut = 7.5;
out_control->write_steps = 0;
out_control->energy_update_freq = 0;
strcpy(out_control->traj_title, "default_title");
out_control->atom_info = 0;
out_control->bond_info = 0;
out_control->angle_info = 0;
/* read control parameters file */
try {
LAMMPS_NS::TextFileReader reader(control_file, "ReaxFF control");
@ -102,12 +93,9 @@ namespace ReaxFF {
if (!values.has_next())
throw parser_error(fmt::format("No value(s) for control parameter: {}\n",keyword));
if (ignored_keywords.find(keyword) != ignored_keywords.end()) {
logmesg(lmp,fmt::format("Ignoring inactive control parameter: {}\n",keyword));
} else if (keyword == "simulation_name") {
strcpy(control->sim_name, values.next_string().c_str());
} else if (keyword == "energy_update_freq") {
out_control->energy_update_freq = values.next_int();
if (inactive_keywords.find(keyword) != inactive_keywords.end()) {
error->warning(FLERR,fmt::format("Ignoring inactive control "
"parameter: {}",keyword));
} else if (keyword == "nbrhood_cutoff") {
control->bond_cut = values.next_double();
} else if (keyword == "bond_graph_cutoff") {
@ -121,21 +109,12 @@ namespace ReaxFF {
} else if (keyword == "tabulate_long_range") {
control->tabulate = values.next_int();
} else if (keyword == "write_freq") {
out_control->write_steps = values.next_int();
} else if (keyword == "traj_title") {
strcpy(out_control->traj_title, values.next_string().c_str());
} else if (keyword == "atom_info") {
out_control->atom_info += values.next_int() * 4;
} else if (keyword == "atom_velocities") {
out_control->atom_info += values.next_int() * 2;
} else if (keyword == "atom_forces") {
out_control->atom_info += values.next_int() * 1;
} else if (keyword == "bond_info") {
out_control->bond_info = values.next_int();
} else if (keyword == "angle_info") {
out_control->angle_info = values.next_int();
if (values.next_int() > 0)
error->warning(FLERR,"Support for writing native trajectories has "
"been removed after LAMMPS version 8 April 2021");
} else {
throw parser_error(fmt::format("Unknown parameter {} in control file", keyword));
throw parser_error(fmt::format("Unknown parameter {} in "
"control file", keyword));
}
}
} catch (LAMMPS_NS::EOFException &) {

27
src/USER-REAXC/reaxc_forces.cpp
View File

@ -34,11 +34,13 @@
namespace ReaxFF {
static void Compute_Bonded_Forces(reax_system *system, control_params *control,
simulation_data *data, storage *workspace,
reax_list **lists, output_controls *out_control)
static void Compute_Bonded_Forces(reax_system *system,
control_params *control,
simulation_data *data,
storage *workspace,
reax_list **lists)
{
BO(system, control, data, workspace, lists, out_control);
BO(system, workspace, lists);
Bonds(system, data, workspace, lists);
Atom_Energy(system, control, data, workspace, lists);
Valence_Angles(system, control, data, workspace, lists);
@ -47,8 +49,10 @@ namespace ReaxFF {
Hydrogen_Bonds(system, control, data, workspace, lists);
}
static void Compute_NonBonded_Forces(reax_system *system, control_params *control,
simulation_data *data, storage *workspace,
static void Compute_NonBonded_Forces(reax_system *system,
control_params *control,
simulation_data *data,
storage *workspace,
reax_list **lists)
{
@ -60,7 +64,7 @@ namespace ReaxFF {
}
static void Compute_Total_Force(reax_system *system, control_params *control,
storage *workspace, reax_list **lists)
storage *workspace, reax_list **lists)
{
int i, pj;
reax_list *bonds = (*lists) + BONDS;
@ -113,10 +117,6 @@ namespace ReaxFF {
system->my_atoms[i].num_hbonds =
(int)(MAX(Num_Entries(Hindex, hbonds)*saferzone, system->minhbonds));
//if(Num_Entries(i, hbonds) >=
//(Start_Index(i+1,hbonds)-Start_Index(i,hbonds))*0.90/*DANGER_ZONE*/) {
// workspace->realloc.hbonds = 1;
if (Hindex < numH-1)
comp = Start_Index(Hindex+1, hbonds);
else comp = hbonds->num_intrs;
@ -374,14 +374,13 @@ namespace ReaxFF {
void Compute_Forces(reax_system *system, control_params *control,
simulation_data *data, storage *workspace,
reax_list **lists, output_controls *out_control)
reax_list **lists)
{
Init_Forces_noQEq(system, control, data, workspace, lists);
/********* bonded interactions ************/
Compute_Bonded_Forces(system, control, data, workspace,
lists, out_control);
Compute_Bonded_Forces(system, control, data, workspace, lists);
/********* nonbonded interactions ************/
Compute_NonBonded_Forces(system, control, data, workspace, lists);

4
src/USER-REAXC/reaxc_init_md.cpp
View File

@ -163,14 +163,12 @@ namespace ReaxFF {
void Initialize(reax_system *system, control_params *control,
simulation_data *data, storage *workspace,
reax_list **lists, output_controls *out_control,
MPI_Comm world)
reax_list **lists, MPI_Comm world)
{
Init_System(system, control);
Init_Simulation_Data(data);
Init_Workspace(system, control, workspace);
Init_Lists(system, control, lists);
Init_Output_Files(system, control, out_control, world);
if (control->tabulate)
Init_Lookup_Tables(system, control, workspace, world);
}

99
src/USER-REAXC/reaxc_io_tools.cpp
View File

@ -1,99 +0,0 @@
/*----------------------------------------------------------------------
PuReMD - Purdue ReaxFF Molecular Dynamics Program
Copyright (2010) Purdue University
Hasan Metin Aktulga, hmaktulga@lbl.gov
Joseph Fogarty, jcfogart@mail.usf.edu
Sagar Pandit, pandit@usf.edu
Ananth Y Grama, ayg@cs.purdue.edu
Please cite the related publication:
H. M. Aktulga, J. C. Fogarty, S. A. Pandit, A. Y. Grama,
"Parallel Reactive Molecular Dynamics: Numerical Methods and
Algorithmic Techniques", Parallel Computing, in press.
This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License as
published by the Free Software Foundation; either version 2 of
the License, or (at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
See the GNU General Public License for more details:
<https://www.gnu.org/licenses/>.
----------------------------------------------------------------------*/
#include "reaxff_api.h"
namespace ReaxFF {
void Collect_System_Energy(reax_system *system, simulation_data *data,
MPI_Comm comm)
{
double my_en[13], sys_en[13];
my_en[0] = data->my_en.e_bond;
my_en[1] = data->my_en.e_ov;
my_en[2] = data->my_en.e_un;
my_en[3] = data->my_en.e_lp;
my_en[4] = data->my_en.e_ang;
my_en[5] = data->my_en.e_pen;
my_en[6] = data->my_en.e_coa;
my_en[7] = data->my_en.e_hb;
my_en[8] = data->my_en.e_tor;
my_en[9] = data->my_en.e_con;
my_en[10] = data->my_en.e_vdW;
my_en[11] = data->my_en.e_ele;
my_en[12] = data->my_en.e_pol;
MPI_Reduce( my_en, sys_en, 13, MPI_DOUBLE, MPI_SUM, MASTER_NODE, comm );
if (system->my_rank == MASTER_NODE) {
data->sys_en.e_bond = sys_en[0];
data->sys_en.e_ov = sys_en[1];
data->sys_en.e_un = sys_en[2];
data->sys_en.e_lp = sys_en[3];
data->sys_en.e_ang = sys_en[4];
data->sys_en.e_pen = sys_en[5];
data->sys_en.e_coa = sys_en[6];
data->sys_en.e_hb = sys_en[7];
data->sys_en.e_tor = sys_en[8];
data->sys_en.e_con = sys_en[9];
data->sys_en.e_vdW = sys_en[10];
data->sys_en.e_ele = sys_en[11];
data->sys_en.e_pol = sys_en[12];
}
}
void Init_Output_Files(reax_system *system, control_params *control,
output_controls *out_control, MPI_Comm world)
{
if (out_control->write_steps > 0)
Init_Traj(system, control, out_control, world);
}
/************************ close output files ************************/
void Close_Output_Files(reax_system *system, output_controls *out_control)
{
if (out_control->write_steps > 0)
End_Traj(system->my_rank, out_control);
}
void Output_Results(reax_system *system, control_params *control,
simulation_data *data, reax_list **lists,
output_controls *out_control, MPI_Comm world)
{
if ((out_control->energy_update_freq > 0 &&
data->step%out_control->energy_update_freq == 0) ||
(out_control->write_steps > 0 &&
data->step%out_control->write_steps == 0)) {
/* update system-wide energies */
Collect_System_Energy(system, data, world);
/* write current frame */
if (out_control->write_steps > 0 && data->step % out_control->write_steps == 0) {
Append_Frame(system, control, data, lists, out_control, world);
}
}
}
}

25
src/USER-REAXC/reaxc_reset_tools.cpp
View File

@ -49,28 +49,9 @@ namespace ReaxFF {
}
}
static void Reset_Energies(energy_data *en)
{
en->e_bond = 0;
en->e_ov = 0;
en->e_un = 0;
en->e_lp = 0;
en->e_ang = 0;
en->e_pen = 0;
en->e_coa = 0;
en->e_hb = 0;
en->e_tor = 0;
en->e_con = 0;
en->e_vdW = 0;
en->e_ele = 0;
en->e_pol = 0;
}
void Reset_Simulation_Data(simulation_data* data)
{
Reset_Energies(&data->my_en);
Reset_Energies(&data->sys_en);
memset(&data->my_en,0,sizeof(energy_data));
}
void Reset_Workspace(reax_system *system, storage *workspace)
@ -136,8 +117,8 @@ namespace ReaxFF {
}
void Reset(reax_system *system, control_params *control, simulation_data *data,
storage *workspace, reax_list **lists)
void Reset(reax_system *system, control_params *control,
simulation_data *data, storage *workspace, reax_list **lists)
{
Reset_Atoms(system, control);
Reset_Simulation_Data(data);

569
src/USER-REAXC/reaxc_traj.cpp
View File

@ -1,569 +0,0 @@
/*----------------------------------------------------------------------
PuReMD - Purdue ReaxFF Molecular Dynamics Program
Copyright (2010) Purdue University
Hasan Metin Aktulga, hmaktulga@lbl.gov
Joseph Fogarty, jcfogart@mail.usf.edu
Sagar Pandit, pandit@usf.edu
Ananth Y Grama, ayg@cs.purdue.edu
Please cite the related publication:
H. M. Aktulga, J. C. Fogarty, S. A. Pandit, A. Y. Grama,
"Parallel Reactive Molecular Dynamics: Numerical Methods and
Algorithmic Techniques", Parallel Computing, in press.
This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License as
published by the Free Software Foundation; either version 2 of
the License, or (at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
See the GNU General Public License for more details:
<https://www.gnu.org/licenses/>.
----------------------------------------------------------------------*/
#include "reaxff_api.h"
#include <cstdio>
#include <cstring>
#include "error.h"
#define MAX_TRJ_LINE_LEN 120
#define MAX_TRJ_BUFFER_SIZE (MAX_TRJ_LINE_LEN * 100)
#define NUM_HEADER_LINES 37
#define HEADER_LINE_LEN 62
#define INIT_DESC_LEN 32
#define ATOM_BASIC_LEN 50
#define ATOM_wV_LEN 80
#define ATOM_wF_LEN 80
#define ATOM_FULL_LEN 110
#define BOND_BASIC_LEN 39
#define BOND_FULL_LEN 69
#define ANGLE_BASIC_LEN 38
namespace ReaxFF {
enum ATOM_LINE_OPTS { OPT_NOATOM = 0, OPT_ATOM_BASIC = 4, OPT_ATOM_wF = 5, OPT_ATOM_wV = 6, OPT_ATOM_FULL = 7, NR_OPT_ATOM = 8 };
enum BOND_LINE_OPTS { OPT_NOBOND, OPT_BOND_BASIC, OPT_BOND_FULL, NR_OPT_BOND };
enum ANGLE_LINE_OPTS { OPT_NOANGLE, OPT_ANGLE_BASIC, NR_OPT_ANGLE };
std::string fmtline(const char *key, double val)
{
return fmt::format("{:<37}{:<24.3f}\n",key,val);
}
std::string fmtline(const char *key, double val1, double val2, double val3)
{
return fmt::format("{:<32}{:>9.3f},{:>9.3f},{:>9.3f}\n",key,val1,val2,val3);
}
template <typename T>
std::string fmtline(const char *key, T val)
{
return fmt::format("{:<37}{:<24}\n",key,val);
}
template <typename T>
std::string fmtline(const char *key, T val1, T val2)
{
return fmt::format("{:<36}{:<12},{:<12}\n",key,val1,val2);
}
static void Reallocate_Output_Buffer(LAMMPS_NS::Error *error_ptr, output_controls *out_control, int req_space)
{
if (out_control->buffer_len > 0)
free(out_control->buffer);
out_control->buffer_len = (int)(req_space*REAX_SAFE_ZONE);
out_control->buffer = (char*) malloc(out_control->buffer_len*sizeof(char));
if (!out_control->buffer)
error_ptr->one(FLERR,fmt::format("Insufficient memory for required buffer "
"size {}", (req_space*REAX_SAFE_ZONE)));
}
static void Write_Skip_Line(output_controls *out_control, int my_rank, int skip, int num_section)
{
if (my_rank == MASTER_NODE)
fmt::print(out_control->strj,fmtline("chars_to_skip_section:",skip,num_section));
}
static void Write_Header(reax_system *system, control_params *control, output_controls *out_control)
{
char atom_formats[8][40] = { "none", "invalid", "invalid", "invalid",
"xyz_q",
"xyz_q_fxfyfz",
"xyz_q_vxvyvz",
"detailed_atom_info" };
char bond_formats[3][30] = { "none",
"basic_bond_info",
"detailed_bond_info" };
char angle_formats[2][30] = { "none", "basic_angle_info" };
/* only the master node writes into trajectory header */
if (system->my_rank == MASTER_NODE) {
std::string buffer("");
/* to skip the header */
buffer += fmtline("chars_to_skip_header:",(NUM_HEADER_LINES-1) * HEADER_LINE_LEN);
/* general simulation info */
buffer += fmtline("simulation_name:", out_control->traj_title);
buffer += fmtline("number_of_atoms:", system->bigN);
buffer += fmtline("ensemble_type:", "(unknown)");
buffer += fmtline("number_of_steps:", 0);
buffer += fmtline("timestep_length_(in_fs):", 0.0);
/* restart info */
buffer += fmtline("is_this_a_restart?:", "no");
buffer += fmtline("write_restart_files?:", "no");
buffer += fmtline("frequency_to_write_restarts:", 0);
/* preferences */
buffer += fmtline("tabulate_long_range_intrs?:",(control->tabulate ? "yes" : "no"));
buffer += fmtline("table_size:", control->tabulate);
buffer += fmtline("restrict_bonds?:", "no");
buffer += fmtline("bond_restriction_length:", 0);
buffer += fmtline("reposition_atoms?:", "fit to periodic box");
buffer += fmtline("remove_CoM_velocity?:", 0);
buffer += fmtline("bonded_intr_dist_cutoff:",control->bond_cut);
buffer += fmtline("nonbonded_intr_dist_cutoff:",control->nonb_cut);
buffer += fmtline("hbond_dist_cutoff:",control->hbond_cut);
buffer += fmtline("reax_bond_threshold:",control->bo_cut);
buffer += fmtline("physical_bond_threshold:", control->bg_cut);
buffer += fmtline("valence_angle_threshold:",control->thb_cut);
buffer += fmtline("QEq_tolerance:", 0.0);
/* temperature controls */
buffer += fmtline("initial_temperature:", 0.0);
buffer += fmtline("target_temperature:", 0.0);
buffer += fmtline("thermal_inertia:", 0.0);
buffer += fmtline("temperature_regime:", "(unknown)");
buffer += fmtline("temperature_change_rate_(K/ps):", 0.0);
/* pressure controls */
buffer += fmtline("target_pressure_(GPa):", 0.0, 0.0, 0.0);
buffer += fmtline("virial_inertia:", 0.0, 0.0, 0.0);
/* trajectory */
buffer += fmtline("energy_dumping_freq:", out_control->energy_update_freq);
buffer += fmtline("trajectory_dumping_freq:",out_control->write_steps);
buffer += fmtline("compress_trajectory_output?:", "no");
buffer += fmtline("trajectory_format:", "regular");
buffer += fmtline("atom_info:", atom_formats[out_control->atom_info]);
buffer += fmtline("bond_info:", bond_formats[out_control->bond_info]);
buffer += fmtline("angle_info:", angle_formats[out_control->angle_info]);
/* analysis */
buffer += fmtline("molecular_analysis:", "no");
buffer += fmtline("molecular_analysis_frequency:", 0);
/* dump out the buffer */
fputs(buffer.c_str(),out_control->strj);
}
}
static void Write_Init_Desc(reax_system *system, output_controls *out_control, MPI_Comm world)
{
int i, me, np, cnt, buffer_len, buffer_req;
reax_atom *p_atom;
MPI_Status status;
me = system->my_rank;
np = system->wsize;
/* skip info */
Write_Skip_Line(out_control, me, system->bigN*INIT_DESC_LEN, system->bigN);
if (me == MASTER_NODE)
buffer_req = system->bigN * INIT_DESC_LEN + 1;
else buffer_req = system->n * INIT_DESC_LEN + 1;
if (buffer_req > out_control->buffer_len * DANGER_ZONE)
Reallocate_Output_Buffer(system->error_ptr, out_control, buffer_req);
out_control->buffer[0] = 0;
for (i = 0; i < system->n; ++i) {
p_atom = &(system->my_atoms[i]);
auto buffer = fmt::format("{:9}{:3}{:9}{:10.3f}\n",
p_atom->orig_id, p_atom->type, p_atom->name,
system->reax_param.sbp[p_atom->type].mass);
strncpy(out_control->buffer + i*INIT_DESC_LEN, buffer.c_str(), INIT_DESC_LEN+1);
}
if (me != MASTER_NODE) {
MPI_Send(out_control->buffer, buffer_req-1, MPI_CHAR, MASTER_NODE,
np * INIT_DESCS + me, world);
} else {
buffer_len = system->n * INIT_DESC_LEN;
for (i = 0; i < np; ++i)
if (i != MASTER_NODE) {
MPI_Recv(out_control->buffer + buffer_len, buffer_req - buffer_len,
MPI_CHAR, i, np*INIT_DESCS+i, world, &status);
MPI_Get_count(&status, MPI_CHAR, &cnt);
buffer_len += cnt;
}
out_control->buffer[buffer_len] = 0;
fputs(out_control->buffer,out_control->strj);
}
}
void Init_Traj(reax_system *system, control_params *control,
output_controls *out_control, MPI_Comm world)
{
int atom_line_len[NR_OPT_ATOM] = { 0, 0, 0, 0, ATOM_BASIC_LEN, ATOM_wV_LEN, ATOM_wF_LEN, ATOM_FULL_LEN};
int bond_line_len[NR_OPT_BOND] = { 0, BOND_BASIC_LEN, BOND_FULL_LEN };
int angle_line_len[NR_OPT_ANGLE] = { 0, ANGLE_BASIC_LEN };
/* generate trajectory name */
auto fname = std::string(control->sim_name) + ".trj";
/* how should I write atoms? */
out_control->atom_line_len = atom_line_len[out_control->atom_info];
out_control->write_atoms = (out_control->atom_line_len ? 1 : 0);
/* bonds? */
out_control->bond_line_len = bond_line_len[out_control->bond_info];
out_control->write_bonds = (out_control->bond_line_len ? 1 : 0);
/* angles? */
out_control->angle_line_len = angle_line_len[out_control->angle_info];
out_control->write_angles = (out_control->angle_line_len ? 1 : 0);
/* allocate line & buffer space */
out_control->buffer_len = 0;
out_control->buffer = nullptr;
/* write trajectory header and atom info, if applicable */
if (system->my_rank == MASTER_NODE)
out_control->strj = fopen(fname.c_str(), "w");
Write_Header(system, control, out_control);
Write_Init_Desc(system, out_control, world);
}
static void Write_Frame_Header(reax_system *system, simulation_data *data, output_controls *out_control)
{
const int me = system->my_rank;
/* frame header lengths */
constexpr int num_frm_hdr_lines = 22;
/* only the master node writes into trajectory header */
if (me == MASTER_NODE) {
std::string buffer("");
/* skip info */
buffer += fmtline("chars_to_skip_frame_header:", (num_frm_hdr_lines-1)*HEADER_LINE_LEN);
/* step & time */
buffer += fmtline("step:", data->step);
buffer += fmtline("time_in_ps:", 0.0);
/* box info */
buffer += fmtline("volume:", 0.0);
buffer += fmtline("box_dimensions:", 0.0, 0.0, 0.0);
buffer += fmtline("coordinate_angles:", 90.0, 90.0, 90.0);
/* system T and P */
buffer += fmtline("temperature:", 0.0);
buffer += fmtline("pressure:", 0.0);
/* energies */
double epot = data->sys_en.e_bond + data->sys_en.e_ov + data->sys_en.e_un
+ data->sys_en.e_lp + data->sys_en.e_ang + data->sys_en.e_pen
+ data->sys_en.e_coa + data->sys_en.e_hb + data->sys_en.e_tor
+ data->sys_en.e_con + data->sys_en.e_vdW
+ data->sys_en.e_ele + data->my_en.e_pol;
buffer += fmtline("total_energy:", epot);
buffer += fmtline("total_kinetic:", 0.0);
buffer += fmtline("total_potential:", epot);
buffer += fmtline("bond_energy:", data->sys_en.e_bond);
buffer += fmtline("atom_energy:", data->sys_en.e_ov + data->sys_en.e_un);
buffer += fmtline("lone_pair_energy:", data->sys_en.e_lp);
buffer += fmtline("valence_angle_energy:", data->sys_en.e_ang + data->sys_en.e_pen);
buffer += fmtline("3-body_conjugation:", data->sys_en.e_coa);
buffer += fmtline("hydrogen_bond_energy:", data->sys_en.e_hb);
buffer += fmtline("torsion_angle_energy:", data->sys_en.e_tor);
buffer += fmtline("4-body_conjugation:", data->sys_en.e_con);
buffer += fmtline("vdWaals_energy:", data->sys_en.e_vdW);
buffer += fmtline("electrostatics_energy:", data->sys_en.e_ele);
buffer += fmtline("polarization_energy:", data->sys_en.e_pol);
/* dump out the buffer */
fputs(buffer.c_str(),out_control->strj);
}
}
static void Write_Atoms(reax_system *system, output_controls *out_control,
MPI_Comm world)
{
int i, me, np, line_len, buffer_len, buffer_req, cnt;
MPI_Status status;
reax_atom *p_atom;
me = system->my_rank;
np = system->wsize;
line_len = out_control->atom_line_len;
Write_Skip_Line(out_control, me, system->bigN*line_len, system->bigN);
if (me == MASTER_NODE)
buffer_req = system->bigN * line_len + 1;
else buffer_req = system->n * line_len + 1;
if (buffer_req > out_control->buffer_len * DANGER_ZONE)
Reallocate_Output_Buffer(system->error_ptr, out_control, buffer_req);
/* fill in buffer */
out_control->buffer[0] = 0;
for (i = 0; i < system->n; ++i) {
std::string buffer("");
p_atom = &(system->my_atoms[i]);
buffer += fmt::format("{:9}{:10.3f}{:10.3f}{:10.3f}",p_atom->orig_id,
p_atom->x[0], p_atom->x[1],p_atom->x[2]);
switch (out_control->atom_info) {
case OPT_ATOM_BASIC:
buffer += fmt::format("{:10.3f}\n",p_atom->q);
break;
case OPT_ATOM_wF:
buffer += fmt::format("{:10.3f}{:10.3f}{:10.3f}{:10.3f}\n",
p_atom->f[0], p_atom->f[1], p_atom->f[2], p_atom->q);
break;
case OPT_ATOM_wV:
buffer += fmt::format("{:10.3f}{:10.3f}{:10.3f}{:10.3f}\n",
p_atom->v[0], p_atom->v[1], p_atom->v[2], p_atom->q);
break;
case OPT_ATOM_FULL:
buffer += fmt::format("{:10.3f}{:10.3f}{:10.3f}{:10.3f}{:10.3f}{:10.3f}\n",
p_atom->v[0], p_atom->v[1], p_atom->v[2],
p_atom->f[0], p_atom->f[1], p_atom->f[2], p_atom->q);
break;
default:
system->error_ptr->one(FLERR,"Write_traj_atoms: unknown atom trajectory format");
}
strncpy(out_control->buffer + i*line_len, buffer.c_str(), line_len+1);
}
if (me != MASTER_NODE) {
MPI_Send(out_control->buffer, buffer_req-1, MPI_CHAR, MASTER_NODE,
np*ATOM_LINES+me, world);
} else {
buffer_len = system->n * line_len;
for (i = 0; i < np; ++i)
if (i != MASTER_NODE) {
MPI_Recv(out_control->buffer + buffer_len, buffer_req - buffer_len,
MPI_CHAR, i, np*ATOM_LINES+i, world, &status);
MPI_Get_count(&status, MPI_CHAR, &cnt);
buffer_len += cnt;
}
out_control->buffer[buffer_len] = 0;
fputs(out_control->buffer, out_control->strj);
}
}
static void Write_Bonds(reax_system *system, control_params *control, reax_list *bonds,
output_controls *out_control, MPI_Comm world)
{
int i, j, pj, me, np;
int my_bonds, num_bonds;
int line_len, buffer_len, buffer_req, cnt;
MPI_Status status;
bond_data *bo_ij;
me = system->my_rank;
np = system->wsize;
line_len = out_control->bond_line_len;
/* count the number of bonds I will write */
my_bonds = 0;
for (i=0; i < system->n; ++i)
for (pj = Start_Index(i, bonds); pj < End_Index(i, bonds); ++pj) {
j = bonds->select.bond_list[pj].nbr;
if (system->my_atoms[i].orig_id <= system->my_atoms[j].orig_id &&
bonds->select.bond_list[pj].bo_data.BO >= control->bg_cut)
++my_bonds;
}
/* allreduce - total number of bonds */
MPI_Allreduce(&my_bonds, &num_bonds, 1, MPI_INT, MPI_SUM, world);
Write_Skip_Line(out_control, me, num_bonds*line_len, num_bonds);
if (me == MASTER_NODE)
buffer_req = num_bonds * line_len + 1;
else buffer_req = my_bonds * line_len + 1;
if (buffer_req > out_control->buffer_len * DANGER_ZONE)
Reallocate_Output_Buffer(system->error_ptr, out_control, buffer_req);
/* fill in the buffer */
out_control->buffer[0] = 0;
my_bonds = 0;
for (i=0; i < system->n; ++i) {
for (pj = Start_Index(i, bonds); pj < End_Index(i, bonds); ++pj) {
bo_ij = &(bonds->select.bond_list[pj]);
j = bo_ij->nbr;
if (system->my_atoms[i].orig_id <= system->my_atoms[j].orig_id &&
bo_ij->bo_data.BO >= control->bg_cut) {
auto buffer = fmt::format("{:9}{:9}{:10.3f}{:10.3f}", system->my_atoms[i].orig_id,
system->my_atoms[j].orig_id,bo_ij->d,bo_ij->bo_data.BO);
switch (out_control->bond_info) {
case OPT_BOND_BASIC:
buffer += "\n";
break;
case OPT_BOND_FULL:
buffer += fmt::format("{:10.3f}{:10.3f}{:10.3f}\n", bo_ij->bo_data.BO_s,
bo_ij->bo_data.BO_pi, bo_ij->bo_data.BO_pi2);
break;
default:
system->error_ptr->one(FLERR, "Write_traj_bonds: FATAL! invalid bond_info option");
}
strncpy(out_control->buffer + my_bonds*line_len, buffer.c_str(), line_len+1);
++my_bonds;
}
}
}
if (me != MASTER_NODE) {
MPI_Send(out_control->buffer, buffer_req-1, MPI_CHAR, MASTER_NODE, np*BOND_LINES+me, world);
} else {
buffer_len = my_bonds * line_len;
for (i = 0; i < np; ++i)
if (i != MASTER_NODE) {
MPI_Recv(out_control->buffer + buffer_len, buffer_req - buffer_len,
MPI_CHAR, i, np*BOND_LINES+i, world, &status);
MPI_Get_count(&status, MPI_CHAR, &cnt);
buffer_len += cnt;
}
out_control->buffer[buffer_len] = 0;
fputs(out_control->buffer,out_control->strj);
}
}
static void Write_Angles(reax_system *system, control_params *control,
reax_list *bonds, reax_list *thb_intrs,
output_controls *out_control, MPI_Comm world)
{
int i, j, k, pi, pk, me, np;
int my_angles, num_angles;
int line_len, buffer_len, buffer_req, cnt;
bond_data *bo_ij, *bo_jk;
three_body_interaction_data *angle_ijk;
MPI_Status status;
me = system->my_rank;
np = system->wsize;
line_len = out_control->angle_line_len;
/* count the number of valence angles I will output */
my_angles = 0;
for (j = 0; j < system->n; ++j)
for (pi = Start_Index(j, bonds); pi < End_Index(j, bonds); ++pi) {
bo_ij = &(bonds->select.bond_list[pi]);
i = bo_ij->nbr;
if (bo_ij->bo_data.BO >= control->bg_cut) // physical j&i bond
for (pk = Start_Index(pi, thb_intrs);
pk < End_Index(pi, thb_intrs); ++pk) {
angle_ijk = &(thb_intrs->select.three_body_list[pk]);
k = angle_ijk->thb;
bo_jk = &(bonds->select.bond_list[ angle_ijk->pthb ]);
if (system->my_atoms[i].orig_id < system->my_atoms[k].orig_id &&
bo_jk->bo_data.BO >= control->bg_cut) // physical j&k bond
++my_angles;
}
}
/* total number of valences */
MPI_Allreduce(&my_angles, &num_angles, 1, MPI_INT, MPI_SUM, world);
Write_Skip_Line(out_control, me, num_angles*line_len, num_angles);
if (me == MASTER_NODE)
buffer_req = num_angles * line_len + 1;
else buffer_req = my_angles * line_len + 1;
if (buffer_req > out_control->buffer_len * DANGER_ZONE)
Reallocate_Output_Buffer(system->error_ptr, out_control, buffer_req);
/* fill in the buffer */
my_angles = 0;
out_control->buffer[0] = 0;
for (j = 0; j < system->n; ++j)
for (pi = Start_Index(j, bonds); pi < End_Index(j, bonds); ++pi) {
bo_ij = &(bonds->select.bond_list[pi]);
i = bo_ij->nbr;
if (bo_ij->bo_data.BO >= control->bg_cut) // physical j&i bond
for (pk = Start_Index(pi, thb_intrs);
pk < End_Index(pi, thb_intrs); ++pk) {
angle_ijk = &(thb_intrs->select.three_body_list[pk]);
k = angle_ijk->thb;
bo_jk = &(bonds->select.bond_list[ angle_ijk->pthb ]);
if (system->my_atoms[i].orig_id < system->my_atoms[k].orig_id &&
bo_jk->bo_data.BO >= control->bg_cut) { // physical j&k bond
auto buffer = fmt::format("{:9}{:9}{:9}{:10.3f}\n",
system->my_atoms[i].orig_id, system->my_atoms[j].orig_id,
system->my_atoms[k].orig_id, RAD2DEG(angle_ijk->theta));
strncpy(out_control->buffer + my_angles*line_len, buffer.c_str(), line_len+1);
++my_angles;
}
}
}
if (me != MASTER_NODE) {
MPI_Send(out_control->buffer, buffer_req-1, MPI_CHAR, MASTER_NODE,
np*ANGLE_LINES+me, world);
} else {
buffer_len = my_angles * line_len;
for (i = 0; i < np; ++i)
if (i != MASTER_NODE) {
MPI_Recv(out_control->buffer + buffer_len, buffer_req - buffer_len,
MPI_CHAR, i, np*ANGLE_LINES+i, world, &status);
MPI_Get_count(&status, MPI_CHAR, &cnt);
buffer_len += cnt;
}
out_control->buffer[buffer_len] = 0;
fputs(out_control->buffer,out_control->strj);
}
}
void Append_Frame(reax_system *system, control_params *control,
simulation_data *data, reax_list **lists,
output_controls *out_control, MPI_Comm world)
{
Write_Frame_Header(system, data, out_control);
if (out_control->write_atoms)
Write_Atoms(system, out_control, world);
if (out_control->write_bonds)
Write_Bonds(system, control, (*lists + BONDS), out_control, world);
if (out_control->write_angles)
Write_Angles(system, control, (*lists + BONDS), (*lists + THREE_BODIES),
out_control, world);
}
void End_Traj(int my_rank, output_controls *out_control)
{
if (my_rank == MASTER_NODE)
fclose(out_control->strj);
free(out_control->buffer);
}
}

39
src/USER-REAXC/reaxff_api.h
View File

@ -31,7 +31,6 @@ namespace ReaxFF
struct API {
control_params *control;
reax_system *system;
output_controls *out_control;
simulation_data *data;
storage *workspace;
reax_list *lists;
@ -39,7 +38,7 @@ namespace ReaxFF
// exported Functions
// allocate X
// allocate
extern void Allocate_Workspace(control_params *, storage *, int);
extern void DeAllocate_System(reax_system *);
@ -50,8 +49,7 @@ namespace ReaxFF
// bond orders
extern void BO(reax_system *, control_params *, simulation_data *,
storage *, reax_list **, output_controls *);
extern void BO(reax_system *, storage *, reax_list **);
extern int BOp(storage *, reax_list *, double, int, int, far_neighbor_data *,
single_body_parameters *, single_body_parameters *,
two_body_parameters *);
@ -64,7 +62,7 @@ namespace ReaxFF
// control
extern void Read_Control_File(const char *, control_params *, output_controls *);
extern void Read_Control_File(const char *, control_params *);
// ffield
@ -73,8 +71,8 @@ namespace ReaxFF
// forces
extern void Compute_Forces(reax_system *, control_params *, simulation_data *,
storage *, reax_list **, output_controls *);
extern void Compute_Forces(reax_system *, control_params *,
simulation_data *, storage *, reax_list **);
extern void Estimate_Storages(reax_system *, control_params *, reax_list **,
int *, int *, int *, int *);
@ -88,16 +86,7 @@ namespace ReaxFF
extern void Init_Simulation_Data(simulation_data *);
extern void Init_Workspace(reax_system *, control_params *, storage *);
extern void Initialize(reax_system *, control_params *, simulation_data *,
storage *, reax_list **, output_controls *, MPI_Comm);
// io tools
extern void Init_Output_Files(reax_system *, control_params *,
output_controls *, MPI_Comm);
extern void Close_Output_Files(reax_system *, output_controls *);
extern void Output_Results(reax_system *, control_params *, simulation_data *,
reax_list **, output_controls *, MPI_Comm);
extern void Collect_System_Energy(reax_system *, simulation_data *, MPI_Comm);
storage *, reax_list **, MPI_Comm);
// lists
@ -163,24 +152,16 @@ namespace ReaxFF
rvec, double, three_body_interaction_data *,
three_body_interaction_data *, rvec, rvec,
rvec, rvec);
extern void Torsion_Angles(reax_system *, control_params *, simulation_data *,
storage *, reax_list **);
// traj
extern void Init_Traj(reax_system *, control_params *,
output_controls *, MPI_Comm);
extern void End_Traj(int, output_controls *);
extern void Append_Frame(reax_system *, control_params *, simulation_data *,
reax_list **, output_controls *, MPI_Comm);
extern void Torsion_Angles(reax_system *, control_params *,
simulation_data *, storage *, reax_list **);
// valence angles
extern void Calculate_Theta(rvec, double, rvec, double, double *, double *);
extern void Calculate_dCos_Theta(rvec, double, rvec, double,
rvec *, rvec *, rvec *);
extern void Valence_Angles(reax_system *, control_params *, simulation_data *,
storage *, reax_list **);
extern void Valence_Angles(reax_system *, control_params *,
simulation_data *, storage *, reax_list **);
// vector

24
src/USER-REAXC/reaxff_types.h
View File

@ -228,7 +228,6 @@ namespace ReaxFF
/* system control parameters */
struct control_params
{
char sim_name[REAX_MAX_STR];
int nthreads;
double bond_cut;
@ -270,9 +269,7 @@ namespace ReaxFF
struct simulation_data
{
rc_bigint step;
energy_data my_en; // per MPI rank energies
energy_data sys_en; // global energies
};
struct three_body_interaction_data
@ -399,27 +396,6 @@ namespace ReaxFF
class LAMMPS_NS::Error *error_ptr;
};
struct output_controls
{
FILE *strj;
int atom_line_len;
int bond_line_len;
int angle_line_len;
int write_atoms;
int write_bonds;
int write_angles;
int buffer_len;
char *buffer;
int write_steps;
char traj_title[81];
int atom_info;
int bond_info;
int angle_info;
int energy_update_freq;
};
struct LR_lookup_table
{
double xmin, xmax;