Commit 9eeb97b0 authored by sjplimp's avatar sjplimp Committed by GitHub
Browse files

Merge pull request #544 from akohlmey/tip4p-triclinic 

Correct handling of triclinic box support in pppm/tip4p and pppm/tip4p/omp
parents 883b7aaa 73b948dc

src/KSPACE/pppm.cpp

+4 −0
Original line number Diff line number Diff line
@@ -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");

src/KSPACE/pppm_tip4p.cpp

+85 −48
Original line number Diff line number Diff line
@@ -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,46 +495,87 @@ 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;



  double delx1 = x[iH1][0] - x[i][0];

  double dely1 = x[iH1][1] - x[i][1];

  double delz1 = x[iH1][2] - x[i][2];

    // 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 ...



    double delx1 = xh1[0] - xo[0];

    double dely1 = xh1[1] - xo[1];

    double delz1 = xh1[2] - xo[2];



    double delx2 = xh2[0] - xo[0];

    double dely2 = xh2[1] - xo[1];

    double delz2 = xh2[2] - xo[2];



    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);



    // ... and convert M to lamda space for PPPM



  double delx2 = x[iH2][0] - x[i][0];

  double dely2 = x[iH2][1] - x[i][1];

  double delz2 = x[iH2][2] - x[i][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);

}



/* ----------------------------------------------------------------------

   triclinic version of find_M()

   calculation of M coords done in real space

------------------------------------------------------------------------- */

    domain->x2lamda(xM,xM);



void PPPMTIP4P::find_M_triclinic(int i, int iH1, int iH2, double *xM)

{

  double **x = atom->x;

  } else {



  // convert from lamda to real space

  // may be possible to do this more efficiently in lamda space

    // set iH1,iH2 to index of closest image to O



  domain->lamda2x(x[i],x[i]);

  domain->lamda2x(x[iH1],x[iH1]);

  domain->lamda2x(x[iH2],x[iH2]);

    iH1 = domain->closest_image(i,iH1);

    iH2 = domain->closest_image(i,iH2);



    double delx1 = x[iH1][0] - x[i][0];

    double dely1 = x[iH1][1] - x[i][1];
@@ -545,11 +588,5 @@ void PPPMTIP4P::find_M_triclinic(int i, int iH1, int iH2, double *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);



  // convert back to lamda space



  domain->x2lamda(x[i],x[i]);

  domain->x2lamda(x[iH1],x[iH1]);

  domain->x2lamda(x[iH2],x[iH2]);

  domain->x2lamda(xM,xM);

  }

}

src/KSPACE/pppm_tip4p.h

+0 −1
Original line number Diff line number Diff line
@@ -39,7 +39,6 @@ class PPPMTIP4P : public PPPM { 


 private:

  void find_M(int, int &, int &, double *);

  void find_M_triclinic(int, int, int, double *);

};



}

src/USER-OMP/pppm_tip4p_omp.cpp

+94 −14
Original line number Diff line number Diff line
@@ -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,6 +744,83 @@ 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");



  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 ...



    double delx1 = xh1[0] - xo[0];

    double dely1 = xh1[1] - xo[1];

    double delz1 = xh1[2] - xo[2];



    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);



    // ... and convert M to lamda space for PPPM



    domain->x2lamda((double *)&xM,(double *)&xM);



  } else {



    // set iH1,iH2 to index of closest image to O



    iH1 = domain->closest_image(i,iH1);
@@ -761,6 +840,7 @@ void PPPMTIP4POMP::find_M_thr(int i, int &iH1, int &iH2, dbl3_t &xM) 
    xM.y = x[i].y + alpha * 0.5 * (dely1 + dely2);

    xM.z = x[i].z + alpha * 0.5 * (delz1 + delz2);

  }

}





/* ----------------------------------------------------------------------

CRA Git | Maintained and supported by SUSTech CRA and CCSE