git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@368 f3b2605a-c512-4ea7-a41b-209d697bcdaa

src/KSPACE/ewald.cpp

@@ -30,6 +30,8 @@
 
 using namespace LAMMPS_NS;
 
+#define SMALL 0.00001
+
 #define MIN(a,b) ((a) < (b) ? (a) : (b))
 #define MAX(a,b) ((a) > (b) ? (a) : (b))
 
@@ -78,8 +80,12 @@ void Ewald::init()
 
   // error check
 
+  if (domain->triclinic) error->all("Cannot use Ewald with triclinic box");
   if (force->dimension == 2) error->all("Cannot use Ewald with 2d simulation");
 
+  if (atom->q == NULL)
+    error->all("Must use charged atom style with kspace style");
+
   if (slabflag == 0 && domain->nonperiodic > 0)
     error->all("Cannot use nonperiodic boundaries with Ewald");
   if (slabflag == 1) {
@@ -99,20 +105,28 @@ void Ewald::init()
   double cutoff;
   pair->extract_long(&cutoff);
 
-  // compute qsum & qsqsum
+  // compute qsum & qsqsum and warn if not charge-neutral
+
+  qsum = qsqsum = 0.0;
+  for (int i = 0; i < atom->nlocal; i++) {
+    qsum += atom->q[i];
+    qsqsum += atom->q[i]*atom->q[i];
+  }
 
   double tmp;
-
-  qsum = 0.0;
-  for (int i = 0; i < atom->nlocal; i++) qsum += atom->q[i];
   MPI_Allreduce(&qsum,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
   qsum = tmp;
-
-  qsqsum = 0.0;
-  for (int i = 0; i < atom->nlocal; i++) qsqsum += atom->q[i]*atom->q[i];
   MPI_Allreduce(&qsqsum,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
   qsqsum = tmp;
 
+  if (qsqsum == 0.0)
+    error->all("Cannot use kspace solver on system with no charge");
+  if (fabs(qsum) > SMALL && comm->me == 0) {
+    char str[128];
+    sprintf(str,"System is not charge neutral, net charge = %g",qsum);
+    error->warning(str);
+  }
+
   // setup K-space resolution
 
   g_ewald = (1.35 - 0.15*log(precision))/cutoff;

src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp

@@ -250,17 +250,17 @@ void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag)
 	      delx1 = x[i][0] - x2[0];
 	      dely1 = x[i][1] - x2[1];
 	      delz1 = x[i][2] - x2[2];
-	      domain->minimum_image(&delx1,&dely1,&delz1);
+	      domain->minimum_image(delx1,dely1,delz1);
 
 	      delx2 = x[iH1][0] - x2[0];
 	      dely2 = x[iH1][1] - x2[1];
 	      delz2 = x[iH1][2] - x2[2];
-	      domain->minimum_image(&delx2,&dely2,&delz2);
+	      domain->minimum_image(delx2,dely2,delz2);
 
 	      delx3 = x[iH2][0] - x2[0];
 	      dely3 = x[iH2][1] - x2[1];
 	      delz3 = x[iH2][2] - x2[2];
-	      domain->minimum_image(&delx3,&dely3,&delz3);
+	      domain->minimum_image(delx3,dely3,delz3);
 
 	      tvirial[0] += 0.5 * (delx1 * fO[0] + (delx2 + delx3) * fH[0]);
 	      tvirial[1] += 0.5 * (dely1 * fO[1] + (dely2 + dely3) * fH[1]);
@@ -329,17 +329,17 @@ void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag)
 	      delx1 = x[j][0] - x1[0];
 	      dely1 = x[j][1] - x1[1];
 	      delz1 = x[j][2] - x1[2];
-	      domain->minimum_image(&delx1,&dely1,&delz1);
+	      domain->minimum_image(delx1,dely1,delz1);
 
 	      delx2 = x[jH1][0] - x1[0];
 	      dely2 = x[jH1][1] - x1[1];
 	      delz2 = x[jH1][2] - x1[2];
-	      domain->minimum_image(&delx2,&dely2,&delz2);
+	      domain->minimum_image(delx2,dely2,delz2);
 
 	      delx3 = x[jH2][0] - x1[0];
 	      dely3 = x[jH2][1] - x1[1];
 	      delz3 = x[jH2][2] - x1[2];
-	      domain->minimum_image(&delx3,&dely3,&delz3);
+	      domain->minimum_image(delx3,dely3,delz3);
 
 	      tvirial[0] += 0.5 * (delx1 * fO[0] + (delx2 + delx3) * fH[0]);
 	      tvirial[1] += 0.5 * (dely1 * fO[1] + (dely2 + dely3) * fH[1]);
@@ -518,12 +518,12 @@ void PairLJCutCoulLongTIP4P::find_M(int i, int &iH1, int &iH2, double *xM)
   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);
+  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);
+  domain->minimum_image(delx2,dely2,delz2);
 
   xM[0] = x[i][0] + alpha * (delx1 + delx2);
   xM[1] = x[i][1] + alpha * (dely1 + dely2);

src/KSPACE/pppm.cpp

@@ -36,15 +36,15 @@
 
 using namespace LAMMPS_NS;
 
-#define MIN(a,b) ((a) < (b) ? (a) : (b))
-#define MAX(a,b) ((a) > (b) ? (a) : (b))
-
 #define MAXORDER 7
 #define OFFSET 4096
 #define SMALL 0.00001
 #define LARGE 10000.0
 #define EPS_HOC 1.0e-7
 
+#define MIN(a,b) ((a) < (b) ? (a) : (b))
+#define MAX(a,b) ((a) > (b) ? (a) : (b))
+
 /* ---------------------------------------------------------------------- */
 
 PPPM::PPPM(LAMMPS *lmp, int narg, char **arg) : KSpace(lmp, narg, arg)
@@ -105,8 +105,13 @@ void PPPM::init()
 
   // error check
 
+  if (domain->triclinic)
+    error->all("Cannot (yet) use PPPM with triclinic box");
   if (force->dimension == 2) error->all("Cannot use PPPM with 2d simulation");
 
+  if (atom->q == NULL)
+    error->all("Must use charged atom style with kspace style");
+
   if (slabflag == 0 && domain->nonperiodic > 0)
     error->all("Cannot use nonperiodic boundaries with PPPM");
   if (slabflag == 1) {
@@ -169,6 +174,8 @@ void PPPM::init()
   MPI_Allreduce(&qsqsum,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
   qsqsum = tmp;
 
+  if (qsqsum == 0.0)
+    error->all("Cannot use kspace solver on system with no charge");
   if (fabs(qsum) > SMALL && me == 0) {
     char str[128];
     sprintf(str,"System is not charge neutral, net charge = %g",qsum);
@@ -214,36 +221,60 @@ void PPPM::init()
   // effectively nlo_in,nhi_in + ghost cells
   // nlo,nhi = global coords of grid pt to "lower left" of smallest/largest
   //           position a particle in my box can be at
-  // particle position bound = subbox + skin/2.0 + qdist
+  // dist[3] = particle position bound = subbox + skin/2.0 + qdist
   //   qdist = offset due to TIP4P fictitious charge
+  //   convert to triclinic if necessary
   // nlo_out,nhi_out = nlo,nhi + stencil size for particle mapping
   // for slab PPPM, assign z grid as if it were not extended
 
-  int nlo,nhi;
+  triclinic = domain->triclinic;
-  double cuthalf = 0.5 * neighbor->skin + qdist;
+  double *prd,*sublo,*subhi;
 
-  double xprd = domain->xprd;
+  if (triclinic == 0) {
-  double yprd = domain->yprd;
+    prd = domain->prd;
-  double zprd = domain->zprd;
+    boxlo = domain->boxlo;
+    sublo = domain->sublo;
+    subhi = domain->subhi;
+  } else {
+    prd = domain->prd_lamda;
+    boxlo = domain->boxlo_lamda;
+    sublo = domain->sublo_lamda;
+    subhi = domain->subhi_lamda;
+  }
+
+  double xprd = prd[0];
+  double yprd = prd[1];
+  double zprd = prd[2];
   double zprd_slab = zprd*slab_volfactor;
 
-  nlo = static_cast<int> ((domain->subxlo-cuthalf-domain->boxxlo) * 
+  double dist[3];
+  double cuthalf = 0.5 * neighbor->skin + qdist;
+  if (triclinic == 0) dist[0] = dist[1] = dist[2] = cuthalf;
+  else {
+    dist[0] = cuthalf/domain->prd[0];
+    dist[1] = cuthalf/domain->prd[1];
+    dist[2] = cuthalf/domain->prd[2];
+  }
+
+  int nlo,nhi;
+  
+  nlo = static_cast<int> ((sublo[0]-dist[0]-boxlo[0]) * 
 			  nx_pppm/xprd + shift) - OFFSET;
-  nhi = static_cast<int> ((domain->subxhi+cuthalf-domain->boxxlo) * 
+  nhi = static_cast<int> ((subhi[0]+dist[0]-boxlo[0]) * 
 			  nx_pppm/xprd + shift) - OFFSET;
   nxlo_out = nlo + nlower;
   nxhi_out = nhi + nupper;
 
-  nlo = static_cast<int> ((domain->subylo-cuthalf-domain->boxylo) * 
+  nlo = static_cast<int> ((sublo[1]-dist[1]-boxlo[1]) * 
 			  ny_pppm/yprd + shift) - OFFSET;
-  nhi = static_cast<int> ((domain->subyhi+cuthalf-domain->boxylo) * 
+  nhi = static_cast<int> ((subhi[1]+dist[1]-boxlo[1]) * 
 			  ny_pppm/yprd + shift) - OFFSET;
   nylo_out = nlo + nlower;
   nyhi_out = nhi + nupper;
 
-  nlo = static_cast<int> ((domain->subzlo-cuthalf-domain->boxzlo) * 
+  nlo = static_cast<int> ((sublo[2]-dist[2]-boxlo[2]) * 
 			  nz_pppm/zprd_slab + shift) - OFFSET;
-  nhi = static_cast<int> ((domain->subzhi+cuthalf-domain->boxzlo) * 
+  nhi = static_cast<int> ((subhi[2]+dist[2]-boxlo[2]) * 
 			  nz_pppm/zprd_slab + shift) - OFFSET;
   nzlo_out = nlo + nlower;
   nzhi_out = nhi + nupper;
@@ -426,18 +457,19 @@ void PPPM::init()
 void PPPM::setup()
 {
   int i,j,k,l,m,n;
+  double *prd;
 
   // volume-dependent factors
+  // adjust z dimension for 2d slab PPPM
+  // z dimension for 3d PPPM is zprd since slab_volfactor = 1.0
 
-  double xprd = domain->xprd;
+  if (triclinic == 0) prd = domain->prd;
-  double yprd = domain->yprd;
+  else prd = domain->prd_lamda;
-  double zprd = domain->zprd;
-    
-  // adjustment of z dimension for 2d slab PPPM
-  // 3d PPPM just uses zprd since slab_volfactor = 1.0
 
+  double xprd = prd[0];
+  double yprd = prd[1];
+  double zprd = prd[2];
   double zprd_slab = zprd*slab_volfactor;
-
   volume = xprd * yprd * zprd_slab;
     
   delxinv = nx_pppm/xprd;
@@ -579,6 +611,14 @@ void PPPM::compute(int eflag, int vflag)
 {
   int i;
 
+  // convert atoms from box to lamda coords
+  
+  if (triclinic == 0) boxlo = domain->boxlo;
+  else {
+    boxlo = domain->boxlo_lamda;
+    domain->x2lamda(atom->nlocal);
+  }
+
   // extend size of nlocal-dependent arrays if necessary
 
   if (atom->nlocal > nmax) {
@@ -641,6 +681,10 @@ void PPPM::compute(int eflag, int vflag)
   // 2d slab correction
 
   if (slabflag) slabcorr(eflag);
+
+  // convert atoms back from lamda to box coords
+  
+  if (triclinic) domain->lamda2x(atom->nlocal);
 }
 
 /* ----------------------------------------------------------------------
@@ -781,7 +825,8 @@ void PPPM::set_grid()
   double q2 = qsqsum / force->dielectric;
   double natoms = atom->natoms;
 
-  // adjustment of z dimension for 2d slab PPPM
+  // use xprd,yprd,zprd even if triclinic so grid size is the same
+  // adjust z dimension for 2d slab PPPM
   // 3d PPPM just uses zprd since slab_volfactor = 1.0
 
   double xprd = domain->xprd;
@@ -1387,9 +1432,6 @@ void PPPM::particle_map()
 
   double **x = atom->x;
   int nlocal = atom->nlocal;
-  double boxxlo = domain->boxxlo;
-  double boxylo = domain->boxylo;
-  double boxzlo = domain->boxzlo;
 
   int flag = 0;
   for (int i = 0; i < nlocal; i++) {
@@ -1398,9 +1440,9 @@ void PPPM::particle_map()
     // current particle coord can be outside global and local box
     // add/subtract OFFSET to avoid int(-0.75) = 0 when want it to be -1
 
-    nx = static_cast<int> ((x[i][0]-boxxlo)*delxinv+shift) - OFFSET;
+    nx = static_cast<int> ((x[i][0]-boxlo[0])*delxinv+shift) - OFFSET;
-    ny = static_cast<int> ((x[i][1]-boxylo)*delyinv+shift) - OFFSET;
+    ny = static_cast<int> ((x[i][1]-boxlo[1])*delyinv+shift) - OFFSET;
-    nz = static_cast<int> ((x[i][2]-boxzlo)*delzinv+shift) - OFFSET;
+    nz = static_cast<int> ((x[i][2]-boxlo[2])*delzinv+shift) - OFFSET;
 
     part2grid[i][0] = nx;
     part2grid[i][1] = ny;
@@ -1443,18 +1485,15 @@ void PPPM::make_rho()
   double *q = atom->q;
   double **x = atom->x;
   int nlocal = atom->nlocal;
-  double boxxlo = domain->boxxlo;
-  double boxylo = domain->boxylo;
-  double boxzlo = domain->boxzlo;
 
   for (int i = 0; i < nlocal; i++) {
 
     nx = part2grid[i][0];
     ny = part2grid[i][1];
     nz = part2grid[i][2];
-    dx = nx+shiftone - (x[i][0]-boxxlo)*delxinv;
+    dx = nx+shiftone - (x[i][0]-boxlo[0])*delxinv;
-    dy = ny+shiftone - (x[i][1]-boxylo)*delyinv;
+    dy = ny+shiftone - (x[i][1]-boxlo[1])*delyinv;
-    dz = nz+shiftone - (x[i][2]-boxzlo)*delzinv;
+    dz = nz+shiftone - (x[i][2]-boxlo[2])*delzinv;
 
     compute_rho1d(dx,dy,dz);
 
@@ -1614,18 +1653,15 @@ void PPPM::fieldforce()
   double **x = atom->x;
   double **f = atom->f;
   int nlocal = atom->nlocal;
-  double boxxlo = domain->boxxlo;
-  double boxylo = domain->boxylo;
-  double boxzlo = domain->boxzlo;
 
   for (i = 0; i < nlocal; i++) {
 
     nx = part2grid[i][0];
     ny = part2grid[i][1];
     nz = part2grid[i][2];
-    dx = nx+shiftone - (x[i][0]-boxxlo)*delxinv;
+    dx = nx+shiftone - (x[i][0]-boxlo[0])*delxinv;
-    dy = ny+shiftone - (x[i][1]-boxylo)*delyinv;
+    dy = ny+shiftone - (x[i][1]-boxlo[1])*delyinv;
-    dz = nz+shiftone - (x[i][2]-boxzlo)*delzinv;
+    dz = nz+shiftone - (x[i][2]-boxlo[2])*delzinv;
 
     compute_rho1d(dx,dy,dz);
 

src/KSPACE/pppm.h

@@ -66,6 +66,8 @@ class PPPM : public KSpace {
   int **part2grid;             // storage for particle -> grid mapping
   int nmax;
 
+  int triclinic;               // domain settings, orthog or triclinic
+  double *boxlo;
                                // TIP4P settings
   int typeH,typeO;             // atom types of TIP4P water H and O atoms
   double qdist;                // distance from O site to negative charge

src/KSPACE/pppm_tip4p.cpp

@@ -45,9 +45,6 @@ void PPPMTIP4P::particle_map()
   int *type = atom->type;
   double **x = atom->x;
   int nlocal = atom->nlocal;
-  double boxxlo = domain->boxxlo;
-  double boxylo = domain->boxylo;
-  double boxzlo = domain->boxzlo;
 
   int flag = 0;
   for (int i = 0; i < nlocal; i++) {
@@ -60,9 +57,9 @@ void PPPMTIP4P::particle_map()
     // current particle coord can be outside global and local box
     // add/subtract OFFSET to avoid int(-0.75) = 0 when want it to be -1
 
-    nx = static_cast<int> ((xi[0]-boxxlo)*delxinv+shift) - OFFSET;
+    nx = static_cast<int> ((xi[0]-boxlo[0])*delxinv+shift) - OFFSET;
-    ny = static_cast<int> ((xi[1]-boxylo)*delyinv+shift) - OFFSET;
+    ny = static_cast<int> ((xi[1]-boxlo[1])*delyinv+shift) - OFFSET;
-    nz = static_cast<int> ((xi[2]-boxzlo)*delzinv+shift) - OFFSET;
+    nz = static_cast<int> ((xi[2]-boxlo[2])*delzinv+shift) - OFFSET;
 
     part2grid[i][0] = nx;
     part2grid[i][1] = ny;
@@ -107,9 +104,6 @@ void PPPMTIP4P::make_rho()
   double *q = atom->q;
   double **x = atom->x;
   int nlocal = atom->nlocal;
-  double boxxlo = domain->boxxlo;
-  double boxylo = domain->boxylo;
-  double boxzlo = domain->boxzlo;
 
   for (int i = 0; i < nlocal; i++) {
     if (type[i] == typeO) {
@@ -120,9 +114,9 @@ void PPPMTIP4P::make_rho()
     nx = part2grid[i][0];
     ny = part2grid[i][1];
     nz = part2grid[i][2];
-    dx = nx+shiftone - (xi[0]-boxxlo)*delxinv;
+    dx = nx+shiftone - (xi[0]-boxlo[0])*delxinv;
-    dy = ny+shiftone - (xi[1]-boxylo)*delyinv;
+    dy = ny+shiftone - (xi[1]-boxlo[1])*delyinv;
-    dz = nz+shiftone - (xi[2]-boxzlo)*delzinv;
+    dz = nz+shiftone - (xi[2]-boxlo[2])*delzinv;
 
     compute_rho1d(dx,dy,dz);
 
@@ -167,9 +161,6 @@ void PPPMTIP4P::fieldforce()
   double **f = atom->f;
   int *type = atom->type;
   int nlocal = atom->nlocal;
-  double boxxlo = domain->boxxlo;
-  double boxylo = domain->boxylo;
-  double boxzlo = domain->boxzlo;
 
   for (i = 0; i < nlocal; i++) {
     if (type[i] == typeO) {
@@ -180,9 +171,9 @@ void PPPMTIP4P::fieldforce()
     nx = part2grid[i][0];
     ny = part2grid[i][1];
     nz = part2grid[i][2];
-    dx = nx+shiftone - (xi[0]-boxxlo)*delxinv;
+    dx = nx+shiftone - (xi[0]-boxlo[0])*delxinv;
-    dy = ny+shiftone - (xi[1]-boxylo)*delyinv;
+    dy = ny+shiftone - (xi[1]-boxlo[1])*delyinv;
-    dz = nz+shiftone - (xi[2]-boxzlo)*delzinv;
+    dz = nz+shiftone - (xi[2]-boxlo[2])*delzinv;
 
     compute_rho1d(dx,dy,dz);
 
@@ -251,12 +242,12 @@ void PPPMTIP4P::find_M(int i, int &iH1, int &iH2, double *xM)
   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);
+  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);
+  domain->minimum_image(delx2,dely2,delz2);
 
   xM[0] = x[i][0] + alpha * (delx1 + delx2);
   xM[1] = x[i][1] + alpha * (dely1 + dely2);