diff --git a/src/KSPACE/pppm.cpp b/src/KSPACE/pppm.cpp
index 232cabe4ea9bed54f8463058bd5f78cefee89b72..05169fca1c05d808e2ff1ea05ff1fb9508be4445 100644
--- a/src/KSPACE/pppm.cpp
+++ b/src/KSPACE/pppm.cpp
@@ -186,6 +186,10 @@ void PPPM::init()
// error check
triclinic_check();
+
+ if (triclinic != domain->triclinic)
+ error->all(FLERR,"Must redefine kspace_style after changing to triclinic box");
+
if (domain->triclinic && differentiation_flag == 1)
error->all(FLERR,"Cannot (yet) use PPPM with triclinic box "
"and kspace_modify diff ad");
diff --git a/src/KSPACE/pppm_tip4p.cpp b/src/KSPACE/pppm_tip4p.cpp
index 259def7c9aebd1901c2577b93a9f3b621ee8f8e1..e807f7e2e98a8ecbddca91fdc6ddb0d29c3ff46d 100644
--- a/src/KSPACE/pppm_tip4p.cpp
+++ b/src/KSPACE/pppm_tip4p.cpp
@@ -332,9 +332,9 @@ void PPPMTIP4P::fieldforce_ad()
const double qfactor = qqrd2e * scale;
- s1 = x[i][0]*hx_inv;
- s2 = x[i][1]*hy_inv;
- s3 = x[i][2]*hz_inv;
+ s1 = xi[0]*hx_inv;
+ s2 = xi[1]*hy_inv;
+ s3 = xi[2]*hz_inv;
sf = sf_coeff[0]*sin(2*MY_PI*s1);
sf += sf_coeff[1]*sin(4*MY_PI*s1);
sf *= 2.0*q[i]*q[i];
@@ -486,6 +486,8 @@ void PPPMTIP4P::fieldforce_peratom()
void PPPMTIP4P::find_M(int i, int &iH1, int &iH2, double *xM)
{
+ double **x = atom->x;
+
iH1 = atom->map(atom->tag[i] + 1);
iH2 = atom->map(atom->tag[i] + 2);
@@ -493,63 +495,98 @@ 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);
-
if (triclinic) {
- find_M_triclinic(i,iH1,iH2,xM);
- return;
- }
- double **x = atom->x;
+ // need to use custom code to find the closest image for triclinic,
+ // since local atoms are in lambda coordinates, but ghosts are not.
+
+ int *sametag = atom->sametag;
+ double xo[3],xh1[3],xh2[3];
+
+ domain->lamda2x(x[i],xo);
+ domain->lamda2x(x[iH1],xh1);
+ domain->lamda2x(x[iH2],xh2);
+
+ double delx = xo[0] - xh1[0];
+ double dely = xo[1] - xh1[1];
+ double delz = xo[2] - xh1[2];
+ double rsqmin = delx*delx + dely*dely + delz*delz;
+ double rsq;
+ int closest = iH1;
+
+ while (sametag[iH1] >= 0) {
+ iH1 = sametag[iH1];
+ delx = xo[0] - x[iH1][0];
+ dely = xo[1] - x[iH1][1];
+ delz = xo[2] - x[iH1][2];
+ rsq = delx*delx + dely*dely + delz*delz;
+ if (rsq < rsqmin) {
+ rsqmin = rsq;
+ closest = iH1;
+ xh1[0] = x[iH1][0];
+ xh1[1] = x[iH1][1];
+ xh1[2] = x[iH1][2];
+ }
+ }
+ iH1 = closest;
+
+ closest = iH2;
+ delx = xo[0] - xh2[0];
+ dely = xo[1] - xh2[1];
+ delz = xo[2] - xh2[2];
+ rsqmin = delx*delx + dely*dely + delz*delz;
+
+ while (sametag[iH2] >= 0) {
+ iH2 = sametag[iH2];
+ delx = xo[0] - x[iH2][0];
+ dely = xo[1] - x[iH2][1];
+ delz = xo[2] - x[iH2][2];
+ rsq = delx*delx + dely*dely + delz*delz;
+ if (rsq < rsqmin) {
+ rsqmin = rsq;
+ closest = iH2;
+ xh2[0] = x[iH2][0];
+ xh2[1] = x[iH2][1];
+ xh2[2] = x[iH2][2];
+ }
+ }
+ iH2 = closest;
- double delx1 = x[iH1][0] - x[i][0];
- double dely1 = x[iH1][1] - x[i][1];
- double delz1 = x[iH1][2] - x[i][2];
+ // finally compute M in real coordinates ...
- double delx2 = x[iH2][0] - x[i][0];
- double dely2 = x[iH2][1] - x[i][1];
- double delz2 = x[iH2][2] - x[i][2];
+ double delx1 = xh1[0] - xo[0];
+ double dely1 = xh1[1] - xo[1];
+ double delz1 = xh1[2] - xo[2];
- xM[0] = x[i][0] + alpha * 0.5 * (delx1 + delx2);
- xM[1] = x[i][1] + alpha * 0.5 * (dely1 + dely2);
- xM[2] = x[i][2] + alpha * 0.5 * (delz1 + delz2);
-}
+ double delx2 = xh2[0] - xo[0];
+ double dely2 = xh2[1] - xo[1];
+ double delz2 = xh2[2] - xo[2];
-/* ----------------------------------------------------------------------
- triclinic version of find_M()
- calculation of M coords done in real space
-------------------------------------------------------------------------- */
+ xM[0] = xo[0] + alpha * 0.5 * (delx1 + delx2);
+ xM[1] = xo[1] + alpha * 0.5 * (dely1 + dely2);
+ xM[2] = xo[2] + alpha * 0.5 * (delz1 + delz2);
-void PPPMTIP4P::find_M_triclinic(int i, int iH1, int iH2, double *xM)
-{
- double **x = atom->x;
+ // ... and convert M to lamda space for PPPM
- // convert from lamda to real space
- // may be possible to do this more efficiently in lamda space
+ domain->x2lamda(xM,xM);
- domain->lamda2x(x[i],x[i]);
- domain->lamda2x(x[iH1],x[iH1]);
- domain->lamda2x(x[iH2],x[iH2]);
+ } else {
- double delx1 = x[iH1][0] - x[i][0];
- double dely1 = x[iH1][1] - x[i][1];
- double delz1 = x[iH1][2] - x[i][2];
+ // set iH1,iH2 to index of closest image to O
- double delx2 = x[iH2][0] - x[i][0];
- double dely2 = x[iH2][1] - x[i][1];
- double delz2 = x[iH2][2] - x[i][2];
+ iH1 = domain->closest_image(i,iH1);
+ iH2 = domain->closest_image(i,iH2);
- xM[0] = x[i][0] + alpha * 0.5 * (delx1 + delx2);
- xM[1] = x[i][1] + alpha * 0.5 * (dely1 + dely2);
- xM[2] = x[i][2] + alpha * 0.5 * (delz1 + delz2);
+ double delx1 = x[iH1][0] - x[i][0];
+ double dely1 = x[iH1][1] - x[i][1];
+ double delz1 = x[iH1][2] - x[i][2];
- // convert back to lamda space
+ 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->x2lamda(x[i],x[i]);
- domain->x2lamda(x[iH1],x[iH1]);
- domain->x2lamda(x[iH2],x[iH2]);
- domain->x2lamda(xM,xM);
+ xM[0] = x[i][0] + alpha * 0.5 * (delx1 + delx2);
+ xM[1] = x[i][1] + alpha * 0.5 * (dely1 + dely2);
+ xM[2] = x[i][2] + alpha * 0.5 * (delz1 + delz2);
+ }
}
diff --git a/src/KSPACE/pppm_tip4p.h b/src/KSPACE/pppm_tip4p.h
index 31e69f785ad322a790d504fb4dd11a2bbc8ca7e6..88a5707f2da95e5691e8079675bc08010814f70a 100644
--- a/src/KSPACE/pppm_tip4p.h
+++ b/src/KSPACE/pppm_tip4p.h
@@ -39,7 +39,6 @@ class PPPMTIP4P : public PPPM {
private:
void find_M(int, int &, int &, double *);
- void find_M_triclinic(int, int, int, double *);
};
}
diff --git a/src/USER-OMP/pppm_tip4p_omp.cpp b/src/USER-OMP/pppm_tip4p_omp.cpp
index fdc58bce10437c8b83f340f06297bf3da32085e1..21e3081c65b3e84eb2e0c0b31e74bf82e0af8b42 100644
--- a/src/USER-OMP/pppm_tip4p_omp.cpp
+++ b/src/USER-OMP/pppm_tip4p_omp.cpp
@@ -48,7 +48,7 @@ using namespace MathSpecial;
PPPMTIP4POMP::PPPMTIP4POMP(LAMMPS *lmp, int narg, char **arg) :
PPPMTIP4P(lmp, narg, arg), ThrOMP(lmp, THR_KSPACE)
{
- triclinic_support = 0;
+ triclinic_support = 1;
suffix_flag |= Suffix::OMP;
}
@@ -735,6 +735,8 @@ void PPPMTIP4POMP::fieldforce_ad()
void PPPMTIP4POMP::find_M_thr(int i, int &iH1, int &iH2, dbl3_t &xM)
{
+ double **x = atom->x;
+
iH1 = atom->map(atom->tag[i] + 1);
iH2 = atom->map(atom->tag[i] + 2);
@@ -742,24 +744,102 @@ 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
+ if (triclinic) {
+
+ // need to use custom code to find the closest image for triclinic,
+ // since local atoms are in lambda coordinates, but ghosts are not.
+
+ int *sametag = atom->sametag;
+ double xo[3],xh1[3],xh2[3];
+
+ domain->lamda2x(x[i],xo);
+ domain->lamda2x(x[iH1],xh1);
+ domain->lamda2x(x[iH2],xh2);
+
+ double delx = xo[0] - xh1[0];
+ double dely = xo[1] - xh1[1];
+ double delz = xo[2] - xh1[2];
+ double rsqmin = delx*delx + dely*dely + delz*delz;
+ double rsq;
+ int closest = iH1;
+
+ while (sametag[iH1] >= 0) {
+ iH1 = sametag[iH1];
+ delx = xo[0] - x[iH1][0];
+ dely = xo[1] - x[iH1][1];
+ delz = xo[2] - x[iH1][2];
+ rsq = delx*delx + dely*dely + delz*delz;
+ if (rsq < rsqmin) {
+ rsqmin = rsq;
+ closest = iH1;
+ xh1[0] = x[iH1][0];
+ xh1[1] = x[iH1][1];
+ xh1[2] = x[iH1][2];
+ }
+ }
+ iH1 = closest;
+
+ closest = iH2;
+ delx = xo[0] - xh2[0];
+ dely = xo[1] - xh2[1];
+ delz = xo[2] - xh2[2];
+ rsqmin = delx*delx + dely*dely + delz*delz;
+
+ while (sametag[iH2] >= 0) {
+ iH2 = sametag[iH2];
+ delx = xo[0] - x[iH2][0];
+ dely = xo[1] - x[iH2][1];
+ delz = xo[2] - x[iH2][2];
+ rsq = delx*delx + dely*dely + delz*delz;
+ if (rsq < rsqmin) {
+ rsqmin = rsq;
+ closest = iH2;
+ xh2[0] = x[iH2][0];
+ xh2[1] = x[iH2][1];
+ xh2[2] = x[iH2][2];
+ }
+ }
+ iH2 = closest;
+
+ // finally compute M in real coordinates ...
- iH1 = domain->closest_image(i,iH1);
- iH2 = domain->closest_image(i,iH2);
+ double delx1 = xh1[0] - xo[0];
+ double dely1 = xh1[1] - xo[1];
+ double delz1 = xh1[2] - xo[2];
- const dbl3_t * _noalias const x = (dbl3_t *) atom->x[0];
+ double delx2 = xh2[0] - xo[0];
+ double dely2 = xh2[1] - xo[1];
+ double delz2 = xh2[2] - xo[2];
+
+ xM.x = xo[0] + alpha * 0.5 * (delx1 + delx2);
+ xM.y = xo[1] + alpha * 0.5 * (dely1 + dely2);
+ xM.z = xo[2] + alpha * 0.5 * (delz1 + delz2);
- double delx1 = x[iH1].x - x[i].x;
- double dely1 = x[iH1].y - x[i].y;
- double delz1 = x[iH1].z - x[i].z;
+ // ... and convert M to lamda space for PPPM
- 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->x2lamda((double *)&xM,(double *)&xM);
- xM.x = x[i].x + alpha * 0.5 * (delx1 + delx2);
- xM.y = x[i].y + alpha * 0.5 * (dely1 + dely2);
- xM.z = x[i].z + alpha * 0.5 * (delz1 + delz2);
+ } else {
+
+ // 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;
+
+ double delx2 = x[iH2].x - x[i].x;
+ double dely2 = x[iH2].y - x[i].y;
+ double delz2 = x[iH2].z - x[i].z;
+
+ xM.x = x[i].x + alpha * 0.5 * (delx1 + delx2);
+ xM.y = x[i].y + alpha * 0.5 * (dely1 + dely2);
+ xM.z = x[i].z + alpha * 0.5 * (delz1 + delz2);
+ }
}