Updated tests, to compare to Gromacs

o Added C12 shift term, to match GMX.
o Added tests to compare energies and forces directly to GMX, with and without exclusions.
o Modified set/get LJPME parameters to look similar to electrostatic equivalents.

olla/include/openmm/kernels.h

@@ -600,12 +600,12 @@ public:
     /**
      * Get the parameters being used for the dispersion terms in LJPME.
      *
-     * @param dalpha   the separation parameter
-     * @param dnx      the number of grid points along the X axis
-     * @param dny      the number of grid points along the Y axis
-     * @param dnz      the number of grid points along the Z axis
+     * @param alpha   the separation parameter
+     * @param nx      the number of grid points along the X axis
+     * @param ny      the number of grid points along the Y axis
+     * @param nz      the number of grid points along the Z axis
      */
-    virtual void getLJPMEParameters(double& dalpha, int& dnx, int& dny, int& dnz) const = 0;
+    virtual void getLJPMEParameters(double& alpha, int& nx, int& ny, int& nz) const = 0;
 };
 
 /**

openmmapi/include/openmm/NonbondedForce.h

@@ -216,12 +216,12 @@ public:
      * Get the parameters to use for dispersion term in LJ-PME calculations.  If alpha is 0 (the default),
      * these parameters are ignored and instead their values are chosen based on the Ewald error tolerance.
      *
-     * @param[out] dalpha   the separation parameter
-     * @param[out] dnx      the number of dispersion grid points along the X axis
-     * @param[out] dny      the number of dispersion grid points along the Y axis
-     * @param[out] dnz      the number of dispersion grid points along the Z axis
+     * @param[out] alpha   the separation parameter
+     * @param[out] nx      the number of dispersion grid points along the X axis
+     * @param[out] ny      the number of dispersion grid points along the Y axis
+     * @param[out] nz      the number of dispersion grid points along the Z axis
      */
-    void getLJPMEParameters(double& dalpha, int& dnx, int& dny, int& dnz) const;
+    void getLJPMEParameters(double& alpha, int& nx, int& ny, int& nz) const;
     /**
      * Set the parameters to use for PME calculations.  If alpha is 0 (the default), these parameters are
      * ignored and instead their values are chosen based on the Ewald error tolerance.
@@ -236,12 +236,12 @@ public:
      * Set the parameters to use for the dispersion term in LJPME calculations.  If alpha is 0 (the default),
      * these parameters are ignored and instead their values are chosen based on the Ewald error tolerance.
      *
-     * @param dalpha   the separation parameter
-     * @param dnx      the number of grid points along the X axis
-     * @param dny      the number of grid points along the Y axis
-     * @param dnz      the number of grid points along the Z axis
+     * @param alpha   the separation parameter
+     * @param nx      the number of grid points along the X axis
+     * @param ny      the number of grid points along the Y axis
+     * @param nz      the number of grid points along the Z axis
      */
-    void setLJPMEParameters(double dalpha, int dnx, int dny, int dnz);
+    void setLJPMEParameters(double alpha, int nx, int ny, int nz);
     /**
      * Get the parameters being used for PME in a particular Context.  Because some platforms have restrictions
      * on the allowed grid sizes, the values that are actually used may be slightly different from those
@@ -262,12 +262,12 @@ public:
      * See the manual for details.
      *
      * @param context      the Context for which to get the parameters
-     * @param[out] dalpha   the separation parameter
-     * @param[out] dnx      the number of grid points along the X axis
-     * @param[out] dny      the number of grid points along the Y axis
-     * @param[out] dnz      the number of grid points along the Z axis
+     * @param[out] alpha   the separation parameter
+     * @param[out] nx      the number of grid points along the X axis
+     * @param[out] ny      the number of grid points along the Y axis
+     * @param[out] nz      the number of grid points along the Z axis
      */
-    void getLJPMEParametersInContext(const Context& context, double& dalpha, int& dnx, int& dny, int& dnz) const;
+    void getLJPMEParametersInContext(const Context& context, double& alpha, int& nx, int& ny, int& nz) const;
     /**
      * Add the nonbonded force parameters for a particle.  This should be called once for each particle
      * in the System.  When it is called for the i'th time, it specifies the parameters for the i'th particle.

openmmapi/src/NonbondedForce.cpp

@@ -106,11 +106,11 @@ void NonbondedForce::getPMEParameters(double& alpha, int& nx, int& ny, int& nz)
     nz = this->nz;
 }
 
-void NonbondedForce::getLJPMEParameters(double& dalpha, int& dnx, int& dny, int& dnz) const {
-    dalpha = this->dalpha;
-    dnx = this->dnx;
-    dny = this->dny;
-    dnz = this->dnz;
+void NonbondedForce::getLJPMEParameters(double& alpha, int& nx, int& ny, int& nz) const {
+    alpha = this->dalpha;
+    nx = this->dnx;
+    ny = this->dny;
+    nz = this->dnz;
 }
 
 void NonbondedForce::setPMEParameters(double alpha, int nx, int ny, int nz) {
@@ -120,19 +120,19 @@ void NonbondedForce::setPMEParameters(double alpha, int nx, int ny, int nz) {
     this->nz = nz;
 }
 
-void NonbondedForce::setLJPMEParameters(double dalpha, int dnx, int dny, int dnz) {
-    this->dalpha = dalpha;
-    this->dnx = dnx;
-    this->dny = dny;
-    this->dnz = dnz;
+void NonbondedForce::setLJPMEParameters(double alpha, int nx, int ny, int nz) {
+    this->dalpha = alpha;
+    this->dnx = nx;
+    this->dny = ny;
+    this->dnz = nz;
 }
 
 void NonbondedForce::getPMEParametersInContext(const Context& context, double& alpha, int& nx, int& ny, int& nz) const {
     dynamic_cast<const NonbondedForceImpl&>(getImplInContext(context)).getPMEParameters(alpha, nx, ny, nz);
 }
 
-void NonbondedForce::getLJPMEParametersInContext(const Context& context, double& dalpha, int& dnx, int& dny, int& dnz) const {
-    dynamic_cast<const NonbondedForceImpl&>(getImplInContext(context)).getLJPMEParameters(dalpha, dnx, dny, dnz);
+void NonbondedForce::getLJPMEParametersInContext(const Context& context, double& alpha, int& nx, int& ny, int& nz) const {
+    dynamic_cast<const NonbondedForceImpl&>(getImplInContext(context)).getLJPMEParameters(alpha, nx, ny, nz);
 }
 
 int NonbondedForce::addParticle(double charge, double sigma, double epsilon) {

platforms/cpu/include/CpuKernels.h

@@ -261,12 +261,12 @@ public:
     /**
      * Get the parameters being used for the dispersion term in LJPME.
      *
-     * @param dalpha   the separation parameter
-     * @param dnx      the number of grid points along the X axis
-     * @param dny      the number of grid points along the Y axis
-     * @param dnz      the number of grid points along the Z axis
+     * @param alpha   the separation parameter
+     * @param nx      the number of grid points along the X axis
+     * @param ny      the number of grid points along the Y axis
+     * @param nz      the number of grid points along the Z axis
      */
-    void getLJPMEParameters(double& dalpha, int& dnx, int& dny, int& dnz) const;
+    void getLJPMEParameters(double& alpha, int& nx, int& ny, int& nz) const;
 private:
     class PmeIO;
     CpuPlatform::PlatformData& data;

platforms/cpu/src/CpuKernels.cpp

@@ -751,16 +751,16 @@ void CpuCalcNonbondedForceKernel::getPMEParameters(double& alpha, int& nx, int&
     }
 }
 
-void CpuCalcNonbondedForceKernel::getLJPMEParameters(double& dalpha, int& dnx, int& dny, int& dnz) const {
+void CpuCalcNonbondedForceKernel::getLJPMEParameters(double& alpha, int& nx, int& ny, int& nz) const {
     if (nonbondedMethod != LJPME)
         throw OpenMMException("getPMEParametersInContext: This Context is not using PME");
     if (useOptimizedPme)
-        optimizedDispersionPme.getAs<const CalcPmeReciprocalForceKernel>().getPMEParameters(dalpha, dnx, dny, dnz);
+        optimizedDispersionPme.getAs<const CalcPmeReciprocalForceKernel>().getPMEParameters(alpha, nx, ny, nz);
     else {
-        dalpha = ewaldDispersionAlpha;
-        dnx = dispersionGridSize[0];
-        dny = dispersionGridSize[1];
-        dnz = dispersionGridSize[2];
+        alpha = ewaldDispersionAlpha;
+        nx = dispersionGridSize[0];
+        ny = dispersionGridSize[1];
+        nz = dispersionGridSize[2];
     }
 }
 

platforms/opencl/include/OpenCLKernels.h

@@ -622,7 +622,7 @@ public:
      * @param ny      the number of grid points along the Y axis
      * @param nz      the number of grid points along the Z axis
      */
-    void getLJPMEParameters(double& dalpha, int& dnx, int& dny, int& dnz) const;
+    void getLJPMEParameters(double& alpha, int& nx, int& ny, int& nz) const;
 private:
     class SortTrait : public OpenCLSort::SortTrait {
         int getDataSize() const {return 8;}

platforms/opencl/include/OpenCLParallelKernels.h

@@ -441,12 +441,12 @@ public:
     /**
      * Get the parameters being used for the dispersion term in LJPME.
      *
-     * @param dalpha   the separation parameter
-     * @param dnx      the number of grid points along the X axis
-     * @param dny      the number of grid points along the Y axis
-     * @param dnz      the number of grid points along the Z axis
+     * @param alpha   the separation parameter
+     * @param nx      the number of grid points along the X axis
+     * @param ny      the number of grid points along the Y axis
+     * @param nz      the number of grid points along the Z axis
      */
-    void getLJPMEParameters(double& dalpha, int& dnx, int& dny, int& dnz) const;
+    void getLJPMEParameters(double& alpha, int& nx, int& ny, int& nz) const;
 private:
     class Task;
     OpenCLPlatform::PlatformData& data;

platforms/opencl/src/OpenCLKernels.cpp

@@ -2205,16 +2205,16 @@ void OpenCLCalcNonbondedForceKernel::getPMEParameters(double& alpha, int& nx, in
     }
 }
 
-void OpenCLCalcNonbondedForceKernel::getLJPMEParameters(double& dalpha, int& dnx, int& dny, int& dnz) const {
+void OpenCLCalcNonbondedForceKernel::getLJPMEParameters(double& alpha, int& nx, int& ny, int& nz) const {
     if (nonbondedMethod != LJPME)
         throw OpenMMException("getPMEParametersInContext: This Context is not using PME");
     if (cl.getPlatformData().useCpuPme)
-        cpuPme.getAs<CalcPmeReciprocalForceKernel>().getPMEParameters(dalpha, dnx, dny, dnz);
+        cpuPme.getAs<CalcPmeReciprocalForceKernel>().getPMEParameters(alpha, nx, ny, nz);
     else {
-        dalpha = this->alpha;
-        dnx = gridSizeX;
-        dny = gridSizeY;
-        dnz = gridSizeZ;
+        alpha = this->alpha;
+        nx = gridSizeX;
+        ny = gridSizeY;
+        nz = gridSizeZ;
     }
 }
 

platforms/opencl/src/OpenCLParallelKernels.cpp

@@ -583,8 +583,8 @@ void OpenCLParallelCalcNonbondedForceKernel::getPMEParameters(double& alpha, int
     dynamic_cast<const OpenCLCalcNonbondedForceKernel&>(kernels[0].getImpl()).getPMEParameters(alpha, nx, ny, nz);
 }
 
-void OpenCLParallelCalcNonbondedForceKernel::getLJPMEParameters(double& dalpha, int& dnx, int& dny, int& dnz) const {
-    dynamic_cast<const OpenCLCalcNonbondedForceKernel&>(kernels[0].getImpl()).getLJPMEParameters(dalpha, dnx, dny, dnz);
+void OpenCLParallelCalcNonbondedForceKernel::getLJPMEParameters(double& alpha, int& nx, int& ny, int& nz) const {
+    dynamic_cast<const OpenCLCalcNonbondedForceKernel&>(kernels[0].getImpl()).getLJPMEParameters(alpha, nx, ny, nz);
 }
 
 class OpenCLParallelCalcCustomNonbondedForceKernel::Task : public OpenCLContext::WorkTask {

platforms/reference/include/ReferenceKernels.h

@@ -607,12 +607,12 @@ public:
     /**
      * Get the dispersion parameters being used for the dispersion term in LJPME.
      *
-     * @param dalpha   the separation parameter
-     * @param dnx      the number of grid points along the X axis
-     * @param dny      the number of grid points along the Y axis
-     * @param dnz      the number of grid points along the Z axis
+     * @param alpha   the separation parameter
+     * @param nx      the number of grid points along the X axis
+     * @param ny      the number of grid points along the Y axis
+     * @param nz      the number of grid points along the Z axis
      */
-    void getLJPMEParameters(double& dalpha, int& dnx, int& dny, int& dnz) const;
+    void getLJPMEParameters(double& alpha, int& nx, int& ny, int& nz) const;
 private:
     int numParticles, num14;
     int **bonded14IndexArray;

platforms/reference/src/ReferenceKernels.cpp

@@ -1080,13 +1080,13 @@ void ReferenceCalcNonbondedForceKernel::getPMEParameters(double& alpha, int& nx,
     nz = gridSize[2];
 }
 
-void ReferenceCalcNonbondedForceKernel::getLJPMEParameters(double& dalpha, int& dnx, int& dny, int& dnz) const {
+void ReferenceCalcNonbondedForceKernel::getLJPMEParameters(double& alpha, int& nx, int& ny, int& nz) const {
     if (nonbondedMethod != LJPME)
         throw OpenMMException("getPMEParametersInContext: This Context is not using LJPME");
-    dalpha = ewaldDispersionAlpha;
-    dnx = dispersionGridSize[0];
-    dny = dispersionGridSize[1];
-    dnz = dispersionGridSize[2];
+    alpha = ewaldDispersionAlpha;
+    nx = dispersionGridSize[0];
+    ny = dispersionGridSize[1];
+    nz = dispersionGridSize[2];
 }
 
 ReferenceCalcCustomNonbondedForceKernel::~ReferenceCalcCustomNonbondedForceKernel() {

platforms/reference/src/SimTKReference/ReferenceLJCoulombIxn.cpp

@@ -166,18 +166,18 @@ void ReferenceLJCoulombIxn::setUsePME(RealOpenMM alpha, int meshSize[3]) {
 
 /**---------------------------------------------------------------------------------------
 
-     Set the force to use Particle-Mesh Ewald (PME) summation for dispersion and terms.
+     Set the force to use Particle-Mesh Ewald (PME) summation for dispersion terms.
 
-     @param dalpha  the dispersion Ewald separation parameter
-     @param dgridSize the dimensions of the dispersion mesh
+     @param alpha  the dispersion Ewald separation parameter
+     @param gridSize the dimensions of the dispersion mesh
 
      --------------------------------------------------------------------------------------- */
 
-void ReferenceLJCoulombIxn::setUseLJPME(RealOpenMM dalpha, int dmeshSize[3]) {
-    alphaDispersionEwald = dalpha;
-    dispersionMeshDim[0] = dmeshSize[0];
-    dispersionMeshDim[1] = dmeshSize[1];
-    dispersionMeshDim[2] = dmeshSize[2];
+void ReferenceLJCoulombIxn::setUseLJPME(RealOpenMM alpha, int meshSize[3]) {
+    alphaDispersionEwald = alpha;
+    dispersionMeshDim[0] = meshSize[0];
+    dispersionMeshDim[1] = meshSize[1];
+    dispersionMeshDim[2] = meshSize[2];
     ljpme = true;
 }
 
@@ -455,13 +455,14 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
             // and multiplicative functional forms at the cutoff.
             RealOpenMM emult = c6i*c6j*inverseR2*inverseR2*inverseR2*(1.0 - EXP(-dar2) * (1.0 + dar2 + 0.5*dar4));
             dEdR += 6.0*c6i*c6j*inverseR2*inverseR2*inverseR2*inverseR2*(1.0 - EXP(-dar2) * (1.0 + dar2 + 0.5*dar4 + dar6/6.0));
-            // The additive part of the potential shift
+
             RealOpenMM inverseCut2 = 1.0/(cutoffDistance*cutoffDistance);
             RealOpenMM inverseCut6 = inverseCut2*inverseCut2*inverseCut2;
             sig2 = atomParameters[ii][SigIndex] +  atomParameters[jj][SigIndex];
             sig2 *= sig2;
             sig6 = sig2*sig2*sig2;
-            RealOpenMM potentialshift = eps*sig6*inverseCut6;
+            // The additive part of the potential shift
+            RealOpenMM potentialshift = eps*(one-sig6*inverseCut6)*sig6*inverseCut6;
             dalphaR   = alphaDispersionEwald * cutoffDistance;
             dar2 = dalphaR*dalphaR;
             dar4 = dar2*dar2;