Revised method for selecting the PME mesh spacing

platforms/cuda/src/CudaKernels.cpp

@@ -322,12 +322,9 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
                 method = EWALD;
             }
             else {
-                int gridSizeX = kmaxx*3;
-                int gridSizeY = kmaxy*3;
-                int gridSizeZ = kmaxz*3;
-                gridSizeX = ((gridSizeX+3)/4)*4;
-                gridSizeY = ((gridSizeY+3)/4)*4;
-                gridSizeZ = ((gridSizeZ+3)/4)*4;
+                int gridSizeX = -0.5*kmaxx*std::log(ewaldErrorTol);
+                int gridSizeY = -0.5*kmaxy*std::log(ewaldErrorTol);
+                int gridSizeZ = -0.5*kmaxz*std::log(ewaldErrorTol);
                 gpuSetPMEParameters(gpu, (float) alpha, gridSizeX, gridSizeY, gridSizeZ);
                 method = PARTICLE_MESH_EWALD;
             }

platforms/reference/src/ReferenceKernels.cpp

@@ -407,12 +407,19 @@ void ReferenceCalcNonbondedForceKernel::initialize(const System& system, const N
         kmax[0] = (int)std::ceil(-(mx/pi)*std::log(ewaldErrorTol));
         kmax[1] = (int)std::ceil(-(my/pi)*std::log(ewaldErrorTol));
         kmax[2] = (int)std::ceil(-(mz/pi)*std::log(ewaldErrorTol));
-        if (kmax[0]%2 == 0)
-            kmax[0]++;
-        if (kmax[1]%2 == 0)
-            kmax[1]++;
-        if (kmax[2]%2 == 0)
-            kmax[2]++;
+        if (nonbondedMethod == Ewald) {
+            if (kmax[0]%2 == 0)
+                kmax[0]++;
+            if (kmax[1]%2 == 0)
+                kmax[1]++;
+            if (kmax[2]%2 == 0)
+                kmax[2]++;
+        }
+        else {
+            gridSize[0] = -0.5*kmax[0]*std::log(ewaldErrorTol);
+            gridSize[1] = -0.5*kmax[1]*std::log(ewaldErrorTol);
+            gridSize[2] = -0.5*kmax[2]*std::log(ewaldErrorTol);
+        }
     }
     rfDielectric = (RealOpenMM)force.getReactionFieldDielectric();
 }
@@ -433,7 +440,7 @@ void ReferenceCalcNonbondedForceKernel::executeForces(ContextImpl& context) {
     if (ewald)
         clj.setUseEwald(ewaldAlpha, kmax[0], kmax[1], kmax[2]);
     if (pme)
-        clj.setUsePME(ewaldAlpha);
+        clj.setUsePME(ewaldAlpha, gridSize);
     clj.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, forceData, 0, 0);
     ReferenceBondForce refBondForce;
     ReferenceLJCoulomb14 nonbonded14;
@@ -459,7 +466,7 @@ double ReferenceCalcNonbondedForceKernel::executeEnergy(ContextImpl& context) {
     if (ewald)
         clj.setUseEwald(ewaldAlpha, kmax[0], kmax[1], kmax[2]);
     if (pme)
-        clj.setUsePME(ewaldAlpha);
+        clj.setUsePME(ewaldAlpha, gridSize);
     clj.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, forceData, 0, &energy);
     ReferenceBondForce refBondForce;
     ReferenceLJCoulomb14 nonbonded14;

platforms/reference/src/ReferenceKernels.h

@@ -269,7 +269,7 @@ private:
     int **exclusionArray, **bonded14IndexArray;
     RealOpenMM **particleParamArray, **bonded14ParamArray;
     RealOpenMM nonbondedCutoff, periodicBoxSize[3], rfDielectric, ewaldAlpha;
-    int kmax[3];
+    int kmax[3], gridSize[3];
     std::vector<std::set<int> > exclusions;
     NonbondedMethod nonbondedMethod;
     NeighborList* neighborList;

platforms/reference/src/SimTKReference/PME.cpp

@@ -627,7 +627,7 @@ int
 pme_init(pme_t *       ppme,
 		 RealOpenMM        ewaldcoeff,
 		 int           natoms,
-		 int           ngrid[3],
+		 const int           ngrid[3],
 		 int           pme_order,
          RealOpenMM        epsilon_r)
 {

platforms/reference/src/SimTKReference/PME.h

@@ -55,7 +55,7 @@ int
 pme_init(pme_t *       ppme,
 		 RealOpenMM        ewaldcoeff,
 		 int           natoms,
-		 int           ngrid[3],
+		 const int           ngrid[3],
 		 int           pme_order,
 	         RealOpenMM        epsilon_r);
 

platforms/reference/src/SimTKReference/ReferenceLJCoulombIxn.cpp

@@ -144,11 +144,15 @@ ReferenceLJCoulombIxn::~ReferenceLJCoulombIxn( ){
      Set the force to use Particle-Mesh Ewald (PME) summation.
 
      @param alpha  the Ewald separation parameter
+     @param gridSize the dimensions of the mesh
 
      --------------------------------------------------------------------------------------- */
 
-  void ReferenceLJCoulombIxn::setUsePME(RealOpenMM alpha) {
+  void ReferenceLJCoulombIxn::setUsePME(RealOpenMM alpha, int meshSize[3]) {
       alphaEwald = alpha;
+      meshDim[0] = meshSize[0];
+      meshDim[1] = meshSize[1];
+      meshDim[2] = meshSize[2];
       pme = true;
   }
 
@@ -273,19 +277,9 @@ int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** at
 
   if (pme) {
 	pme_t          pmedata; /* abstract handle for PME data */
-	int            ngrid[3];
 	RealOpenMM virial[3][3];
 
-	/* PME grid dimensions.
-	 * We typically want to set this as the spacing rather than absolute dimensions, but
-	 * to be able to reproduce results from other programs (e.g. Gromacs) we need to be
-	 * able to set exact grid dimenisions occasionally.
-	 */
-	ngrid[0] = 16;
-	ngrid[1] = 16;
-	ngrid[2] = 16;
-
-	pme_init(&pmedata,alphaEwald,numberOfAtoms,ngrid,4,1);
+	pme_init(&pmedata,alphaEwald,numberOfAtoms,meshDim,4,1);
 
 	pme_exec(pmedata,atomCoordinates,forces,atomParameters,periodicBoxSize,&recipEnergy,virial);
 

platforms/reference/src/SimTKReference/ReferenceLJCoulombIxn.h

@@ -44,6 +44,7 @@ class ReferenceLJCoulombIxn : public ReferencePairIxn {
       RealOpenMM krf, crf;
       RealOpenMM alphaEwald;
       int numRx, numRy, numRz;
+      int meshDim[3];
 
       // parameter indices
 
@@ -136,11 +137,12 @@ class ReferenceLJCoulombIxn : public ReferencePairIxn {
       
          Set the force to use Particle-Mesh Ewald (PME) summation.
       
-         @param alpha  the Ewald separation parameter
+         @param alpha    the Ewald separation parameter
+         @param gridSize the dimensions of the mesh
       
          --------------------------------------------------------------------------------------- */
       
-      void setUsePME(RealOpenMM alpha);
+      void setUsePME(RealOpenMM alpha, int meshSize[3]);
 
 
       /**---------------------------------------------------------------------------------------