From f6a819580c6b0fc145b01ba8ef0fcf924f30e11f Mon Sep 17 00:00:00 2001
From: Steve Plimpton <sjplimp@sandia.gov>
Date: Fri, 6 Jan 2017 09:57:27 -0700
Subject: [PATCH] pair TIP4P bug fix for cutoffs >> box size
---
src/KOKKOS/modify_kokkos.cpp | 70 +++++++++++++++++-
src/KOKKOS/modify_kokkos.h | 3 +
src/KSPACE/pair_lj_cut_tip4p_long.cpp | 40 ++++++----
src/KSPACE/pair_lj_long_tip4p_long.cpp | 40 ++++++----
src/KSPACE/pair_tip4p_long.cpp | 40 ++++++----
src/KSPACE/pppm_disp_tip4p.cpp | 7 +-
src/KSPACE/pppm_tip4p.cpp | 7 +-
src/MOLECULE/fix_cmap.cpp | 18 ++++-
src/MOLECULE/fix_cmap.h | 1 +
src/MOLECULE/pair_lj_cut_tip4p_cut.cpp | 42 ++++++-----
src/MOLECULE/pair_tip4p_cut.cpp | 40 ++++++----
src/OPT/pair_lj_cut_tip4p_long_opt.cpp | 38 ++++++----
src/QEQ/fix_qeq.cpp | 7 --
src/QEQ/fix_qeq.h | 1 -
src/REPLICA/compute_event_displace.cpp | 5 +-
src/REPLICA/fix_neb.cpp | 2 +-
src/USER-ATC/fix_atc.cpp | 23 ------
src/USER-ATC/fix_atc.h | 2 -
src/USER-DPD/fix_shardlow.cpp | 7 --
src/USER-DPD/fix_shardlow.h | 1 -
src/USER-FEP/pair_lj_cut_tip4p_long_soft.cpp | 40 ++++++----
src/USER-FEP/pair_tip4p_long_soft.cpp | 40 ++++++----
src/USER-INTEL/fix_intel.cpp | 7 ++
src/USER-INTEL/fix_intel.h | 2 +
src/USER-OMP/pair_lj_cut_tip4p_cut_omp.cpp | 47 ++++++------
src/USER-OMP/pair_lj_cut_tip4p_long_omp.cpp | 46 ++++++------
.../pair_lj_cut_tip4p_long_soft_omp.cpp | 46 ++++++------
src/USER-OMP/pair_lj_long_tip4p_long_omp.cpp | 64 +++++++++-------
src/USER-OMP/pair_tip4p_cut_omp.cpp | 47 ++++++------
src/USER-OMP/pair_tip4p_long_omp.cpp | 48 ++++++------
src/USER-OMP/pair_tip4p_long_soft_omp.cpp | 48 ++++++------
src/USER-OMP/pppm_disp_tip4p_omp.cpp | 7 +-
src/USER-OMP/pppm_tip4p_omp.cpp | 7 +-
src/USER-REAXC/fix_qeq_reax.cpp | 7 --
src/USER-REAXC/fix_qeq_reax.h | 1 -
src/domain.cpp | 40 ++++++----
src/finish.cpp | 73 +++++++++++++++++--
src/fix.h | 40 +++++-----
src/fix_restrain.cpp | 3 +-
src/min.cpp | 13 +++-
src/modify.cpp | 38 ++++++++--
src/modify.h | 6 +-
src/respa.cpp | 10 ++-
src/run.cpp | 1 +
src/verlet.cpp | 4 +-
45 files changed, 663 insertions(+), 416 deletions(-)
diff --git a/src/KOKKOS/modify_kokkos.cpp b/src/KOKKOS/modify_kokkos.cpp
index 9a035ac20c..ec3831dff8 100644
--- a/src/KOKKOS/modify_kokkos.cpp
+++ b/src/KOKKOS/modify_kokkos.cpp
@@ -81,7 +81,7 @@ void ModifyKokkos::setup_pre_exchange()
atomKK->sync(fix[list_min_pre_exchange[i]]->execution_space,
fix[list_min_pre_exchange[i]]->datamask_read);
if (!fix[list_min_pre_exchange[i]]->kokkosable) lmp->kokkos->auto_sync = 1;
- fix[list_min_pre_exchange[i]]->min_setup_pre_exchange();
+ fix[list_min_pre_exchange[i]]->setup_pre_exchange();
lmp->kokkos->auto_sync = 0;
atomKK->modified(fix[list_min_pre_exchange[i]]->execution_space,
fix[list_min_pre_exchange[i]]->datamask_modify);
@@ -110,7 +110,7 @@ void ModifyKokkos::setup_pre_neighbor()
atomKK->sync(fix[list_min_pre_neighbor[i]]->execution_space,
fix[list_min_pre_neighbor[i]]->datamask_read);
if (!fix[list_min_pre_neighbor[i]]->kokkosable) lmp->kokkos->auto_sync = 1;
- fix[list_min_pre_neighbor[i]]->min_setup_pre_neighbor();
+ fix[list_min_pre_neighbor[i]]->setup_pre_neighbor();
lmp->kokkos->auto_sync = 0;
atomKK->modified(fix[list_min_pre_neighbor[i]]->execution_space,
fix[list_min_pre_neighbor[i]]->datamask_modify);
@@ -139,13 +139,42 @@ void ModifyKokkos::setup_pre_force(int vflag)
atomKK->sync(fix[list_min_pre_force[i]]->execution_space,
fix[list_min_pre_force[i]]->datamask_read);
if (!fix[list_min_pre_force[i]]->kokkosable) lmp->kokkos->auto_sync = 1;
- fix[list_min_pre_force[i]]->min_setup_pre_force(vflag);
+ fix[list_min_pre_force[i]]->setup_pre_force(vflag);
lmp->kokkos->auto_sync = 0;
atomKK->modified(fix[list_min_pre_force[i]]->execution_space,
fix[list_min_pre_force[i]]->datamask_modify);
}
}
+/* ----------------------------------------------------------------------
+ setup pre_reverse call, only for fixes that define pre_reverse
+ called from Verlet, RESPA, Min
+------------------------------------------------------------------------- */
+
+void ModifyKokkos::setup_pre_reverse(int eflag, int vflag)
+{
+ if (update->whichflag == 1)
+ for (int i = 0; i < n_pre_reverse; i++) {
+ atomKK->sync(fix[list_pre_reverse[i]]->execution_space,
+ fix[list_pre_reverse[i]]->datamask_read);
+ if (!fix[list_pre_reverse[i]]->kokkosable) lmp->kokkos->auto_sync = 1;
+ fix[list_pre_reverse[i]]->setup_pre_reverse(eflag,vflag);
+ lmp->kokkos->auto_sync = 0;
+ atomKK->modified(fix[list_pre_reverse[i]]->execution_space,
+ fix[list_pre_reverse[i]]->datamask_modify);
+ }
+ else if (update->whichflag == 2)
+ for (int i = 0; i < n_min_pre_reverse; i++) {
+ atomKK->sync(fix[list_min_pre_reverse[i]]->execution_space,
+ fix[list_min_pre_reverse[i]]->datamask_read);
+ if (!fix[list_min_pre_reverse[i]]->kokkosable) lmp->kokkos->auto_sync = 1;
+ fix[list_min_pre_reverse[i]]->setup_pre_reverse(eflag,vflag);
+ lmp->kokkos->auto_sync = 0;
+ atomKK->modified(fix[list_min_pre_reverse[i]]->execution_space,
+ fix[list_min_pre_reverse[i]]->datamask_modify);
+ }
+}
+
/* ----------------------------------------------------------------------
1st half of integrate call, only for relevant fixes
------------------------------------------------------------------------- */
@@ -231,6 +260,23 @@ void ModifyKokkos::pre_force(int vflag)
}
}
+/* ----------------------------------------------------------------------
+ pre_reverse call, only for relevant fixes
+------------------------------------------------------------------------- */
+
+void ModifyKokkos::pre_reverse(int eflag, int vflag)
+{
+ for (int i = 0; i < n_pre_reverse; i++) {
+ atomKK->sync(fix[list_pre_reverse[i]]->execution_space,
+ fix[list_pre_reverse[i]]->datamask_read);
+ if (!fix[list_pre_reverse[i]]->kokkosable) lmp->kokkos->auto_sync = 1;
+ fix[list_pre_reverse[i]]->pre_reverse(eflag,vflag);
+ lmp->kokkos->auto_sync = 0;
+ atomKK->modified(fix[list_pre_reverse[i]]->execution_space,
+ fix[list_pre_reverse[i]]->datamask_modify);
+ }
+}
+
/* ----------------------------------------------------------------------
post_force call, only for relevant fixes
------------------------------------------------------------------------- */
@@ -476,6 +522,23 @@ void ModifyKokkos::min_pre_force(int vflag)
}
}
+/* ----------------------------------------------------------------------
+ minimizer pre-reverse call, only for relevant fixes
+------------------------------------------------------------------------- */
+
+void ModifyKokkos::min_pre_reverse(int eflag, int vflag)
+{
+ for (int i = 0; i < n_min_pre_reverse; i++) {
+ atomKK->sync(fix[list_min_pre_reverse[i]]->execution_space,
+ fix[list_min_pre_reverse[i]]->datamask_read);
+ if (!fix[list_min_pre_reverse[i]]->kokkosable) lmp->kokkos->auto_sync = 1;
+ fix[list_min_pre_reverse[i]]->min_pre_reverse(eflag,vflag);
+ lmp->kokkos->auto_sync = 0;
+ atomKK->modified(fix[list_min_pre_reverse[i]]->execution_space,
+ fix[list_min_pre_reverse[i]]->datamask_modify);
+ }
+}
+
/* ----------------------------------------------------------------------
minimizer force adjustment call, only for relevant fixes
------------------------------------------------------------------------- */
@@ -658,4 +721,3 @@ int ModifyKokkos::min_reset_ref()
}
return itmpall;
}
-
diff --git a/src/KOKKOS/modify_kokkos.h b/src/KOKKOS/modify_kokkos.h
index b8e5f8de30..32dfb2fd97 100644
--- a/src/KOKKOS/modify_kokkos.h
+++ b/src/KOKKOS/modify_kokkos.h
@@ -26,12 +26,14 @@ class ModifyKokkos : public Modify {
void setup_pre_exchange();
void setup_pre_neighbor();
void setup_pre_force(int);
+ void setup_pre_reverse(int, int);
void initial_integrate(int);
void post_integrate();
void pre_decide();
void pre_exchange();
void pre_neighbor();
void pre_force(int);
+ void pre_reverse(int,int);
void post_force(int);
void final_integrate();
void end_of_step();
@@ -48,6 +50,7 @@ class ModifyKokkos : public Modify {
void min_pre_exchange();
void min_pre_neighbor();
void min_pre_force(int);
+ void min_pre_reverse(int,int);
void min_post_force(int);
double min_energy(double *);
diff --git a/src/KSPACE/pair_lj_cut_tip4p_long.cpp b/src/KSPACE/pair_lj_cut_tip4p_long.cpp
index 6e855e41f0..146d4e6f37 100644
--- a/src/KSPACE/pair_lj_cut_tip4p_long.cpp
+++ b/src/KSPACE/pair_lj_cut_tip4p_long.cpp
@@ -88,8 +88,8 @@ void PairLJCutTIP4PLong::compute(int eflag, int vflag)
double r,r2inv,r6inv,forcecoul,forcelj,cforce;
double factor_coul,factor_lj;
double grij,expm2,prefactor,t,erfc;
- double fO[3],fH[3],fd[3],v[6],xH1[3],xH2[3];
- double *x1,*x2;
+ double fO[3],fH[3],fd[3],v[6];
+ double *x1,*x2,*xH1,*xH2;
int *ilist,*jlist,*numneigh,**firstneigh;
double rsq;
@@ -146,14 +146,20 @@ void PairLJCutTIP4PLong::compute(int eflag, int vflag)
if (itype == typeO) {
if (hneigh[i][0] < 0) {
- hneigh[i][0] = iH1 = atom->map(tag[i] + 1);
- hneigh[i][1] = iH2 = atom->map(tag[i] + 2);
- hneigh[i][2] = 1;
+ iH1 = atom->map(tag[i] + 1);
+ iH2 = atom->map(tag[i] + 2);
if (iH1 == -1 || iH2 == -1)
error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+ // set iH1,iH2 to closest image to O
+ iH1 = domain->closest_image(i,iH1);
+ iH2 = domain->closest_image(i,iH2);
compute_newsite(x[i],x[iH1],x[iH2],newsite[i]);
+ hneigh[i][0] = iH1;
+ hneigh[i][1] = iH2;
+ hneigh[i][2] = 1;
+
} else {
iH1 = hneigh[i][0];
iH2 = hneigh[i][1];
@@ -216,14 +222,20 @@ void PairLJCutTIP4PLong::compute(int eflag, int vflag)
if (jtype == typeO) {
if (hneigh[j][0] < 0) {
- hneigh[j][0] = jH1 = atom->map(tag[j] + 1);
- hneigh[j][1] = jH2 = atom->map(tag[j] + 2);
- hneigh[j][2] = 1;
+ jH1 = atom->map(tag[j] + 1);
+ jH2 = atom->map(tag[j] + 2);
if (jH1 == -1 || jH2 == -1)
error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[jH1] != typeH || atom->type[jH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+ // set jH1,jH2 to closest image to O
+ jH1 = domain->closest_image(j,jH1);
+ jH2 = domain->closest_image(j,jH2);
compute_newsite(x[j],x[jH1],x[jH2],newsite[j]);
+ hneigh[j][0] = jH1;
+ hneigh[j][1] = jH2;
+ hneigh[j][2] = 1;
+
} else {
jH1 = hneigh[j][0];
jH2 = hneigh[j][1];
@@ -327,9 +339,8 @@ void PairLJCutTIP4PLong::compute(int eflag, int vflag)
f[iH2][2] += fH[2];
if (vflag) {
- domain->closest_image(x[i],x[iH1],xH1);
- domain->closest_image(x[i],x[iH2],xH2);
-
+ xH1 = x[iH1];
+ xH2 = x[iH2];
v[0] = x[i][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0];
v[1] = x[i][1]*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1];
v[2] = x[i][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2];
@@ -385,9 +396,8 @@ void PairLJCutTIP4PLong::compute(int eflag, int vflag)
f[jH2][2] += fH[2];
if (vflag) {
- domain->closest_image(x[j],x[jH1],xH1);
- domain->closest_image(x[j],x[jH2],xH2);
-
+ xH1 = x[jH1];
+ xH2 = x[jH2];
v[0] += x[j][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0];
v[1] += x[j][1]*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1];
v[2] += x[j][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2];
@@ -565,12 +575,10 @@ void PairLJCutTIP4PLong::compute_newsite(double *xO, double *xH1,
double delx1 = xH1[0] - xO[0];
double dely1 = xH1[1] - xO[1];
double delz1 = xH1[2] - xO[2];
- domain->minimum_image(delx1,dely1,delz1);
double delx2 = xH2[0] - xO[0];
double dely2 = xH2[1] - xO[1];
double delz2 = xH2[2] - xO[2];
- domain->minimum_image(delx2,dely2,delz2);
xM[0] = xO[0] + alpha * 0.5 * (delx1 + delx2);
xM[1] = xO[1] + alpha * 0.5 * (dely1 + dely2);
diff --git a/src/KSPACE/pair_lj_long_tip4p_long.cpp b/src/KSPACE/pair_lj_long_tip4p_long.cpp
index b8ce8ece4c..c3d95c37a6 100644
--- a/src/KSPACE/pair_lj_long_tip4p_long.cpp
+++ b/src/KSPACE/pair_lj_long_tip4p_long.cpp
@@ -86,8 +86,8 @@ void PairLJLongTIP4PLong::compute(int eflag, int vflag)
double r,r2inv,forcecoul,forcelj,cforce;
double factor_coul;
double grij,expm2,prefactor,t,erfc;
- double fO[3],fH[3],fd[3],v[6],xH1[3],xH2[3];// f1[3];
- double *x1,*x2;
+ double fO[3],fH[3],fd[3],v[6];
+ double *x1,*x2,*xH1,*xH2;
int *ilist,*jlist,*numneigh,**firstneigh;
double rsq;
@@ -145,14 +145,20 @@ void PairLJLongTIP4PLong::compute(int eflag, int vflag)
itype = type[i];
if (itype == typeO) {
if (hneigh[i][0] < 0) {
- hneigh[i][0] = iH1 = atom->map(tag[i] + 1);
- hneigh[i][1] = iH2 = atom->map(tag[i] + 2);
- hneigh[i][2] = 1;
+ iH1 = atom->map(tag[i] + 1);
+ iH2 = atom->map(tag[i] + 2);
if (iH1 == -1 || iH2 == -1)
error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+ // set iH1,iH2 to closest image to O
+ iH1 = domain->closest_image(i,iH1);
+ iH2 = domain->closest_image(i,iH2);
compute_newsite(x[i],x[iH1],x[iH2],newsite[i]);
+ hneigh[i][0] = iH1;
+ hneigh[i][1] = iH2;
+ hneigh[i][2] = 1;
+
} else {
iH1 = hneigh[i][0];
iH2 = hneigh[i][1];
@@ -253,14 +259,20 @@ void PairLJLongTIP4PLong::compute(int eflag, int vflag)
if (itype == typeO || jtype == typeO) {
if (jtype == typeO) {
if (hneigh[j][0] < 0) {
- hneigh[j][0] = jH1 = atom->map(tag[j] + 1);
- hneigh[j][1] = jH2 = atom->map(tag[j] + 2);
- hneigh[j][2] = 1;
+ jH1 = atom->map(tag[j] + 1);
+ jH2 = atom->map(tag[j] + 2);
if (jH1 == -1 || jH2 == -1)
error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[jH1] != typeH || atom->type[jH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+ // set jH1,jH2 to closest image to O
+ jH1 = domain->closest_image(j,jH1);
+ jH2 = domain->closest_image(j,jH2);
compute_newsite(x[j],x[jH1],x[jH2],newsite[j]);
+ hneigh[j][0] = jH1;
+ hneigh[j][1] = jH2;
+ hneigh[j][2] = 1;
+
} else {
jH1 = hneigh[j][0];
jH2 = hneigh[j][1];
@@ -365,9 +377,8 @@ void PairLJLongTIP4PLong::compute(int eflag, int vflag)
f[iH2][2] += fH[2];
if (vflag) {
- domain->closest_image(x[i],x[iH1],xH1);
- domain->closest_image(x[i],x[iH2],xH2);
-
+ xH1 = x[iH1];
+ xH2 = x[iH2];
v[0] = x[i][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0];
v[1] = x[i][1]*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1];
v[2] = x[i][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2];
@@ -423,9 +434,8 @@ void PairLJLongTIP4PLong::compute(int eflag, int vflag)
f[jH2][2] += fH[2];
if (vflag) {
- domain->closest_image(x[j],x[jH1],xH1);
- domain->closest_image(x[j],x[jH2],xH2);
-
+ xH1 = x[jH1];
+ xH2 = x[jH2];
v[0] += x[j][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0];
v[1] += x[j][1]*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1];
v[2] += x[j][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2];
@@ -1548,12 +1558,10 @@ void PairLJLongTIP4PLong::compute_newsite(double *xO, double *xH1,
double delx1 = xH1[0] - xO[0];
double dely1 = xH1[1] - xO[1];
double delz1 = xH1[2] - xO[2];
- domain->minimum_image(delx1,dely1,delz1);
double delx2 = xH2[0] - xO[0];
double dely2 = xH2[1] - xO[1];
double delz2 = xH2[2] - xO[2];
- domain->minimum_image(delx2,dely2,delz2);
xM[0] = xO[0] + alpha * 0.5 * (delx1 + delx2);
xM[1] = xO[1] + alpha * 0.5 * (dely1 + dely2);
diff --git a/src/KSPACE/pair_tip4p_long.cpp b/src/KSPACE/pair_tip4p_long.cpp
index 0753a3f34d..1dd57ce822 100644
--- a/src/KSPACE/pair_tip4p_long.cpp
+++ b/src/KSPACE/pair_tip4p_long.cpp
@@ -84,8 +84,8 @@ void PairTIP4PLong::compute(int eflag, int vflag)
double r,r2inv,forcecoul,cforce;
double factor_coul;
double grij,expm2,prefactor,t,erfc;
- double fO[3],fH[3],fd[3],v[6],xH1[3],xH2[3];
- double *x1,*x2;
+ double fO[3],fH[3],fd[3],v[6];
+ double *x1,*x2,*xH1,*xH2;
int *ilist,*jlist,*numneigh,**firstneigh;
double rsq;
@@ -140,14 +140,20 @@ void PairTIP4PLong::compute(int eflag, int vflag)
if (itype == typeO) {
if (hneigh[i][0] < 0) {
- hneigh[i][0] = iH1 = atom->map(tag[i] + 1);
- hneigh[i][1] = iH2 = atom->map(tag[i] + 2);
- hneigh[i][2] = 1;
+ iH1 = atom->map(tag[i] + 1);
+ iH2 = atom->map(tag[i] + 2);
if (iH1 == -1 || iH2 == -1)
error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+ // set iH1,iH2 to closest image to O
+ iH1 = domain->closest_image(i,iH1);
+ iH2 = domain->closest_image(i,iH2);
compute_newsite(x[i],x[iH1],x[iH2],newsite[i]);
+ hneigh[i][0] = iH1;
+ hneigh[i][1] = iH2;
+ hneigh[i][2] = 1;
+
} else {
iH1 = hneigh[i][0];
iH2 = hneigh[i][1];
@@ -184,14 +190,20 @@ void PairTIP4PLong::compute(int eflag, int vflag)
if (jtype == typeO) {
if (hneigh[j][0] < 0) {
- hneigh[j][0] = jH1 = atom->map(tag[j] + 1);
- hneigh[j][1] = jH2 = atom->map(tag[j] + 2);
- hneigh[j][2] = 1;
+ jH1 = atom->map(tag[j] + 1);
+ jH2 = atom->map(tag[j] + 2);
if (jH1 == -1 || jH2 == -1)
error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[jH1] != typeH || atom->type[jH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+ // set jH1,jH2 to closest image to O
+ jH1 = domain->closest_image(j,jH1);
+ jH2 = domain->closest_image(j,jH2);
compute_newsite(x[j],x[jH1],x[jH2],newsite[j]);
+ hneigh[j][0] = jH1;
+ hneigh[j][1] = jH2;
+ hneigh[j][2] = 1;
+
} else {
jH1 = hneigh[j][0];
jH2 = hneigh[j][1];
@@ -295,9 +307,8 @@ void PairTIP4PLong::compute(int eflag, int vflag)
f[iH2][2] += fH[2];
if (vflag) {
- domain->closest_image(x[i],x[iH1],xH1);
- domain->closest_image(x[i],x[iH2],xH2);
-
+ xH1 = x[iH1];
+ xH2 = x[iH2];
v[0] = x[i][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0];
v[1] = x[i][1]*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1];
v[2] = x[i][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2];
@@ -353,9 +364,8 @@ void PairTIP4PLong::compute(int eflag, int vflag)
f[jH2][2] += fH[2];
if (vflag) {
- domain->closest_image(x[j],x[jH1],xH1);
- domain->closest_image(x[j],x[jH2],xH2);
-
+ xH1 = x[jH1];
+ xH2 = x[jH2];
v[0] += x[j][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0];
v[1] += x[j][1]*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1];
v[2] += x[j][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2];
@@ -488,12 +498,10 @@ void PairTIP4PLong::compute_newsite(double *xO, double *xH1,
double delx1 = xH1[0] - xO[0];
double dely1 = xH1[1] - xO[1];
double delz1 = xH1[2] - xO[2];
- domain->minimum_image(delx1,dely1,delz1);
double delx2 = xH2[0] - xO[0];
double dely2 = xH2[1] - xO[1];
double delz2 = xH2[2] - xO[2];
- domain->minimum_image(delx2,dely2,delz2);
xM[0] = xO[0] + alpha * 0.5 * (delx1 + delx2);
xM[1] = xO[1] + alpha * 0.5 * (dely1 + dely2);
diff --git a/src/KSPACE/pppm_disp_tip4p.cpp b/src/KSPACE/pppm_disp_tip4p.cpp
index d5a85abb4c..cf3a3943b9 100644
--- a/src/KSPACE/pppm_disp_tip4p.cpp
+++ b/src/KSPACE/pppm_disp_tip4p.cpp
@@ -502,17 +502,20 @@ void PPPMDispTIP4P::find_M(int i, int &iH1, int &iH2, double *xM)
if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+ // set iH1,iH2 to index of closest image to O
+
+ iH1 = domain->closest_image(i,iH1);
+ iH2 = domain->closest_image(i,iH2);
+
double **x = atom->x;
double delx1 = x[iH1][0] - x[i][0];
double dely1 = x[iH1][1] - x[i][1];
double delz1 = x[iH1][2] - x[i][2];
- domain->minimum_image(delx1,dely1,delz1);
double delx2 = x[iH2][0] - x[i][0];
double dely2 = x[iH2][1] - x[i][1];
double delz2 = x[iH2][2] - x[i][2];
- domain->minimum_image(delx2,dely2,delz2);
xM[0] = x[i][0] + alpha * 0.5 * (delx1 + delx2);
xM[1] = x[i][1] + alpha * 0.5 * (dely1 + dely2);
diff --git a/src/KSPACE/pppm_tip4p.cpp b/src/KSPACE/pppm_tip4p.cpp
index 34d0616fe9..fa5820cb34 100644
--- a/src/KSPACE/pppm_tip4p.cpp
+++ b/src/KSPACE/pppm_tip4p.cpp
@@ -493,17 +493,20 @@ void PPPMTIP4P::find_M(int i, int &iH1, int &iH2, double *xM)
if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+ // set iH1,iH2 to index of closest image to O
+
+ iH1 = domain->closest_image(i,iH1);
+ iH2 = domain->closest_image(i,iH2);
+
double **x = atom->x;
double delx1 = x[iH1][0] - x[i][0];
double dely1 = x[iH1][1] - x[i][1];
double delz1 = x[iH1][2] - x[i][2];
- domain->minimum_image(delx1,dely1,delz1);
double delx2 = x[iH2][0] - x[i][0];
double dely2 = x[iH2][1] - x[i][1];
double delz2 = x[iH2][2] - x[i][2];
- domain->minimum_image(delx2,dely2,delz2);
xM[0] = x[i][0] + alpha * 0.5 * (delx1 + delx2);
xM[1] = x[i][1] + alpha * 0.5 * (dely1 + dely2);
diff --git a/src/MOLECULE/fix_cmap.cpp b/src/MOLECULE/fix_cmap.cpp
index da7d337b9b..26aaef60d7 100644
--- a/src/MOLECULE/fix_cmap.cpp
+++ b/src/MOLECULE/fix_cmap.cpp
@@ -62,10 +62,13 @@ using namespace MathConst;
/* ---------------------------------------------------------------------- */
-FixCMAP::FixCMAP(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg),
- crosstermlist(NULL), num_crossterm(NULL), crossterm_type(NULL), crossterm_atom1(NULL),
- crossterm_atom2(NULL), crossterm_atom3(NULL), crossterm_atom4(NULL), crossterm_atom5(NULL),
- g_axis(NULL), cmapgrid(NULL), d1cmapgrid(NULL), d2cmapgrid(NULL), d12cmapgrid(NULL)
+FixCMAP::FixCMAP(LAMMPS *lmp, int narg, char **arg) :
+ Fix(lmp, narg, arg),
+ crosstermlist(NULL), num_crossterm(NULL), crossterm_type(NULL),
+ crossterm_atom1(NULL), crossterm_atom2(NULL), crossterm_atom3(NULL),
+ crossterm_atom4(NULL), crossterm_atom5(NULL),
+ g_axis(NULL), cmapgrid(NULL), d1cmapgrid(NULL), d2cmapgrid(NULL),
+ d12cmapgrid(NULL)
{
if (narg != 4) error->all(FLERR,"Illegal fix cmap command");
@@ -208,6 +211,13 @@ void FixCMAP::setup_pre_neighbor()
/* --------------------------------------------------------------------- */
+void FixCMAP::setup_pre_reverse(int eflag, int vflag)
+{
+ pre_reverse(eflag,vflag);
+}
+
+/* --------------------------------------------------------------------- */
+
void FixCMAP::min_setup(int vflag)
{
pre_neighbor();
diff --git a/src/MOLECULE/fix_cmap.h b/src/MOLECULE/fix_cmap.h
index 3573d213e8..d384d61d7d 100644
--- a/src/MOLECULE/fix_cmap.h
+++ b/src/MOLECULE/fix_cmap.h
@@ -31,6 +31,7 @@ class FixCMAP : public Fix {
void init();
void setup(int);
void setup_pre_neighbor();
+ void setup_pre_reverse(int, int);
void min_setup(int);
void pre_neighbor();
void pre_reverse(int, int);
diff --git a/src/MOLECULE/pair_lj_cut_tip4p_cut.cpp b/src/MOLECULE/pair_lj_cut_tip4p_cut.cpp
index a9e00e80b4..15f5d52961 100644
--- a/src/MOLECULE/pair_lj_cut_tip4p_cut.cpp
+++ b/src/MOLECULE/pair_lj_cut_tip4p_cut.cpp
@@ -87,8 +87,8 @@ void PairLJCutTIP4PCut::compute(int eflag, int vflag)
int n,vlist[6];
int iH1,iH2,jH1,jH2;
double cforce;
- double fO[3],fH[3],fd[3],v[6],xH1[3],xH2[3];
- double *x1,*x2;
+ double fO[3],fH[3],fd[3],v[6];
+ double *x1,*x2,*xH1,*xH2;
evdwl = ecoul = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
@@ -139,14 +139,20 @@ void PairLJCutTIP4PCut::compute(int eflag, int vflag)
if (itype == typeO) {
if (hneigh[i][0] < 0) {
- hneigh[i][0] = iH1 = atom->map(tag[i] + 1);
- hneigh[i][1] = iH2 = atom->map(tag[i] + 2);
- hneigh[i][2] = 1;
+ iH1 = atom->map(tag[i] + 1);
+ iH2 = atom->map(tag[i] + 2);
if (iH1 == -1 || iH2 == -1)
error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+ // set iH1,iH2 to index of closest image to O
+ iH1 = domain->closest_image(i,iH1);
+ iH2 = domain->closest_image(i,iH2);
compute_newsite(x[i],x[iH1],x[iH2],newsite[i]);
+ hneigh[i][0] = iH1;
+ hneigh[i][1] = iH2;
+ hneigh[i][2] = 1;
+
} else {
iH1 = hneigh[i][0];
iH2 = hneigh[i][1];
@@ -209,14 +215,20 @@ void PairLJCutTIP4PCut::compute(int eflag, int vflag)
if (jtype == typeO) {
if (hneigh[j][0] < 0) {
- hneigh[j][0] = jH1 = atom->map(tag[j] + 1);
- hneigh[j][1] = jH2 = atom->map(tag[j] + 2);
- hneigh[j][2] = 1;
+ jH1 = atom->map(tag[j] + 1);
+ jH2 = atom->map(tag[j] + 2);
if (jH1 == -1 || jH2 == -1)
error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[jH1] != typeH || atom->type[jH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+ // set jH1,jH2 to closest image to O
+ jH1 = domain->closest_image(j,jH1);
+ jH2 = domain->closest_image(j,jH2);
compute_newsite(x[j],x[jH1],x[jH2],newsite[j]);
+ hneigh[j][0] = jH1;
+ hneigh[j][1] = jH2;
+ hneigh[j][2] = 1;
+
} else {
jH1 = hneigh[j][0];
jH2 = hneigh[j][1];
@@ -294,10 +306,9 @@ void PairLJCutTIP4PCut::compute(int eflag, int vflag)
f[iH2][1] += fH[1];
f[iH2][2] += fH[2];
- if(vflag) {
- domain->closest_image(x[i],x[iH1],xH1);
- domain->closest_image(x[i],x[iH2],xH2);
-
+ if (vflag) {
+ xH1 = x[iH1];
+ xH2 = x[iH2];
v[0] = x[i][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0];
v[1] = x[i][1]*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1];
v[2] = x[i][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2];
@@ -353,9 +364,8 @@ void PairLJCutTIP4PCut::compute(int eflag, int vflag)
f[jH2][2] += fH[2];
if (vflag) {
- domain->closest_image(x[j],x[jH1],xH1);
- domain->closest_image(x[j],x[jH2],xH2);
-
+ xH1 = x[jH1];
+ xH2 = x[jH2];
v[0] += x[j][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0];
v[1] += x[j][1]*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1];
v[2] += x[j][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2];
@@ -711,12 +721,10 @@ void PairLJCutTIP4PCut::compute_newsite(double *xO, double *xH1,
double delx1 = xH1[0] - xO[0];
double dely1 = xH1[1] - xO[1];
double delz1 = xH1[2] - xO[2];
- domain->minimum_image(delx1,dely1,delz1);
double delx2 = xH2[0] - xO[0];
double dely2 = xH2[1] - xO[1];
double delz2 = xH2[2] - xO[2];
- domain->minimum_image(delx2,dely2,delz2);
xM[0] = xO[0] + alpha * 0.5 * (delx1 + delx2);
xM[1] = xO[1] + alpha * 0.5 * (dely1 + dely2);
diff --git a/src/MOLECULE/pair_tip4p_cut.cpp b/src/MOLECULE/pair_tip4p_cut.cpp
index dd48637f23..954bd35435 100644
--- a/src/MOLECULE/pair_tip4p_cut.cpp
+++ b/src/MOLECULE/pair_tip4p_cut.cpp
@@ -75,8 +75,8 @@ void PairTIP4PCut::compute(int eflag, int vflag)
int n,vlist[6];
int iH1,iH2,jH1,jH2;
double cforce;
- double fO[3],fH[3],fd[3],v[6],xH1[3],xH2[3];
- double *x1,*x2;
+ double fO[3],fH[3],fd[3],v[6];
+ double *x1,*x2,*xH1,*xH2;
ecoul = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
@@ -125,14 +125,20 @@ void PairTIP4PCut::compute(int eflag, int vflag)
if (itype == typeO) {
if (hneigh[i][0] < 0) {
- hneigh[i][0] = iH1 = atom->map(tag[i] + 1);
- hneigh[i][1] = iH2 = atom->map(tag[i] + 2);
- hneigh[i][2] = 1;
+ iH1 = atom->map(tag[i] + 1);
+ iH2 = atom->map(tag[i] + 2);
if (iH1 == -1 || iH2 == -1)
error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+ // set iH1,iH2 to closest image to O
+ iH1 = domain->closest_image(i,iH1);
+ iH2 = domain->closest_image(i,iH2);
compute_newsite(x[i],x[iH1],x[iH2],newsite[i]);
+ hneigh[i][0] = iH1;
+ hneigh[i][1] = iH2;
+ hneigh[i][2] = 1;
+
} else {
iH1 = hneigh[i][0];
iH2 = hneigh[i][1];
@@ -169,14 +175,20 @@ void PairTIP4PCut::compute(int eflag, int vflag)
if (jtype == typeO) {
if (hneigh[j][0] < 0) {
- hneigh[j][0] = jH1 = atom->map(tag[j] + 1);
- hneigh[j][1] = jH2 = atom->map(tag[j] + 2);
- hneigh[j][2] = 1;
+ jH1 = atom->map(tag[j] + 1);
+ jH2 = atom->map(tag[j] + 2);
if (jH1 == -1 || jH2 == -1)
error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[jH1] != typeH || atom->type[jH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+ // set jH1,jH2 to closest image to O
+ jH1 = domain->closest_image(j,jH1);
+ jH2 = domain->closest_image(j,jH2);
compute_newsite(x[j],x[jH1],x[jH2],newsite[j]);
+ hneigh[j][0] = jH1;
+ hneigh[j][1] = jH2;
+ hneigh[j][2] = 1;
+
} else {
jH1 = hneigh[j][0];
jH2 = hneigh[j][1];
@@ -255,9 +267,8 @@ void PairTIP4PCut::compute(int eflag, int vflag)
f[iH2][2] += fH[2];
if(vflag) {
- domain->closest_image(x[i],x[iH1],xH1);
- domain->closest_image(x[i],x[iH2],xH2);
-
+ xH1 = x[iH1];
+ xH2 = x[iH2];
v[0] = x[i][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0];
v[1] = x[i][1]*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1];
v[2] = x[i][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2];
@@ -313,9 +324,8 @@ void PairTIP4PCut::compute(int eflag, int vflag)
f[jH2][2] += fH[2];
if (vflag) {
- domain->closest_image(x[j],x[jH1],xH1);
- domain->closest_image(x[j],x[jH2],xH2);
-
+ xH1 = x[jH1];
+ xH2 = x[jH2];
v[0] += x[j][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0];
v[1] += x[j][1]*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1];
v[2] += x[j][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2];
@@ -526,12 +536,10 @@ void PairTIP4PCut::compute_newsite(double *xO, double *xH1,
double delx1 = xH1[0] - xO[0];
double dely1 = xH1[1] - xO[1];
double delz1 = xH1[2] - xO[2];
- domain->minimum_image(delx1,dely1,delz1);
double delx2 = xH2[0] - xO[0];
double dely2 = xH2[1] - xO[1];
double delz2 = xH2[2] - xO[2];
- domain->minimum_image(delx2,dely2,delz2);
xM[0] = xO[0] + alpha * 0.5 * (delx1 + delx2);
xM[1] = xO[1] + alpha * 0.5 * (dely1 + dely2);
diff --git a/src/OPT/pair_lj_cut_tip4p_long_opt.cpp b/src/OPT/pair_lj_cut_tip4p_long_opt.cpp
index 2cca32baec..db9402b521 100644
--- a/src/OPT/pair_lj_cut_tip4p_long_opt.cpp
+++ b/src/OPT/pair_lj_cut_tip4p_long_opt.cpp
@@ -111,9 +111,9 @@ void PairLJCutTIP4PLongOpt::eval()
double r,rsq,r2inv,r6inv,forcecoul,forcelj,cforce;
double factor_coul,factor_lj;
double grij,expm2,prefactor,t,erfc;
- double v[6],xH1[3],xH2[3];
+ double v[6];
double fdx,fdy,fdz,fOx,fOy,fOz,fHx,fHy,fHz;
- const double *x1,*x2;
+ const double *x1,*x2,*xH1,*xH2;
int *ilist,*jlist,*numneigh,**firstneigh;
int i,j,ii,jj,inum,jnum,itype,jtype,itable,key;
@@ -155,14 +155,21 @@ void PairLJCutTIP4PLongOpt::eval()
if (itype == typeO) {
if (hneigh[i][0] < 0) {
- hneigh[i][0] = iH1 = atom->map(tag[i] + 1);
- hneigh[i][1] = iH2 = atom->map(tag[i] + 2);
+ iH1 = atom->map(tag[i] + 1);
+ iH2 = atom->map(tag[i] + 2);
hneigh[i][2] = 1;
if (iH1 == -1 || iH2 == -1)
error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+ // set iH1,iH2 to closest image to O
+ iH1 = domain->closest_image(i,iH1);
+ iH2 = domain->closest_image(i,iH2);
compute_newsite_opt(x[i],x[iH1],x[iH2],newsite[i]);
+ hneigh[i][0] = iH1;
+ hneigh[i][1] = iH2;
+ hneigh[i][2] = 1;
+
} else {
iH1 = hneigh[i][0];
iH2 = hneigh[i][1];
@@ -226,14 +233,21 @@ void PairLJCutTIP4PLongOpt::eval()
if (jtype == typeO) {
if (hneigh[j][0] < 0) {
- hneigh[j][0] = jH1 = atom->map(tag[j] + 1);
- hneigh[j][1] = jH2 = atom->map(tag[j] + 2);
+ jH1 = atom->map(tag[j] + 1);
+ jH2 = atom->map(tag[j] + 2);
hneigh[j][2] = 1;
if (jH1 == -1 || jH2 == -1)
error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[jH1] != typeH || atom->type[jH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+ // set jH1,jH2 to closest image to O
+ jH1 = domain->closest_image(j,jH1);
+ jH2 = domain->closest_image(j,jH2);
compute_newsite_opt(x[j],x[jH1],x[jH2],newsite[j]);
+ hneigh[j][0] = jH1;
+ hneigh[j][1] = jH2;
+ hneigh[j][2] = 1;
+
} else {
jH1 = hneigh[j][0];
jH2 = hneigh[j][1];
@@ -339,9 +353,8 @@ void PairLJCutTIP4PLongOpt::eval()
f[iH2][2] += fHz;
if (VFLAG) {
- domain->closest_image(x[i],x[iH1],xH1);
- domain->closest_image(x[i],x[iH2],xH2);
-
+ xH1 = x[iH1];
+ xH2 = x[iH2];
v[0] = x[i][0]*fOx + xH1[0]*fHx + xH2[0]*fHx;
v[1] = x[i][1]*fOy + xH1[1]*fHy + xH2[1]*fHy;
v[2] = x[i][2]*fOz + xH1[2]*fHz + xH2[2]*fHz;
@@ -399,9 +412,8 @@ void PairLJCutTIP4PLongOpt::eval()
f[jH2][2] += fHz;
if (VFLAG) {
- domain->closest_image(x[j],x[jH1],xH1);
- domain->closest_image(x[j],x[jH2],xH2);
-
+ xH1 = x[jH1];
+ xH2 = x[jH2];
v[0] += x[j][0]*fOx + xH1[0]*fHx + xH2[0]*fHx;
v[1] += x[j][1]*fOy + xH1[1]*fHy + xH2[1]*fHy;
v[2] += x[j][2]*fOz + xH1[2]*fHz + xH2[2]*fHz;
@@ -448,12 +460,10 @@ void PairLJCutTIP4PLongOpt::compute_newsite_opt(const double * xO,
double delx1 = xH1[0] - xO[0];
double dely1 = xH1[1] - xO[1];
double delz1 = xH1[2] - xO[2];
- domain->minimum_image(delx1,dely1,delz1);
double delx2 = xH2[0] - xO[0];
double dely2 = xH2[1] - xO[1];
double delz2 = xH2[2] - xO[2];
- domain->minimum_image(delx2,dely2,delz2);
const double prefac = alpha * 0.5;
xM[0] = xO[0] + prefac * (delx1 + delx2);
diff --git a/src/QEQ/fix_qeq.cpp b/src/QEQ/fix_qeq.cpp
index 01d51e12ed..7e8db7632c 100644
--- a/src/QEQ/fix_qeq.cpp
+++ b/src/QEQ/fix_qeq.cpp
@@ -310,13 +310,6 @@ void FixQEq::setup_pre_force_respa(int vflag, int ilevel)
/* ---------------------------------------------------------------------- */
-void FixQEq::min_setup_pre_force(int vflag)
-{
- setup_pre_force(vflag);
-}
-
-/* ---------------------------------------------------------------------- */
-
void FixQEq::init_storage()
{
nlocal = atom->nlocal;
diff --git a/src/QEQ/fix_qeq.h b/src/QEQ/fix_qeq.h
index 152ffd4e32..ed87412284 100644
--- a/src/QEQ/fix_qeq.h
+++ b/src/QEQ/fix_qeq.h
@@ -33,7 +33,6 @@ class FixQEq : public Fix {
void setup_pre_force(int);
void setup_pre_force_respa(int, int);
void pre_force_respa(int, int, int);
- void min_setup_pre_force(int);
void min_pre_force(int);
// derived child classes must provide these functions
diff --git a/src/REPLICA/compute_event_displace.cpp b/src/REPLICA/compute_event_displace.cpp
index cdd0e5de40..1431fc202e 100644
--- a/src/REPLICA/compute_event_displace.cpp
+++ b/src/REPLICA/compute_event_displace.cpp
@@ -48,7 +48,7 @@ ComputeEventDisplace::ComputeEventDisplace(LAMMPS *lmp, int narg, char **arg) :
error->all(FLERR,"Distance must be > 0 for compute event/displace");
displace_distsq = displace_dist * displace_dist;
- // fix event ID will be set later by PRD
+ // fix event ID will be set later by accelerated dynamics method
id_event = NULL;
}
@@ -75,7 +75,8 @@ void ComputeEventDisplace::init()
fix_event = (FixEvent*) modify->fix[ifix];
if (strcmp(fix_event->style,"EVENT/PRD") != 0 &&
- strcmp(fix_event->style,"EVENT/TAD") != 0)
+ strcmp(fix_event->style,"EVENT/TAD") != 0 &&
+ strcmp(fix_event->style,"EVENT/HYPER") != 0)
error->all(FLERR,"Compute event/displace has invalid fix event assigned");
}
diff --git a/src/REPLICA/fix_neb.cpp b/src/REPLICA/fix_neb.cpp
index 939e4d7fac..249339191a 100644
--- a/src/REPLICA/fix_neb.cpp
+++ b/src/REPLICA/fix_neb.cpp
@@ -395,7 +395,7 @@ void FixNEB::inter_replica_comm()
// -----------------------------------------------------
// single proc per replica
- // all atoms are NEB atoms and no atom sorting
+ // all atoms are NEB atoms and no atom sorting is enabled
// direct comm of x -> xprev and x -> xnext
if (cmode == SINGLE_PROC_DIRECT) {
diff --git a/src/USER-ATC/fix_atc.cpp b/src/USER-ATC/fix_atc.cpp
index 0dc10e1773..96ad93481d 100644
--- a/src/USER-ATC/fix_atc.cpp
+++ b/src/USER-ATC/fix_atc.cpp
@@ -632,16 +632,6 @@ void FixATC::min_pre_exchange()
throw;
}
}
-void FixATC::min_setup_pre_exchange()
-{
- try {
- atc_->setup_pre_exchange();
- }
- catch (ATC::ATC_Error& atcError) {
- ATC::LammpsInterface::instance()->print_msg(atcError.error_description());
- throw;
- }
-}
double FixATC::memory_usage()
{
@@ -896,19 +886,6 @@ void FixATC::setup_pre_neighbor()
}
}
/* ---------------------------------------------------------------------- */
-void FixATC::min_setup_pre_neighbor()
-{
- if (atc_->is_initialized()) {
- try {
- atc_->pre_neighbor();
- }
- catch (ATC::ATC_Error& atcError) {
- ATC::LammpsInterface::instance()->print_msg(atcError.error_description());
- throw;
- }
- }
-}
-/* ---------------------------------------------------------------------- */
void FixATC::min_pre_force(int vflag)
{
try {
diff --git a/src/USER-ATC/fix_atc.h b/src/USER-ATC/fix_atc.h
index b87c585e73..2512c658db 100644
--- a/src/USER-ATC/fix_atc.h
+++ b/src/USER-ATC/fix_atc.h
@@ -52,7 +52,6 @@ namespace LAMMPS_NS {
and is called before domain->pbc() and comm->exchange(). */
void setup_pre_exchange();
void pre_exchange();
- void min_setup_pre_exchange();
void min_pre_exchange();
double memory_usage();
@@ -84,7 +83,6 @@ namespace LAMMPS_NS {
neighbor->build(). */
void pre_neighbor();
void setup_pre_neighbor();
- void min_setup_pre_neighbor();
/** pre/post_force is used to modify fix-specific data
and is before/after the various force computations. */
diff --git a/src/USER-DPD/fix_shardlow.cpp b/src/USER-DPD/fix_shardlow.cpp
index 28c5382237..1b3d9f1c34 100644
--- a/src/USER-DPD/fix_shardlow.cpp
+++ b/src/USER-DPD/fix_shardlow.cpp
@@ -178,13 +178,6 @@ void FixShardlow::min_pre_exchange()
/* ---------------------------------------------------------------------- */
-void FixShardlow::min_setup_pre_exchange()
-{
- memset(atom->ssaAIR, 0, sizeof(int)*atom->nlocal);
-}
-
-/* ---------------------------------------------------------------------- */
-
void FixShardlow::setup(int vflag)
{
bool fixShardlow = false;
diff --git a/src/USER-DPD/fix_shardlow.h b/src/USER-DPD/fix_shardlow.h
index 45a7b030b9..139a6de02c 100644
--- a/src/USER-DPD/fix_shardlow.h
+++ b/src/USER-DPD/fix_shardlow.h
@@ -37,7 +37,6 @@ class FixShardlow : public Fix {
virtual void initial_integrate(int);
void setup_pre_exchange();
void pre_exchange();
- void min_setup_pre_exchange();
void min_pre_exchange();
void grow_arrays(int);
diff --git a/src/USER-FEP/pair_lj_cut_tip4p_long_soft.cpp b/src/USER-FEP/pair_lj_cut_tip4p_long_soft.cpp
index 2edad28bf3..5beed08b72 100644
--- a/src/USER-FEP/pair_lj_cut_tip4p_long_soft.cpp
+++ b/src/USER-FEP/pair_lj_cut_tip4p_long_soft.cpp
@@ -86,8 +86,8 @@ void PairLJCutTIP4PLongSoft::compute(int eflag, int vflag)
double factor_coul,factor_lj;
double grij,expm2,prefactor,t,erfc;
double denc, denlj, r4sig6;
- double fO[3],fH[3],fd[3],v[6],xH1[3],xH2[3];
- double *x1,*x2;
+ double fO[3],fH[3],fd[3],v[6];
+ double *x1,*x2,*xH1,*xH2;
int *ilist,*jlist,*numneigh,**firstneigh;
double rsq;
@@ -143,14 +143,20 @@ void PairLJCutTIP4PLongSoft::compute(int eflag, int vflag)
if (itype == typeO) {
if (hneigh[i][0] < 0) {
- hneigh[i][0] = iH1 = atom->map(atom->tag[i] + 1);
- hneigh[i][1] = iH2 = atom->map(atom->tag[i] + 2);
- hneigh[i][2] = 1;
+ iH1 = atom->map(atom->tag[i] + 1);
+ iH2 = atom->map(atom->tag[i] + 2);
if (iH1 == -1 || iH2 == -1)
error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+ // set iH1,iH2 to closest image to O
+ iH1 = domain->closest_image(i,iH1);
+ iH2 = domain->closest_image(i,iH2);
compute_newsite(x[i],x[iH1],x[iH2],newsite[i]);
+ hneigh[i][0] = iH1;
+ hneigh[i][1] = iH2;
+ hneigh[i][2] = 1;
+
} else {
iH1 = hneigh[i][0];
iH2 = hneigh[i][1];
@@ -216,14 +222,20 @@ void PairLJCutTIP4PLongSoft::compute(int eflag, int vflag)
if (jtype == typeO) {
if (hneigh[j][0] < 0) {
- hneigh[j][0] = jH1 = atom->map(atom->tag[j] + 1);
- hneigh[j][1] = jH2 = atom->map(atom->tag[j] + 2);
- hneigh[j][2] = 1;
+ jH1 = atom->map(atom->tag[j] + 1);
+ jH2 = atom->map(atom->tag[j] + 2);
if (jH1 == -1 || jH2 == -1)
error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[jH1] != typeH || atom->type[jH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+ // set jH1,jH2 to closest image to O
+ jH1 = domain->closest_image(j,jH1);
+ jH2 = domain->closest_image(j,jH2);
compute_newsite(x[j],x[jH1],x[jH2],newsite[j]);
+ hneigh[j][0] = jH1;
+ hneigh[j][1] = jH2;
+ hneigh[j][2] = 1;
+
} else {
jH1 = hneigh[j][0];
jH2 = hneigh[j][1];
@@ -314,9 +326,8 @@ void PairLJCutTIP4PLongSoft::compute(int eflag, int vflag)
f[iH2][2] += fH[2];
if (vflag) {
- domain->closest_image(x[i],x[iH1],xH1);
- domain->closest_image(x[i],x[iH2],xH2);
-
+ xH1 = x[iH1];
+ xH2 = x[iH2];
v[0] = x[i][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0];
v[1] = x[i][1]*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1];
v[2] = x[i][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2];
@@ -372,9 +383,8 @@ void PairLJCutTIP4PLongSoft::compute(int eflag, int vflag)
f[jH2][2] += fH[2];
if (vflag) {
- domain->closest_image(x[j],x[jH1],xH1);
- domain->closest_image(x[j],x[jH2],xH2);
-
+ xH1 = x[jH1];
+ xH2 = x[jH2];
v[0] += x[j][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0];
v[1] += x[j][1]*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1];
v[2] += x[j][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2];
@@ -549,12 +559,10 @@ void PairLJCutTIP4PLongSoft::compute_newsite(double *xO, double *xH1,
double delx1 = xH1[0] - xO[0];
double dely1 = xH1[1] - xO[1];
double delz1 = xH1[2] - xO[2];
- domain->minimum_image(delx1,dely1,delz1);
double delx2 = xH2[0] - xO[0];
double dely2 = xH2[1] - xO[1];
double delz2 = xH2[2] - xO[2];
- domain->minimum_image(delx2,dely2,delz2);
xM[0] = xO[0] + alpha * 0.5 * (delx1 + delx2);
xM[1] = xO[1] + alpha * 0.5 * (dely1 + dely2);
diff --git a/src/USER-FEP/pair_tip4p_long_soft.cpp b/src/USER-FEP/pair_tip4p_long_soft.cpp
index 027882cfba..7ad654d95c 100644
--- a/src/USER-FEP/pair_tip4p_long_soft.cpp
+++ b/src/USER-FEP/pair_tip4p_long_soft.cpp
@@ -85,8 +85,8 @@ void PairTIP4PLongSoft::compute(int eflag, int vflag)
double factor_coul;
double grij,expm2,prefactor,t,erfc;
double denc;
- double fO[3],fH[3],fd[3],v[6],xH1[3],xH2[3];
- double *x1,*x2;
+ double fO[3],fH[3],fd[3],v[6];
+ double *x1,*x2,*xH1,*xH2;
int *ilist,*jlist,*numneigh,**firstneigh;
double rsq;
@@ -140,14 +140,20 @@ void PairTIP4PLongSoft::compute(int eflag, int vflag)
if (itype == typeO) {
if (hneigh[i][0] < 0) {
- hneigh[i][0] = iH1 = atom->map(atom->tag[i] + 1);
- hneigh[i][1] = iH2 = atom->map(atom->tag[i] + 2);
- hneigh[i][2] = 1;
+ iH1 = atom->map(atom->tag[i] + 1);
+ iH2 = atom->map(atom->tag[i] + 2);
if (iH1 == -1 || iH2 == -1)
error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+ // set iH1,iH2 to closest image to O
+ iH1 = domain->closest_image(i,iH1);
+ iH2 = domain->closest_image(i,iH2);
compute_newsite(x[i],x[iH1],x[iH2],newsite[i]);
+ hneigh[i][0] = iH1;
+ hneigh[i][1] = iH2;
+ hneigh[i][2] = 1;
+
} else {
iH1 = hneigh[i][0];
iH2 = hneigh[i][1];
@@ -184,14 +190,20 @@ void PairTIP4PLongSoft::compute(int eflag, int vflag)
if (jtype == typeO) {
if (hneigh[j][0] < 0) {
- hneigh[j][0] = jH1 = atom->map(atom->tag[j] + 1);
- hneigh[j][1] = jH2 = atom->map(atom->tag[j] + 2);
- hneigh[j][2] = 1;
+ jH1 = atom->map(atom->tag[j] + 1);
+ jH2 = atom->map(atom->tag[j] + 2);
if (jH1 == -1 || jH2 == -1)
error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[jH1] != typeH || atom->type[jH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+ // set jH1,jH2 to closest image to O
+ jH1 = domain->closest_image(j,jH1);
+ jH2 = domain->closest_image(j,jH2);
compute_newsite(x[j],x[jH1],x[jH2],newsite[j]);
+ hneigh[j][0] = jH1;
+ hneigh[j][1] = jH2;
+ hneigh[j][2] = 1;
+
} else {
jH1 = hneigh[j][0];
jH2 = hneigh[j][1];
@@ -282,9 +294,8 @@ void PairTIP4PLongSoft::compute(int eflag, int vflag)
f[iH2][2] += fH[2];
if (vflag) {
- domain->closest_image(x[i],x[iH1],xH1);
- domain->closest_image(x[i],x[iH2],xH2);
-
+ xH1 = x[iH1];
+ xH2 = x[iH2];
v[0] = x[i][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0];
v[1] = x[i][1]*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1];
v[2] = x[i][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2];
@@ -340,9 +351,8 @@ void PairTIP4PLongSoft::compute(int eflag, int vflag)
f[jH2][2] += fH[2];
if (vflag) {
- domain->closest_image(x[j],x[jH1],xH1);
- domain->closest_image(x[j],x[jH2],xH2);
-
+ xH1 = x[jH1];
+ xH2 = x[jH2];
v[0] += x[j][0]*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0];
v[1] += x[j][1]*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1];
v[2] += x[j][2]*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2];
@@ -473,12 +483,10 @@ void PairTIP4PLongSoft::compute_newsite(double *xO, double *xH1,
double delx1 = xH1[0] - xO[0];
double dely1 = xH1[1] - xO[1];
double delz1 = xH1[2] - xO[2];
- domain->minimum_image(delx1,dely1,delz1);
double delx2 = xH2[0] - xO[0];
double dely2 = xH2[1] - xO[1];
double delz2 = xH2[2] - xO[2];
- domain->minimum_image(delx2,dely2,delz2);
xM[0] = xO[0] + alpha * 0.5 * (delx1 + delx2);
xM[1] = xO[1] + alpha * 0.5 * (dely1 + dely2);
diff --git a/src/USER-INTEL/fix_intel.cpp b/src/USER-INTEL/fix_intel.cpp
index 7c5dd0b34c..edd33eb72b 100644
--- a/src/USER-INTEL/fix_intel.cpp
+++ b/src/USER-INTEL/fix_intel.cpp
@@ -356,6 +356,13 @@ void FixIntel::setup(int vflag)
/* ---------------------------------------------------------------------- */
+void FixIntel::setup_pre_reverse(int eflag, int vflag)
+{
+ pre_reverse(eflag,vflag);
+}
+
+/* ---------------------------------------------------------------------- */
+
void FixIntel::pair_init_check(const bool cdmessage)
{
#ifdef INTEL_VMASK
diff --git a/src/USER-INTEL/fix_intel.h b/src/USER-INTEL/fix_intel.h
index 6ffcb4ccab..f4c02b37b5 100644
--- a/src/USER-INTEL/fix_intel.h
+++ b/src/USER-INTEL/fix_intel.h
@@ -43,6 +43,8 @@ class FixIntel : public Fix {
virtual int setmask();
virtual void init();
virtual void setup(int);
+ void setup_pre_reverse(int eflag = 0, int vflag = 0);
+
void pair_init_check(const bool cdmessage=false);
void bond_init_check();
void kspace_init_check();
diff --git a/src/USER-OMP/pair_lj_cut_tip4p_cut_omp.cpp b/src/USER-OMP/pair_lj_cut_tip4p_cut_omp.cpp
index 3b028fb67b..3b4d616149 100644
--- a/src/USER-OMP/pair_lj_cut_tip4p_cut_omp.cpp
+++ b/src/USER-OMP/pair_lj_cut_tip4p_cut_omp.cpp
@@ -127,9 +127,9 @@ void PairLJCutTIP4PCutOMP::eval(int iifrom, int iito, ThrData * const thr)
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul;
double r,rsq,r2inv,r6inv,forcecoul,forcelj,cforce;
double factor_coul,factor_lj;
- double v[6],xH1[3],xH2[3];
+ double v[6];
double fdx,fdy,fdz,fOx,fOy,fOz,fHx,fHy,fHz;
- dbl3_t x1,x2;
+ dbl3_t x1,x2,xH1,xH2;
int *ilist,*jlist,*numneigh,**firstneigh;
int i,j,ii,jj,jnum,itype,jtype,key;
@@ -179,10 +179,14 @@ void PairLJCutTIP4PCutOMP::eval(int iifrom, int iito, ThrData * const thr)
error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+ // set iH1,iH2 to index of closest image to O
+ iH1 = domain->closest_image(i,iH1);
+ iH2 = domain->closest_image(i,iH2);
compute_newsite_thr(x[i],x[iH1],x[iH2],newsite_thr[i]);
hneigh_thr[i].t = 1;
hneigh_thr[i].b = iH2;
hneigh_thr[i].a = iH1;
+
} else {
iH1 = hneigh_thr[i].a;
iH2 = hneigh_thr[i].b;
@@ -256,6 +260,9 @@ void PairLJCutTIP4PCutOMP::eval(int iifrom, int iito, ThrData * const thr)
error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[jH1] != typeH || atom->type[jH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+ // set jH1,jH2 to closest image to O
+ jH1 = domain->closest_image(j,jH1);
+ jH2 = domain->closest_image(j,jH2);
compute_newsite_thr(x[j],x[jH1],x[jH2],newsite_thr[j]);
hneigh_thr[j].t = 1;
hneigh_thr[j].b = jH2;
@@ -341,15 +348,14 @@ void PairLJCutTIP4PCutOMP::eval(int iifrom, int iito, ThrData * const thr)
f[iH2].z += fHz;
if (VFLAG) {
- domain->closest_image(&x[i].x,&x[iH1].x,xH1);
- domain->closest_image(&x[i].x,&x[iH2].x,xH2);
-
- v[0] = x[i].x*fOx + xH1[0]*fHx + xH2[0]*fHx;
- v[1] = x[i].y*fOy + xH1[1]*fHy + xH2[1]*fHy;
- v[2] = x[i].z*fOz + xH1[2]*fHz + xH2[2]*fHz;
- v[3] = x[i].x*fOy + xH1[0]*fHy + xH2[0]*fHy;
- v[4] = x[i].x*fOz + xH1[0]*fHz + xH2[0]*fHz;
- v[5] = x[i].y*fOz + xH1[1]*fHz + xH2[1]*fHz;
+ xH1 = x[iH1];
+ xH2 = x[iH2];
+ v[0] = x[i].x*fOx + xH1.x*fHx + xH2.x*fHx;
+ v[1] = x[i].y*fOy + xH1.y*fHy + xH2.y*fHy;
+ v[2] = x[i].z*fOz + xH1.z*fHz + xH2.z*fHz;
+ v[3] = x[i].x*fOy + xH1.x*fHy + xH2.x*fHy;
+ v[4] = x[i].x*fOz + xH1.x*fHz + xH2.x*fHz;
+ v[5] = x[i].y*fOz + xH1.y*fHz + xH2.y*fHz;
}
if (EVFLAG) {
vlist[n++] = i;
@@ -401,15 +407,14 @@ void PairLJCutTIP4PCutOMP::eval(int iifrom, int iito, ThrData * const thr)
f[jH2].z += fHz;
if (VFLAG) {
- domain->closest_image(&x[j].x,&x[jH1].x,xH1);
- domain->closest_image(&x[j].x,&x[jH2].x,xH2);
-
- v[0] += x[j].x*fOx + xH1[0]*fHx + xH2[0]*fHx;
- v[1] += x[j].y*fOy + xH1[1]*fHy + xH2[1]*fHy;
- v[2] += x[j].z*fOz + xH1[2]*fHz + xH2[2]*fHz;
- v[3] += x[j].x*fOy + xH1[0]*fHy + xH2[0]*fHy;
- v[4] += x[j].x*fOz + xH1[0]*fHz + xH2[0]*fHz;
- v[5] += x[j].y*fOz + xH1[1]*fHz + xH2[1]*fHz;
+ xH1 = x[jH1];
+ xH2 = x[jH2];
+ v[0] += x[j].x*fOx + xH1.x*fHx + xH2.x*fHx;
+ v[1] += x[j].y*fOy + xH1.y*fHy + xH2.y*fHy;
+ v[2] += x[j].z*fOz + xH1.z*fHz + xH2.z*fHz;
+ v[3] += x[j].x*fOy + xH1.x*fHy + xH2.x*fHy;
+ v[4] += x[j].x*fOz + xH1.x*fHz + xH2.x*fHz;
+ v[5] += x[j].y*fOz + xH1.y*fHz + xH2.y*fHz;
}
if (EVFLAG) {
vlist[n++] = j;
@@ -446,12 +451,10 @@ void PairLJCutTIP4PCutOMP::compute_newsite_thr(const dbl3_t &xO,
double delx1 = xH1.x - xO.x;
double dely1 = xH1.y - xO.y;
double delz1 = xH1.z - xO.z;
- domain->minimum_image(delx1,dely1,delz1);
double delx2 = xH2.x - xO.x;
double dely2 = xH2.y - xO.y;
double delz2 = xH2.z - xO.z;
- domain->minimum_image(delx2,dely2,delz2);
const double prefac = alpha * 0.5;
xM.x = xO.x + prefac * (delx1 + delx2);
diff --git a/src/USER-OMP/pair_lj_cut_tip4p_long_omp.cpp b/src/USER-OMP/pair_lj_cut_tip4p_long_omp.cpp
index 70be3e5228..2f1b98fc4f 100644
--- a/src/USER-OMP/pair_lj_cut_tip4p_long_omp.cpp
+++ b/src/USER-OMP/pair_lj_cut_tip4p_long_omp.cpp
@@ -141,9 +141,9 @@ void PairLJCutTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
double r,rsq,r2inv,r6inv,forcecoul,forcelj,cforce;
double factor_coul,factor_lj;
double grij,expm2,prefactor,t,erfc;
- double v[6],xH1[3],xH2[3];
+ double v[6];
double fdx,fdy,fdz,fOx,fOy,fOz,fHx,fHy,fHz;
- dbl3_t x1,x2;
+ dbl3_t x1,x2,xH1,xH2;
int *ilist,*jlist,*numneigh,**firstneigh;
int i,j,ii,jj,jnum,itype,jtype,itable, key;
@@ -193,6 +193,9 @@ void PairLJCutTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+ // set iH1,iH2 to index of closest image to O
+ iH1 = domain->closest_image(i,iH1);
+ iH2 = domain->closest_image(i,iH2);
compute_newsite_thr(x[i],x[iH1],x[iH2],newsite_thr[i]);
hneigh_thr[i].t = 1;
hneigh_thr[i].b = iH2;
@@ -270,6 +273,9 @@ void PairLJCutTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[jH1] != typeH || atom->type[jH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+ // set jH1,jH2 to closest image to O
+ jH1 = domain->closest_image(j,jH1);
+ jH2 = domain->closest_image(j,jH2);
compute_newsite_thr(x[j],x[jH1],x[jH2],newsite_thr[j]);
hneigh_thr[j].t = 1;
hneigh_thr[j].b = jH2;
@@ -379,15 +385,14 @@ void PairLJCutTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
f[iH2].z += fHz;
if (VFLAG) {
- domain->closest_image(&x[i].x,&x[iH1].x,xH1);
- domain->closest_image(&x[i].x,&x[iH2].x,xH2);
-
- v[0] = x[i].x*fOx + xH1[0]*fHx + xH2[0]*fHx;
- v[1] = x[i].y*fOy + xH1[1]*fHy + xH2[1]*fHy;
- v[2] = x[i].z*fOz + xH1[2]*fHz + xH2[2]*fHz;
- v[3] = x[i].x*fOy + xH1[0]*fHy + xH2[0]*fHy;
- v[4] = x[i].x*fOz + xH1[0]*fHz + xH2[0]*fHz;
- v[5] = x[i].y*fOz + xH1[1]*fHz + xH2[1]*fHz;
+ xH1 = x[iH1];
+ xH2 = x[iH2];
+ v[0] = x[i].x*fOx + xH1.x*fHx + xH2.x*fHx;
+ v[1] = x[i].y*fOy + xH1.y*fHy + xH2.y*fHy;
+ v[2] = x[i].z*fOz + xH1.z*fHz + xH2.z*fHz;
+ v[3] = x[i].x*fOy + xH1.x*fHy + xH2.x*fHy;
+ v[4] = x[i].x*fOz + xH1.x*fHz + xH2.x*fHz;
+ v[5] = x[i].y*fOz + xH1.y*fHz + xH2.y*fHz;
}
if (EVFLAG) {
vlist[n++] = i;
@@ -439,15 +444,14 @@ void PairLJCutTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
f[jH2].z += fHz;
if (VFLAG) {
- domain->closest_image(&x[j].x,&x[jH1].x,xH1);
- domain->closest_image(&x[j].x,&x[jH2].x,xH2);
-
- v[0] += x[j].x*fOx + xH1[0]*fHx + xH2[0]*fHx;
- v[1] += x[j].y*fOy + xH1[1]*fHy + xH2[1]*fHy;
- v[2] += x[j].z*fOz + xH1[2]*fHz + xH2[2]*fHz;
- v[3] += x[j].x*fOy + xH1[0]*fHy + xH2[0]*fHy;
- v[4] += x[j].x*fOz + xH1[0]*fHz + xH2[0]*fHz;
- v[5] += x[j].y*fOz + xH1[1]*fHz + xH2[1]*fHz;
+ xH1 = x[jH1];
+ xH2 = x[jH2];
+ v[0] += x[j].x*fOx + xH1.x*fHx + xH2.x*fHx;
+ v[1] += x[j].y*fOy + xH1.y*fHy + xH2.y*fHy;
+ v[2] += x[j].z*fOz + xH1.z*fHz + xH2.z*fHz;
+ v[3] += x[j].x*fOy + xH1.x*fHy + xH2.x*fHy;
+ v[4] += x[j].x*fOz + xH1.x*fHz + xH2.x*fHz;
+ v[5] += x[j].y*fOz + xH1.y*fHz + xH2.y*fHz;
}
if (EVFLAG) {
vlist[n++] = j;
@@ -489,12 +493,10 @@ void PairLJCutTIP4PLongOMP::compute_newsite_thr(const dbl3_t &xO,
double delx1 = xH1.x - xO.x;
double dely1 = xH1.y - xO.y;
double delz1 = xH1.z - xO.z;
- domain->minimum_image(delx1,dely1,delz1);
double delx2 = xH2.x - xO.x;
double dely2 = xH2.y - xO.y;
double delz2 = xH2.z - xO.z;
- domain->minimum_image(delx2,dely2,delz2);
const double prefac = alpha * 0.5;
xM.x = xO.x + prefac * (delx1 + delx2);
diff --git a/src/USER-OMP/pair_lj_cut_tip4p_long_soft_omp.cpp b/src/USER-OMP/pair_lj_cut_tip4p_long_soft_omp.cpp
index 2d4be0d884..b31c1c109b 100644
--- a/src/USER-OMP/pair_lj_cut_tip4p_long_soft_omp.cpp
+++ b/src/USER-OMP/pair_lj_cut_tip4p_long_soft_omp.cpp
@@ -129,9 +129,9 @@ void PairLJCutTIP4PLongSoftOMP::eval(int iifrom, int iito, ThrData * const thr)
double factor_coul,factor_lj;
double grij,expm2,prefactor,t,erfc;
double denc, denlj, r4sig6;
- double v[6],xH1[3],xH2[3];
+ double v[6];
double fdx,fdy,fdz,fOx,fOy,fOz,fHx,fHy,fHz;
- dbl3_t x1,x2;
+ dbl3_t x1,x2,xH1,xH2;
int *ilist,*jlist,*numneigh,**firstneigh;
int i,j,ii,jj,jnum,itype,jtype,key;
@@ -181,6 +181,9 @@ void PairLJCutTIP4PLongSoftOMP::eval(int iifrom, int iito, ThrData * const thr)
error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+ // set iH1,iH2 to index of closest image to O
+ iH1 = domain->closest_image(i,iH1);
+ iH2 = domain->closest_image(i,iH2);
compute_newsite_thr(x[i],x[iH1],x[iH2],newsite_thr[i]);
hneigh_thr[i].t = 1;
hneigh_thr[i].b = iH2;
@@ -261,6 +264,9 @@ void PairLJCutTIP4PLongSoftOMP::eval(int iifrom, int iito, ThrData * const thr)
error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[jH1] != typeH || atom->type[jH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+ // set jH1,jH2 to closest image to O
+ jH1 = domain->closest_image(j,jH1);
+ jH2 = domain->closest_image(j,jH2);
compute_newsite_thr(x[j],x[jH1],x[jH2],newsite_thr[j]);
hneigh_thr[j].t = 1;
hneigh_thr[j].b = jH2;
@@ -357,15 +363,14 @@ void PairLJCutTIP4PLongSoftOMP::eval(int iifrom, int iito, ThrData * const thr)
f[iH2].z += fHz;
if (VFLAG) {
- domain->closest_image(&x[i].x,&x[iH1].x,xH1);
- domain->closest_image(&x[i].x,&x[iH2].x,xH2);
-
- v[0] = x[i].x*fOx + xH1[0]*fHx + xH2[0]*fHx;
- v[1] = x[i].y*fOy + xH1[1]*fHy + xH2[1]*fHy;
- v[2] = x[i].z*fOz + xH1[2]*fHz + xH2[2]*fHz;
- v[3] = x[i].x*fOy + xH1[0]*fHy + xH2[0]*fHy;
- v[4] = x[i].x*fOz + xH1[0]*fHz + xH2[0]*fHz;
- v[5] = x[i].y*fOz + xH1[1]*fHz + xH2[1]*fHz;
+ xH1 = x[iH1];
+ xH2 = x[iH2];
+ v[0] = x[i].x*fOx + xH1.x*fHx + xH2.x*fHx;
+ v[1] = x[i].y*fOy + xH1.y*fHy + xH2.y*fHy;
+ v[2] = x[i].z*fOz + xH1.z*fHz + xH2.z*fHz;
+ v[3] = x[i].x*fOy + xH1.x*fHy + xH2.x*fHy;
+ v[4] = x[i].x*fOz + xH1.x*fHz + xH2.x*fHz;
+ v[5] = x[i].y*fOz + xH1.y*fHz + xH2.y*fHz;
}
if (EVFLAG) {
vlist[n++] = i;
@@ -417,15 +422,14 @@ void PairLJCutTIP4PLongSoftOMP::eval(int iifrom, int iito, ThrData * const thr)
f[jH2].z += fHz;
if (VFLAG) {
- domain->closest_image(&x[j].x,&x[jH1].x,xH1);
- domain->closest_image(&x[j].x,&x[jH2].x,xH2);
-
- v[0] += x[j].x*fOx + xH1[0]*fHx + xH2[0]*fHx;
- v[1] += x[j].y*fOy + xH1[1]*fHy + xH2[1]*fHy;
- v[2] += x[j].z*fOz + xH1[2]*fHz + xH2[2]*fHz;
- v[3] += x[j].x*fOy + xH1[0]*fHy + xH2[0]*fHy;
- v[4] += x[j].x*fOz + xH1[0]*fHz + xH2[0]*fHz;
- v[5] += x[j].y*fOz + xH1[1]*fHz + xH2[1]*fHz;
+ xH1 = x[jH1];
+ xH2 = x[jH2];
+ v[0] += x[j].x*fOx + xH1.x*fHx + xH2.x*fHx;
+ v[1] += x[j].y*fOy + xH1.y*fHy + xH2.y*fHy;
+ v[2] += x[j].z*fOz + xH1.z*fHz + xH2.z*fHz;
+ v[3] += x[j].x*fOy + xH1.x*fHy + xH2.x*fHy;
+ v[4] += x[j].x*fOz + xH1.x*fHz + xH2.x*fHz;
+ v[5] += x[j].y*fOz + xH1.y*fHz + xH2.y*fHz;
}
if (EVFLAG) {
vlist[n++] = j;
@@ -463,12 +467,10 @@ void PairLJCutTIP4PLongSoftOMP::compute_newsite_thr(const dbl3_t &xO,
double delx1 = xH1.x - xO.x;
double dely1 = xH1.y - xO.y;
double delz1 = xH1.z - xO.z;
- domain->minimum_image(delx1,dely1,delz1);
double delx2 = xH2.x - xO.x;
double dely2 = xH2.y - xO.y;
double delz2 = xH2.z - xO.z;
- domain->minimum_image(delx2,dely2,delz2);
const double prefac = alpha * 0.5;
xM.x = xO.x + prefac * (delx1 + delx2);
diff --git a/src/USER-OMP/pair_lj_long_tip4p_long_omp.cpp b/src/USER-OMP/pair_lj_long_tip4p_long_omp.cpp
index f04e6380c8..5072496cc6 100644
--- a/src/USER-OMP/pair_lj_long_tip4p_long_omp.cpp
+++ b/src/USER-OMP/pair_lj_long_tip4p_long_omp.cpp
@@ -737,8 +737,8 @@ void PairLJLongTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
double r,r2inv,forcecoul,forcelj,cforce;
double factor_coul;
double grij,expm2,prefactor,t,erfc;
- double fO[3],fH[3],fd[3],v[6],xH1[3],xH2[3];
- dbl3_t x1,x2;
+ double fO[3],fH[3],fd[3],v[6];
+ dbl3_t x1,x2,xH1,xH2;
int *ilist,*jlist,*numneigh,**firstneigh;
double rsq;
@@ -763,14 +763,20 @@ void PairLJLongTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
itype = type[i];
if (itype == typeO) {
if (hneigh_thr[i].a < 0) {
- hneigh_thr[i].a = iH1 = atom->map(atom->tag[i] + 1);
- hneigh_thr[i].b = iH2 = atom->map(atom->tag[i] + 2);
- hneigh_thr[i].t = 1;
+ iH1 = atom->map(atom->tag[i] + 1);
+ iH2 = atom->map(atom->tag[i] + 2);
if (iH1 == -1 || iH2 == -1)
error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+ // set iH1,iH2 to index of closest image to O
+ iH1 = domain->closest_image(i,iH1);
+ iH2 = domain->closest_image(i,iH2);
compute_newsite_thr(x[i],x[iH1],x[iH2],newsite_thr[i]);
+ hneigh_thr[i].a = iH1;
+ hneigh_thr[i].b = iH2;
+ hneigh_thr[i].t = 1;
+
} else {
iH1 = hneigh_thr[i].a;
iH2 = hneigh_thr[i].b;
@@ -871,14 +877,20 @@ void PairLJLongTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
if (itype == typeO || jtype == typeO) {
if (jtype == typeO) {
if (hneigh_thr[j].a < 0) {
- hneigh_thr[j].a = jH1 = atom->map(atom->tag[j] + 1);
- hneigh_thr[j].b = jH2 = atom->map(atom->tag[j] + 2);
- hneigh_thr[j].t = 1;
+ jH1 = atom->map(atom->tag[j] + 1);
+ jH2 = atom->map(atom->tag[j] + 2);
if (jH1 == -1 || jH2 == -1)
error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[jH1] != typeH || atom->type[jH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+ // set jH1,jH2 to closest image to O
+ jH1 = domain->closest_image(j,jH1);
+ jH2 = domain->closest_image(j,jH2);
compute_newsite_thr(x[j],x[jH1],x[jH2],newsite_thr[j]);
+ hneigh_thr[j].a = jH1;
+ hneigh_thr[j].b = jH2;
+ hneigh_thr[j].t = 1;
+
} else {
jH1 = hneigh_thr[j].a;
jH2 = hneigh_thr[j].b;
@@ -983,15 +995,14 @@ void PairLJLongTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
f[iH2][2] += fH[2];
if (EVFLAG) {
- domain->closest_image(&x[i].x,&x[iH1].x,xH1);
- domain->closest_image(&x[i].x,&x[iH2].x,xH2);
-
- v[0] = x[i].x*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0];
- v[1] = x[i].y*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1];
- v[2] = x[i].z*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2];
- v[3] = x[i].x*fO[1] + xH1[0]*fH[1] + xH2[0]*fH[1];
- v[4] = x[i].x*fO[2] + xH1[0]*fH[2] + xH2[0]*fH[2];
- v[5] = x[i].y*fO[2] + xH1[1]*fH[2] + xH2[1]*fH[2];
+ xH1 = x[iH1];
+ xH2 = x[iH2];
+ v[0] = x[i].x*fO[0] + xH1.x*fH[0] + xH2.x*fH[0];
+ v[1] = x[i].y*fO[1] + xH1.y*fH[1] + xH2.y*fH[1];
+ v[2] = x[i].z*fO[2] + xH1.z*fH[2] + xH2.z*fH[2];
+ v[3] = x[i].x*fO[1] + xH1.x*fH[1] + xH2.x*fH[1];
+ v[4] = x[i].x*fO[2] + xH1.x*fH[2] + xH2.x*fH[2];
+ v[5] = x[i].y*fO[2] + xH1.y*fH[2] + xH2.y*fH[2];
}
vlist[n++] = i;
vlist[n++] = iH1;
@@ -1041,15 +1052,14 @@ void PairLJLongTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
f[jH2][2] += fH[2];
if (EVFLAG) {
- domain->closest_image(&x[j].x,&x[jH1].x,xH1);
- domain->closest_image(&x[j].x,&x[jH2].x,xH2);
-
- v[0] += x[j].x*fO[0] + xH1[0]*fH[0] + xH2[0]*fH[0];
- v[1] += x[j].y*fO[1] + xH1[1]*fH[1] + xH2[1]*fH[1];
- v[2] += x[j].z*fO[2] + xH1[2]*fH[2] + xH2[2]*fH[2];
- v[3] += x[j].x*fO[1] + xH1[0]*fH[1] + xH2[0]*fH[1];
- v[4] += x[j].x*fO[2] + xH1[0]*fH[2] + xH2[0]*fH[2];
- v[5] += x[j].y*fO[2] + xH1[1]*fH[2] + xH2[1]*fH[2];
+ xH1 = x[jH1];
+ xH2 = x[jH2];
+ v[0] += x[j].x*fO[0] + xH1.x*fH[0] + xH2.x*fH[0];
+ v[1] += x[j].y*fO[1] + xH1.y*fH[1] + xH2.y*fH[1];
+ v[2] += x[j].z*fO[2] + xH1.z*fH[2] + xH2.z*fH[2];
+ v[3] += x[j].x*fO[1] + xH1.x*fH[1] + xH2.x*fH[1];
+ v[4] += x[j].x*fO[2] + xH1.x*fH[2] + xH2.x*fH[2];
+ v[5] += x[j].y*fO[2] + xH1.y*fH[2] + xH2.y*fH[2];
}
vlist[n++] = j;
vlist[n++] = jH1;
@@ -1981,12 +1991,10 @@ void PairLJLongTIP4PLongOMP::compute_newsite_thr(const dbl3_t &xO,
double delx1 = xH1.x - xO.x;
double dely1 = xH1.y - xO.y;
double delz1 = xH1.z - xO.z;
- domain->minimum_image(delx1,dely1,delz1);
double delx2 = xH2.x - xO.x;
double dely2 = xH2.y - xO.y;
double delz2 = xH2.z - xO.z;
- domain->minimum_image(delx2,dely2,delz2);
const double prefac = alpha * 0.5;
xM.x = xO.x + prefac * (delx1 + delx2);
diff --git a/src/USER-OMP/pair_tip4p_cut_omp.cpp b/src/USER-OMP/pair_tip4p_cut_omp.cpp
index 2b125fde7f..3edd2b688f 100644
--- a/src/USER-OMP/pair_tip4p_cut_omp.cpp
+++ b/src/USER-OMP/pair_tip4p_cut_omp.cpp
@@ -126,9 +126,9 @@ void PairTIP4PCutOMP::eval(int iifrom, int iito, ThrData * const thr)
double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul;
double r,rsq,r2inv,forcecoul,cforce;
double factor_coul;
- double v[6],xH1[3],xH2[3];
+ double v[6];
double fdx,fdy,fdz,fOx,fOy,fOz,fHx,fHy,fHz;
- dbl3_t x1,x2;
+ dbl3_t x1,x2,xH1,xH2;
int *ilist,*jlist,*numneigh,**firstneigh;
int i,j,ii,jj,jnum,itype,jtype,key;
@@ -176,10 +176,14 @@ void PairTIP4PCutOMP::eval(int iifrom, int iito, ThrData * const thr)
error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+ // set iH1,iH2 to index of closest image to O
+ iH1 = domain->closest_image(i,iH1);
+ iH2 = domain->closest_image(i,iH2);
compute_newsite_thr(x[i],x[iH1],x[iH2],newsite_thr[i]);
hneigh_thr[i].t = 1;
hneigh_thr[i].b = iH2;
hneigh_thr[i].a = iH1;
+
} else {
iH1 = hneigh_thr[i].a;
iH2 = hneigh_thr[i].b;
@@ -227,6 +231,9 @@ void PairTIP4PCutOMP::eval(int iifrom, int iito, ThrData * const thr)
error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[jH1] != typeH || atom->type[jH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+ // set jH1,jH2 to closest image to O
+ jH1 = domain->closest_image(j,jH1);
+ jH2 = domain->closest_image(j,jH2);
compute_newsite_thr(x[j],x[jH1],x[jH2],newsite_thr[j]);
hneigh_thr[j].t = 1;
hneigh_thr[j].b = jH2;
@@ -312,15 +319,14 @@ void PairTIP4PCutOMP::eval(int iifrom, int iito, ThrData * const thr)
f[iH2].z += fHz;
if (VFLAG) {
- domain->closest_image(&x[i].x,&x[iH1].x,xH1);
- domain->closest_image(&x[i].x,&x[iH2].x,xH2);
-
- v[0] = x[i].x*fOx + xH1[0]*fHx + xH2[0]*fHx;
- v[1] = x[i].y*fOy + xH1[1]*fHy + xH2[1]*fHy;
- v[2] = x[i].z*fOz + xH1[2]*fHz + xH2[2]*fHz;
- v[3] = x[i].x*fOy + xH1[0]*fHy + xH2[0]*fHy;
- v[4] = x[i].x*fOz + xH1[0]*fHz + xH2[0]*fHz;
- v[5] = x[i].y*fOz + xH1[1]*fHz + xH2[1]*fHz;
+ xH1 = x[iH1];
+ xH2 = x[iH2];
+ v[0] = x[i].x*fOx + xH1.x*fHx + xH2.x*fHx;
+ v[1] = x[i].y*fOy + xH1.y*fHy + xH2.y*fHy;
+ v[2] = x[i].z*fOz + xH1.z*fHz + xH2.z*fHz;
+ v[3] = x[i].x*fOy + xH1.x*fHy + xH2.x*fHy;
+ v[4] = x[i].x*fOz + xH1.x*fHz + xH2.x*fHz;
+ v[5] = x[i].y*fOz + xH1.y*fHz + xH2.y*fHz;
}
if (EVFLAG) {
vlist[n++] = i;
@@ -372,15 +378,14 @@ void PairTIP4PCutOMP::eval(int iifrom, int iito, ThrData * const thr)
f[jH2].z += fHz;
if (VFLAG) {
- domain->closest_image(&x[j].x,&x[jH1].x,xH1);
- domain->closest_image(&x[j].x,&x[jH2].x,xH2);
-
- v[0] += x[j].x*fOx + xH1[0]*fHx + xH2[0]*fHx;
- v[1] += x[j].y*fOy + xH1[1]*fHy + xH2[1]*fHy;
- v[2] += x[j].z*fOz + xH1[2]*fHz + xH2[2]*fHz;
- v[3] += x[j].x*fOy + xH1[0]*fHy + xH2[0]*fHy;
- v[4] += x[j].x*fOz + xH1[0]*fHz + xH2[0]*fHz;
- v[5] += x[j].y*fOz + xH1[1]*fHz + xH2[1]*fHz;
+ xH1 = x[jH1];
+ xH2 = x[jH2];
+ v[0] += x[j].x*fOx + xH1.x*fHx + xH2.x*fHx;
+ v[1] += x[j].y*fOy + xH1.y*fHy + xH2.y*fHy;
+ v[2] += x[j].z*fOz + xH1.z*fHz + xH2.z*fHz;
+ v[3] += x[j].x*fOy + xH1.x*fHy + xH2.x*fHy;
+ v[4] += x[j].x*fOz + xH1.x*fHz + xH2.x*fHz;
+ v[5] += x[j].y*fOz + xH1.y*fHz + xH2.y*fHz;
}
if (EVFLAG) {
vlist[n++] = j;
@@ -417,12 +422,10 @@ void PairTIP4PCutOMP::compute_newsite_thr(const dbl3_t &xO,
double delx1 = xH1.x - xO.x;
double dely1 = xH1.y - xO.y;
double delz1 = xH1.z - xO.z;
- domain->minimum_image(delx1,dely1,delz1);
double delx2 = xH2.x - xO.x;
double dely2 = xH2.y - xO.y;
double delz2 = xH2.z - xO.z;
- domain->minimum_image(delx2,dely2,delz2);
const double prefac = alpha * 0.5;
xM.x = xO.x + prefac * (delx1 + delx2);
diff --git a/src/USER-OMP/pair_tip4p_long_omp.cpp b/src/USER-OMP/pair_tip4p_long_omp.cpp
index 96ebc64a90..bcd8659a54 100644
--- a/src/USER-OMP/pair_tip4p_long_omp.cpp
+++ b/src/USER-OMP/pair_tip4p_long_omp.cpp
@@ -141,9 +141,9 @@ void PairTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
double r,rsq,r2inv,forcecoul,cforce;
double factor_coul;
double grij,expm2,prefactor,t,erfc;
- double v[6],xH1[3],xH2[3];
+ double v[6];
double fdx,fdy,fdz,fOx,fOy,fOz,fHx,fHy,fHz;
- dbl3_t x1,x2;
+ dbl3_t x1,x2,xH1,xH2;
int *ilist,*jlist,*numneigh,**firstneigh;
int i,j,ii,jj,jnum,itype,jtype,itable, key;
@@ -191,10 +191,14 @@ void PairTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+ // set iH1,iH2 to index of closest image to O
+ iH1 = domain->closest_image(i,iH1);
+ iH2 = domain->closest_image(i,iH2);
compute_newsite_thr(x[i],x[iH1],x[iH2],newsite_thr[i]);
hneigh_thr[i].t = 1;
hneigh_thr[i].b = iH2;
hneigh_thr[i].a = iH1;
+
} else {
iH1 = hneigh_thr[i].a;
iH2 = hneigh_thr[i].b;
@@ -242,10 +246,14 @@ void PairTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[jH1] != typeH || atom->type[jH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+ // set jH1,jH2 to closest image to O
+ jH1 = domain->closest_image(j,jH1);
+ jH2 = domain->closest_image(j,jH2);
compute_newsite_thr(x[j],x[jH1],x[jH2],newsite_thr[j]);
hneigh_thr[j].t = 1;
hneigh_thr[j].b = jH2;
hneigh_thr[j].a = jH1;
+
} else {
jH1 = hneigh_thr[j].a;
jH2 = hneigh_thr[j].b;
@@ -351,15 +359,14 @@ void PairTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
f[iH2].z += fHz;
if (VFLAG) {
- domain->closest_image(&x[i].x,&x[iH1].x,xH1);
- domain->closest_image(&x[i].x,&x[iH2].x,xH2);
-
- v[0] = x[i].x*fOx + xH1[0]*fHx + xH2[0]*fHx;
- v[1] = x[i].y*fOy + xH1[1]*fHy + xH2[1]*fHy;
- v[2] = x[i].z*fOz + xH1[2]*fHz + xH2[2]*fHz;
- v[3] = x[i].x*fOy + xH1[0]*fHy + xH2[0]*fHy;
- v[4] = x[i].x*fOz + xH1[0]*fHz + xH2[0]*fHz;
- v[5] = x[i].y*fOz + xH1[1]*fHz + xH2[1]*fHz;
+ xH1 = x[iH1];
+ xH2 = x[iH2];
+ v[0] = x[i].x*fOx + xH1.x*fHx + xH2.x*fHx;
+ v[1] = x[i].y*fOy + xH1.y*fHy + xH2.y*fHy;
+ v[2] = x[i].z*fOz + xH1.z*fHz + xH2.z*fHz;
+ v[3] = x[i].x*fOy + xH1.x*fHy + xH2.x*fHy;
+ v[4] = x[i].x*fOz + xH1.x*fHz + xH2.x*fHz;
+ v[5] = x[i].y*fOz + xH1.y*fHz + xH2.y*fHz;
}
if (EVFLAG) {
vlist[n++] = i;
@@ -411,15 +418,14 @@ void PairTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
f[jH2].z += fHz;
if (VFLAG) {
- domain->closest_image(&x[j].x,&x[jH1].x,xH1);
- domain->closest_image(&x[j].x,&x[jH2].x,xH2);
-
- v[0] += x[j].x*fOx + xH1[0]*fHx + xH2[0]*fHx;
- v[1] += x[j].y*fOy + xH1[1]*fHy + xH2[1]*fHy;
- v[2] += x[j].z*fOz + xH1[2]*fHz + xH2[2]*fHz;
- v[3] += x[j].x*fOy + xH1[0]*fHy + xH2[0]*fHy;
- v[4] += x[j].x*fOz + xH1[0]*fHz + xH2[0]*fHz;
- v[5] += x[j].y*fOz + xH1[1]*fHz + xH2[1]*fHz;
+ xH1 = x[jH1];
+ xH2 = x[jH2];
+ v[0] += x[j].x*fOx + xH1.x*fHx + xH2.x*fHx;
+ v[1] += x[j].y*fOy + xH1.y*fHy + xH2.y*fHy;
+ v[2] += x[j].z*fOz + xH1.z*fHz + xH2.z*fHz;
+ v[3] += x[j].x*fOy + xH1.x*fHy + xH2.x*fHy;
+ v[4] += x[j].x*fOz + xH1.x*fHz + xH2.x*fHz;
+ v[5] += x[j].y*fOz + xH1.y*fHz + xH2.y*fHz;
}
if (EVFLAG) {
vlist[n++] = j;
@@ -461,12 +467,10 @@ void PairTIP4PLongOMP::compute_newsite_thr(const dbl3_t &xO,
double delx1 = xH1.x - xO.x;
double dely1 = xH1.y - xO.y;
double delz1 = xH1.z - xO.z;
- domain->minimum_image(delx1,dely1,delz1);
double delx2 = xH2.x - xO.x;
double dely2 = xH2.y - xO.y;
double delz2 = xH2.z - xO.z;
- domain->minimum_image(delx2,dely2,delz2);
const double prefac = alpha * 0.5;
xM.x = xO.x + prefac * (delx1 + delx2);
diff --git a/src/USER-OMP/pair_tip4p_long_soft_omp.cpp b/src/USER-OMP/pair_tip4p_long_soft_omp.cpp
index ac3bcfc393..2f16a28b8c 100644
--- a/src/USER-OMP/pair_tip4p_long_soft_omp.cpp
+++ b/src/USER-OMP/pair_tip4p_long_soft_omp.cpp
@@ -128,9 +128,9 @@ void PairTIP4PLongSoftOMP::eval(int iifrom, int iito, ThrData * const thr)
double r,rsq,forcecoul,cforce;
double factor_coul,denc;
double grij,expm2,prefactor,t,erfc;
- double v[6],xH1[3],xH2[3];
+ double v[6];
double fdx,fdy,fdz,fOx,fOy,fOz,fHx,fHy,fHz;
- dbl3_t x1,x2;
+ dbl3_t x1,x2,xH1,xH2;
int *ilist,*jlist,*numneigh,**firstneigh;
int i,j,ii,jj,jnum,itype,jtype,key;
@@ -178,10 +178,14 @@ void PairTIP4PLongSoftOMP::eval(int iifrom, int iito, ThrData * const thr)
error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+ // set iH1,iH2 to index of closest image to O
+ iH1 = domain->closest_image(i,iH1);
+ iH2 = domain->closest_image(i,iH2);
compute_newsite_thr(x[i],x[iH1],x[iH2],newsite_thr[i]);
hneigh_thr[i].t = 1;
hneigh_thr[i].b = iH2;
hneigh_thr[i].a = iH1;
+
} else {
iH1 = hneigh_thr[i].a;
iH2 = hneigh_thr[i].b;
@@ -229,10 +233,14 @@ void PairTIP4PLongSoftOMP::eval(int iifrom, int iito, ThrData * const thr)
error->one(FLERR,"TIP4P hydrogen is missing");
if (atom->type[jH1] != typeH || atom->type[jH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+ // set jH1,jH2 to closest image to O
+ jH1 = domain->closest_image(j,jH1);
+ jH2 = domain->closest_image(j,jH2);
compute_newsite_thr(x[j],x[jH1],x[jH2],newsite_thr[j]);
hneigh_thr[j].t = 1;
hneigh_thr[j].b = jH2;
hneigh_thr[j].a = jH1;
+
} else {
jH1 = hneigh_thr[j].a;
jH2 = hneigh_thr[j].b;
@@ -326,15 +334,14 @@ void PairTIP4PLongSoftOMP::eval(int iifrom, int iito, ThrData * const thr)
f[iH2].z += fHz;
if (VFLAG) {
- domain->closest_image(&x[i].x,&x[iH1].x,xH1);
- domain->closest_image(&x[i].x,&x[iH2].x,xH2);
-
- v[0] = x[i].x*fOx + xH1[0]*fHx + xH2[0]*fHx;
- v[1] = x[i].y*fOy + xH1[1]*fHy + xH2[1]*fHy;
- v[2] = x[i].z*fOz + xH1[2]*fHz + xH2[2]*fHz;
- v[3] = x[i].x*fOy + xH1[0]*fHy + xH2[0]*fHy;
- v[4] = x[i].x*fOz + xH1[0]*fHz + xH2[0]*fHz;
- v[5] = x[i].y*fOz + xH1[1]*fHz + xH2[1]*fHz;
+ xH1 = x[iH1];
+ xH2 = x[iH2];
+ v[0] = x[i].x*fOx + xH1.x*fHx + xH2.x*fHx;
+ v[1] = x[i].y*fOy + xH1.y*fHy + xH2.y*fHy;
+ v[2] = x[i].z*fOz + xH1.z*fHz + xH2.z*fHz;
+ v[3] = x[i].x*fOy + xH1.x*fHy + xH2.x*fHy;
+ v[4] = x[i].x*fOz + xH1.x*fHz + xH2.x*fHz;
+ v[5] = x[i].y*fOz + xH1.y*fHz + xH2.y*fHz;
}
if (EVFLAG) {
vlist[n++] = i;
@@ -386,15 +393,14 @@ void PairTIP4PLongSoftOMP::eval(int iifrom, int iito, ThrData * const thr)
f[jH2].z += fHz;
if (VFLAG) {
- domain->closest_image(&x[j].x,&x[jH1].x,xH1);
- domain->closest_image(&x[j].x,&x[jH2].x,xH2);
-
- v[0] += x[j].x*fOx + xH1[0]*fHx + xH2[0]*fHx;
- v[1] += x[j].y*fOy + xH1[1]*fHy + xH2[1]*fHy;
- v[2] += x[j].z*fOz + xH1[2]*fHz + xH2[2]*fHz;
- v[3] += x[j].x*fOy + xH1[0]*fHy + xH2[0]*fHy;
- v[4] += x[j].x*fOz + xH1[0]*fHz + xH2[0]*fHz;
- v[5] += x[j].y*fOz + xH1[1]*fHz + xH2[1]*fHz;
+ xH1 = x[jH1];
+ xH2 = x[jH2];
+ v[0] += x[j].x*fOx + xH1.x*fHx + xH2.x*fHx;
+ v[1] += x[j].y*fOy + xH1.y*fHy + xH2.y*fHy;
+ v[2] += x[j].z*fOz + xH1.z*fHz + xH2.z*fHz;
+ v[3] += x[j].x*fOy + xH1.x*fHy + xH2.x*fHy;
+ v[4] += x[j].x*fOz + xH1.x*fHz + xH2.x*fHz;
+ v[5] += x[j].y*fOz + xH1.y*fHz + xH2.y*fHz;
}
if (EVFLAG) {
vlist[n++] = j;
@@ -432,12 +438,10 @@ void PairTIP4PLongSoftOMP::compute_newsite_thr(const dbl3_t &xO,
double delx1 = xH1.x - xO.x;
double dely1 = xH1.y - xO.y;
double delz1 = xH1.z - xO.z;
- domain->minimum_image(delx1,dely1,delz1);
double delx2 = xH2.x - xO.x;
double dely2 = xH2.y - xO.y;
double delz2 = xH2.z - xO.z;
- domain->minimum_image(delx2,dely2,delz2);
const double prefac = alpha * 0.5;
xM.x = xO.x + prefac * (delx1 + delx2);
diff --git a/src/USER-OMP/pppm_disp_tip4p_omp.cpp b/src/USER-OMP/pppm_disp_tip4p_omp.cpp
index 6606c9602c..8872c849f3 100644
--- a/src/USER-OMP/pppm_disp_tip4p_omp.cpp
+++ b/src/USER-OMP/pppm_disp_tip4p_omp.cpp
@@ -1852,17 +1852,20 @@ void PPPMDispTIP4POMP::find_M_thr(int i, int &iH1, int &iH2, dbl3_t &xM)
if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+ // set iH1,iH2 to index of closest image to O
+
+ iH1 = domain->closest_image(i,iH1);
+ iH2 = domain->closest_image(i,iH2);
+
const dbl3_t * _noalias const x = (dbl3_t *) atom->x[0];
double delx1 = x[iH1].x - x[i].x;
double dely1 = x[iH1].y - x[i].y;
double delz1 = x[iH1].z - x[i].z;
- domain->minimum_image(delx1,dely1,delz1);
double delx2 = x[iH2].x - x[i].x;
double dely2 = x[iH2].y - x[i].y;
double delz2 = x[iH2].z - x[i].z;
- domain->minimum_image(delx2,dely2,delz2);
xM.x = x[i].x + alpha * 0.5 * (delx1 + delx2);
xM.y = x[i].y + alpha * 0.5 * (dely1 + dely2);
diff --git a/src/USER-OMP/pppm_tip4p_omp.cpp b/src/USER-OMP/pppm_tip4p_omp.cpp
index 21da813123..fdc58bce10 100644
--- a/src/USER-OMP/pppm_tip4p_omp.cpp
+++ b/src/USER-OMP/pppm_tip4p_omp.cpp
@@ -742,17 +742,20 @@ void PPPMTIP4POMP::find_M_thr(int i, int &iH1, int &iH2, dbl3_t &xM)
if (atom->type[iH1] != typeH || atom->type[iH2] != typeH)
error->one(FLERR,"TIP4P hydrogen has incorrect atom type");
+ // set iH1,iH2 to index of closest image to O
+
+ iH1 = domain->closest_image(i,iH1);
+ iH2 = domain->closest_image(i,iH2);
+
const dbl3_t * _noalias const x = (dbl3_t *) atom->x[0];
double delx1 = x[iH1].x - x[i].x;
double dely1 = x[iH1].y - x[i].y;
double delz1 = x[iH1].z - x[i].z;
- domain->minimum_image(delx1,dely1,delz1);
double delx2 = x[iH2].x - x[i].x;
double dely2 = x[iH2].y - x[i].y;
double delz2 = x[iH2].z - x[i].z;
- domain->minimum_image(delx2,dely2,delz2);
xM.x = x[i].x + alpha * 0.5 * (delx1 + delx2);
xM.y = x[i].y + alpha * 0.5 * (dely1 + dely2);
diff --git a/src/USER-REAXC/fix_qeq_reax.cpp b/src/USER-REAXC/fix_qeq_reax.cpp
index 5fa0ef6d79..26cf03f60a 100644
--- a/src/USER-REAXC/fix_qeq_reax.cpp
+++ b/src/USER-REAXC/fix_qeq_reax.cpp
@@ -440,13 +440,6 @@ void FixQEqReax::setup_pre_force_respa(int vflag, int ilevel)
/* ---------------------------------------------------------------------- */
-void FixQEqReax::min_setup_pre_force(int vflag)
-{
- setup_pre_force(vflag);
-}
-
-/* ---------------------------------------------------------------------- */
-
void FixQEqReax::init_storage()
{
int NN;
diff --git a/src/USER-REAXC/fix_qeq_reax.h b/src/USER-REAXC/fix_qeq_reax.h
index a0d29d6561..7c3e8a8f96 100644
--- a/src/USER-REAXC/fix_qeq_reax.h
+++ b/src/USER-REAXC/fix_qeq_reax.h
@@ -48,7 +48,6 @@ class FixQEqReax : public Fix {
void setup_pre_force_respa(int, int);
void pre_force_respa(int, int, int);
- void min_setup_pre_force(int);
void min_pre_force(int);
int matvecs;
diff --git a/src/domain.cpp b/src/domain.cpp
index 54183f6f2c..bc184e79b5 100644
--- a/src/domain.cpp
+++ b/src/domain.cpp
@@ -938,28 +938,31 @@ void Domain::subbox_too_small_check(double thresh)
}
/* ----------------------------------------------------------------------
- minimum image convention
+ minimum image convention in periodic dimensions
use 1/2 of box size as test
for triclinic, also add/subtract tilt factors in other dims as needed
+ changed "if" to "while" to enable distance to
+ far-away ghost atom returned by atom->map() to be wrapped back into box
+ could be problem for looking up atom IDs when cutoff > boxsize
------------------------------------------------------------------------- */
void Domain::minimum_image(double &dx, double &dy, double &dz)
{
if (triclinic == 0) {
if (xperiodic) {
- if (fabs(dx) > xprd_half) {
+ while (fabs(dx) > xprd_half) {
if (dx < 0.0) dx += xprd;
else dx -= xprd;
}
}
if (yperiodic) {
- if (fabs(dy) > yprd_half) {
+ while (fabs(dy) > yprd_half) {
if (dy < 0.0) dy += yprd;
else dy -= yprd;
}
}
if (zperiodic) {
- if (fabs(dz) > zprd_half) {
+ while (fabs(dz) > zprd_half) {
if (dz < 0.0) dz += zprd;
else dz -= zprd;
}
@@ -967,7 +970,7 @@ void Domain::minimum_image(double &dx, double &dy, double &dz)
} else {
if (zperiodic) {
- if (fabs(dz) > zprd_half) {
+ while (fabs(dz) > zprd_half) {
if (dz < 0.0) {
dz += zprd;
dy += yz;
@@ -980,7 +983,7 @@ void Domain::minimum_image(double &dx, double &dy, double &dz)
}
}
if (yperiodic) {
- if (fabs(dy) > yprd_half) {
+ while (fabs(dy) > yprd_half) {
if (dy < 0.0) {
dy += yprd;
dx += xy;
@@ -991,7 +994,7 @@ void Domain::minimum_image(double &dx, double &dy, double &dz)
}
}
if (xperiodic) {
- if (fabs(dx) > xprd_half) {
+ while (fabs(dx) > xprd_half) {
if (dx < 0.0) dx += xprd;
else dx -= xprd;
}
@@ -1000,28 +1003,31 @@ void Domain::minimum_image(double &dx, double &dy, double &dz)
}
/* ----------------------------------------------------------------------
- minimum image convention
+ minimum image convention in periodic dimensions
use 1/2 of box size as test
for triclinic, also add/subtract tilt factors in other dims as needed
+ changed "if" to "while" to enable distance to
+ far-away ghost atom returned by atom->map() to be wrapped back into box
+ could be problem for looking up atom IDs when cutoff > boxsize
------------------------------------------------------------------------- */
void Domain::minimum_image(double *delta)
{
if (triclinic == 0) {
if (xperiodic) {
- if (fabs(delta[0]) > xprd_half) {
+ while (fabs(delta[0]) > xprd_half) {
if (delta[0] < 0.0) delta[0] += xprd;
else delta[0] -= xprd;
}
}
if (yperiodic) {
- if (fabs(delta[1]) > yprd_half) {
+ while (fabs(delta[1]) > yprd_half) {
if (delta[1] < 0.0) delta[1] += yprd;
else delta[1] -= yprd;
}
}
if (zperiodic) {
- if (fabs(delta[2]) > zprd_half) {
+ while (fabs(delta[2]) > zprd_half) {
if (delta[2] < 0.0) delta[2] += zprd;
else delta[2] -= zprd;
}
@@ -1029,7 +1035,7 @@ void Domain::minimum_image(double *delta)
} else {
if (zperiodic) {
- if (fabs(delta[2]) > zprd_half) {
+ while (fabs(delta[2]) > zprd_half) {
if (delta[2] < 0.0) {
delta[2] += zprd;
delta[1] += yz;
@@ -1042,7 +1048,7 @@ void Domain::minimum_image(double *delta)
}
}
if (yperiodic) {
- if (fabs(delta[1]) > yprd_half) {
+ while (fabs(delta[1]) > yprd_half) {
if (delta[1] < 0.0) {
delta[1] += yprd;
delta[0] += xy;
@@ -1053,7 +1059,7 @@ void Domain::minimum_image(double *delta)
}
}
if (xperiodic) {
- if (fabs(delta[0]) > xprd_half) {
+ while (fabs(delta[0]) > xprd_half) {
if (delta[0] < 0.0) delta[0] += xprd;
else delta[0] -= xprd;
}
@@ -1092,12 +1098,16 @@ int Domain::closest_image(int i, int j)
closest = j;
}
}
+
return closest;
}
/* ----------------------------------------------------------------------
find and return Xj image = periodic image of Xj that is closest to Xi
for triclinic, add/subtract tilt factors in other dims as needed
+ not currently used (Jan 2017):
+ used to be called by pair TIP4P styles but no longer,
+ due to use of other closest_image() method
------------------------------------------------------------------------- */
void Domain::closest_image(const double * const xi, const double * const xj,
@@ -1984,7 +1994,7 @@ void Domain::subbox_corners()
/* ----------------------------------------------------------------------
compute 8 corner pts of any triclinic box with lo/hi in lamda coords
- 8 output conners are ordered with x changing fastest, then y, finally z
+ 8 output corners are ordered with x changing fastest, then y, finally z
could be more efficient if just coded with xy,yz,xz explicitly
------------------------------------------------------------------------- */
diff --git a/src/finish.cpp b/src/finish.cpp
index 850bc0b019..9222d4e868 100644
--- a/src/finish.cpp
+++ b/src/finish.cpp
@@ -65,7 +65,8 @@ void Finish::end(int flag)
{
int i,m,nneigh,nneighfull;
int histo[10];
- int minflag,prdflag,tadflag,timeflag,fftflag,histoflag,neighflag;
+ int minflag,prdflag,tadflag,hyperflag;
+ int timeflag,fftflag,histoflag,neighflag;
double time,tmp,ave,max,min;
double time_loop,time_other,cpu_loop;
@@ -86,9 +87,11 @@ void Finish::end(int flag)
// flag = 1 = dynamics or minimization
// flag = 2 = PRD
// flag = 3 = TAD
+ // flag = 4 = HYPER
// turn off neighflag for Kspace partition of verlet/split integrator
- minflag = prdflag = tadflag = timeflag = fftflag = histoflag = neighflag = 0;
+ minflag = prdflag = tadflag = hyperflag = 0;
+ timeflag = fftflag = histoflag = neighflag = 0;
time_loop = cpu_loop = time_other = 0.0;
if (flag == 1) {
@@ -103,6 +106,7 @@ void Finish::end(int flag)
}
if (flag == 2) prdflag = timeflag = histoflag = neighflag = 1;
if (flag == 3) tadflag = histoflag = neighflag = 1;
+ if (flag == 4) hyperflag = timeflag = histoflag = neighflag = 1;
// loop stats
@@ -161,8 +165,10 @@ void Finish::end(int flag)
if (lmp->kokkos) {
const char fmt2[] =
"%.1f%% CPU use with %d MPI tasks x %d OpenMP threads\n";
- if (screen) fprintf(screen,fmt2,cpu_loop,nprocs,lmp->kokkos->num_threads);
- if (logfile) fprintf(logfile,fmt2,cpu_loop,nprocs,lmp->kokkos->num_threads);
+ if (screen) fprintf(screen,fmt2,cpu_loop,nprocs,
+ lmp->kokkos->num_threads);
+ if (logfile) fprintf(logfile,fmt2,cpu_loop,nprocs,
+ lmp->kokkos->num_threads);
} else {
#if defined(_OPENMP)
const char fmt2[] =
@@ -409,12 +415,61 @@ void Finish::end(int flag)
}
}
+ // HYPER stats using PAIR,BOND,KSPACE for dynamics,quench
+
+ if (hyperflag) {
+ if (me == 0) {
+ if (screen) fprintf(screen,"\nHyper stats:\n");
+ if (logfile) fprintf(logfile,"\nHyper stats:\n");
+ }
+
+ time = timer->get_wall(Timer::DYNAMICS);
+ MPI_Allreduce(&time,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
+ time = tmp/nprocs;
+ if (me == 0) {
+ if (screen)
+ fprintf(screen," Dynamics time (%%) = %g (%g)\n",
+ time,time/time_loop*100.0);
+ if (logfile)
+ fprintf(logfile," Dynamics time (%%) = %g (%g)\n",
+ time,time/time_loop*100.0);
+ }
+
+ time = timer->get_wall(Timer::QUENCH);
+ MPI_Allreduce(&time,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
+ time = tmp/nprocs;
+ if (me == 0) {
+ if (screen)
+ fprintf(screen," Quench time (%%) = %g (%g)\n",
+ time,time/time_loop*100.0);
+ if (logfile)
+ fprintf(logfile," Quench time (%%) = %g (%g)\n",
+ time,time/time_loop*100.0);
+ }
+
+ time = time_other;
+ MPI_Allreduce(&time,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
+ time = tmp/nprocs;
+ if (me == 0) {
+ if (screen)
+ fprintf(screen," Other time (%%) = %g (%g)\n",
+ time,time/time_loop*100.0);
+ if (logfile)
+ fprintf(logfile," Other time (%%) = %g (%g)\n",
+ time,time/time_loop*100.0);
+ }
+ }
+
+ // further timing breakdowns
+
if (timeflag && timer->has_normal()) {
if (timer->has_full()) {
const char hdr[] = "\nMPI task timing breakdown:\n"
- "Section | min time | avg time | max time |%varavg| %CPU | %total\n"
- "-----------------------------------------------------------------------\n";
+ "Section | min time | avg time | max time "
+ "|%varavg| %CPU | %total\n"
+ "-----------------------------------------------"
+ "------------------------\n";
if (me == 0) {
if (screen) fputs(hdr,screen);
if (logfile) fputs(hdr,logfile);
@@ -516,8 +571,10 @@ void Finish::end(int flag)
if (lmp->kokkos && lmp->kokkos->ngpu > 0)
if (const char* env_clb = getenv("CUDA_LAUNCH_BLOCKING"))
if (!(strcmp(env_clb,"1") == 0)) {
- error->warning(FLERR,"Timing breakdown may not be accurate since GPU/CPU overlap is enabled. "
- "Using 'export CUDA_LAUNCH_BLOCKING=1' will give an accurate timing breakdown but will reduce performance");
+ error->warning(FLERR,"Timing breakdown may not be accurate "
+ "since GPU/CPU overlap is enabled\n"
+ "Using 'export CUDA_LAUNCH_BLOCKING=1' will give an "
+ "accurate timing breakdown but will reduce performance");
}
// FFT timing statistics
diff --git a/src/fix.h b/src/fix.h
index 65ababf287..bd3189afe7 100644
--- a/src/fix.h
+++ b/src/fix.h
@@ -113,6 +113,7 @@ class Fix : protected Pointers {
virtual void setup_pre_exchange() {}
virtual void setup_pre_neighbor() {}
virtual void setup_pre_force(int) {}
+ virtual void setup_pre_reverse(int, int) {}
virtual void min_setup(int) {}
virtual void initial_integrate(int) {}
virtual void post_integrate() {}
@@ -151,12 +152,10 @@ class Fix : protected Pointers {
virtual void post_force_respa(int, int, int) {}
virtual void final_integrate_respa(int, int) {}
- virtual void min_setup_pre_exchange() {}
- virtual void min_setup_pre_neighbor() {}
- virtual void min_setup_pre_force(int) {}
virtual void min_pre_exchange() {}
virtual void min_pre_neighbor() {}
virtual void min_pre_force(int) {}
+ virtual void min_pre_reverse(int,int) {}
virtual void min_post_force(int) {}
virtual double min_energy(double *) {return 0.0;}
@@ -241,23 +240,24 @@ namespace FixConst {
static const int PRE_EXCHANGE = 1<<2;
static const int PRE_NEIGHBOR = 1<<3;
static const int PRE_FORCE = 1<<4;
- static const int POST_FORCE = 1<<5;
- static const int FINAL_INTEGRATE = 1<<6;
- static const int END_OF_STEP = 1<<7;
- static const int THERMO_ENERGY = 1<<8;
- static const int INITIAL_INTEGRATE_RESPA = 1<<9;
- static const int POST_INTEGRATE_RESPA = 1<<10;
- static const int PRE_FORCE_RESPA = 1<<11;
- static const int POST_FORCE_RESPA = 1<<12;
- static const int FINAL_INTEGRATE_RESPA = 1<<13;
- static const int MIN_PRE_EXCHANGE = 1<<14;
- static const int MIN_PRE_NEIGHBOR = 1<<15;
- static const int MIN_PRE_FORCE = 1<<16;
- 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 PRE_REVERSE = 1<<20;
- static const int FIX_CONST_LAST = 1<<21;
+ static const int PRE_REVERSE = 1<<5;
+ static const int POST_FORCE = 1<<6;
+ static const int FINAL_INTEGRATE = 1<<7;
+ static const int END_OF_STEP = 1<<8;
+ static const int POST_RUN = 1<<9;
+ static const int THERMO_ENERGY = 1<<10;
+ static const int INITIAL_INTEGRATE_RESPA = 1<<11;
+ static const int POST_INTEGRATE_RESPA = 1<<12;
+ static const int PRE_FORCE_RESPA = 1<<13;
+ static const int POST_FORCE_RESPA = 1<<14;
+ static const int FINAL_INTEGRATE_RESPA = 1<<15;
+ static const int MIN_PRE_EXCHANGE = 1<<16;
+ static const int MIN_PRE_NEIGHBOR = 1<<17;
+ static const int MIN_PRE_FORCE = 1<<18;
+ static const int MIN_PRE_REVERSE = 1<<19;
+ static const int MIN_POST_FORCE = 1<<20;
+ static const int MIN_ENERGY = 1<<21;
+ static const int FIX_CONST_LAST = 1<<22;
}
}
diff --git a/src/fix_restrain.cpp b/src/fix_restrain.cpp
index 067f1014d5..f276be76af 100644
--- a/src/fix_restrain.cpp
+++ b/src/fix_restrain.cpp
@@ -45,7 +45,8 @@ enum{BOND,ANGLE,DIHEDRAL};
FixRestrain::FixRestrain(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg),
- rstyle(NULL), ids(NULL), kstart(NULL), kstop(NULL), target(NULL), cos_target(NULL), sin_target(NULL)
+ rstyle(NULL), ids(NULL), kstart(NULL), kstop(NULL), target(NULL),
+ cos_target(NULL), sin_target(NULL)
{
if (narg < 4) error->all(FLERR,"Illegal fix restrain command");
diff --git a/src/min.cpp b/src/min.cpp
index 6e568aee0b..9207c6bdc2 100644
--- a/src/min.cpp
+++ b/src/min.cpp
@@ -284,7 +284,7 @@ void Min::setup()
else force->kspace->compute_dummy(eflag,vflag);
}
- modify->pre_reverse(eflag,vflag);
+ modify->setup_pre_reverse(eflag,vflag);
if (force->newton) comm->reverse_comm();
// update per-atom minimization variables stored by pair styles
@@ -365,7 +365,7 @@ void Min::setup_minimal(int flag)
else force->kspace->compute_dummy(eflag,vflag);
}
- modify->pre_reverse(eflag,vflag);
+ modify->setup_pre_reverse(eflag,vflag);
if (force->newton) comm->reverse_comm();
// update per-atom minimization variables stored by pair styles
@@ -493,6 +493,11 @@ double Min::energy_force(int resetflag)
comm->borders();
if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
timer->stamp(Timer::COMM);
+ if (modify->n_min_pre_neighbor) {
+ timer->stamp();
+ modify->min_pre_neighbor();
+ timer->stamp(Timer::MODIFY);
+ }
neighbor->build();
timer->stamp(Timer::NEIGH);
}
@@ -525,8 +530,8 @@ double Min::energy_force(int resetflag)
timer->stamp(Timer::KSPACE);
}
- if (modify->n_pre_reverse) {
- modify->pre_reverse(eflag,vflag);
+ if (modify->n_min_pre_reverse) {
+ modify->min_pre_reverse(eflag,vflag);
timer->stamp(Timer::MODIFY);
}
diff --git a/src/modify.cpp b/src/modify.cpp
index 2e5bd78504..14e55ab414 100644
--- a/src/modify.cpp
+++ b/src/modify.cpp
@@ -47,7 +47,8 @@ Modify::Modify(LAMMPS *lmp) : Pointers(lmp)
n_thermo_energy_atom = 0;
n_initial_integrate_respa = n_post_integrate_respa = 0;
n_pre_force_respa = n_post_force_respa = n_final_integrate_respa = 0;
- n_min_pre_exchange = n_min_pre_force = n_min_post_force = n_min_energy = 0;
+ n_min_pre_exchange = n_min_pre_force = n_min_pre_reverse = 0;
+ n_min_post_force = n_min_energy = 0;
fix = NULL;
fmask = NULL;
@@ -60,7 +61,7 @@ Modify::Modify(LAMMPS *lmp) : Pointers(lmp)
list_pre_force_respa = list_post_force_respa = NULL;
list_final_integrate_respa = NULL;
list_min_pre_exchange = list_min_pre_neighbor = NULL;
- list_min_pre_force = list_min_post_force = NULL;
+ list_min_pre_force = list_min_pre_reverse = list_min_post_force = NULL;
list_min_energy = NULL;
end_of_step_every = NULL;
@@ -136,6 +137,7 @@ Modify::~Modify()
delete [] list_min_pre_exchange;
delete [] list_min_pre_neighbor;
delete [] list_min_pre_force;
+ delete [] list_min_pre_reverse;
delete [] list_min_post_force;
delete [] list_min_energy;
@@ -188,6 +190,7 @@ void Modify::init()
list_init(MIN_PRE_EXCHANGE,n_min_pre_exchange,list_min_pre_exchange);
list_init(MIN_PRE_NEIGHBOR,n_min_pre_neighbor,list_min_pre_neighbor);
list_init(MIN_PRE_FORCE,n_min_pre_force,list_min_pre_force);
+ list_init(MIN_PRE_REVERSE,n_min_pre_reverse,list_min_pre_reverse);
list_init(MIN_POST_FORCE,n_min_post_force,list_min_post_force);
list_init(MIN_ENERGY,n_min_energy,list_min_energy);
@@ -307,7 +310,7 @@ void Modify::setup_pre_exchange()
fix[list_pre_exchange[i]]->setup_pre_exchange();
else if (update->whichflag == 2)
for (int i = 0; i < n_min_pre_exchange; i++)
- fix[list_min_pre_exchange[i]]->min_setup_pre_exchange();
+ fix[list_min_pre_exchange[i]]->setup_pre_exchange();
}
/* ----------------------------------------------------------------------
@@ -322,7 +325,7 @@ void Modify::setup_pre_neighbor()
fix[list_pre_neighbor[i]]->setup_pre_neighbor();
else if (update->whichflag == 2)
for (int i = 0; i < n_min_pre_neighbor; i++)
- fix[list_min_pre_neighbor[i]]->min_setup_pre_neighbor();
+ fix[list_min_pre_neighbor[i]]->setup_pre_neighbor();
}
/* ----------------------------------------------------------------------
@@ -337,7 +340,22 @@ void Modify::setup_pre_force(int vflag)
fix[list_pre_force[i]]->setup_pre_force(vflag);
else if (update->whichflag == 2)
for (int i = 0; i < n_min_pre_force; i++)
- fix[list_min_pre_force[i]]->min_setup_pre_force(vflag);
+ fix[list_min_pre_force[i]]->setup_pre_force(vflag);
+}
+
+/* ----------------------------------------------------------------------
+ setup pre_reverse call, only for fixes that define pre_reverse
+ called from Verlet, RESPA, Min
+------------------------------------------------------------------------- */
+
+void Modify::setup_pre_reverse(int eflag, int vflag)
+{
+ if (update->whichflag == 1)
+ for (int i = 0; i < n_pre_reverse; i++)
+ fix[list_pre_reverse[i]]->setup_pre_reverse(eflag,vflag);
+ else if (update->whichflag == 2)
+ for (int i = 0; i < n_min_pre_reverse; i++)
+ fix[list_min_pre_reverse[i]]->setup_pre_reverse(eflag,vflag);
}
/* ----------------------------------------------------------------------
@@ -580,6 +598,16 @@ void Modify::min_pre_force(int vflag)
fix[list_min_pre_force[i]]->min_pre_force(vflag);
}
+/* ----------------------------------------------------------------------
+ minimizer pre-reverse call, only for relevant fixes
+------------------------------------------------------------------------- */
+
+void Modify::min_pre_reverse(int eflag, int vflag)
+{
+ for (int i = 0; i < n_min_pre_reverse; i++)
+ fix[list_min_pre_reverse[i]]->min_pre_reverse(eflag,vflag);
+}
+
/* ----------------------------------------------------------------------
minimizer force adjustment call, only for relevant fixes
------------------------------------------------------------------------- */
diff --git a/src/modify.h b/src/modify.h
index f6f59a6a9c..d6793a754e 100644
--- a/src/modify.h
+++ b/src/modify.h
@@ -35,7 +35,7 @@ class Modify : protected Pointers {
int n_initial_integrate_respa,n_post_integrate_respa;
int n_pre_force_respa,n_post_force_respa,n_final_integrate_respa;
int n_min_pre_exchange,n_min_pre_neighbor;
- int n_min_pre_force,n_min_post_force,n_min_energy;
+ int n_min_pre_force,n_min_pre_reverse,n_min_post_force,n_min_energy;
int restart_pbc_any; // 1 if any fix sets restart_pbc
int nfix_restart_global; // stored fix global info from restart file
@@ -54,6 +54,7 @@ class Modify : protected Pointers {
virtual void setup_pre_exchange();
virtual void setup_pre_neighbor();
virtual void setup_pre_force(int);
+ virtual void setup_pre_reverse(int, int);
virtual void initial_integrate(int);
virtual void post_integrate();
virtual void pre_exchange();
@@ -78,6 +79,7 @@ class Modify : protected Pointers {
virtual void min_pre_exchange();
virtual void min_pre_neighbor();
virtual void min_pre_force(int);
+ virtual void min_pre_reverse(int,int);
virtual void min_post_force(int);
virtual double min_energy(double *);
@@ -124,7 +126,7 @@ class Modify : protected Pointers {
int *list_pre_force_respa,*list_post_force_respa;
int *list_final_integrate_respa;
int *list_min_pre_exchange,*list_min_pre_neighbor;
- int *list_min_pre_force,*list_min_post_force;
+ int *list_min_pre_force,*list_min_pre_reverse,*list_min_post_force;
int *list_min_energy;
int *end_of_step_every;
diff --git a/src/respa.cpp b/src/respa.cpp
index db0f11d5b3..7646115fa9 100644
--- a/src/respa.cpp
+++ b/src/respa.cpp
@@ -44,8 +44,10 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
-Respa::Respa(LAMMPS *lmp, int narg, char **arg) : Integrate(lmp, narg, arg),
-step(NULL), loop(NULL), hybrid_level(NULL), hybrid_compute(NULL), newton(NULL), fix_respa(NULL)
+Respa::Respa(LAMMPS *lmp, int narg, char **arg) :
+ Integrate(lmp, narg, arg),
+ step(NULL), loop(NULL), hybrid_level(NULL), hybrid_compute(NULL),
+ newton(NULL), fix_respa(NULL)
{
nhybrid_styles = 0;
if (narg < 1) error->all(FLERR,"Illegal run_style respa command");
@@ -473,7 +475,7 @@ void Respa::setup()
if (kspace_compute_flag) force->kspace->compute(eflag,vflag);
}
- modify->pre_reverse(eflag,vflag);
+ modify->setup_pre_reverse(eflag,vflag);
if (newton[ilevel]) comm->reverse_comm();
copy_f_flevel(ilevel);
}
@@ -549,7 +551,7 @@ void Respa::setup_minimal(int flag)
if (kspace_compute_flag) force->kspace->compute(eflag,vflag);
}
- modify->pre_reverse(eflag,vflag);
+ modify->setup_pre_reverse(eflag,vflag);
if (newton[ilevel]) comm->reverse_comm();
copy_f_flevel(ilevel);
}
diff --git a/src/run.cpp b/src/run.cpp
index f238e7852b..3317545342 100644
--- a/src/run.cpp
+++ b/src/run.cpp
@@ -43,6 +43,7 @@ void Run::command(int narg, char **arg)
error->all(FLERR,"Run command before simulation box is defined");
// ignore run command, if walltime limit was already reached
+
if (timer->is_timeout()) return;
bigint nsteps_input = force->bnumeric(FLERR,arg[0]);
diff --git a/src/verlet.cpp b/src/verlet.cpp
index 035cab79ef..915648040e 100644
--- a/src/verlet.cpp
+++ b/src/verlet.cpp
@@ -144,7 +144,7 @@ void Verlet::setup()
else force->kspace->compute_dummy(eflag,vflag);
}
- modify->pre_reverse(eflag,vflag);
+ modify->setup_pre_reverse(eflag,vflag);
if (force->newton) comm->reverse_comm();
modify->setup(vflag);
@@ -205,7 +205,7 @@ void Verlet::setup_minimal(int flag)
else force->kspace->compute_dummy(eflag,vflag);
}
- modify->pre_reverse(eflag,vflag);
+ modify->setup_pre_reverse(eflag,vflag);
if (force->newton) comm->reverse_comm();
modify->setup(vflag);
--
GitLab