Adding a method for estimating coulomb table precision and factoring that into overall kspace precision estimate.

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

src/KSPACE/ewald.cpp

@@ -176,6 +176,8 @@ void Ewald::init()
   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);
+  double tpr = estimate_table_accuracy(spr);
+  double accuracy = sqrt(lpr*lpr + spr*spr + tpr*tpr);
 
   // stats
 
@@ -183,18 +185,18 @@ void Ewald::init()
     if (screen) {
       fprintf(screen,"  G vector (1/distance) = %g\n",g_ewald);
       fprintf(screen,"  estimated absolute RMS force accuracy = %g\n",
-              MAX(lpr,spr));
+              accuracy);
       fprintf(screen,"  estimated relative force accuracy = %g\n",
-              MAX(lpr,spr)/two_charge_force);
+              accuracy/two_charge_force);
       fprintf(screen,"  KSpace vectors: actual max1d max3d = %d %d %d\n",
               kcount,kmax,kmax3d);
     }
     if (logfile) {
       fprintf(logfile,"  G vector (1/distnace) = %g\n",g_ewald);
       fprintf(logfile,"  estimated absolute RMS force accuracy = %g\n",
-              MAX(lpr,spr));
+              accuracy);
       fprintf(logfile,"  estimated relative force accuracy = %g\n",
-              MAX(lpr,spr)/two_charge_force);
+              accuracy/two_charge_force);
       fprintf(logfile,"  KSpace vectors: actual max1d max3d = %d %d %d\n",
               kcount,kmax,kmax3d);
     }

src/KSPACE/pppm.cpp

@@ -1091,6 +1091,8 @@ void PPPM::set_grid()
   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);
+  double tpr = estimate_table_accuracy(spr);
+  double accuracy = sqrt(lpr*lpr + spr*spr + tpr*tpr);
 
   // free local memory
 
@@ -1109,9 +1111,9 @@ void PPPM::set_grid()
       fprintf(screen,"  grid = %d %d %d\n",nx_pppm,ny_pppm,nz_pppm);
       fprintf(screen,"  stencil order = %d\n",order);
       fprintf(screen,"  estimated absolute RMS force accuracy = %g\n",
-              MAX(lpr,spr));
+              accuracy);
       fprintf(screen,"  estimated relative force accuracy = %g\n",
-              MAX(lpr,spr)/two_charge_force);
+              accuracy/two_charge_force);
       fprintf(screen,"  using %s precision FFTs\n",fft_prec);
     }
     if (logfile) {
@@ -1119,9 +1121,9 @@ void PPPM::set_grid()
       fprintf(logfile,"  grid = %d %d %d\n",nx_pppm,ny_pppm,nz_pppm);
       fprintf(logfile,"  stencil order = %d\n",order);
       fprintf(logfile,"  estimated absolute RMS force accuracy = %g\n",
-              MAX(lpr,spr));
+              accuracy);
       fprintf(logfile,"  estimated relative force accuracy = %g\n",
-              MAX(lpr,spr)/two_charge_force);
+              accuracy/two_charge_force);
       fprintf(logfile,"  using %s precision FFTs\n",fft_prec);
     }
   }

src/kspace.cpp

@@ -17,6 +17,7 @@
 #include "atom.h"
 #include "comm.h"
 #include "force.h"
+#include "pair.h"
 #include "memory.h"
 #include "error.h"
 #include "suffix.h"
@@ -124,6 +125,37 @@ void KSpace::ev_setup(int eflag, int vflag)
   }
 }
 
+/* ----------------------------------------------------------------------
+   estimate the accuracy of the short-range coulomb tables
+------------------------------------------------------------------------- */
+
+double KSpace::estimate_table_accuracy(double spr)
+{
+  double table_accuracy = 0.0;
+  int nctb = force->pair->ncoultablebits;
+  if (nctb) {
+    double empirical_precision[17];
+    empirical_precision[6] = 2.12E-04;
+    empirical_precision[7] = 4.97E-05;
+    empirical_precision[8] = 1.24E-05;
+    empirical_precision[9] = 3.04E-06;
+    empirical_precision[10] = 8.51E-07;
+    empirical_precision[11] = 1.85E-07;
+    empirical_precision[12] = 5.87E-08;
+    empirical_precision[13] = 2.81E-08;
+    empirical_precision[14] = 2.20E-08;
+    empirical_precision[15] = 2.13E-08;
+    empirical_precision[16] = 2.11E-08;
+    if (nctb <= 6) table_accuracy = empirical_precision[6];
+    else if (nctb <= 16) table_accuracy = empirical_precision[nctb];
+    else table_accuracy = empirical_precision[16];
+    table_accuracy *= two_charge_force;
+    if (table_accuracy > spr)
+      error->warning(FLERR,"For better accuracy use 'pair_modify table 0'");
+  }
+  return table_accuracy;
+}
+
 /* ----------------------------------------------------------------------
    modify parameters of the KSpace style
 ------------------------------------------------------------------------- */

src/kspace.h

@@ -70,6 +70,7 @@ class KSpace : protected Pointers {
   int maxeatom,maxvatom;
 
   void ev_setup(int, int);
+  double estimate_table_accuracy(double);
 };
 
 }
@@ -93,4 +94,9 @@ W: Kspace_modify slab param < 2.0 may cause unphysical behavior
 The kspace_modify slab parameter should be larger to insure periodic
 grids padded with empty space do not overlap.
 
+W: For better accuracy use 'pair_modify table 0'
+
+The user-specified force accuracy cannot be achieved unless the table
+feature is disabled by using 'pair_modify table 0'.
+
 */

src/pair.h

@@ -57,6 +57,8 @@ class Pair : protected Pointers {
   int evflag;                    // energy,virial settings
   int eflag_either,eflag_global,eflag_atom;
   int vflag_either,vflag_global,vflag_atom;
+  
+  int ncoultablebits;                  // size of Coulomb table
 
   int nextra;                    // # of extra quantities pair style calculates
   double *pvector;               // vector of extra pair quantities
@@ -146,7 +148,6 @@ class Pair : protected Pointers {
 
                                        // pair_modify settings
   int offset_flag,mix_flag;            // flags for offset and mixing
-  int ncoultablebits;                  // size of Coulomb table
   double tabinner;                     // inner cutoff for Coulomb table
 
   // custom data type for accessing Coulomb tables