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

src/KSPACE/pppm.cpp

@@ -795,7 +795,7 @@ void PPPM::deallocate()
 
 void PPPM::set_grid()
 {
-  // see JCP 109, pg. 7698 for derivation of coefficients
+  // see JCP 109, pg 7698 for derivation of coefficients
   // higher order coefficients may be computed if needed
 
   double **acons = memory->create_2d_double_array(8,7,"pppm:acons");
@@ -846,119 +846,91 @@ void PPPM::set_grid()
   // fluid-occupied volume used to estimate real-space error
   // zprd used rather than zprd_slab
 
+  double hx,hy,hz;
+
   if (!gewaldflag)
     g_ewald = sqrt(-log(precision*sqrt(natoms*cutoff*xprd*yprd*zprd) / 
 			(2.0*q2))) / cutoff;
 
   // set optimal nx_pppm,ny_pppm,nz_pppm based on order and precision
   // nz_pppm uses extended zprd_slab instead of zprd
-
+  // h = 1/g_ewald is upper bound on h such that h*g_ewald <= 1
-  double h,h1,h2,err,er1,er2,lpr;
+  // reduce it until precision target is met
-  int ncount;
 
   if (!gridflag) {
-    h = 1.0;
+    double err;
-    h1 = 2.0;
+    hx = hy = hz = 1/g_ewald;  
 
-    ncount = 0;
+    nx_pppm = static_cast<int> (xprd/hx + 1);
-    err = LARGE;
+    ny_pppm = static_cast<int> (yprd/hy + 1);
-    er1 = 0.0;
+    nz_pppm = static_cast<int> (zprd_slab/hz + 1);
-    while (fabs(err) > SMALL) {
+
-      lpr = rms(h,xprd,natoms,q2,acons);
+    err = rms(hx,xprd,natoms,q2,acons);
-      err = log(lpr) - log(precision);
+    while (err > precision) {
-      er2 = er1;
+      err = rms(hx,xprd,natoms,q2,acons);
-      er1 = err;
+      nx_pppm++;
-      h2 = h1;
+      hx = xprd/nx_pppm;
-      h1 = h;
-      if ((er1 - er2) == 0.0) h = h1 + er1;
-      else h = h1 + er1*(h2 - h1)/(er1 - er2);
-      ncount++;
-      if (ncount > LARGE) error->all("Cannot compute PPPM X grid spacing");
     }
-    nx_pppm = static_cast<int> (xprd/h + 1);
 
-    ncount = 0;
+    err = rms(hy,yprd,natoms,q2,acons);
-    err = LARGE;
+    while (err > precision) {
-    er1 = 0.0;
+      err = rms(hy,yprd,natoms,q2,acons);
-    while (fabs(err) > SMALL) {
+      ny_pppm++;
-      lpr = rms(h,yprd,natoms,q2,acons);
+      hy = yprd/ny_pppm;
-      err = log(lpr) - log(precision);
-      er2 = er1;
-      er1 = err;
-      h2 = h1;
-      h1 = h;
-      if ((er1 - er2) == 0.0) h = h1 + er1;
-      else h = h1 + er1*(h2 - h1)/(er1 - er2);
-      ncount++;
-      if (ncount > LARGE) error->all("Cannot compute PPPM Y grid spacing");
     }
-    ny_pppm = static_cast<int> (yprd/h + 1);
 
-    ncount = 0;
+    err = rms(hz,zprd_slab,natoms,q2,acons);
-    err = LARGE;
+    while (err > precision) {
-    er1 = 0.0;
+      err = rms(hz,zprd_slab,natoms,q2,acons);
-    while (fabs(err) > SMALL) {
+      nz_pppm++;
-      lpr = rms(h,zprd_slab,natoms,q2,acons);
+      hz = zprd_slab/nz_pppm;
-      err = log(lpr) - log(precision);
-      er2 = er1;
-      er1 = err;
-      h2 = h1;
-      h1 = h;
-      if ((er1 - er2) == 0.0) h = h1 + er1;
-      else h = h1 + er1*(h2 - h1)/(er1 - er2);
-      ncount++;
-      if (ncount > LARGE) error->all("Cannot compute PPPM Z grid spacing");
     }
-    nz_pppm = static_cast<int> (zprd_slab/h + 1);
-
   }
 
-  // convert grid into sizes that are factorable
+  // boost grid size until it is factorable
 
   while (!factorable(nx_pppm)) nx_pppm++;
   while (!factorable(ny_pppm)) ny_pppm++;
   while (!factorable(nz_pppm)) nz_pppm++;
 
-  // adjust g_ewald for new grid
+  // adjust g_ewald for new grid size
 
-  double dx = xprd/nx_pppm;
+  hx = xprd/nx_pppm;
-  double dy = yprd/ny_pppm;
+  hy = yprd/ny_pppm;
-  double dz = zprd_slab/nz_pppm;
+  hz = zprd_slab/nz_pppm;
 
   if (!gewaldflag) {
-    double gew1,gew2,lprx,lpry,lprz,spr;
+    double gew1,gew2,dgew,f,fmid,hmin,rtb;
+    int ncount;
 
-    gew1 = g_ewald + 0.01;
+    gew1 = 0.0;
+    g_ewald = gew1;
+    f = diffpr(hx,hy,hz,q2,acons);
+
+    hmin = MIN(hx,MIN(hy,hz));
+    gew2 = 10/hmin;
+    g_ewald = gew2;
+    fmid = diffpr(hx,hy,hz,q2,acons);
+
+    if (f*fmid >= 0.0) error->all("Cannot compute PPPM G");
+    rtb = f < 0.0 ? (dgew=gew2-gew1,gew1) : (dgew=gew1-gew2,gew2);
     ncount = 0;
-    err = LARGE;
+    while (fabs(dgew) > SMALL && fmid != 0.0) {
-    er1 = 0.0;
+      dgew *= 0.5;
-
+      g_ewald = rtb + dgew;
-    while (fabs(err) > SMALL) {
+      fmid = diffpr(hx,hy,hz,q2,acons);      
-      lprx = rms(dx,xprd,natoms,q2,acons);
+      if (fmid <= 0.0) rtb = g_ewald;
-      lpry = rms(dy,yprd,natoms,q2,acons);
-      lprz = rms(dz,zprd_slab,natoms,q2,acons);
-      lpr = sqrt(lprx*lprx + lpry*lpry + lprz*lprz) / sqrt(3.0);
-      spr = 2.0*q2 * exp(-g_ewald*g_ewald*cutoff*cutoff) / 
-	sqrt(natoms*cutoff*xprd*yprd*zprd_slab);
-      
-      err = log(lpr) - log(spr);
-      er2 = er1;
-      er1 = err;
-      gew2 = gew1;
-      gew1 = g_ewald;
-      if ((er1 - er2) == 0.0) g_ewald = gew1 + er1;
-      else g_ewald = gew1 + er1*(gew2 - gew1)/(er1 - er2);
       ncount++;
       if (ncount > LARGE) error->all("Cannot compute PPPM G");
     }
   }
 
-  // compute final RMS precision
+  // final RMS precision
 
-  double lprx = rms(dx,xprd,natoms,q2,acons);
+  double lprx = rms(hx,xprd,natoms,q2,acons);
-  double lpry = rms(dy,yprd,natoms,q2,acons);
+  double lpry = rms(hy,yprd,natoms,q2,acons);
-  double lprz = rms(dz,zprd_slab,natoms,q2,acons);
+  double lprz = rms(hz,zprd_slab,natoms,q2,acons);
-  lpr = sqrt(lprx*lprx + lpry*lpry + lprz*lprz) / sqrt(3.0);
+  double lpr = sqrt(lprx*lprx + lpry*lpry + lprz*lprz) / sqrt(3.0);
   double spr = 2.0*q2 * exp(-g_ewald*g_ewald*cutoff*cutoff) / 
     sqrt(natoms*cutoff*xprd*yprd*zprd_slab);
 
@@ -1019,6 +991,28 @@ double PPPM::rms(double h, double prd, double natoms,
   return value;
 }
 
+/* ----------------------------------------------------------------------
+   compute difference in real-space and kspace RMS precision
+------------------------------------------------------------------------- */
+
+double PPPM::diffpr(double hx, double hy, double hz, double q2, double **acons)
+{
+  double lprx,lpry,lprz,kspace_prec,real_prec;
+  double xprd = domain->xprd;
+  double yprd = domain->yprd;
+  double zprd = domain->zprd;
+  double natoms = atom->natoms;
+
+  lprx = rms(hx,xprd,natoms,q2,acons);
+  lpry = rms(hy,yprd,natoms,q2,acons);
+  lprz = rms(hz,zprd*slab_volfactor,natoms,q2,acons);
+  kspace_prec = sqrt(lprx*lprx + lpry*lpry + lprz*lprz) / sqrt(3.0);
+  real_prec = 2.0*q2 * exp(-g_ewald*g_ewald*cutoff*cutoff) / 
+   sqrt(natoms*cutoff*xprd*yprd*zprd);
+  double value = kspace_prec - real_prec;
+  return value;
+}
+
 /* ----------------------------------------------------------------------
    denominator for Hockney-Eastwood Green's function
      of x,y,z = sin(kx*deltax/2), etc

src/KSPACE/pppm.h

@@ -78,6 +78,7 @@ class PPPM : public KSpace {
   void deallocate();
   int factorable(int);
   double rms(double, double, double, double, double **);
+  double diffpr(double, double, double, double, double **);
   void compute_gf_denom();
   double gf_denom(double, double, double);
   virtual void particle_map();