From cee848f948653c39f50b5686a1926f373e30c458 Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Thu, 22 Oct 2015 22:06:49 +0000
Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14164
f3b2605a-c512-4ea7-a41b-209d697bcdaa
---
src/CORESHELL/pair_lj_cut_coul_long_cs.cpp | 3 +-
src/DIPOLE/pair_lj_cut_dipole_long.cpp | 4 +-
src/KOKKOS/verlet_kokkos.cpp | 18 +-
src/KSPACE/pair_lj_cut_coul_long.cpp | 4 +-
src/KSPACE/pair_lj_cut_coul_msm.cpp | 4 +-
src/MANYBODY/pair_adp.cpp | 1 +
src/MANYBODY/pair_airebo.cpp | 1 +
src/MANYBODY/pair_comb.cpp | 2 +
src/MANYBODY/pair_comb3.cpp | 2 +
src/MANYBODY/pair_eam.cpp | 2 +
src/MANYBODY/pair_eim.cpp | 1 +
src/MANYBODY/pair_nb3b_harmonic.cpp | 1 +
src/MANYBODY/pair_polymorphic.cpp | 2 +-
src/MANYBODY/pair_sw.cpp | 1 +
src/MANYBODY/pair_tersoff.cpp | 1 +
src/PYTHON/python.cpp | 2 +
src/REPLICA/verlet_split.cpp | 13 +-
src/USER-AWPMD/pair_awpmd_cut.h | 2 +-
src/USER-CUDA/verlet_cuda.cpp | 23 +-
src/USER-DIFFRACTION/compute_saed.h | 1 -
src/USER-DRUDE/pair_lj_cut_thole_long.cpp | 4 +-
src/USER-DRUDE/pair_thole.cpp | 2 +-
.../pair_lj_charmm_coul_long_soft.cpp | 2 +
src/USER-FEP/pair_lj_cut_coul_long_soft.cpp | 2 +
src/USER-FEP/pair_lj_cut_soft.cpp | 1 +
src/USER-H5MD/dump_h5md.cpp | 2 +-
src/USER-H5MD/dump_h5md.h | 1 -
src/USER-INTEL/verlet_intel.cpp | 22 +-
src/USER-MISC/fix_imd.cpp | 20 +-
src/USER-MISC/fix_pimd.cpp | 7 +-
src/USER-MISC/fix_srp.cpp | 12 +-
src/USER-MISC/fix_ti_rs.cpp | 26 +-
src/USER-MISC/fix_ti_rs.h | 6 +-
src/USER-MISC/fix_ti_spring.cpp | 23 +-
src/USER-MISC/fix_ti_spring.h | 6 +-
src/USER-MISC/pair_tersoff_table.cpp | 11 -
src/USER-MOLFILE/dump_molfile.cpp | 1 +
src/USER-OMP/fix_omp.cpp | 40 +-
src/USER-OMP/pair_tersoff_table_omp.cpp | 11 -
src/USER-OMP/respa_omp.cpp | 17 +
src/USER-REAXC/reaxc_bond_orders.cpp | 6 -
src/USER-REAXC/reaxc_bond_orders.h | 17 -
src/USER-REAXC/reaxc_ffield.cpp | 2 +
src/USER-REAXC/reaxc_forces.cpp | 219 +--------
src/USER-REAXC/reaxc_io_tools.cpp | 428 ------------------
src/USER-REAXC/reaxc_io_tools.h | 74 ---
src/USER-REAXC/reaxc_lookup.cpp | 24 -
src/USER-REAXC/reaxc_reset_tools.cpp | 23 -
src/USER-REAXC/reaxc_reset_tools.h | 5 -
src/USER-REAXC/reaxc_system_props.cpp | 215 ---------
src/USER-REAXC/reaxc_system_props.h | 12 -
src/USER-REAXC/reaxc_tool_box.cpp | 120 -----
src/USER-REAXC/reaxc_tool_box.h | 24 -
src/USER-REAXC/reaxc_vector.cpp | 76 ----
src/USER-REAXC/reaxc_vector.h | 9 -
src/USER-SPH/atom_vec_meso.cpp | 26 +-
src/fix.h | 4 +-
src/info.cpp | 1 +
src/input.cpp | 69 ++-
src/lammps.cpp | 1 +
src/min_cg.cpp | 1 +
src/min_fire.cpp | 1 +
src/min_quickmin.cpp | 1 +
src/min_sd.cpp | 1 +
src/modify.cpp | 19 +-
src/modify.h | 5 +-
src/read_restart.h | 2 -
src/respa.cpp | 8 +
src/set.cpp | 16 +-
src/verlet.cpp | 8 +
70 files changed, 331 insertions(+), 1390 deletions(-)
diff --git a/src/CORESHELL/pair_lj_cut_coul_long_cs.cpp b/src/CORESHELL/pair_lj_cut_coul_long_cs.cpp
index 6f068f4be7..fe9ba28696 100644
--- a/src/CORESHELL/pair_lj_cut_coul_long_cs.cpp
+++ b/src/CORESHELL/pair_lj_cut_coul_long_cs.cpp
@@ -555,9 +555,10 @@ void PairLJCutCoulLongCS::compute_outer(int eflag, int vflag)
if (rsq <= cut_in_off_sq) {
r6inv = r2inv*r2inv*r2inv;
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
- } else if (rsq <= cut_in_on_sq)
+ } else if (rsq <= cut_in_on_sq) {
r6inv = r2inv*r2inv*r2inv;
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
+ }
fpair = (forcecoul + factor_lj*forcelj) * r2inv;
}
diff --git a/src/DIPOLE/pair_lj_cut_dipole_long.cpp b/src/DIPOLE/pair_lj_cut_dipole_long.cpp
index 60190ce905..824bbf734a 100755
--- a/src/DIPOLE/pair_lj_cut_dipole_long.cpp
+++ b/src/DIPOLE/pair_lj_cut_dipole_long.cpp
@@ -291,10 +291,10 @@ void PairLJCutDipoleLong::compute(int eflag, int vflag)
}
if (eflag) {
- if (rsq < cut_coulsq) {
+ if (rsq < cut_coulsq && factor_coul > 0.0) {
ecoul = qqrd2e*(b0*g0 + b1*g1 + b2*g2);
if (factor_coul < 1.0) {
- ecoul *= factor_coul;
+ ecoul *= factor_coul;
ecoul += (1-factor_coul) * qqrd2e * (d0*g0 + d1*g1 + d2*g2);
}
} else ecoul = 0.0;
diff --git a/src/KOKKOS/verlet_kokkos.cpp b/src/KOKKOS/verlet_kokkos.cpp
index b63ca98ba4..cdd3cecf8e 100644
--- a/src/KOKKOS/verlet_kokkos.cpp
+++ b/src/KOKKOS/verlet_kokkos.cpp
@@ -53,8 +53,24 @@ VerletKokkos::VerletKokkos(LAMMPS *lmp, int narg, char **arg) :
void VerletKokkos::setup()
{
+ if (comm->me == 0 && screen) {
+ fprintf(screen,"Setting up Verlet run ...\n");
+ fprintf(screen," Unit style : %s\n", update->unit_style);
+ fprintf(screen," Current step: " BIGINT_FORMAT "\n", update->ntimestep);
+ fprintf(screen," Time step : %g\n", update->dt);
+ if (update->max_wall > 0) {
+ char outtime[128];
+ double totalclock = update->max_wall;
+ int seconds = fmod(totalclock,60.0);
+ totalclock = (totalclock - seconds) / 60.0;
+ int minutes = fmod(totalclock,60.0);
+ int hours = (totalclock - minutes) / 60.0;
+ sprintf(outtime," Max walltime: "
+ "%d:%02d:%02d\n", hours, minutes, seconds);
+ fputs(outtime,screen);
+ }
+ }
- if (comm->me == 0 && screen) fprintf(screen,"Setting up run ...\n");
update->setupflag = 1;
// setup domain, communication and neighboring
diff --git a/src/KSPACE/pair_lj_cut_coul_long.cpp b/src/KSPACE/pair_lj_cut_coul_long.cpp
index aa19c7f533..c113aea83a 100644
--- a/src/KSPACE/pair_lj_cut_coul_long.cpp
+++ b/src/KSPACE/pair_lj_cut_coul_long.cpp
@@ -550,10 +550,10 @@ void PairLJCutCoulLong::compute_outer(int eflag, int vflag)
if (rsq <= cut_in_off_sq) {
r6inv = r2inv*r2inv*r2inv;
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
- } else if (rsq <= cut_in_on_sq)
+ } else if (rsq <= cut_in_on_sq) {
r6inv = r2inv*r2inv*r2inv;
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
-
+ }
fpair = (forcecoul + factor_lj*forcelj) * r2inv;
}
diff --git a/src/KSPACE/pair_lj_cut_coul_msm.cpp b/src/KSPACE/pair_lj_cut_coul_msm.cpp
index a503b54aa3..f9fc9be3fc 100644
--- a/src/KSPACE/pair_lj_cut_coul_msm.cpp
+++ b/src/KSPACE/pair_lj_cut_coul_msm.cpp
@@ -403,10 +403,10 @@ void PairLJCutCoulMSM::compute_outer(int eflag, int vflag)
if (rsq <= cut_in_off_sq) {
r6inv = r2inv*r2inv*r2inv;
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
- } else if (rsq <= cut_in_on_sq)
+ } else if (rsq <= cut_in_on_sq) {
r6inv = r2inv*r2inv*r2inv;
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
-
+ }
fpair = (forcecoul + factor_lj*forcelj) * r2inv;
}
diff --git a/src/MANYBODY/pair_adp.cpp b/src/MANYBODY/pair_adp.cpp
index 369bcc8c9e..5cb861b563 100644
--- a/src/MANYBODY/pair_adp.cpp
+++ b/src/MANYBODY/pair_adp.cpp
@@ -44,6 +44,7 @@ PairADP::PairADP(LAMMPS *lmp) : Pair(lmp)
fp = NULL;
mu = NULL;
lambda = NULL;
+ map = NULL;
setfl = NULL;
diff --git a/src/MANYBODY/pair_airebo.cpp b/src/MANYBODY/pair_airebo.cpp
index 58bbeffbfc..f3a7be5a2d 100644
--- a/src/MANYBODY/pair_airebo.cpp
+++ b/src/MANYBODY/pair_airebo.cpp
@@ -59,6 +59,7 @@ PairAIREBO::PairAIREBO(LAMMPS *lmp) : Pair(lmp)
pgsize = oneatom = 0;
nC = nH = NULL;
+ map = NULL;
manybody_flag = 1;
}
diff --git a/src/MANYBODY/pair_comb.cpp b/src/MANYBODY/pair_comb.cpp
index da9c4a49aa..26ff69d564 100644
--- a/src/MANYBODY/pair_comb.cpp
+++ b/src/MANYBODY/pair_comb.cpp
@@ -56,6 +56,8 @@ PairComb::PairComb(LAMMPS *lmp) : Pair(lmp)
nmax = 0;
NCo = NULL;
bbij = NULL;
+ map = NULL;
+ esm = NULL;
nelements = 0;
elements = NULL;
diff --git a/src/MANYBODY/pair_comb3.cpp b/src/MANYBODY/pair_comb3.cpp
index 95e6239bce..2d1a782070 100644
--- a/src/MANYBODY/pair_comb3.cpp
+++ b/src/MANYBODY/pair_comb3.cpp
@@ -54,6 +54,8 @@ PairComb3::PairComb3(LAMMPS *lmp) : Pair(lmp)
nmax = 0;
NCo = NULL;
bbij = NULL;
+ map = NULL;
+ esm = NULL;
nelements = 0;
elements = NULL;
diff --git a/src/MANYBODY/pair_eam.cpp b/src/MANYBODY/pair_eam.cpp
index 82f69c76b6..5b832d0ab4 100644
--- a/src/MANYBODY/pair_eam.cpp
+++ b/src/MANYBODY/pair_eam.cpp
@@ -42,6 +42,8 @@ PairEAM::PairEAM(LAMMPS *lmp) : Pair(lmp)
nmax = 0;
rho = NULL;
fp = NULL;
+ map = NULL;
+ type2frho = NULL;
nfuncfl = 0;
funcfl = NULL;
diff --git a/src/MANYBODY/pair_eim.cpp b/src/MANYBODY/pair_eim.cpp
index 9206229e9d..196e899534 100644
--- a/src/MANYBODY/pair_eim.cpp
+++ b/src/MANYBODY/pair_eim.cpp
@@ -45,6 +45,7 @@ PairEIM::PairEIM(LAMMPS *lmp) : Pair(lmp)
nmax = 0;
rho = NULL;
fp = NULL;
+ map = NULL;
nelements = 0;
elements = NULL;
diff --git a/src/MANYBODY/pair_nb3b_harmonic.cpp b/src/MANYBODY/pair_nb3b_harmonic.cpp
index ef1432827d..a98bacd274 100644
--- a/src/MANYBODY/pair_nb3b_harmonic.cpp
+++ b/src/MANYBODY/pair_nb3b_harmonic.cpp
@@ -52,6 +52,7 @@ PairNb3bHarmonic::PairNb3bHarmonic(LAMMPS *lmp) : Pair(lmp)
nparams = maxparam = 0;
params = NULL;
elem2param = NULL;
+ map = NULL;
}
/* ----------------------------------------------------------------------
diff --git a/src/MANYBODY/pair_polymorphic.cpp b/src/MANYBODY/pair_polymorphic.cpp
index 67342d91d6..b9d3755e6c 100755
--- a/src/MANYBODY/pair_polymorphic.cpp
+++ b/src/MANYBODY/pair_polymorphic.cpp
@@ -51,7 +51,7 @@ PairPolymorphic::PairPolymorphic(LAMMPS *lmp) : Pair(lmp)
tripletParameters = NULL;
elem2param = NULL;
elem3param = NULL;
-
+ type_map = NULL;
}
/* ----------------------------------------------------------------------
diff --git a/src/MANYBODY/pair_sw.cpp b/src/MANYBODY/pair_sw.cpp
index 27394c7a67..5cdd8f778e 100755
--- a/src/MANYBODY/pair_sw.cpp
+++ b/src/MANYBODY/pair_sw.cpp
@@ -50,6 +50,7 @@ PairSW::PairSW(LAMMPS *lmp) : Pair(lmp)
nparams = maxparam = 0;
params = NULL;
elem2param = NULL;
+ map = NULL;
}
/* ----------------------------------------------------------------------
diff --git a/src/MANYBODY/pair_tersoff.cpp b/src/MANYBODY/pair_tersoff.cpp
index f725df0cc6..44712e691e 100755
--- a/src/MANYBODY/pair_tersoff.cpp
+++ b/src/MANYBODY/pair_tersoff.cpp
@@ -51,6 +51,7 @@ PairTersoff::PairTersoff(LAMMPS *lmp) : Pair(lmp)
nparams = maxparam = 0;
params = NULL;
elem2param = NULL;
+ map = NULL;
}
/* ----------------------------------------------------------------------
diff --git a/src/PYTHON/python.cpp b/src/PYTHON/python.cpp
index 8c34098126..607181a1ec 100644
--- a/src/PYTHON/python.cpp
+++ b/src/PYTHON/python.cpp
@@ -117,6 +117,7 @@ void Python::command(int narg, char **arg)
iarg += 2;
} else if (strcmp(arg[iarg],"file") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Invalid python command");
+ delete[] pyfile;
int n = strlen(arg[iarg+1]) + 1;
pyfile = new char[n];
strcpy(pyfile,arg[iarg+1]);
@@ -173,6 +174,7 @@ void Python::command(int narg, char **arg)
if (fp == NULL) error->all(FLERR,"Could not open Python file");
int err = PyRun_SimpleFile(fp,pyfile);
if (err) error->all(FLERR,"Could not process Python file");
+ fclose(fp);
} else if (herestr) {
int err = PyRun_SimpleString(herestr);
if (err) error->all(FLERR,"Could not process Python string");
diff --git a/src/REPLICA/verlet_split.cpp b/src/REPLICA/verlet_split.cpp
index 408821fe22..dbb2660eb5 100644
--- a/src/REPLICA/verlet_split.cpp
+++ b/src/REPLICA/verlet_split.cpp
@@ -241,7 +241,8 @@ void VerletSplit::init()
void VerletSplit::setup()
{
- if (comm->me == 0 && screen) fprintf(screen,"Setting up run ...\n");
+ if (comm->me == 0 && screen)
+ fprintf(screen,"Setting up Verlet/split run ...\n");
if (!master) force->kspace->setup();
else Verlet::setup();
@@ -298,6 +299,7 @@ void VerletSplit::run(int n)
int n_pre_exchange = modify->n_pre_exchange;
int n_pre_neighbor = modify->n_pre_neighbor;
int n_pre_force = modify->n_pre_force;
+ int n_pre_reverse = modify->n_pre_reverse;
int n_post_force = modify->n_post_force;
int n_end_of_step = modify->n_end_of_step;
@@ -374,6 +376,10 @@ void VerletSplit::run(int n)
timer->stamp(Timer::BOND);
}
+ if (n_pre_reverse) {
+ modify->pre_reverse(eflag,vflag);
+ timer->stamp(Timer::MODIFY);
+ }
if (force->newton) {
comm->reverse_comm();
timer->stamp(Timer::COMM);
@@ -391,6 +397,11 @@ void VerletSplit::run(int n)
timer->stamp(Timer::KSPACE);
}
+ if (n_pre_reverse) {
+ modify->pre_reverse(eflag,vflag);
+ timer->stamp(Timer::MODIFY);
+ }
+
// TIP4P PPPM puts forces on ghost atoms, so must reverse_comm()
if (tip4p_flag && force->newton) {
diff --git a/src/USER-AWPMD/pair_awpmd_cut.h b/src/USER-AWPMD/pair_awpmd_cut.h
index b07db3abbb..fb8b9c65dc 100644
--- a/src/USER-AWPMD/pair_awpmd_cut.h
+++ b/src/USER-AWPMD/pair_awpmd_cut.h
@@ -69,7 +69,7 @@ class PairAWPMDCut : public Pair {
void virial_eradius_compute();
- AWPMD_split *wpmd; // solver oblect
+ AWPMD_split *wpmd; // solver object
double ermscale; // scale of width mass for motion
double width_pbc; // setting for width pbc
double half_box_length; // calculated by coeff function
diff --git a/src/USER-CUDA/verlet_cuda.cpp b/src/USER-CUDA/verlet_cuda.cpp
index eba4a0b09c..2a216b5b46 100644
--- a/src/USER-CUDA/verlet_cuda.cpp
+++ b/src/USER-CUDA/verlet_cuda.cpp
@@ -118,7 +118,23 @@ void VerletCuda::setup()
cuda->allocate();
- if(comm->me == 0 && screen) fprintf(screen, "Setting up run ...\n");
+ if (comm->me == 0 && screen) {
+ fprintf(screen,"Setting up Verlet run ...\n");
+ fprintf(screen," Unit style : %s\n", update->unit_style);
+ fprintf(screen," Current step: " BIGINT_FORMAT "\n", update->ntimestep);
+ fprintf(screen," Time step : %g\n", update->dt);
+ if (update->max_wall > 0) {
+ char outtime[128];
+ double totalclock = update->max_wall;
+ int seconds = fmod(totalclock,60.0);
+ totalclock = (totalclock - seconds) / 60.0;
+ int minutes = fmod(totalclock,60.0);
+ int hours = (totalclock - minutes) / 60.0;
+ sprintf(outtime," Max walltime: "
+ "%d:%02d:%02d\n", hours, minutes, seconds);
+ fputs(outtime,screen);
+ }
+ }
// setup domain, communication and neighboring
// acquire ghosts
@@ -639,6 +655,11 @@ void VerletCuda::run(int n)
int firstreneigh = 1;
for(int i = 0; i < n; i++) {
+ if (update->time_expired()) {
+ update->nsteps = i;
+ break;
+ }
+
if(atom->nlocal == 0)
error->warning(FLERR, "# CUDA: There are currently no atoms on one of the MPI processes. This is currently prone to encountering errors with USER-CUDA package. Please use the 'processors' keyword to use a more balanced processor layout.");
diff --git a/src/USER-DIFFRACTION/compute_saed.h b/src/USER-DIFFRACTION/compute_saed.h
index ba221125a0..0118b64811 100644
--- a/src/USER-DIFFRACTION/compute_saed.h
+++ b/src/USER-DIFFRACTION/compute_saed.h
@@ -42,7 +42,6 @@ class ComputeSAED : public Compute {
double prd_inv[3]; // Inverse spacing of unit cell
bool echo; // echo compute_array progress
bool manual; // Turn on manual recpiprocal map
- double *f;
int nRows; // Number of relp explored
double Zone[3]; // Zone axis to view SAED
diff --git a/src/USER-DRUDE/pair_lj_cut_thole_long.cpp b/src/USER-DRUDE/pair_lj_cut_thole_long.cpp
index a31a3b29f7..d517eb7932 100644
--- a/src/USER-DRUDE/pair_lj_cut_thole_long.cpp
+++ b/src/USER-DRUDE/pair_lj_cut_thole_long.cpp
@@ -637,7 +637,7 @@ double PairLJCutTholeLong::single(int i, int j, int itype, int jtype,
if (drudetype[type[i]] != NOPOL_TYPE && drudetype[type[j]] != NOPOL_TYPE) {
di = atom->map(drudeid[i]);
di_closest = domain->closest_image(i, di);
- if (di_closest != dj){
+ if (j != di_closest){
if (drudetype[i] == CORE_TYPE) dqi = -atom->q[di];
else if (drudetype[i] == DRUDE_TYPE) dqi = atom->q[i];
if (drudetype[j] == CORE_TYPE) {
@@ -672,7 +672,7 @@ double PairLJCutTholeLong::single(int i, int j, int itype, int jtype,
}
if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor;
if (drudetype[type[i]] != NOPOL_TYPE && drudetype[type[j]] != NOPOL_TYPE &&
- di_closest != dj)
+ di_closest != j)
phicoul += factor_e * dcoul;
eng += phicoul;
}
diff --git a/src/USER-DRUDE/pair_thole.cpp b/src/USER-DRUDE/pair_thole.cpp
index f9a16267d0..5840dbc51b 100644
--- a/src/USER-DRUDE/pair_thole.cpp
+++ b/src/USER-DRUDE/pair_thole.cpp
@@ -364,7 +364,7 @@ double PairThole::single(int i, int j, int itype, int jtype,
double rsq, double factor_coul, double factor_lj,
double &fforce)
{
- double r2inv,rinv,r,forcecoul,phicoul;
+ double r2inv,rinv,r,phicoul;
double qi,qj,factor_f,factor_e,dcoul,asr,exp_asr;
int di, dj;
diff --git a/src/USER-FEP/pair_lj_charmm_coul_long_soft.cpp b/src/USER-FEP/pair_lj_charmm_coul_long_soft.cpp
index 81464e3606..1f4a2529f5 100644
--- a/src/USER-FEP/pair_lj_charmm_coul_long_soft.cpp
+++ b/src/USER-FEP/pair_lj_charmm_coul_long_soft.cpp
@@ -533,6 +533,8 @@ void PairLJCharmmCoulLongSoft::compute_outer(int eflag, int vflag)
} else ecoul = 0.0;
if (rsq < cut_ljsq) {
+ r4sig6 = rsq*rsq / lj2[itype][jtype];
+ denlj = lj3[itype][jtype] + rsq*r4sig6;
evdwl = lj1[itype][jtype] * 4.0 * epsilon[itype][jtype] *
(1.0/(denlj*denlj) - 1.0/denlj);
if (rsq > cut_lj_innersq) {
diff --git a/src/USER-FEP/pair_lj_cut_coul_long_soft.cpp b/src/USER-FEP/pair_lj_cut_coul_long_soft.cpp
index daeeffd28c..6bb9538023 100644
--- a/src/USER-FEP/pair_lj_cut_coul_long_soft.cpp
+++ b/src/USER-FEP/pair_lj_cut_coul_long_soft.cpp
@@ -498,6 +498,8 @@ void PairLJCutCoulLongSoft::compute_outer(int eflag, int vflag)
} else ecoul = 0.0;
if (rsq < cut_ljsq[itype][jtype]) {
+ r4sig6 = rsq*rsq / lj2[itype][jtype];
+ denlj = lj3[itype][jtype] + rsq*r4sig6;
evdwl = lj1[itype][jtype] * 4.0 * epsilon[itype][jtype] *
(1.0/(denlj*denlj) - 1.0/denlj) - offset[itype][jtype];
evdwl *= factor_lj;
diff --git a/src/USER-FEP/pair_lj_cut_soft.cpp b/src/USER-FEP/pair_lj_cut_soft.cpp
index 1979508224..95b74d0d1f 100644
--- a/src/USER-FEP/pair_lj_cut_soft.cpp
+++ b/src/USER-FEP/pair_lj_cut_soft.cpp
@@ -393,6 +393,7 @@ void PairLJCutSoft::compute_outer(int eflag, int vflag)
}
if (eflag) {
+ r4sig6 = rsq*rsq / lj2[itype][jtype];
denlj = lj3[itype][jtype] + rsq*r4sig6;
evdwl = lj1[itype][jtype] * 4.0 * epsilon[itype][jtype] *
(1.0/(denlj*denlj) - 1.0/denlj) - offset[itype][jtype];
diff --git a/src/USER-H5MD/dump_h5md.cpp b/src/USER-H5MD/dump_h5md.cpp
index 1c35c0ca52..4af1e2ed4c 100644
--- a/src/USER-H5MD/dump_h5md.cpp
+++ b/src/USER-H5MD/dump_h5md.cpp
@@ -36,7 +36,7 @@ using namespace LAMMPS_NS;
/** Scan common options for the dump elements
*/
-int element_args(int narg, char **arg, int *every)
+static int element_args(int narg, char **arg, int *every)
{
int iarg=0;
while (iarg<narg) {
diff --git a/src/USER-H5MD/dump_h5md.h b/src/USER-H5MD/dump_h5md.h
index 464733419b..496c30441e 100644
--- a/src/USER-H5MD/dump_h5md.h
+++ b/src/USER-H5MD/dump_h5md.h
@@ -33,7 +33,6 @@ class DumpH5MD : public Dump {
private:
int natoms,ntotal;
- int nevery_save;
int unwrap_flag; // 1 if atom coords are unwrapped, 0 if no
h5md_file datafile;
int datafile_from_dump;
diff --git a/src/USER-INTEL/verlet_intel.cpp b/src/USER-INTEL/verlet_intel.cpp
index 039e3bc36e..4eba2d4516 100644
--- a/src/USER-INTEL/verlet_intel.cpp
+++ b/src/USER-INTEL/verlet_intel.cpp
@@ -97,7 +97,23 @@ void VerletIntel::init()
void VerletIntel::setup()
{
- if (comm->me == 0 && screen) fprintf(screen,"Setting up run ...\n");
+ if (comm->me == 0 && screen) {
+ fprintf(screen,"Setting up Verlet run ...\n");
+ fprintf(screen," Unit style : %s\n", update->unit_style);
+ fprintf(screen," Current step: " BIGINT_FORMAT "\n", update->ntimestep);
+ fprintf(screen," Time step : %g\n", update->dt);
+ if (update->max_wall > 0) {
+ char outtime[128];
+ double totalclock = update->max_wall;
+ int seconds = fmod(totalclock,60.0);
+ totalclock = (totalclock - seconds) / 60.0;
+ int minutes = fmod(totalclock,60.0);
+ int hours = (totalclock - minutes) / 60.0;
+ sprintf(outtime," Max walltime: "
+ "%d:%02d:%02d\n", hours, minutes, seconds);
+ fputs(outtime,screen);
+ }
+ }
update->setupflag = 1;
@@ -266,6 +282,10 @@ void VerletIntel::run(int n)
else sortflag = 0;
for (int i = 0; i < n; i++) {
+ if (update->time_expired()) {
+ update->nsteps = i;
+ return;
+ }
ntimestep = ++update->ntimestep;
ev_set(ntimestep);
diff --git a/src/USER-MISC/fix_imd.cpp b/src/USER-MISC/fix_imd.cpp
index a5c411bf88..b5bd8caf0b 100644
--- a/src/USER-MISC/fix_imd.cpp
+++ b/src/USER-MISC/fix_imd.cpp
@@ -1205,6 +1205,7 @@ void * imdsock_create(void) {
s = (imdsocket *) malloc(sizeof(imdsocket));
if (s != NULL)
memset(s, 0, sizeof(imdsocket));
+ else return NULL;
if ((s->sd = socket(PF_INET, SOCK_STREAM, 0)) == -1) {
printf("Failed to open socket.");
@@ -1344,19 +1345,6 @@ int imdsock_selwrite(void *v, int sec) {
/*************************************************************************/
/* start of imd API code. */
-/* Only works with aligned 4-byte quantities, will cause a bus error */
-/* on some platforms if used on unaligned data. */
-void swap4_aligned(void *v, long ndata) {
- int *data = (int *) v;
- long i;
- int *N;
- for (i=0; i<ndata; i++) {
- N = data + i;
- *N=(((*N>>24)&0xff) | ((*N&0xff)<<24) |
- ((*N>>8)&0xff00) | ((*N&0xff00)<<8));
- }
-}
-
/** structure used to perform byte swapping operations */
typedef union {
@@ -1431,12 +1419,6 @@ static int32 imd_writen(void *s, const char *ptr, int32 n) {
return n;
}
-int imd_disconnect(void *s) {
- IMDheader header;
- imd_fill_header(&header, IMD_DISCONNECT, 0);
- return (imd_writen(s, (char *)&header, IMDHEADERSIZE) != IMDHEADERSIZE);
-}
-
int imd_handshake(void *s) {
IMDheader header;
imd_fill_header(&header, IMD_HANDSHAKE, 1);
diff --git a/src/USER-MISC/fix_pimd.cpp b/src/USER-MISC/fix_pimd.cpp
index 5ade11ce26..75662a8207 100644
--- a/src/USER-MISC/fix_pimd.cpp
+++ b/src/USER-MISC/fix_pimd.cpp
@@ -659,8 +659,11 @@ void FixPIMD::comm_exec(double **ptr)
{
char error_line[256];
- sprintf(error_line, "Atom %d is missing at world [%d] rank [%d] required by rank [%d] (%d, %d, %d).\n",
- tag_send[i], universe->iworld, comm->me, plan_recv[iplan], atom->tag[0], atom->tag[1], atom->tag[2]);
+ sprintf(error_line, "Atom " TAGINT_FORMAT " is missing at world [%d] "
+ "rank [%d] required by rank [%d] (" TAGINT_FORMAT ", "
+ TAGINT_FORMAT ", " TAGINT_FORMAT ").\n",tag_send[i],
+ universe->iworld, comm->me, plan_recv[iplan],
+ atom->tag[0], atom->tag[1], atom->tag[2]);
error->universe_one(FLERR,error_line);
}
diff --git a/src/USER-MISC/fix_srp.cpp b/src/USER-MISC/fix_srp.cpp
index 613204a55e..08224d2950 100644
--- a/src/USER-MISC/fix_srp.cpp
+++ b/src/USER-MISC/fix_srp.cpp
@@ -152,7 +152,8 @@ void FixSRP::setup_pre_force(int zz)
double delx, dely, delz, rmax, rsq, rsqmax;
double **x = atom->x;
bigint nall = atom->nlocal + atom->nghost;
- double xold[nall][3];
+ double **xold;
+ memory->create(xold,nall,3,"fix_srp:xold");
// make a copy of all coordinates and tags
// need this because create_atom overwrites ghost atoms
@@ -163,7 +164,8 @@ void FixSRP::setup_pre_force(int zz)
}
tagint *tag = atom->tag;
- tagint tagold[nall];
+ tagint *tagold;
+ memory->create(tagold,nall,"fix_srp:tagold");
for(int i = 0; i < nall; i++){
tagold[i]=tag[i];
@@ -212,7 +214,11 @@ void FixSRP::setup_pre_force(int zz)
array[atom->nlocal-1][1] = (double)tagold[j];
nadd++;
}
- }
+ }
+
+ // free temporary storage
+ memory->destroy(xold);
+ memory->destroy(tagold);
// new total # of atoms
bigint nblocal = atom->nlocal;
diff --git a/src/USER-MISC/fix_ti_rs.cpp b/src/USER-MISC/fix_ti_rs.cpp
index f80fc4cb17..d5036eb520 100755
--- a/src/USER-MISC/fix_ti_rs.cpp
+++ b/src/USER-MISC/fix_ti_rs.cpp
@@ -25,6 +25,7 @@
#include "update.h"
#include "respa.h"
#include "error.h"
+#include "force.h"
using namespace LAMMPS_NS;
using namespace FixConst;
@@ -45,18 +46,18 @@ FixTIRS::FixTIRS(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
extvector = 1;
// Time variables.
- t_switch = atoi(arg[5]);
- t_equil = atoi(arg[6]);
+ t_switch = force->bnumeric(FLERR,arg[5]);
+ t_equil = force->bnumeric(FLERR,arg[6]);
t0 = update->ntimestep;
- if (t_switch < 0.0) error->all(FLERR,"Illegal fix ti/rs command");
- if (t_equil < 0.0) error->all(FLERR,"Illegal fix ti/rs command");
+ if (t_switch <= 0) error->all(FLERR,"Illegal fix ti/rs command");
+ if (t_equil <= 0) error->all(FLERR,"Illegal fix ti/rs command");
// Coupling parameter limits and initialization.
- l_initial = atof(arg[3]);
- l_final = atof(arg[4]);
+ l_initial = force->numeric(FLERR,arg[3]);
+ l_final = force->numeric(FLERR,arg[4]);
sf = 1;
if (narg > 7) {
- if (strcmp(arg[7], "function") == 0) sf = atoi(arg[8]);
+ if (strcmp(arg[7], "function") == 0) sf = force->inumeric(FLERR,arg[8]);
else error->all(FLERR,"Illegal fix ti/rs switching function");
if ((sf<1) || (sf>3))
error->all(FLERR,"Illegal fix ti/rs switching function");
@@ -151,16 +152,17 @@ void FixTIRS::min_post_force(int vflag)
void FixTIRS::initial_integrate(int vflag)
{
// Update the coupling parameter value.
- double t = update->ntimestep - (t0+t_equil);
+ const bigint t = update->ntimestep - (t0+t_equil);
+ const double r_switch = 1.0/t_switch;
if( (t >= 0) && (t <= t_switch) ) {
- lambda = switch_func(t/t_switch);
- dlambda = dswitch_func(t/t_switch);
+ lambda = switch_func(t*r_switch);
+ dlambda = dswitch_func(t*r_switch);
}
if( (t >= t_equil+t_switch) && (t <= (t_equil+2*t_switch)) ) {
- lambda = switch_func(1.0 - (t - t_switch - t_equil)/t_switch);
- dlambda = - dswitch_func(1.0 - (t - t_switch - t_equil)/t_switch);
+ lambda = switch_func(1.0 - (t - t_switch - t_equil)*r_switch);
+ dlambda = - dswitch_func(1.0 - (t - t_switch - t_equil)*r_switch);
}
}
diff --git a/src/USER-MISC/fix_ti_rs.h b/src/USER-MISC/fix_ti_rs.h
index 58466501a9..6b5124f711 100755
--- a/src/USER-MISC/fix_ti_rs.h
+++ b/src/USER-MISC/fix_ti_rs.h
@@ -53,9 +53,9 @@ class FixTIRS : public Fix {
double l_initial; // Lambda initial value.
double l_final; // Lambda final value.
double linfo[2]; // Current lambda status.
- int t_switch; // Total switching steps.
- int t_equil; // Equilibration time.
- int t0; // Initial time.
+ bigint t_switch; // Total switching steps.
+ bigint t_equil; // Equilibration time.
+ bigint t0; // Initial time.
int sf; // Switching function option.
int nlevels_respa;
};
diff --git a/src/USER-MISC/fix_ti_spring.cpp b/src/USER-MISC/fix_ti_spring.cpp
index 529e650978..3d10312720 100755
--- a/src/USER-MISC/fix_ti_spring.cpp
+++ b/src/USER-MISC/fix_ti_spring.cpp
@@ -73,16 +73,16 @@ FixTISpring::FixTISpring(LAMMPS *lmp, int narg, char **arg) :
}
// Time variables.
- t_switch = atoi(arg[4]); // Switching time.
- t_equil = atoi(arg[5]); // Equilibration time.
+ t_switch = force->bnumeric(FLERR,arg[4]); // Number of steps for switching.
+ t_equil = force->bnumeric(FLERR,arg[5]); // Number of steps for equilibration.
t0 = update->ntimestep; // Initial time.
- if (t_switch <= 0.0) error->all(FLERR,"Illegal fix ti/spring command");
- if (t_equil <= 0.0) error->all(FLERR,"Illegal fix ti/spring command");
+ if (t_switch <= 0) error->all(FLERR,"Illegal fix ti/spring command");
+ if (t_equil <= 0) error->all(FLERR,"Illegal fix ti/spring command");
// Coupling parameter initialization.
sf = 1;
if (narg > 6) {
- if (strcmp(arg[6], "function") == 0) sf = atoi(arg[7]);
+ if (strcmp(arg[6], "function") == 0) sf = force->inumeric(FLERR,arg[7]);
else error->all(FLERR,"Illegal fix ti/spring switching function");
if ((sf!=1) && (sf!=2))
error->all(FLERR,"Illegal fix ti/spring switching function");
@@ -151,7 +151,7 @@ void FixTISpring::min_setup(int vflag)
void FixTISpring::post_force(int vflag)
{
// If on the first equilibration do not calculate forces.
- int t = update->ntimestep - t0;
+ bigint t = update->ntimestep - t0;
if(t < t_equil) return;
double **x = atom->x;
@@ -199,16 +199,17 @@ void FixTISpring::min_post_force(int vflag)
void FixTISpring::initial_integrate(int vflag)
{
// Update the coupling parameter value.
- double t = update->ntimestep - (t0+t_equil);
+ const bigint t = update->ntimestep - (t0+t_equil);
+ const double r_switch = 1.0/t_switch;
if( (t >= 0) && (t <= t_switch) ) {
- lambda = switch_func(t/t_switch);
- dlambda = dswitch_func(t/t_switch);
+ lambda = switch_func(t*r_switch);
+ dlambda = dswitch_func(t*r_switch);
}
if( (t >= t_equil+t_switch) && (t <= (t_equil+2*t_switch)) ) {
- lambda = switch_func(1.0 - (t - t_switch - t_equil)/t_switch );
- dlambda = - dswitch_func(1.0 - (t - t_switch - t_equil)/t_switch );
+ lambda = switch_func(1.0 - (t - t_switch - t_equil)*r_switch);
+ dlambda = - dswitch_func(1.0 - (t - t_switch - t_equil)*r_switch);
}
}
diff --git a/src/USER-MISC/fix_ti_spring.h b/src/USER-MISC/fix_ti_spring.h
index ce67968d56..a5ff1fb677 100755
--- a/src/USER-MISC/fix_ti_spring.h
+++ b/src/USER-MISC/fix_ti_spring.h
@@ -65,9 +65,9 @@ class FixTISpring : public Fix {
double lambda; // Coupling parameter.
double dlambda; // Lambda variation with t.
double linfo[2]; // Current lambda status.
- int t_switch; // Total switching steps.
- int t_equil; // Equilibration time.
- int t0; // Initial time.
+ bigint t_switch; // Total switching steps.
+ bigint t_equil; // Equilibration time.
+ bigint t0; // Initial time.
int sf; // Switching function option.
int nlevels_respa;
};
diff --git a/src/USER-MISC/pair_tersoff_table.cpp b/src/USER-MISC/pair_tersoff_table.cpp
index 77002a5ba3..c5865c0a99 100644
--- a/src/USER-MISC/pair_tersoff_table.cpp
+++ b/src/USER-MISC/pair_tersoff_table.cpp
@@ -306,10 +306,6 @@ void PairTersoffTable::compute(int eflag, int vflag)
if (r_ik > params[ikparam].cutsq) continue;
- r_ik = sqrt(r_ik);
-
- invR_ik = 1.0 / r_ik;
-
gtetaFunctionIJK = preGtetaFunction[neighbor_j][neighbor_k];
cutoffFunctionIK = preCutoffFunction[neighbor_k];
@@ -335,13 +331,6 @@ void PairTersoffTable::compute(int eflag, int vflag)
if (r_ik > params[ikparam].cutsq) continue;
- r_ik = sqrt(r_ik);
- invR_ik = 1.0 / r_ik;
-
- directorCos_ik_x = invR_ik * dr_ik[0];
- directorCos_ik_y = invR_ik * dr_ik[1];
- directorCos_ik_z = invR_ik * dr_ik[2];
-
gtetaFunctionIJK = preGtetaFunction[neighbor_j][neighbor_k];
cutoffFunctionIK = preCutoffFunction[neighbor_k];
diff --git a/src/USER-MOLFILE/dump_molfile.cpp b/src/USER-MOLFILE/dump_molfile.cpp
index add7bf695b..54e86f3bc5 100644
--- a/src/USER-MOLFILE/dump_molfile.cpp
+++ b/src/USER-MOLFILE/dump_molfile.cpp
@@ -387,6 +387,7 @@ void DumpMolfile::write_data(int n, double *mybuf)
if (need_structure) {
mf->property(MFI::P_NAME,types,typenames);
+ mf->property(MFI::P_TYPE,types,typenames);
if (atom->molecule_flag)
mf->property(MFI::P_RESI,molids);
diff --git a/src/USER-OMP/fix_omp.cpp b/src/USER-OMP/fix_omp.cpp
index 2b3e2de299..7cb63bcc16 100644
--- a/src/USER-OMP/fix_omp.cpp
+++ b/src/USER-OMP/fix_omp.cpp
@@ -147,13 +147,9 @@ FixOMP::FixOMP(LAMMPS *lmp, int narg, char **arg)
FixOMP::~FixOMP()
{
-#if defined(_OPENMP)
-#pragma omp parallel default(none)
-#endif
- {
- const int tid = get_tid();
- delete thr[tid];
- }
+ for (int i=0; i < _nthr; ++i)
+ delete thr[i];
+
delete[] thr;
}
@@ -189,8 +185,30 @@ void FixOMP::init()
error->all(FLERR,"USER-OMP package does not (yet) work with "
"atom_style template");
+ // adjust number of data objects when the number of OpenMP
+ // threads has been changed somehow
+ const int nthreads = comm->nthreads;
+ if (_nthr != nthreads) {
+ if (screen) fprintf(screen,"Re-init USER-OMP for %d OpenMP thread(s)\n", nthreads);
+ if (logfile) fprintf(logfile,"Re-init USER-OMP for %d OpenMP thread(s)\n", nthreads);
+
+ for (int i=0; i < _nthr; ++i)
+ delete thr[i];
+
+ thr = new ThrData *[nthreads];
+ _nthr = nthreads;
+#if defined(_OPENMP)
+#pragma omp parallel default(none)
+#endif
+ {
+ const int tid = get_tid();
+ Timer *t = new Timer(lmp);
+ thr[tid] = new ThrData(tid,t);
+ }
+ }
+
// reset per thread timer
- for (int i=0; i < comm->nthreads; ++i) {
+ for (int i=0; i < nthreads; ++i) {
thr[i]->_timer_active=1;
thr[i]->timer(Timer::RESET);
thr[i]->_timer_active=-1;
@@ -323,7 +341,7 @@ void FixOMP::set_neighbor_omp()
void FixOMP::setup(int)
{
// we are post the force compute in setup. turn on timers
- for (int i=0; i < comm->nthreads; ++i)
+ for (int i=0; i < _nthr; ++i)
thr[i]->_timer_active=0;
}
@@ -356,8 +374,8 @@ void FixOMP::pre_force(int)
double FixOMP::memory_usage()
{
- double bytes = comm->nthreads * (sizeof(ThrData *) + sizeof(ThrData));
- bytes += comm->nthreads * thr[0]->memory_usage();
+ double bytes = _nthr * (sizeof(ThrData *) + sizeof(ThrData));
+ bytes += _nthr * thr[0]->memory_usage();
return bytes;
}
diff --git a/src/USER-OMP/pair_tersoff_table_omp.cpp b/src/USER-OMP/pair_tersoff_table_omp.cpp
index bd786c7ca9..ece331abd9 100644
--- a/src/USER-OMP/pair_tersoff_table_omp.cpp
+++ b/src/USER-OMP/pair_tersoff_table_omp.cpp
@@ -297,10 +297,6 @@ void PairTersoffTableOMP::eval(int iifrom, int iito, ThrData * const thr)
if (r_ik > params[ikparam].cutsq) continue;
- r_ik = sqrt(r_ik);
-
- invR_ik = 1.0 / r_ik;
-
gtetaFunctionIJK = thrGtetaFunction[tid][neighbor_j][neighbor_k];
cutoffFunctionIK = thrCutoffFunction[tid][neighbor_k];
@@ -326,13 +322,6 @@ void PairTersoffTableOMP::eval(int iifrom, int iito, ThrData * const thr)
if (r_ik > params[ikparam].cutsq) continue;
- r_ik = sqrt(r_ik);
- invR_ik = 1.0 / r_ik;
-
- directorCos_ik_x = invR_ik * dr_ik[0];
- directorCos_ik_y = invR_ik * dr_ik[1];
- directorCos_ik_z = invR_ik * dr_ik[2];
-
gtetaFunctionIJK = thrGtetaFunction[tid][neighbor_j][neighbor_k];
cutoffFunctionIK = thrCutoffFunction[tid][neighbor_k];
diff --git a/src/USER-OMP/respa_omp.cpp b/src/USER-OMP/respa_omp.cpp
index ed08f019fb..f5603adf92 100644
--- a/src/USER-OMP/respa_omp.cpp
+++ b/src/USER-OMP/respa_omp.cpp
@@ -74,6 +74,17 @@ void RespaOMP::setup()
fprintf(screen," Unit style : %s\n", update->unit_style);
fprintf(screen," Current step : " BIGINT_FORMAT "\n", update->ntimestep);
fprintf(screen," OuterTime step: %g\n", update->dt);
+ if (update->max_wall > 0) {
+ char outtime[128];
+ double totalclock = update->max_wall;
+ int seconds = fmod(totalclock,60.0);
+ totalclock = (totalclock - seconds) / 60.0;
+ int minutes = fmod(totalclock,60.0);
+ int hours = (totalclock - minutes) / 60.0;
+ sprintf(outtime," Max walltime: "
+ "%d:%02d:%02d\n", hours, minutes, seconds);
+ fputs(outtime,screen);
+ }
}
update->setupflag = 1;
@@ -150,6 +161,7 @@ void RespaOMP::setup()
fix->did_reduce();
}
+ modify->pre_reverse(eflag,vflag);
if (newton[ilevel]) comm->reverse_comm();
copy_f_flevel(ilevel);
}
@@ -243,6 +255,7 @@ void RespaOMP::setup_minimal(int flag)
fix->did_reduce();
}
+ modify->pre_reverse(eflag,vflag);
if (newton[ilevel]) comm->reverse_comm();
copy_f_flevel(ilevel);
}
@@ -391,6 +404,10 @@ void RespaOMP::recurse(int ilevel)
fix->did_reduce();
}
+ if (modify->n_pre_reverse) {
+ modify->pre_reverse(eflag,vflag);
+ timer->stamp(Timer::MODIFY);
+ }
if (newton[ilevel]) {
comm->reverse_comm();
timer->stamp(Timer::COMM);
diff --git a/src/USER-REAXC/reaxc_bond_orders.cpp b/src/USER-REAXC/reaxc_bond_orders.cpp
index e524c75e87..0b4ca21adf 100644
--- a/src/USER-REAXC/reaxc_bond_orders.cpp
+++ b/src/USER-REAXC/reaxc_bond_orders.cpp
@@ -359,12 +359,6 @@ int BOp( storage *workspace, reax_list *bonds, double bo_cut,
}
-int compare_bonds( const void *p1, const void *p2 )
-{
- return ((bond_data *)p1)->nbr - ((bond_data *)p2)->nbr;
-}
-
-
void BO( reax_system *system, control_params *control, simulation_data *data,
storage *workspace, reax_list **lists, output_controls *out_control )
{
diff --git a/src/USER-REAXC/reaxc_bond_orders.h b/src/USER-REAXC/reaxc_bond_orders.h
index 8f1e182aa2..3631d90c89 100644
--- a/src/USER-REAXC/reaxc_bond_orders.h
+++ b/src/USER-REAXC/reaxc_bond_orders.h
@@ -36,23 +36,6 @@ typedef struct{
double C1dDelta, C2dDelta, C3dDelta;
}dbond_coefficients;
-#ifdef TEST_FORCES
-void Get_dBO( reax_system*, reax_list**, int, int, double, rvec* );
-void Get_dBOpinpi2( reax_system*, reax_list**,
- int, int, double, double, rvec*, rvec* );
-
-void Add_dBO( reax_system*, reax_list**, int, int, double, rvec* );
-void Add_dBOpinpi2( reax_system*, reax_list**,
- int, int, double, double, rvec*, rvec* );
-
-void Add_dBO_to_Forces( reax_system*, reax_list**, int, int, double );
-void Add_dBOpinpi2_to_Forces( reax_system*, reax_list**,
- int, int, double, double );
-
-void Add_dDelta( reax_system*, reax_list**, int, double, rvec* );
-void Add_dDelta_to_Forces( reax_system *, reax_list**, int, double );
-#endif
-
void Add_dBond_to_Forces( reax_system*, int, int, storage*, reax_list** );
void Add_dBond_to_Forces_NPT( int, int, simulation_data*,
storage*, reax_list** );
diff --git a/src/USER-REAXC/reaxc_ffield.cpp b/src/USER-REAXC/reaxc_ffield.cpp
index 34ac121856..6007cb7344 100644
--- a/src/USER-REAXC/reaxc_ffield.cpp
+++ b/src/USER-REAXC/reaxc_ffield.cpp
@@ -60,6 +60,8 @@ char Read_Force_Field( FILE *fp, reax_interaction *reax,
if (n < 1) {
fprintf( stderr, "WARNING: number of globals in ffield file is 0!\n" );
fclose(fp);
+ free(s);
+ free(tmp);
return 1;
}
diff --git a/src/USER-REAXC/reaxc_forces.cpp b/src/USER-REAXC/reaxc_forces.cpp
index 84bb86e5d2..7f11f5565f 100644
--- a/src/USER-REAXC/reaxc_forces.cpp
+++ b/src/USER-REAXC/reaxc_forces.cpp
@@ -171,213 +171,6 @@ void Validate_Lists( reax_system *system, storage *workspace, reax_list **lists,
}
-double Compute_H( double r, double gamma, double *ctap )
-{
- double taper, dr3gamij_1, dr3gamij_3;
-
- taper = ctap[7] * r + ctap[6];
- taper = taper * r + ctap[5];
- taper = taper * r + ctap[4];
- taper = taper * r + ctap[3];
- taper = taper * r + ctap[2];
- taper = taper * r + ctap[1];
- taper = taper * r + ctap[0];
-
- dr3gamij_1 = ( r*r*r + gamma );
- dr3gamij_3 = pow( dr3gamij_1 , 0.33333333333333 );
- return taper * EV_to_KCALpMOL / dr3gamij_3;
-}
-
-
-double Compute_tabH( double r_ij, int ti, int tj )
-{
- int r, tmin, tmax;
- double val, dif, base;
- LR_lookup_table *t;
-
- tmin = MIN( ti, tj );
- tmax = MAX( ti, tj );
- t = &( LR[tmin][tmax] );
-
- /* cubic spline interpolation */
- r = (int)(r_ij * t->inv_dx);
- if( r == 0 ) ++r;
- base = (double)(r+1) * t->dx;
- dif = r_ij - base;
- val = ((t->ele[r].d*dif + t->ele[r].c)*dif + t->ele[r].b)*dif +
- t->ele[r].a;
- val *= EV_to_KCALpMOL / C_ele;
-
- return val;
-}
-
-
-void Init_Forces( reax_system *system, control_params *control,
- simulation_data *data, storage *workspace, reax_list **lists,
- output_controls *out_control, MPI_Comm comm ) {
- int i, j, pj;
- int start_i, end_i;
- int type_i, type_j;
- int Htop, btop_i, btop_j, num_bonds, num_hbonds;
- int ihb, jhb, ihb_top, jhb_top;
- int local, flag, renbr;
- double r_ij, cutoff;
- sparse_matrix *H;
- reax_list *far_nbrs, *bonds, *hbonds;
- single_body_parameters *sbp_i, *sbp_j;
- two_body_parameters *twbp;
- far_neighbor_data *nbr_pj;
- reax_atom *atom_i, *atom_j;
-
- far_nbrs = *lists + FAR_NBRS;
- bonds = *lists + BONDS;
- hbonds = *lists + HBONDS;
-
- for( i = 0; i < system->n; ++i )
- workspace->bond_mark[i] = 0;
- for( i = system->n; i < system->N; ++i ) {
- workspace->bond_mark[i] = 1000; // put ghost atoms to an infinite distance
- }
-
- H = workspace->H;
- H->n = system->n;
- Htop = 0;
- num_bonds = 0;
- num_hbonds = 0;
- btop_i = btop_j = 0;
- renbr = (data->step-data->prev_steps) % control->reneighbor == 0;
-
- for( i = 0; i < system->N; ++i ) {
- atom_i = &(system->my_atoms[i]);
- type_i = atom_i->type;
- if (type_i < 0) continue;
- start_i = Start_Index(i, far_nbrs);
- end_i = End_Index(i, far_nbrs);
- btop_i = End_Index( i, bonds );
- sbp_i = &(system->reax_param.sbp[type_i]);
-
- if( i < system->n ) {
- local = 1;
- cutoff = control->nonb_cut;
- }
- else {
- local = 0;
- cutoff = control->bond_cut;
- }
-
- ihb = -1;
- ihb_top = -1;
- if( local ) {
- H->start[i] = Htop;
- H->entries[Htop].j = i;
- H->entries[Htop].val = sbp_i->eta;
- ++Htop;
-
- if( control->hbond_cut > 0 ) {
- ihb = sbp_i->p_hbond;
- if( ihb == 1 )
- ihb_top = End_Index( atom_i->Hindex, hbonds );
- else ihb_top = -1;
- }
- }
-
- /* update i-j distance - check if j is within cutoff */
- for( pj = start_i; pj < end_i; ++pj ) {
- nbr_pj = &( far_nbrs->select.far_nbr_list[pj] );
- j = nbr_pj->nbr;
- atom_j = &(system->my_atoms[j]);
- if( renbr ) {
- if(nbr_pj->d <= cutoff)
- flag = 1;
- else flag = 0;
- }
- else{
- nbr_pj->dvec[0] = atom_j->x[0] - atom_i->x[0];
- nbr_pj->dvec[1] = atom_j->x[1] - atom_i->x[1];
- nbr_pj->dvec[2] = atom_j->x[2] - atom_i->x[2];
- nbr_pj->d = rvec_Norm_Sqr( nbr_pj->dvec );
- if( nbr_pj->d <= SQR(cutoff) ) {
- nbr_pj->d = sqrt(nbr_pj->d);
- flag = 1;
- }
- else {
- flag = 0;
- }
- }
-
- if( flag ){
- type_j = atom_j->type;
- if (type_j < 0) continue;
- r_ij = nbr_pj->d;
- sbp_j = &(system->reax_param.sbp[type_j]);
- twbp = &(system->reax_param.tbp[type_i][type_j]);
-
- if( local ) {
- /* H matrix entry */
- if( j < system->n || atom_i->orig_id < atom_j->orig_id ) {//tryQEq||1
- H->entries[Htop].j = j;
- if( control->tabulate == 0 )
- H->entries[Htop].val = Compute_H(r_ij,twbp->gamma,workspace->Tap);
- else H->entries[Htop].val = Compute_tabH(r_ij, type_i, type_j);
- ++Htop;
- }
-
- /* hydrogen bond lists */
- if( control->hbond_cut > 0 && (ihb==1 || ihb==2) &&
- nbr_pj->d <= control->hbond_cut ) {
- jhb = sbp_j->p_hbond;
- if( ihb == 1 && jhb == 2 ) {
- hbonds->select.hbond_list[ihb_top].nbr = j;
- hbonds->select.hbond_list[ihb_top].scl = 1;
- hbonds->select.hbond_list[ihb_top].ptr = nbr_pj;
- ++ihb_top;
- ++num_hbonds;
- }
- else if( j < system->n && ihb == 2 && jhb == 1 ) {
- jhb_top = End_Index( atom_j->Hindex, hbonds );
- hbonds->select.hbond_list[jhb_top].nbr = i;
- hbonds->select.hbond_list[jhb_top].scl = -1;
- hbonds->select.hbond_list[jhb_top].ptr = nbr_pj;
- Set_End_Index( atom_j->Hindex, jhb_top+1, hbonds );
- ++num_hbonds;
- }
- }
- }
-
- /* uncorrected bond orders */
- if( //(workspace->bond_mark[i] < 3 || workspace->bond_mark[j] < 3) &&
- nbr_pj->d <= control->bond_cut &&
- BOp( workspace, bonds, control->bo_cut,
- i , btop_i, nbr_pj, sbp_i, sbp_j, twbp ) ) {
- num_bonds += 2;
- ++btop_i;
-
- if( workspace->bond_mark[j] > workspace->bond_mark[i] + 1 )
- workspace->bond_mark[j] = workspace->bond_mark[i] + 1;
- else if( workspace->bond_mark[i] > workspace->bond_mark[j] + 1 ) {
- workspace->bond_mark[i] = workspace->bond_mark[j] + 1;
- }
- }
- }
- }
-
- Set_End_Index( i, btop_i, bonds );
- if( local ) {
- H->end[i] = Htop;
- if( ihb == 1 )
- Set_End_Index( atom_i->Hindex, ihb_top, hbonds );
- }
- }
-
- workspace->realloc.Htop = Htop;
- workspace->realloc.num_bonds = num_bonds;
- workspace->realloc.num_hbonds = num_hbonds;
-
- Validate_Lists( system, workspace, lists, data->step,
- system->n, system->N, system->numH, comm );
-}
-
-
void Init_Forces_noQEq( reax_system *system, control_params *control,
simulation_data *data, storage *workspace,
reax_list **lists, output_controls *out_control,
@@ -647,19 +440,11 @@ void Compute_Forces( reax_system *system, control_params *control,
reax_list **lists, output_controls *out_control,
mpi_datatypes *mpi_data )
{
- MPI_Comm comm;
- int qeq_flag;
-
- comm = mpi_data->world;
- qeq_flag = 0;
+ MPI_Comm comm = mpi_data->world;
- if( qeq_flag )
- Init_Forces( system, control, data, workspace, lists, out_control, comm );
- else
- Init_Forces_noQEq( system, control, data, workspace,
+ Init_Forces_noQEq( system, control, data, workspace,
lists, out_control, comm );
-
/********* bonded interactions ************/
Compute_Bonded_Forces( system, control, data, workspace,
lists, out_control, mpi_data->world );
diff --git a/src/USER-REAXC/reaxc_io_tools.cpp b/src/USER-REAXC/reaxc_io_tools.cpp
index 1a978d0073..0c14dad5d4 100644
--- a/src/USER-REAXC/reaxc_io_tools.cpp
+++ b/src/USER-REAXC/reaxc_io_tools.cpp
@@ -34,8 +34,6 @@
#include "reaxc_traj.h"
#include "reaxc_vector.h"
-print_interaction Print_Interactions[NUM_INTRS];
-
int Init_Output_Files( reax_system *system, control_params *control,
output_controls *out_control, mpi_datatypes *mpi_data,
char *msg )
@@ -109,432 +107,6 @@ int Close_Output_Files( reax_system *system, control_params *control,
}
-
-void Print_Box( simulation_box* box, char *name, FILE *out )
-{
- // int i, j;
-
- fprintf( out, "%s:\n", name );
- fprintf( out, "\tmin[%8.3f %8.3f %8.3f]\n",
- box->min[0], box->min[1], box->min[2] );
- fprintf( out, "\tmax[%8.3f %8.3f %8.3f]\n",
- box->max[0], box->max[1], box->max[2] );
- fprintf( out, "\tdims[%8.3f%8.3f%8.3f]\n",
- box->box_norms[0], box->box_norms[1], box->box_norms[2] );
-
-}
-
-
-
-void Print_Grid( grid* g, FILE *out )
-{
- int x, y, z, gc_type;
- ivec gc_str;
- char gcell_type_text[10][12] =
- { "NO_NBRS", "NEAR_ONLY", "HBOND_ONLY", "FAR_ONLY",
- "NEAR_HBOND", "NEAR_FAR", "HBOND_FAR", "FULL_NBRS", "NATIVE" };
-
- fprintf( out, "\tnumber of grid cells: %d %d %d\n",
- g->ncells[0], g->ncells[1], g->ncells[2] );
- fprintf( out, "\tgcell lengths: %8.3f %8.3f %8.3f\n",
- g->cell_len[0], g->cell_len[1], g->cell_len[2] );
- fprintf( out, "\tinverses of gcell lengths: %8.3f %8.3f %8.3f\n",
- g->inv_len[0], g->inv_len[1], g->inv_len[2] );
- fprintf( out, "\t---------------------------------\n" );
- fprintf( out, "\tnumber of native gcells: %d %d %d\n",
- g->native_cells[0], g->native_cells[1], g->native_cells[2] );
- fprintf( out, "\tnative gcell span: %d-%d %d-%d %d-%d\n",
- g->native_str[0], g->native_end[0],
- g->native_str[1], g->native_end[1],
- g->native_str[2], g->native_end[2] );
- fprintf( out, "\t---------------------------------\n" );
- fprintf( out, "\tvlist gcell stretch: %d %d %d\n",
- g->vlist_span[0], g->vlist_span[1], g->vlist_span[2] );
- fprintf( out, "\tnonbonded nbrs gcell stretch: %d %d %d\n",
- g->nonb_span[0], g->nonb_span[1], g->nonb_span[2] );
- fprintf( out, "\tbonded nbrs gcell stretch: %d %d %d\n",
- g->bond_span[0], g->bond_span[1], g->bond_span[2] );
- fprintf( out, "\t---------------------------------\n" );
- fprintf( out, "\tghost gcell span: %d %d %d\n",
- g->ghost_span[0], g->ghost_span[1], g->ghost_span[2] );
- fprintf( out, "\tnonbonded ghost gcell span: %d %d %d\n",
- g->ghost_nonb_span[0],g->ghost_nonb_span[1],g->ghost_nonb_span[2]);
- fprintf(out, "\thbonded ghost gcell span: %d %d %d\n",
- g->ghost_hbond_span[0],g->ghost_hbond_span[1],g->ghost_hbond_span[2]);
- fprintf( out, "\tbonded ghost gcell span: %d %d %d\n",
- g->ghost_bond_span[0],g->ghost_bond_span[1],g->ghost_bond_span[2]);
- fprintf( out, "\t---------------------------------\n" );
-
- fprintf( stderr, "GCELL MARKS:\n" );
- gc_type = g->cells[0][0][0].type;
- ivec_MakeZero( gc_str );
-
- x = y = z = 0;
- for( x = 0; x < g->ncells[0]; ++x )
- for( y = 0; y < g->ncells[1]; ++y )
- for( z = 0; z < g->ncells[2]; ++z )
- if( g->cells[x][y][z].type != gc_type ){
- fprintf( stderr,
- "\tgcells from(%2d %2d %2d) to (%2d %2d %2d): %d - %s\n",
- gc_str[0], gc_str[1], gc_str[2], x, y, z,
- gc_type, gcell_type_text[gc_type] );
- gc_type = g->cells[x][y][z].type;
- gc_str[0] = x;
- gc_str[1] = y;
- gc_str[2] = z;
- }
- fprintf( stderr, "\tgcells from(%2d %2d %2d) to (%2d %2d %2d): %d - %s\n",
- gc_str[0], gc_str[1], gc_str[2], x, y, z,
- gc_type, gcell_type_text[gc_type] );
- fprintf( out, "-------------------------------------\n" );
-}
-
-void Print_Native_GCells( reax_system *system )
-{
- int i, j, k, l;
- char fname[100];
- FILE *f;
- grid *g;
- grid_cell *gc;
- char gcell_type_text[10][12] =
- { "NO_NBRS", "NEAR_ONLY", "HBOND_ONLY", "FAR_ONLY",
- "NEAR_HBOND", "NEAR_FAR", "HBOND_FAR", "FULL_NBRS", "NATIVE" };
-
- sprintf( fname, "native_gcells.%d", system->my_rank );
- f = fopen( fname, "w" );
- g = &(system->my_grid);
-
- for( i = g->native_str[0]; i < g->native_end[0]; i++ )
- for( j = g->native_str[1]; j < g->native_end[1]; j++ )
- for( k = g->native_str[2]; k < g->native_end[2]; k++ )
- {
- gc = &( g->cells[i][j][k] );
-
- fprintf( f, "p%d gcell(%2d %2d %2d) of type %d(%s)\n",
- system->my_rank, i, j, k,
- gc->type, gcell_type_text[gc->type] );
-
- fprintf( f, "\tatom list start: %d, end: %d\n\t", gc->str, gc->end );
-
- for( l = gc->str; l < gc->end; ++l )
- fprintf( f, TAGINT_FORMAT, system->my_atoms[l].orig_id );
- fprintf( f, "\n" );
- }
-
- fclose(f);
-}
-
-
-
-void Print_All_GCells( reax_system *system )
-{
- int i, j, k, l;
- char fname[100];
- FILE *f;
- grid *g;
- grid_cell *gc;
- char gcell_type_text[10][12] =
- { "NO_NBRS", "NEAR_ONLY", "HBOND_ONLY", "FAR_ONLY",
- "NEAR_HBOND", "NEAR_FAR", "HBOND_FAR", "FULL_NBRS", "NATIVE" };
-
- sprintf( fname, "all_gcells.%d", system->my_rank );
- f = fopen( fname, "w" );
- g = &(system->my_grid);
-
- for( i = 0; i < g->ncells[0]; i++ )
- for( j = 0; j < g->ncells[1]; j++ )
- for( k = 0; k < g->ncells[2]; k++ )
- {
- gc = &( g->cells[i][j][k] );
-
- fprintf( f, "p%d gcell(%2d %2d %2d) of type %d(%s)\n",
- system->my_rank, i, j, k,
- gc->type, gcell_type_text[gc->type] );
-
- fprintf( f, "\tatom list start: %d, end: %d\n\t", gc->str, gc->end );
-
- for( l = gc->str; l < gc->end; ++l )
- fprintf( f, TAGINT_FORMAT, system->my_atoms[l].orig_id );
- fprintf( f, "\n" );
- }
-
- fclose(f);
-}
-
-
-
-void Print_My_Atoms( reax_system *system )
-{
- int i;
- char fname[100];
- FILE *fh;
-
- sprintf( fname, "my_atoms.%d", system->my_rank );
- if( (fh = fopen( fname, "w" )) == NULL )
- {
- fprintf( stderr, "error in opening my_atoms file" );
- MPI_Abort( MPI_COMM_WORLD, FILE_NOT_FOUND );
- }
-
- for( i = 0; i < system->n; ++i )
- fprintf( fh, "p%-2d %-5d %2d %24.15e%24.15e%24.15e\n",
- system->my_rank,
- system->my_atoms[i].orig_id, system->my_atoms[i].type,
- system->my_atoms[i].x[0],
- system->my_atoms[i].x[1],
- system->my_atoms[i].x[2] );
-
- fclose( fh );
-}
-
-
-void Print_My_Ext_Atoms( reax_system *system )
-{
- int i;
- char fname[100];
- FILE *fh;
-
- sprintf( fname, "my_ext_atoms.%d", system->my_rank );
- if( (fh = fopen( fname, "w" )) == NULL )
- {
- fprintf( stderr, "error in opening my_ext_atoms file" );
- MPI_Abort( MPI_COMM_WORLD, FILE_NOT_FOUND );
- }
-
- for( i = 0; i < system->N; ++i )
- fprintf( fh, "p%-2d %-5d imprt%-5d %2d %24.15e%24.15e%24.15e\n",
- system->my_rank, system->my_atoms[i].orig_id,
- system->my_atoms[i].imprt_id, system->my_atoms[i].type,
- system->my_atoms[i].x[0],
- system->my_atoms[i].x[1],
- system->my_atoms[i].x[2] );
-
- fclose( fh );
-}
-
-
-void Print_Far_Neighbors( reax_system *system, reax_list **lists,
- control_params *control )
-{
- char fname[100];
- int i, j, nbr, natoms;
- rc_tagint id_i, id_j;
- FILE *fout;
- reax_list *far_nbrs;
-
- sprintf( fname, "%s.far_nbrs.%d", control->sim_name, system->my_rank );
- fout = fopen( fname, "w" );
- far_nbrs = (*lists) + FAR_NBRS;
- natoms = system->N;
-
- for( i = 0; i < natoms; ++i ) {
- id_i = system->my_atoms[i].orig_id;
-
- for( j = Start_Index(i,far_nbrs); j < End_Index(i,far_nbrs); ++j ) {
- nbr = far_nbrs->select.far_nbr_list[j].nbr;
- id_j = system->my_atoms[nbr].orig_id;
-
- fprintf( fout, "%6d%6d%24.15e%24.15e%24.15e%24.15e\n",
- id_i, id_j, far_nbrs->select.far_nbr_list[j].d,
- far_nbrs->select.far_nbr_list[j].dvec[0],
- far_nbrs->select.far_nbr_list[j].dvec[1],
- far_nbrs->select.far_nbr_list[j].dvec[2] );
-
- fprintf( fout, "%6d%6d%24.15e%24.15e%24.15e%24.15e\n",
- id_j, id_i, far_nbrs->select.far_nbr_list[j].d,
- -far_nbrs->select.far_nbr_list[j].dvec[0],
- -far_nbrs->select.far_nbr_list[j].dvec[1],
- -far_nbrs->select.far_nbr_list[j].dvec[2] );
- }
- }
-
- fclose( fout );
-}
-
-
-void Print_Sparse_Matrix( reax_system *system, sparse_matrix *A )
-{
- int i, j;
-
- for( i = 0; i < A->n; ++i )
- for( j = A->start[i]; j < A->end[i]; ++j )
- fprintf( stderr, "%d %d %.15e\n",
- system->my_atoms[i].orig_id,
- system->my_atoms[A->entries[j].j].orig_id,
- A->entries[j].val );
-}
-
-
-void Print_Sparse_Matrix2( reax_system *system, sparse_matrix *A, char *fname )
-{
- int i, j;
- FILE *f = fopen( fname, "w" );
-
- for( i = 0; i < A->n; ++i )
- for( j = A->start[i]; j < A->end[i]; ++j )
- fprintf( f, "%d %d %.15e\n",
- system->my_atoms[i].orig_id,
- system->my_atoms[A->entries[j].j].orig_id,
- A->entries[j].val );
-
- fclose(f);
-}
-
-
-void Print_Symmetric_Sparse(reax_system *system, sparse_matrix *A, char *fname)
-{
- int i, j;
- reax_atom *ai, *aj;
- FILE *f = fopen( fname, "w" );
-
- for( i = 0; i < A->n; ++i ) {
- ai = &(system->my_atoms[i]);
- for( j = A->start[i]; j < A->end[i]; ++j ) {
- aj = &(system->my_atoms[A->entries[j].j]);
- fprintf( f, "%d %d %.15e\n",
- ai->renumber, aj->renumber, A->entries[j].val );
- if( A->entries[j].j < system->n && ai->renumber != aj->renumber )
- fprintf( f, "%d %d %.15e\n",
- aj->renumber, ai->renumber, A->entries[j].val );
- }
- }
-
- fclose(f);
-}
-
-
-void Print_Linear_System( reax_system *system, control_params *control,
- storage *workspace, int step )
-{
- int i;
- char fname[100];
- reax_atom *ai;
- FILE *out;
-
- // print rhs and init guesses for QEq
- sprintf( fname, "%s.p%dstate%d", control->sim_name, system->my_rank, step );
- out = fopen( fname, "w" );
- for( i = 0; i < system->n; i++ ) {
- ai = &(system->my_atoms[i]);
- fprintf( out, "%6d%2d%24.15e%24.15e%24.15e%24.15e%24.15e%24.15e%24.15e\n",
- ai->renumber, ai->type, ai->x[0], ai->x[1], ai->x[2],
- workspace->s[i], workspace->b_s[i],
- workspace->t[i], workspace->b_t[i] );
- }
- fclose( out );
-
- // print QEq coef matrix
- sprintf( fname, "%s.p%dH%d", control->sim_name, system->my_rank, step );
- Print_Symmetric_Sparse( system, workspace->H, fname );
-
-}
-
-
-void Print_LinSys_Soln( reax_system *system, double *x, double *b_prm, double *b )
-{
- int i;
- char fname[100];
- FILE *fout;
-
- sprintf( fname, "qeq.%d.out", system->my_rank );
- fout = fopen( fname, "w" );
-
- for( i = 0; i < system->n; ++i )
- fprintf( fout, "%6d%10.4f%10.4f%10.4f\n",
- system->my_atoms[i].orig_id, x[i], b_prm[i], b[i] );
-
- fclose( fout );
-}
-
-
-void Print_Charges( reax_system *system )
-{
- int i;
- char fname[100];
- FILE *fout;
-
- sprintf( fname, "q.%d.out", system->my_rank );
- fout = fopen( fname, "w" );
-
- for( i = 0; i < system->n; ++i )
- fprintf( fout, "%6d %10.7f %10.7f %10.7f\n",
- system->my_atoms[i].orig_id,
- system->my_atoms[i].s[0],
- system->my_atoms[i].t[0],
- system->my_atoms[i].q );
-
- fclose( fout );
-}
-
-
-void Print_Bonds( reax_system *system, reax_list *bonds, char *fname )
-{
- int i, j, pj;
- bond_data *pbond;
- bond_order_data *bo_ij;
- FILE *f = fopen( fname, "w" );
-
- for( i = 0; i < system->N; ++i )
- for( pj = Start_Index(i, bonds); pj < End_Index(i, bonds); ++pj ) {
- pbond = &(bonds->select.bond_list[pj]);
- bo_ij = &(pbond->bo_data);
- j = pbond->nbr;
- fprintf( f, "%8d%8d %24.15f %24.15f\n",
- i, j,//system->my_atoms[i].orig_id, system->my_atoms[j].orig_id,
- pbond->d, bo_ij->BO );
- }
-
- fclose(f);
-}
-
-
-int fn_qsort_intcmp( const void *a, const void *b )
-{
- return( *(int *)a - *(int *)b );
-}
-
-void Print_Bond_List2( reax_system *system, reax_list *bonds, char *fname )
-{
- int i,j, nbr, pj;
- rc_tagint id_i, id_j;
- FILE *f = fopen( fname, "w" );
- int temp[500];
- int num=0;
-
- for( i = 0; i < system->n; ++i ) {
- num=0;
- id_i = system->my_atoms[i].orig_id;
- fprintf( f, "%6d:", id_i);
- for( pj = Start_Index(i, bonds); pj < End_Index(i, bonds); ++pj ) {
- nbr = bonds->select.bond_list[pj].nbr;
- id_j = system->my_atoms[nbr].orig_id;
- if( id_i < id_j )
- temp[num++] = id_j;
- }
-
- qsort(&temp, num, sizeof(int), fn_qsort_intcmp);
- for(j=0; j < num; j++)
- fprintf(f, "%6d", temp[j] );
- fprintf(f, "\n");
- }
-}
-
-
-void Print_Total_Force( reax_system *system, simulation_data *data,
- storage *workspace )
-{
- int i;
-
- fprintf( stderr, "step: %d\n", data->step );
- fprintf( stderr, "%6s\t%-38s\n", "atom", "atom.f[0,1,2]");
-
- for( i = 0; i < system->N; ++i )
- fprintf( stderr, "%6d %f %f %f\n",
- //"%6d%24.15e%24.15e%24.15e\n",
- system->my_atoms[i].orig_id,
- workspace->f[i][0], workspace->f[i][1], workspace->f[i][2] );
-}
-
void Output_Results( reax_system *system, control_params *control,
simulation_data *data, reax_list **lists,
output_controls *out_control, mpi_datatypes *mpi_data )
diff --git a/src/USER-REAXC/reaxc_io_tools.h b/src/USER-REAXC/reaxc_io_tools.h
index 12312d344e..a3f22fccc2 100644
--- a/src/USER-REAXC/reaxc_io_tools.h
+++ b/src/USER-REAXC/reaxc_io_tools.h
@@ -33,80 +33,6 @@ int Init_Output_Files( reax_system*, control_params*,
output_controls*, mpi_datatypes*, char* );
int Close_Output_Files( reax_system*, control_params*,
output_controls*, mpi_datatypes* );
-
-void Print_Box( simulation_box*, char*, FILE* );
-
-void Print_Grid( grid*, FILE* );
-void Print_GCell_Exchange_Bounds( int, neighbor_proc* );
-void Print_Native_GCells( reax_system* );
-void Print_All_GCells( reax_system*);
-
-void Print_Init_Atoms( reax_system*, storage* );
-void Print_My_Atoms( reax_system* );
-void Print_My_Ext_Atoms( reax_system* );
-
-void Print_Far_Neighbors( reax_system*, reax_list**, control_params *);
-void Print_Sparse_Matrix( reax_system*, sparse_matrix* );
-void Print_Sparse_Matrix2( reax_system*, sparse_matrix*, char* );
-void Print_Linear_System( reax_system*, control_params*, storage*, int );
-void Print_LinSys_Soln( reax_system*, double*, double*, double* );
-void Print_Charges( reax_system* );
-void Print_Bonds( reax_system*, reax_list*, char* );
-void Print_Bond_List2( reax_system*, reax_list*, char* );
-void Print_Total_Force( reax_system*, simulation_data*, storage* );
void Output_Results( reax_system*, control_params*, simulation_data*,
reax_list**, output_controls*, mpi_datatypes* );
-
-#if defined(DEBUG_FOCUS) || defined(TEST_FORCES) || defined(TEST_ENERGY)
-void Debug_Marker_Bonded( output_controls*, int );
-void Debug_Marker_Nonbonded( output_controls*, int );
-void Print_Near_Neighbors_List( reax_system*, reax_list**, control_params*,
- simulation_data*, output_controls* );
-void Print_Far_Neighbors_List( reax_system*, reax_list**, control_params*,
- simulation_data*, output_controls* );
-void Print_Bond_List( reax_system*, control_params*, simulation_data*,
- reax_list**, output_controls* );
-/*void Dummy_Printer( reax_system*, control_params*, simulation_data*,
- storage*, reax_list**, output_controls* );
-void Print_Bond_Orders( reax_system*, control_params*, simulation_data*,
- storage*, reax_list**, output_controls* );
-void Print_Bond_Forces( reax_system*, control_params*, simulation_data*,
- storage*, reax_list**, output_controls* );
-void Print_LonePair_Forces( reax_system*, control_params*, simulation_data*,
- storage*, reax_list**, output_controls* );
-void Print_OverUnderCoor_Forces( reax_system*, control_params*,
- simulation_data*, storage*, reax_list**,
- output_controls* );
-void Print_Three_Body_Forces( reax_system*, control_params*, simulation_data*,
- storage*, reax_list**, output_controls* );
-void Print_Hydrogen_Bond_Forces( reax_system*, control_params*,
- simulation_data*, storage*, reax_list**,
- output_controls* );
-void Print_Four_Body_Forces( reax_system*, control_params*, simulation_data*,
- storage*, reax_list**, output_controls* );
-void Print_vdW_Coulomb_Forces( reax_system*, control_params*,
- simulation_data*, storage*, reax_list**,
- output_controls* );
-void Print_Total_Force( reax_system*, control_params*, simulation_data*,
- storage*, reax_list**, output_controls* );
-void Compare_Total_Forces( reax_system*, control_params*, simulation_data*,
-storage*, reax_list**, output_controls* );*/
-//void Print_Total_Force( reax_system*, control_params* );
-void Print_Force_Files( reax_system*, control_params*, simulation_data*,
- storage*, reax_list**, output_controls*,
- mpi_datatypes * );
-//void Init_Force_Test_Functions( );
-
-int fn_qsort_intcmp( const void *, const void * );
-
-void Print_Far_Neighbors_List( reax_system*, reax_list**, control_params*,
- simulation_data*, output_controls* );
-
-void Print_Near_Neighbors_List( reax_system*, reax_list**, control_params*,
- simulation_data*, output_controls* );
-
-void Print_Bond_List( reax_system*, control_params*, simulation_data*,
- reax_list**, output_controls*);
-
-#endif
#endif
diff --git a/src/USER-REAXC/reaxc_lookup.cpp b/src/USER-REAXC/reaxc_lookup.cpp
index 7497a5c096..903e54962d 100644
--- a/src/USER-REAXC/reaxc_lookup.cpp
+++ b/src/USER-REAXC/reaxc_lookup.cpp
@@ -150,30 +150,6 @@ void Complete_Cubic_Spline( const double *h, const double *f, double v0, double
}
-void LR_Lookup( LR_lookup_table *t, double r, LR_data *y )
-{
- int i;
- double base, dif;
-
- i = (int)(r * t->inv_dx);
- if( i == 0 ) ++i;
- base = (double)(i+1) * t->dx;
- dif = r - base;
-
- y->e_vdW = ((t->vdW[i].d*dif + t->vdW[i].c)*dif + t->vdW[i].b)*dif +
- t->vdW[i].a;
- y->CEvd = ((t->CEvd[i].d*dif + t->CEvd[i].c)*dif +
- t->CEvd[i].b)*dif + t->CEvd[i].a;
-
- y->e_ele = ((t->ele[i].d*dif + t->ele[i].c)*dif + t->ele[i].b)*dif +
- t->ele[i].a;
- y->CEclmb = ((t->CEclmb[i].d*dif + t->CEclmb[i].c)*dif + t->CEclmb[i].b)*dif +
- t->CEclmb[i].a;
-
- y->H = y->e_ele * EV_to_KCALpMOL / C_ele;
-}
-
-
int Init_Lookup_Tables( reax_system *system, control_params *control,
storage *workspace, mpi_datatypes *mpi_data, char *msg )
{
diff --git a/src/USER-REAXC/reaxc_reset_tools.cpp b/src/USER-REAXC/reaxc_reset_tools.cpp
index 7d9e9deb76..1e6aeab475 100644
--- a/src/USER-REAXC/reaxc_reset_tools.cpp
+++ b/src/USER-REAXC/reaxc_reset_tools.cpp
@@ -119,29 +119,6 @@ void Reset_Workspace( reax_system *system, storage *workspace )
}
-void Reset_Grid( grid *g )
-{
- int i, j, k;
-
- for( i = 0; i < g->ncells[0]; i++ )
- for( j = 0; j < g->ncells[1]; j++ )
- for( k = 0; k < g->ncells[2]; k++ ) {
- g->cells[i][j][k].top = 0;
- g->cells[i][j][k].str = 0;
- g->cells[i][j][k].end = 0;
- }
-}
-
-
-void Reset_Out_Buffers( mpi_out_data *out_buf, int n )
-{
- int i;
-
- for( i = 0; i < n; ++i )
- out_buf[i].cnt = 0;
-}
-
-
void Reset_Neighbor_Lists( reax_system *system, control_params *control,
storage *workspace, reax_list **lists,
MPI_Comm comm )
diff --git a/src/USER-REAXC/reaxc_reset_tools.h b/src/USER-REAXC/reaxc_reset_tools.h
index adc562735c..c2a90072d5 100644
--- a/src/USER-REAXC/reaxc_reset_tools.h
+++ b/src/USER-REAXC/reaxc_reset_tools.h
@@ -33,13 +33,8 @@ void Reset_Pressures( simulation_data* );
void Reset_Simulation_Data( simulation_data*, int );
void Reset_Timing( reax_timing* );
void Reset_Workspace( reax_system*, storage* );
-void Reset_Grid( grid* );
-void Reset_Out_Buffers( mpi_out_data*, int );
void Reset_Neighbor_Lists( reax_system*, control_params*, storage*,
reax_list**, MPI_Comm );
void Reset( reax_system*, control_params*, simulation_data*, storage*,
reax_list**, MPI_Comm );
-#ifdef TEST_FORCES
-void Reset_Test_Forces( reax_system*, storage* );
-#endif
#endif
diff --git a/src/USER-REAXC/reaxc_system_props.cpp b/src/USER-REAXC/reaxc_system_props.cpp
index 554151ba82..6b4551a03f 100644
--- a/src/USER-REAXC/reaxc_system_props.cpp
+++ b/src/USER-REAXC/reaxc_system_props.cpp
@@ -29,55 +29,6 @@
#include "reaxc_tool_box.h"
#include "reaxc_vector.h"
-void Temperature_Control( control_params *control, simulation_data *data )
-{
- double tmp;
-
- if( control->T_mode == 1 ) {// step-wise temperature control
- if((data->step-data->prev_steps) % ((int)(control->T_freq/control->dt))==0){
- if( fabs( control->T - control->T_final ) >= fabs( control->T_rate ) )
- control->T += control->T_rate;
- else control->T = control->T_final;
- }
- }
- else if( control->T_mode == 2 ) { // constant slope control
- tmp = control->T_rate * control->dt / control->T_freq;
-
- if( fabs( control->T - control->T_final ) >= fabs( tmp ) )
- control->T += tmp;
- }
-}
-
-
-
-void Compute_Kinetic_Energy( reax_system* system, simulation_data* data,
- MPI_Comm comm )
-{
- int i;
- rvec p;
- double m;
-
- data->my_en.e_kin = 0.0;
- data->sys_en.e_kin = 0.0;
- data->therm.T = 0;
-
- for( i = 0; i < system->n; i++ ) {
- m = system->reax_param.sbp[system->my_atoms[i].type].mass;
-
- rvec_Scale( p, m, system->my_atoms[i].v );
- data->my_en.e_kin += 0.5 * rvec_Dot( p, system->my_atoms[i].v );
- }
-
- MPI_Allreduce( &data->my_en.e_kin, &data->sys_en.e_kin,
- 1, MPI_DOUBLE, MPI_SUM, comm );
-
- data->therm.T = (2. * data->sys_en.e_kin) / (data->N_f * K_B);
-
- // avoid T being an absolute zero, might cause F.P.E!
- if( fabs(data->therm.T) < ALMOST_ZERO )
- data->therm.T = ALMOST_ZERO;
-}
-
void Compute_System_Energy( reax_system *system, simulation_data *data,
MPI_Comm comm )
@@ -135,169 +86,3 @@ void Compute_System_Energy( reax_system *system, simulation_data *data,
data->sys_en.e_tot = data->sys_en.e_pot + E_CONV * data->sys_en.e_kin;
}
}
-
-
-void Compute_Total_Mass( reax_system *system, simulation_data *data,
- MPI_Comm comm )
-{
- int i;
- double tmp;
-
- tmp = 0;
- for( i = 0; i < system->n; i++ )
- tmp += system->reax_param.sbp[ system->my_atoms[i].type ].mass;
-
- MPI_Allreduce( &tmp, &data->M, 1, MPI_DOUBLE, MPI_SUM, comm );
-
- data->inv_M = 1. / data->M;
-}
-
-
-
-void Compute_Center_of_Mass( reax_system *system, simulation_data *data,
- mpi_datatypes *mpi_data, MPI_Comm comm )
-{
- int i;
- double m, det; //xx, xy, xz, yy, yz, zz;
- double tmp_mat[6], tot_mat[6];
- rvec my_xcm, my_vcm, my_amcm, my_avcm;
- rvec tvec, diff;
- rtensor mat, inv;
-
- rvec_MakeZero( my_xcm ); // position of CoM
- rvec_MakeZero( my_vcm ); // velocity of CoM
- rvec_MakeZero( my_amcm ); // angular momentum of CoM
- rvec_MakeZero( my_avcm ); // angular velocity of CoM
-
- /* Compute the position, vel. and ang. momentum about the centre of mass */
- for( i = 0; i < system->n; ++i ) {
- m = system->reax_param.sbp[ system->my_atoms[i].type ].mass;
-
- rvec_ScaledAdd( my_xcm, m, system->my_atoms[i].x );
- rvec_ScaledAdd( my_vcm, m, system->my_atoms[i].v );
-
- rvec_Cross( tvec, system->my_atoms[i].x, system->my_atoms[i].v );
- rvec_ScaledAdd( my_amcm, m, tvec );
- }
-
- MPI_Allreduce( my_xcm, data->xcm, 3, MPI_DOUBLE, MPI_SUM, comm );
- MPI_Allreduce( my_vcm, data->vcm, 3, MPI_DOUBLE, MPI_SUM, comm );
- MPI_Allreduce( my_amcm, data->amcm, 3, MPI_DOUBLE, MPI_SUM, comm );
-
- rvec_Scale( data->xcm, data->inv_M, data->xcm );
- rvec_Scale( data->vcm, data->inv_M, data->vcm );
- rvec_Cross( tvec, data->xcm, data->vcm );
- rvec_ScaledAdd( data->amcm, -data->M, tvec );
- data->etran_cm = 0.5 * data->M * rvec_Norm_Sqr( data->vcm );
-
- /* Calculate and then invert the inertial tensor */
- for( i = 0; i < 6; ++i )
- tmp_mat[i] = 0;
- //my_xx = my_xy = my_xz = my_yy = my_yz = my_zz = 0;
-
- for( i = 0; i < system->n; ++i ){
- m = system->reax_param.sbp[ system->my_atoms[i].type ].mass;
- rvec_ScaledSum( diff, 1., system->my_atoms[i].x, -1., data->xcm );
-
- tmp_mat[0]/*my_xx*/ += diff[0] * diff[0] * m;
- tmp_mat[1]/*my_xy*/ += diff[0] * diff[1] * m;
- tmp_mat[2]/*my_xz*/ += diff[0] * diff[2] * m;
- tmp_mat[3]/*my_yy*/ += diff[1] * diff[1] * m;
- tmp_mat[4]/*my_yz*/ += diff[1] * diff[2] * m;
- tmp_mat[5]/*my_zz*/ += diff[2] * diff[2] * m;
- }
-
- MPI_Reduce( tmp_mat, tot_mat, 6, MPI_DOUBLE, MPI_SUM, MASTER_NODE, comm );
-
- if( system->my_rank == MASTER_NODE ) {
- mat[0][0] = tot_mat[3] + tot_mat[5]; // yy + zz;
- mat[0][1] = mat[1][0] = -tot_mat[1]; // -xy;
- mat[0][2] = mat[2][0] = -tot_mat[2]; // -xz;
- mat[1][1] = tot_mat[0] + tot_mat[5]; // xx + zz;
- mat[2][1] = mat[1][2] = -tot_mat[4]; // -yz;
- mat[2][2] = tot_mat[0] + tot_mat[3]; // xx + yy;
-
- /* invert the inertial tensor */
- det = ( mat[0][0] * mat[1][1] * mat[2][2] +
- mat[0][1] * mat[1][2] * mat[2][0] +
- mat[0][2] * mat[1][0] * mat[2][1] ) -
- ( mat[0][0] * mat[1][2] * mat[2][1] +
- mat[0][1] * mat[1][0] * mat[2][2] +
- mat[0][2] * mat[1][1] * mat[2][0] );
-
- inv[0][0] = mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1];
- inv[0][1] = mat[0][2] * mat[2][1] - mat[0][1] * mat[2][2];
- inv[0][2] = mat[0][1] * mat[1][2] - mat[0][2] * mat[1][1];
- inv[1][0] = mat[1][2] * mat[2][0] - mat[1][0] * mat[2][2];
- inv[1][1] = mat[0][0] * mat[2][2] - mat[0][2] * mat[2][0];
- inv[1][2] = mat[0][2] * mat[1][0] - mat[0][0] * mat[1][2];
- inv[2][0] = mat[1][0] * mat[2][1] - mat[2][0] * mat[1][1];
- inv[2][1] = mat[2][0] * mat[0][1] - mat[0][0] * mat[2][1];
- inv[2][2] = mat[0][0] * mat[1][1] - mat[1][0] * mat[0][1];
-
- if( det > ALMOST_ZERO )
- rtensor_Scale( inv, 1./det, inv );
- else rtensor_MakeZero( inv );
-
- /* Compute the angular velocity about the centre of mass */
- rtensor_MatVec( data->avcm, inv, data->amcm );
- }
-
- MPI_Bcast( data->avcm, 3, MPI_DOUBLE, MASTER_NODE, comm );
-
- /* Compute the rotational energy */
- data->erot_cm = 0.5 * E_CONV * rvec_Dot( data->avcm, data->amcm );
-
-}
-
-void Compute_Pressure(reax_system* system, control_params *control,
- simulation_data* data, mpi_datatypes *mpi_data)
-{
- int i;
- reax_atom *p_atom;
- rvec tmp, tx, int_press;
- simulation_box *big_box = &(system->big_box);
-
- /* Calculate internal pressure */
- rvec_MakeZero( int_press );
-
- // 0: both int and ext, 1: ext only, 2: int only
- if( control->press_mode == 0 || control->press_mode == 2 ) {
- for( i = 0; i < system->n; ++i ) {
- p_atom = &( system->my_atoms[i] );
-
- /* transform x into unitbox coordinates */
- Transform_to_UnitBox( p_atom->x, big_box, 1, tx );
-
- /* this atom's contribution to internal pressure */
- rvec_Multiply( tmp, p_atom->f, tx );
- rvec_Add( int_press, tmp );
-
- }
- }
-
- /* sum up internal and external pressure */
- MPI_Allreduce( int_press, data->int_press,
- 3, MPI_DOUBLE, MPI_SUM, mpi_data->comm_mesh3D );
- MPI_Allreduce( data->my_ext_press, data->ext_press,
- 3, MPI_DOUBLE, MPI_SUM, mpi_data->comm_mesh3D );
-
- /* kinetic contribution */
- data->kin_press = 2.*(E_CONV*data->sys_en.e_kin) / (3.*big_box->V*P_CONV);
-
- /* Calculate total pressure in each direction */
- data->tot_press[0] = data->kin_press -
- (( data->int_press[0] + data->ext_press[0] ) /
- ( big_box->box_norms[1] * big_box->box_norms[2] * P_CONV ));
-
- data->tot_press[1] = data->kin_press -
- (( data->int_press[1] + data->ext_press[1] ) /
- ( big_box->box_norms[0] * big_box->box_norms[2] * P_CONV ));
-
- data->tot_press[2] = data->kin_press -
- (( data->int_press[2] + data->ext_press[2] ) /
- ( big_box->box_norms[0] * big_box->box_norms[1] * P_CONV ));
-
- data->iso_bar.P =
- ( data->tot_press[0] + data->tot_press[1] + data->tot_press[2] ) / 3.;
-}
diff --git a/src/USER-REAXC/reaxc_system_props.h b/src/USER-REAXC/reaxc_system_props.h
index 544b051dce..161060e184 100644
--- a/src/USER-REAXC/reaxc_system_props.h
+++ b/src/USER-REAXC/reaxc_system_props.h
@@ -29,18 +29,6 @@
#include "reaxc_types.h"
-void Temperature_Control( control_params*, simulation_data* );
-
-void Compute_Kinetic_Energy( reax_system*, simulation_data*, MPI_Comm );
-
void Compute_System_Energy( reax_system*, simulation_data*, MPI_Comm );
-void Compute_Total_Mass( reax_system*, simulation_data*, MPI_Comm );
-
-void Compute_Center_of_Mass( reax_system*, simulation_data*,
- mpi_datatypes*, MPI_Comm );
-
-void Compute_Pressure( reax_system*, control_params*,
- simulation_data*, mpi_datatypes* );
-//void Compute_Pressure( reax_system*, simulation_data* );
#endif
diff --git a/src/USER-REAXC/reaxc_tool_box.cpp b/src/USER-REAXC/reaxc_tool_box.cpp
index 4096c75118..fa1f8b8a75 100644
--- a/src/USER-REAXC/reaxc_tool_box.cpp
+++ b/src/USER-REAXC/reaxc_tool_box.cpp
@@ -27,57 +27,6 @@
#include "pair_reax_c.h"
#include "reaxc_tool_box.h"
-void Transform( rvec x1, simulation_box *box, char flag, rvec x2 )
-{
- int i, j;
- double tmp;
-
- if (flag > 0) {
- for (i=0; i < 3; i++) {
- tmp = 0.0;
- for (j=0; j < 3; j++)
- tmp += box->trans[i][j]*x1[j];
- x2[i] = tmp;
- }
- }
- else {
- for (i=0; i < 3; i++) {
- tmp = 0.0;
- for (j=0; j < 3; j++)
- tmp += box->trans_inv[i][j]*x1[j];
- x2[i] = tmp;
- }
- }
-}
-
-
-void Transform_to_UnitBox( rvec x1, simulation_box *box, char flag, rvec x2 )
-{
- Transform( x1, box, flag, x2 );
-
- x2[0] /= box->box_norms[0];
- x2[1] /= box->box_norms[1];
- x2[2] /= box->box_norms[2];
-}
-
-void Fit_to_Periodic_Box( simulation_box *box, rvec *p )
-{
- int i;
-
- for( i = 0; i < 3; ++i ) {
- if( (*p)[i] < box->min[i] ) {
- /* handle lower coords */
- while( (*p)[i] < box->min[i] )
- (*p)[i] += box->box_norms[i];
- }
- else if( (*p)[i] >= box->max[i] ) {
- /* handle higher coords */
- while( (*p)[i] >= box->max[i] )
- (*p)[i] -= box->box_norms[i];
- }
- }
-}
-
struct timeval tim;
double t_end;
@@ -87,75 +36,6 @@ double Get_Time( )
return( tim.tv_sec + (tim.tv_usec / 1000000.0) );
}
-
-double Get_Timing_Info( double t_start )
-{
- gettimeofday(&tim, NULL );
- t_end = tim.tv_sec + (tim.tv_usec / 1000000.0);
- return (t_end - t_start);
-}
-
-
-void Update_Timing_Info( double *t_start, double *timing )
-{
- gettimeofday(&tim, NULL );
- t_end = tim.tv_sec + (tim.tv_usec / 1000000.0);
- *timing += (t_end - *t_start);
- *t_start = t_end;
-}
-
-int Get_Atom_Type( reax_interaction *reax_param, char *s, MPI_Comm comm )
-{
- int i;
-
- for( i = 0; i < reax_param->num_atom_types; ++i )
- if( !strcmp( reax_param->sbp[i].name, s ) )
- return i;
-
- fprintf( stderr, "Unknown atom type %s. Terminating...\n", s );
- MPI_Abort( comm, UNKNOWN_ATOM_TYPE );
-
- return -1;
-}
-
-
-
-char *Get_Element( reax_system *system, int i )
-{
- return &( system->reax_param.sbp[system->my_atoms[i].type].name[0] );
-}
-
-
-
-char *Get_Atom_Name( reax_system *system, int i )
-{
- return &(system->my_atoms[i].name[0]);
-}
-
-
-
-int Allocate_Tokenizer_Space( char **line, char **backup, char ***tokens )
-{
- int i;
-
- if( (*line = (char*) malloc( sizeof(char) * MAX_LINE )) == NULL )
- return FAILURE;
-
- if( (*backup = (char*) malloc( sizeof(char) * MAX_LINE )) == NULL )
- return FAILURE;
-
- if( (*tokens = (char**) malloc( sizeof(char*) * MAX_TOKENS )) == NULL )
- return FAILURE;
-
- for( i = 0; i < MAX_TOKENS; i++ )
- if( ((*tokens)[i] = (char*) malloc(sizeof(char) * MAX_TOKEN_LEN)) == NULL )
- return FAILURE;
-
- return SUCCESS;
-}
-
-
-
int Tokenize( char* s, char*** tok )
{
char test[MAX_LINE];
diff --git a/src/USER-REAXC/reaxc_tool_box.h b/src/USER-REAXC/reaxc_tool_box.h
index 2a5ea1d6db..90a711fdb7 100644
--- a/src/USER-REAXC/reaxc_tool_box.h
+++ b/src/USER-REAXC/reaxc_tool_box.h
@@ -30,34 +30,10 @@
#include "reaxc_types.h"
#include "reaxc_defs.h"
-/* from comm_tools.h */
-int SumScan( int, int, int, MPI_Comm );
-void SumScanB( int, int, int, int, MPI_Comm, int* );
-
-/* from box.h */
-void Transform_to_UnitBox( rvec, simulation_box*, char, rvec );
-void Fit_to_Periodic_Box( simulation_box*, rvec* );
-void Box_Touch_Point( simulation_box*, ivec, rvec );
-int is_Inside_Box( simulation_box*, rvec );
-int iown_midpoint( simulation_box*, rvec, rvec );
-
-/* from grid.h */
-void GridCell_Closest_Point( grid_cell*, grid_cell*, ivec, ivec, rvec );
-void GridCell_to_Box_Points( grid_cell*, ivec, rvec, rvec );
-double DistSqr_between_Special_Points( rvec, rvec );
-double DistSqr_to_Special_Point( rvec, rvec );
-int Relative_Coord_Encoding( ivec );
-
/* from system_props.h */
double Get_Time( );
-double Get_Timing_Info( double );
-void Update_Timing_Info( double*, double* );
/* from io_tools.h */
-int Get_Atom_Type( reax_interaction*, char*, MPI_Comm );
-char *Get_Element( reax_system*, int );
-char *Get_Atom_Name( reax_system*, int );
-int Allocate_Tokenizer_Space( char**, char**, char*** );
int Tokenize( char*, char*** );
/* from lammps */
diff --git a/src/USER-REAXC/reaxc_vector.cpp b/src/USER-REAXC/reaxc_vector.cpp
index 17a2851b00..ee63e94280 100644
--- a/src/USER-REAXC/reaxc_vector.cpp
+++ b/src/USER-REAXC/reaxc_vector.cpp
@@ -52,14 +52,6 @@ void rvec_ScaledAdd( rvec ret, double c, rvec v )
}
-void rvec_Sum( rvec ret, rvec v1 ,rvec v2 )
-{
- ret[0] = v1[0] + v2[0];
- ret[1] = v1[1] + v2[1];
- ret[2] = v1[2] + v2[2];
-}
-
-
void rvec_ScaledSum( rvec ret, double c1, rvec v1 ,double c2, rvec v2 )
{
ret[0] = c1 * v1[0] + c2 * v2[0];
@@ -74,20 +66,6 @@ double rvec_Dot( rvec v1, rvec v2 )
}
-double rvec_ScaledDot( double c1, rvec v1, double c2, rvec v2 )
-{
- return (c1*c2) * (v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2]);
-}
-
-
-void rvec_Multiply( rvec r, rvec v1, rvec v2 )
-{
- r[0] = v1[0] * v2[0];
- r[1] = v1[1] * v2[1];
- r[2] = v1[2] * v2[2];
-}
-
-
void rvec_iMultiply( rvec r, ivec v1, rvec v2 )
{
r[0] = v1[0] * v2[0];
@@ -96,30 +74,6 @@ void rvec_iMultiply( rvec r, ivec v1, rvec v2 )
}
-void rvec_Divide( rvec r, rvec v1, rvec v2 )
-{
- r[0] = v1[0] / v2[0];
- r[1] = v1[1] / v2[1];
- r[2] = v1[2] / v2[2];
-}
-
-
-void rvec_iDivide( rvec r, rvec v1, ivec v2 )
-{
- r[0] = v1[0] / v2[0];
- r[1] = v1[1] / v2[1];
- r[2] = v1[2] / v2[2];
-}
-
-
-void rvec_Invert( rvec r, rvec v )
-{
- r[0] = 1. / v[0];
- r[1] = 1. / v[1];
- r[2] = 1. / v[2];
-}
-
-
void rvec_Cross( rvec ret, rvec v1, rvec v2 )
{
ret[0] = v1[1] * v2[2] - v1[2] * v2[1];
@@ -128,16 +82,6 @@ void rvec_Cross( rvec ret, rvec v1, rvec v2 )
}
-void rvec_OuterProduct( rtensor r, rvec v1, rvec v2 )
-{
- int i, j;
-
- for( i = 0; i < 3; ++i )
- for( j = 0; j < 3; ++j )
- r[i][j] = v1[i] * v2[j];
-}
-
-
double rvec_Norm_Sqr( rvec v )
{
return SQR(v[0]) + SQR(v[1]) + SQR(v[2]);
@@ -150,16 +94,6 @@ double rvec_Norm( rvec v )
}
-int rvec_isZero( rvec v )
-{
- if( fabs(v[0]) > ALMOST_ZERO ||
- fabs(v[1]) > ALMOST_ZERO ||
- fabs(v[2]) > ALMOST_ZERO )
- return 0;
- return 1;
-}
-
-
void rvec_MakeZero( rvec v )
{
v[0] = v[1] = v[2] = 0.000000000000000e+00;
@@ -187,16 +121,6 @@ void rtensor_MatVec( rvec ret, rtensor m, rvec v )
}
-void rtensor_Scale( rtensor ret, double c, rtensor m )
-{
- int i, j;
-
- for( i = 0; i < 3; ++i )
- for( j = 0; j < 3; ++j )
- ret[i][j] = c * m[i][j];
-}
-
-
void rtensor_MakeZero( rtensor t )
{
t[0][0] = t[0][1] = t[0][2] = 0;
diff --git a/src/USER-REAXC/reaxc_vector.h b/src/USER-REAXC/reaxc_vector.h
index 4aafff87e9..906b200dc9 100644
--- a/src/USER-REAXC/reaxc_vector.h
+++ b/src/USER-REAXC/reaxc_vector.h
@@ -35,26 +35,17 @@ void rvec_Copy( rvec, rvec );
void rvec_Scale( rvec, double, rvec );
void rvec_Add( rvec, rvec );
void rvec_ScaledAdd( rvec, double, rvec );
-void rvec_Sum( rvec, rvec, rvec );
void rvec_ScaledSum( rvec, double, rvec, double, rvec );
double rvec_Dot( rvec, rvec );
-double rvec_ScaledDot( double, rvec, double, rvec );
-void rvec_Multiply( rvec, rvec, rvec );
void rvec_iMultiply( rvec, ivec, rvec );
-void rvec_Divide( rvec, rvec, rvec );
-void rvec_iDivide( rvec, rvec, ivec );
-void rvec_Invert( rvec, rvec );
void rvec_Cross( rvec, rvec, rvec );
-void rvec_OuterProduct( rtensor, rvec, rvec );
double rvec_Norm_Sqr( rvec );
double rvec_Norm( rvec );
-int rvec_isZero( rvec );
void rvec_MakeZero( rvec );
void rvec_Random( rvec );
void rtensor_MakeZero( rtensor );
void rtensor_MatVec( rvec, rtensor, rvec );
-void rtensor_Scale( rtensor, double, rtensor );
void ivec_MakeZero( ivec );
void ivec_Copy( ivec, ivec );
diff --git a/src/USER-SPH/atom_vec_meso.cpp b/src/USER-SPH/atom_vec_meso.cpp
index e4677703db..947139e791 100644
--- a/src/USER-SPH/atom_vec_meso.cpp
+++ b/src/USER-SPH/atom_vec_meso.cpp
@@ -482,7 +482,7 @@ int AtomVecMeso::pack_border_vel(int n, int *list, double *buf, int pbc_flag,
int *pbc)
{
int i,j,m;
- double dx,dy,dz;
+ double dx,dy,dz,dvx,dvy,dvz;
m = 0;
if (pbc_flag == 0) {
@@ -534,6 +534,9 @@ int AtomVecMeso::pack_border_vel(int n, int *list, double *buf, int pbc_flag,
buf[m++] = vest[j][2];
}
} else {
+ dvx = pbc[0] * h_rate[0] + pbc[5] * h_rate[5] + pbc[4] * h_rate[4];
+ dvy = pbc[1] * h_rate[1] + pbc[3] * h_rate[3];
+ dvz = pbc[2] * h_rate[2];
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x[j][0] + dx;
@@ -543,20 +546,23 @@ int AtomVecMeso::pack_border_vel(int n, int *list, double *buf, int pbc_flag,
buf[m++] = ubuf(type[j]).d;
buf[m++] = ubuf(mask[j]).d;
if (mask[i] & deform_groupbit) {
- buf[m++] = v[j][0];
- buf[m++] = v[j][1];
- buf[m++] = v[j][2];
+ buf[m++] = v[j][0] + dvx;
+ buf[m++] = v[j][1] + dvy;
+ buf[m++] = v[j][2] + dvz;
+ buf[m++] = vest[j][0] + dvx;
+ buf[m++] = vest[j][1] + dvy;
+ buf[m++] = vest[j][2] + dvz;
} else {
buf[m++] = v[j][0];
buf[m++] = v[j][1];
buf[m++] = v[j][2];
+ buf[m++] = vest[j][0];
+ buf[m++] = vest[j][1];
+ buf[m++] = vest[j][2];
}
buf[m++] = rho[j];
buf[m++] = e[j];
buf[m++] = cv[j];
- buf[m++] = vest[j][0];
- buf[m++] = vest[j][1];
- buf[m++] = vest[j][2];
}
}
}
@@ -619,12 +625,12 @@ void AtomVecMeso::unpack_border_vel(int n, int first, double *buf) {
v[i][0] = buf[m++];
v[i][1] = buf[m++];
v[i][2] = buf[m++];
- rho[i] = buf[m++];
- e[i] = buf[m++];
- cv[i] = buf[m++];
vest[i][0] = buf[m++];
vest[i][1] = buf[m++];
vest[i][2] = buf[m++];
+ rho[i] = buf[m++];
+ e[i] = buf[m++];
+ cv[i] = buf[m++];
}
if (atom->nextra_border)
diff --git a/src/fix.h b/src/fix.h
index f0836de021..24d209f019 100644
--- a/src/fix.h
+++ b/src/fix.h
@@ -120,6 +120,7 @@ class Fix : protected Pointers {
virtual void pre_exchange() {}
virtual void pre_neighbor() {}
virtual void pre_force(int) {}
+ virtual void pre_reverse(int,int) {}
virtual void post_force(int) {}
virtual void final_integrate() {}
virtual void end_of_step() {}
@@ -251,7 +252,8 @@ namespace FixConst {
static const int MIN_POST_FORCE = 1<<17;
static const int MIN_ENERGY = 1<<18;
static const int POST_RUN = 1<<19;
- static const int FIX_CONST_LAST = 1<<20;
+ static const int PRE_REVERSE = 1<<20;
+ static const int FIX_CONST_LAST = 1<<21;
}
}
diff --git a/src/info.cpp b/src/info.cpp
index c34eb401cd..1eba5a8b92 100644
--- a/src/info.cpp
+++ b/src/info.cpp
@@ -109,6 +109,7 @@ void Info::command(int narg, char **arg)
idx += 2;
} else if ((idx+1 < narg) && (strncmp(arg[idx],"out",3) == 0)
&& (strncmp(arg[idx+1],"log",3) == 0)) {
+ if ((out != screen) && (out != logfile)) fclose(out);
out = logfile;
idx += 2;
} else if ((idx+2 < narg) && (strncmp(arg[idx],"out",3) == 0)
diff --git a/src/input.cpp b/src/input.cpp
index 5b0c27e603..b048389c3c 100644
--- a/src/input.cpp
+++ b/src/input.cpp
@@ -15,6 +15,7 @@
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
+#include "errno.h"
#include "ctype.h"
#include "unistd.h"
#include "sys/stat.h"
@@ -1086,48 +1087,96 @@ void Input::quit()
/* ---------------------------------------------------------------------- */
+char *shell_failed_message(const char* cmd, int errnum)
+{
+ const char *errmsg = strerror(errnum);
+ int len = strlen(cmd)+strlen(errmsg)+64;
+ char *msg = new char[len];
+ sprintf(msg,"shell command '%s' failed with error '%s'", cmd, errmsg);
+ return msg;
+}
+
void Input::shell()
{
+ int rv,err;
+
if (narg < 1) error->all(FLERR,"Illegal shell command");
if (strcmp(arg[0],"cd") == 0) {
if (narg != 2) error->all(FLERR,"Illegal shell cd command");
- chdir(arg[1]);
+ rv = (chdir(arg[1]) < 0) ? errno : 0;
+ MPI_Reduce(&rv,&err,1,MPI_INT,MPI_MAX,0,world);
+ if (me == 0 && err != 0) {
+ char *message = shell_failed_message("cd",err);
+ error->warning(FLERR,message);
+ delete [] message;
+ }
} else if (strcmp(arg[0],"mkdir") == 0) {
if (narg < 2) error->all(FLERR,"Illegal shell mkdir command");
if (me == 0)
for (int i = 1; i < narg; i++) {
#if defined(_WIN32)
- _mkdir(arg[i]);
+ rv = _mkdir(arg[i]);
#else
- mkdir(arg[i], S_IRWXU | S_IRGRP | S_IXGRP);
+ rv = mkdir(arg[i], S_IRWXU | S_IRGRP | S_IXGRP);
#endif
+ if (rv < 0) {
+ char *message = shell_failed_message("mkdir",errno);
+ error->warning(FLERR,message);
+ delete [] message;
+ }
}
} else if (strcmp(arg[0],"mv") == 0) {
if (narg != 3) error->all(FLERR,"Illegal shell mv command");
- if (me == 0) rename(arg[1],arg[2]);
+ rv = (rename(arg[1],arg[2]) < 0) ? errno : 0;
+ MPI_Reduce(&rv,&err,1,MPI_INT,MPI_MAX,0,world);
+ if (me == 0 && err != 0) {
+ char *message = shell_failed_message("mv",err);
+ error->warning(FLERR,message);
+ delete [] message;
+ }
} else if (strcmp(arg[0],"rm") == 0) {
if (narg < 2) error->all(FLERR,"Illegal shell rm command");
if (me == 0)
- for (int i = 1; i < narg; i++) unlink(arg[i]);
+ for (int i = 1; i < narg; i++) {
+ if (unlink(arg[i]) < 0) {
+ char *message = shell_failed_message("rm",errno);
+ error->warning(FLERR,message);
+ delete [] message;
+ }
+ }
} else if (strcmp(arg[0],"rmdir") == 0) {
if (narg < 2) error->all(FLERR,"Illegal shell rmdir command");
if (me == 0)
- for (int i = 1; i < narg; i++) rmdir(arg[i]);
+ for (int i = 1; i < narg; i++) {
+ if (rmdir(arg[i]) < 0) {
+ char *message = shell_failed_message("rmdir",errno);
+ error->warning(FLERR,message);
+ delete [] message;
+ }
+ }
} else if (strcmp(arg[0],"putenv") == 0) {
if (narg < 2) error->all(FLERR,"Illegal shell putenv command");
for (int i = 1; i < narg; i++) {
char *ptr = strdup(arg[i]);
+ rv = 0;
#ifdef _WIN32
- if (ptr != NULL) _putenv(ptr);
+ if (ptr != NULL) rv = _putenv(ptr);
#else
- if (ptr != NULL) putenv(ptr);
+ if (ptr != NULL) rv = putenv(ptr);
#endif
+ rv = (rv < 0) ? errno : 0;
+ MPI_Reduce(&rv,&err,1,MPI_INT,MPI_MAX,0,world);
+ if (me == 0 && err != 0) {
+ char *message = shell_failed_message("putenv",err);
+ error->warning(FLERR,message);
+ delete [] message;
+ }
}
// use work string to concat args back into one string separated by spaces
@@ -1144,7 +1193,9 @@ void Input::shell()
strcat(work,arg[i]);
}
- if (me == 0) system(work);
+ if (me == 0)
+ if (system(work) != 0)
+ error->warning(FLERR,"shell command returned with non-zero status");
}
}
diff --git a/src/lammps.cpp b/src/lammps.cpp
index 4014090aac..95d925b6cb 100644
--- a/src/lammps.cpp
+++ b/src/lammps.cpp
@@ -195,6 +195,7 @@ LAMMPS::LAMMPS(int narg, char **arg, MPI_Comm communicator)
error->universe_all(FLERR,"Invalid command-line argument");
delete [] suffix;
delete [] suffix2;
+ suffix2 = NULL;
suffix_enable = 1;
// hybrid option to set fall-back for suffix2
if (strcmp(arg[iarg+1],"hybrid") == 0) {
diff --git a/src/min_cg.cpp b/src/min_cg.cpp
index 4953370562..af3ef5b7a4 100644
--- a/src/min_cg.cpp
+++ b/src/min_cg.cpp
@@ -66,6 +66,7 @@ int MinCG::iterate(int maxiter)
gg = fnorm_sqr();
for (int iter = 0; iter < maxiter; iter++) {
+
ntimestep = ++update->ntimestep;
niter++;
diff --git a/src/min_fire.cpp b/src/min_fire.cpp
index 8d0debf349..115e91e63c 100644
--- a/src/min_fire.cpp
+++ b/src/min_fire.cpp
@@ -92,6 +92,7 @@ int MinFire::iterate(int maxiter)
alpha_final = 0.0;
for (int iter = 0; iter < maxiter; iter++) {
+
ntimestep = ++update->ntimestep;
niter++;
diff --git a/src/min_quickmin.cpp b/src/min_quickmin.cpp
index 124b5bf575..84abf5487f 100644
--- a/src/min_quickmin.cpp
+++ b/src/min_quickmin.cpp
@@ -87,6 +87,7 @@ int MinQuickMin::iterate(int maxiter)
alpha_final = 0.0;
for (int iter = 0; iter < maxiter; iter++) {
+
ntimestep = ++update->ntimestep;
niter++;
diff --git a/src/min_sd.cpp b/src/min_sd.cpp
index 44936ce32a..1394f379fa 100644
--- a/src/min_sd.cpp
+++ b/src/min_sd.cpp
@@ -57,6 +57,7 @@ int MinSD::iterate(int maxiter)
for (i = 0; i < nextra_global; i++) hextra[i] = fextra[i];
for (int iter = 0; iter < maxiter; iter++) {
+
ntimestep = ++update->ntimestep;
niter++;
diff --git a/src/modify.cpp b/src/modify.cpp
index 0745c1c6ca..9ca992738a 100644
--- a/src/modify.cpp
+++ b/src/modify.cpp
@@ -42,7 +42,7 @@ Modify::Modify(LAMMPS *lmp) : Pointers(lmp)
nfix = maxfix = 0;
n_initial_integrate = n_post_integrate = 0;
n_pre_exchange = n_pre_neighbor = 0;
- n_pre_force = n_post_force = 0;
+ n_pre_force = n_pre_reverse = n_post_force = 0;
n_final_integrate = n_end_of_step = n_thermo_energy = 0;
n_initial_integrate_respa = n_post_integrate_respa = 0;
n_pre_force_respa = n_post_force_respa = n_final_integrate_respa = 0;
@@ -52,7 +52,7 @@ Modify::Modify(LAMMPS *lmp) : Pointers(lmp)
fmask = NULL;
list_initial_integrate = list_post_integrate = NULL;
list_pre_exchange = list_pre_neighbor = NULL;
- list_pre_force = list_post_force = NULL;
+ list_pre_force = list_pre_reverse = list_post_force = NULL;
list_final_integrate = list_end_of_step = NULL;
list_thermo_energy = NULL;
list_initial_integrate_respa = list_post_integrate_respa = NULL;
@@ -119,6 +119,7 @@ Modify::~Modify()
delete [] list_pre_exchange;
delete [] list_pre_neighbor;
delete [] list_pre_force;
+ delete [] list_pre_reverse;
delete [] list_post_force;
delete [] list_final_integrate;
delete [] list_end_of_step;
@@ -162,6 +163,7 @@ void Modify::init()
list_init(PRE_EXCHANGE,n_pre_exchange,list_pre_exchange);
list_init(PRE_NEIGHBOR,n_pre_neighbor,list_pre_neighbor);
list_init(PRE_FORCE,n_pre_force,list_pre_force);
+ list_init(PRE_REVERSE,n_pre_reverse,list_pre_reverse);
list_init(POST_FORCE,n_post_force,list_post_force);
list_init(FINAL_INTEGRATE,n_final_integrate,list_final_integrate);
list_init_end_of_step(END_OF_STEP,n_end_of_step,list_end_of_step);
@@ -382,6 +384,15 @@ void Modify::pre_force(int vflag)
for (int i = 0; i < n_pre_force; i++)
fix[list_pre_force[i]]->pre_force(vflag);
}
+/* ----------------------------------------------------------------------
+ pre_reverse call, only for relevant fixes
+------------------------------------------------------------------------- */
+
+void Modify::pre_reverse(int eflag, int vflag)
+{
+ for (int i = 0; i < n_pre_reverse; i++)
+ fix[list_pre_reverse[i]]->pre_reverse(eflag,vflag);
+}
/* ----------------------------------------------------------------------
post_force call, only for relevant fixes
@@ -789,8 +800,8 @@ void Modify::add_fix(int narg, char **arg, int trysuffix)
strcmp(style_restart_global[i],fix[ifix]->style) == 0) {
fix[ifix]->restart(state_restart_global[i]);
if (comm->me == 0) {
- char *str = (char *) ("Resetting global state of Fix %s Style %s "
- "from restart file info\n");
+ const char *str = (const char *) ("Resetting global state of Fix %s "
+ "Style %s from restart file info\n");
if (screen) fprintf(screen,str,fix[ifix]->id,fix[ifix]->style);
if (logfile) fprintf(logfile,str,fix[ifix]->id,fix[ifix]->style);
}
diff --git a/src/modify.h b/src/modify.h
index 7c11714a74..12a475cf0b 100644
--- a/src/modify.h
+++ b/src/modify.h
@@ -26,7 +26,7 @@ class Modify : protected Pointers {
public:
int nfix,maxfix;
int n_initial_integrate,n_post_integrate,n_pre_exchange,n_pre_neighbor;
- int n_pre_force,n_post_force;
+ int n_pre_force,n_pre_reverse,n_post_force;
int n_final_integrate,n_end_of_step,n_thermo_energy;
int n_initial_integrate_respa,n_post_integrate_respa;
int n_pre_force_respa,n_post_force_respa,n_final_integrate_respa;
@@ -55,6 +55,7 @@ class Modify : protected Pointers {
virtual void pre_exchange();
virtual void pre_neighbor();
virtual void pre_force(int);
+ virtual void pre_reverse(int,int);
virtual void post_force(int);
virtual void final_integrate();
virtual void end_of_step();
@@ -111,7 +112,7 @@ class Modify : protected Pointers {
int *list_initial_integrate,*list_post_integrate;
int *list_pre_exchange,*list_pre_neighbor;
- int *list_pre_force,*list_post_force;
+ int *list_pre_force,*list_pre_reverse,*list_post_force;
int *list_final_integrate,*list_end_of_step,*list_thermo_energy;
int *list_initial_integrate_respa,*list_post_integrate_respa;
int *list_pre_force_respa,*list_post_force_respa;
diff --git a/src/read_restart.h b/src/read_restart.h
index 1d84cf238e..586a13c210 100644
--- a/src/read_restart.h
+++ b/src/read_restart.h
@@ -33,7 +33,6 @@ class ReadRestart : protected Pointers {
private:
int me,nprocs,nprocs_file,multiproc_file;
FILE *fp;
- int nfix_restart_global,nfix_restart_peratom;
int multiproc; // 0 = proc 0 writes for all
// else # of procs writing files
@@ -42,7 +41,6 @@ class ReadRestart : protected Pointers {
int mpiioflag; // 1 for MPIIO output, else 0
class RestartMPIIO *mpiio; // MPIIO for restart file input
- int numChunksAssigned;
bigint assignedChunkSize;
MPI_Offset assignedChunkOffset,headerOffset;
diff --git a/src/respa.cpp b/src/respa.cpp
index 632d2a109f..dbd3effb63 100644
--- a/src/respa.cpp
+++ b/src/respa.cpp
@@ -453,6 +453,7 @@ void Respa::setup()
if (kspace_compute_flag) force->kspace->compute(eflag,vflag);
}
+ modify->pre_reverse(eflag,vflag);
if (newton[ilevel]) comm->reverse_comm();
copy_f_flevel(ilevel);
}
@@ -527,6 +528,8 @@ void Respa::setup_minimal(int flag)
force->kspace->setup();
if (kspace_compute_flag) force->kspace->compute(eflag,vflag);
}
+
+ modify->pre_reverse(eflag,vflag);
if (newton[ilevel]) comm->reverse_comm();
copy_f_flevel(ilevel);
}
@@ -711,6 +714,11 @@ void Respa::recurse(int ilevel)
timer->stamp(Timer::KSPACE);
}
+ if (modify->n_pre_reverse) {
+ modify->pre_reverse(eflag,vflag);
+ timer->stamp(Timer::MODIFY);
+ }
+
if (newton[ilevel]) {
comm->reverse_comm();
timer->stamp(Timer::COMM);
diff --git a/src/set.cpp b/src/set.cpp
index 914d725d81..9919518294 100644
--- a/src/set.cpp
+++ b/src/set.cpp
@@ -565,14 +565,16 @@ void Set::set(int keyword)
// overwrite dvalue, ivalue, xyzw value if variables defined
// else the input script scalar value remains in place
-
- if (varflag1) {
- dvalue = xvalue = vec1[i];
- ivalue = static_cast<int> (dvalue);
+
+ if (varflag) {
+ if (varflag1) {
+ dvalue = xvalue = vec1[i];
+ ivalue = static_cast<int> (dvalue);
+ }
+ if (varflag2) yvalue = vec2[i];
+ if (varflag3) zvalue = vec3[i];
+ if (varflag4) wvalue = vec4[i];
}
- if (varflag2) yvalue = vec2[i];
- if (varflag3) zvalue = vec3[i];
- if (varflag4) wvalue = vec4[i];
// set values in per-atom arrays
// error check here in case atom-style variables generated bogus value
diff --git a/src/verlet.cpp b/src/verlet.cpp
index 345549d914..e5a067777e 100644
--- a/src/verlet.cpp
+++ b/src/verlet.cpp
@@ -139,6 +139,7 @@ void Verlet::setup()
else force->kspace->compute_dummy(eflag,vflag);
}
+ modify->pre_reverse(eflag,vflag);
if (force->newton) comm->reverse_comm();
modify->setup(vflag);
@@ -199,6 +200,7 @@ void Verlet::setup_minimal(int flag)
else force->kspace->compute_dummy(eflag,vflag);
}
+ modify->pre_reverse(eflag,vflag);
if (force->newton) comm->reverse_comm();
modify->setup(vflag);
@@ -218,6 +220,7 @@ void Verlet::run(int n)
int n_pre_exchange = modify->n_pre_exchange;
int n_pre_neighbor = modify->n_pre_neighbor;
int n_pre_force = modify->n_pre_force;
+ int n_pre_reverse = modify->n_pre_reverse;
int n_post_force = modify->n_post_force;
int n_end_of_step = modify->n_end_of_step;
@@ -303,6 +306,11 @@ void Verlet::run(int n)
timer->stamp(Timer::KSPACE);
}
+ if (n_pre_reverse) {
+ modify->pre_reverse(eflag,vflag);
+ timer->stamp(Timer::MODIFY);
+ }
+
// reverse communication of forces
if (force->newton) {
--
GitLab