Initial reference implementation of LJPME.

olla/include/openmm/kernels.h

@@ -555,7 +555,8 @@ public:
         CutoffNonPeriodic = 1,
         CutoffPeriodic = 2,
         Ewald = 3,
-        PME = 4
+        PME = 4,
+        LJPME = 5
     };
     static std::string Name() {
         return "CalcNonbondedForce";
@@ -589,13 +590,22 @@ public:
     virtual void copyParametersToContext(ContextImpl& context, const NonbondedForce& force) = 0;
     /**
      * Get the parameters being used for PME.
-     * 
+     *
      * @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 getPMEParameters(double& alpha, int& nx, int& ny, int& nz) const = 0;
+    /**
+     * 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
+     */
+    virtual void getLJPMEParameters(double& dalpha, int& dnx, int& dny, int& dnz) const = 0;
 };
 
 /**

openmmapi/include/openmm/NonbondedForce.h

@@ -109,7 +109,12 @@ public:
          * Periodic boundary conditions are used, and Particle-Mesh Ewald (PME) summation is used to compute the interaction of each particle
          * with all periodic copies of every other particle.
          */
-        PME = 4
+        PME = 4,
+        /**
+         * Periodic boundary conditions are used, and Particle-Mesh Ewald (PME) summation is used to compute the interaction of each particle
+         * with all periodic copies of every other particle for both electrostatics and dispersion.  No switching is used for either interaction.
+         */
+        LJPME = 5
     };
     /**
      * Create a NonbondedForce.
@@ -207,6 +212,16 @@ public:
      * @param[out] nz      the number of grid points along the Z axis
      */
     void getPMEParameters(double& alpha, int& nx, int& ny, int& nz) const;
+    /**
+     * 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
+     */
+    void getLJPMEParameters(double& dalpha, int& dnx, int& dny, int& dnz) 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.
@@ -217,6 +232,16 @@ public:
      * @param nz      the number of grid points along the Z axis
      */
     void setPMEParameters(double alpha, int nx, int ny, int nz);
+    /**
+     * 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
+     */
+    void setLJPMEParameters(double dalpha, int dnx, int dny, int dnz);
     /**
      * 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
@@ -230,6 +255,19 @@ public:
      * @param[out] nz      the number of grid points along the Z axis
      */
     void getPMEParametersInContext(const Context& context, double& alpha, int& nx, int& ny, int& nz) const;
+    /**
+     * Get the PME parameters being used for the dispersion term for LJPME 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 specified with setPMEParameters(), or the standard values calculated based on the Ewald error tolerance.
+     * 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
+     */
+    void getLJPMEParametersInContext(const Context& context, double& dalpha, int& dnx, int& dny, int& dnz) 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.
@@ -382,9 +420,9 @@ private:
     class ParticleInfo;
     class ExceptionInfo;
     NonbondedMethod nonbondedMethod;
-    double cutoffDistance, switchingDistance, rfDielectric, ewaldErrorTol, alpha;
+    double cutoffDistance, switchingDistance, rfDielectric, ewaldErrorTol, alpha, dalpha;
     bool useSwitchingFunction, useDispersionCorrection;
-    int recipForceGroup, nx, ny, nz;
+    int recipForceGroup, nx, ny, nz, dnx, dny, dnz;
     void addExclusionsToSet(const std::vector<std::set<int> >& bonded12, std::set<int>& exclusions, int baseParticle, int fromParticle, int currentLevel) const;
     std::vector<ParticleInfo> particles;
     std::vector<ExceptionInfo> exceptions;

openmmapi/include/openmm/internal/NonbondedForceImpl.h

@@ -65,6 +65,7 @@ public:
     std::vector<std::string> getKernelNames();
     void updateParametersInContext(ContextImpl& context);
     void getPMEParameters(double& alpha, int& nx, int& ny, int& nz) const;
+    void getLJPMEParameters(double& alpha, int& nx, int& ny, int& nz) const;
     /**
      * This is a utility routine that calculates the values to use for alpha and kmax when using
      * Ewald summation.
@@ -74,7 +75,7 @@ public:
      * This is a utility routine that calculates the values to use for alpha and grid size when using
      * Particle Mesh Ewald.
      */
-    static void calcPMEParameters(const System& system, const NonbondedForce& force, double& alpha, int& xsize, int& ysize, int& zsize);
+    static void calcPMEParameters(const System& system, const NonbondedForce& force, double& alpha, int& xsize, int& ysize, int& zsize, bool LJ);
     /**
      * Compute the coefficient which, when divided by the periodic box volume, gives the
      * long range dispersion correction to the energy.

openmmapi/src/NonbondedForce.cpp

@@ -106,6 +106,13 @@ 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::setPMEParameters(double alpha, int nx, int ny, int nz) {
     this->alpha = alpha;
     this->nx = nx;
@@ -113,10 +120,21 @@ 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::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);
+}
+
 int NonbondedForce::addParticle(double charge, double sigma, double epsilon) {
     particles.push_back(ParticleInfo(charge, sigma, epsilon));
     return particles.size()-1;

openmmapi/src/NonbondedForceImpl.cpp

@@ -151,8 +151,12 @@ void NonbondedForceImpl::calcEwaldParameters(const System& system, const Nonbond
         kmaxz++;
 }
 
-void NonbondedForceImpl::calcPMEParameters(const System& system, const NonbondedForce& force, double& alpha, int& xsize, int& ysize, int& zsize) {
-    force.getPMEParameters(alpha, xsize, ysize, zsize);
+void NonbondedForceImpl::calcPMEParameters(const System& system, const NonbondedForce& force, double& alpha, int& xsize, int& ysize, int& zsize, bool LJ) {
+    if(LJ) {
+        force.getLJPMEParameters(alpha, xsize, ysize, zsize);
+    } else {
+        force.getPMEParameters(alpha, xsize, ysize, zsize);
+    }
     if (alpha == 0.0) {
         Vec3 boxVectors[3];
         system.getDefaultPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
@@ -283,3 +287,7 @@ void NonbondedForceImpl::updateParametersInContext(ContextImpl& context) {
 void NonbondedForceImpl::getPMEParameters(double& alpha, int& nx, int& ny, int& nz) const {
     kernel.getAs<CalcNonbondedForceKernel>().getPMEParameters(alpha, nx, ny, nz);
 }
+
+void NonbondedForceImpl::getLJPMEParameters(double& alpha, int& nx, int& ny, int& nz) const {
+    kernel.getAs<CalcNonbondedForceKernel>().getLJPMEParameters(alpha, nx, ny, nz);
+}

platforms/cpu/include/CpuKernels.h

@@ -251,27 +251,36 @@ public:
     void copyParametersToContext(ContextImpl& context, const NonbondedForce& force);
     /**
      * Get the parameters being used for PME.
-     * 
+     *
      * @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 getPMEParameters(double& alpha, int& nx, int& ny, int& nz) const;
+    /**
+     * 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
+     */
+    void getLJPMEParameters(double& dalpha, int& dnx, int& dny, int& dnz) const;
 private:
     class PmeIO;
     CpuPlatform::PlatformData& data;
     int numParticles, num14;
     int **bonded14IndexArray;
     double **bonded14ParamArray;
-    double nonbondedCutoff, switchingDistance, rfDielectric, ewaldAlpha, ewaldSelfEnergy, dispersionCoefficient;
-    int kmax[3], gridSize[3];
-    bool useSwitchingFunction, useOptimizedPme, hasInitializedPme;
+    double nonbondedCutoff, switchingDistance, rfDielectric, ewaldAlpha, ewaldDispersionAlpha, ewaldSelfEnergy, dispersionCoefficient;
+    int kmax[3], gridSize[3], dispersionGridSize[3];
+    bool useSwitchingFunction, useOptimizedPme, hasInitializedPme, hasInitializedDispersionPme;
     std::vector<std::set<int> > exclusions;
     std::vector<std::pair<float, float> > particleParams;
     NonbondedMethod nonbondedMethod;
     CpuNonbondedForce* nonbonded;
-    Kernel optimizedPme;
+    Kernel optimizedPme, optimizedDispersionPme;
     CpuBondForce bondForce;
 };
 

platforms/cpu/src/CpuKernels.cpp

@@ -528,7 +528,7 @@ CpuNonbondedForce* createCpuNonbondedForceVec4();
 CpuNonbondedForce* createCpuNonbondedForceVec8();
 
 CpuCalcNonbondedForceKernel::CpuCalcNonbondedForceKernel(string name, const Platform& platform, CpuPlatform::PlatformData& data) : CalcNonbondedForceKernel(name, platform),
-        data(data), bonded14IndexArray(NULL), bonded14ParamArray(NULL), hasInitializedPme(false), nonbonded(NULL) {
+        data(data), bonded14IndexArray(NULL), bonded14ParamArray(NULL), hasInitializedPme(false), hasInitializedDispersionPme(false), nonbonded(NULL) {
     if (isVec8Supported())
         nonbonded = createCpuNonbondedForceVec8();
     else
@@ -616,7 +616,7 @@ void CpuCalcNonbondedForceKernel::initialize(const System& system, const Nonbond
     }
     else if (nonbondedMethod == PME) {
         double alpha;
-        NonbondedForceImpl::calcPMEParameters(system, force, alpha, gridSize[0], gridSize[1], gridSize[2]);
+        NonbondedForceImpl::calcPMEParameters(system, force, alpha, gridSize[0], gridSize[1], gridSize[2], false);
         ewaldAlpha = alpha;
     }
     if (nonbondedMethod == Ewald || nonbondedMethod == PME)
@@ -739,7 +739,7 @@ void CpuCalcNonbondedForceKernel::copyParametersToContext(ContextImpl& context,
 }
 
 void CpuCalcNonbondedForceKernel::getPMEParameters(double& alpha, int& nx, int& ny, int& nz) const {
-    if (nonbondedMethod != PME)
+    if (nonbondedMethod != PME && nonbondedMethod != LJPME)
         throw OpenMMException("getPMEParametersInContext: This Context is not using PME");
     if (useOptimizedPme)
         optimizedPme.getAs<const CalcPmeReciprocalForceKernel>().getPMEParameters(alpha, nx, ny, nz);
@@ -751,6 +751,19 @@ void CpuCalcNonbondedForceKernel::getPMEParameters(double& alpha, int& nx, int&
     }
 }
 
+void CpuCalcNonbondedForceKernel::getLJPMEParameters(double& dalpha, int& dnx, int& dny, int& dnz) const {
+    if (nonbondedMethod != LJPME)
+        throw OpenMMException("getPMEParametersInContext: This Context is not using PME");
+    if (useOptimizedPme)
+        optimizedDispersionPme.getAs<const CalcPmeReciprocalForceKernel>().getPMEParameters(dalpha, dnx, dny, dnz);
+    else {
+        dalpha = ewaldDispersionAlpha;
+        dnx = dispersionGridSize[0];
+        dny = dispersionGridSize[1];
+        dnz = dispersionGridSize[2];
+    }
+}
+
 CpuCalcCustomNonbondedForceKernel::CpuCalcCustomNonbondedForceKernel(string name, const Platform& platform, CpuPlatform::PlatformData& data) :
             CalcCustomNonbondedForceKernel(name, platform), data(data), forceCopy(NULL), nonbonded(NULL) {
 }

platforms/opencl/include/OpenCLKernels.h

@@ -607,13 +607,22 @@ public:
     void copyParametersToContext(ContextImpl& context, const NonbondedForce& force);
     /**
      * Get the parameters being used for PME.
-     * 
+     *
      * @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 getPMEParameters(double& alpha, int& nx, int& ny, int& nz) const;
+    /**
+     * Get the parameters being used for the dispersion term in LJPME.
+     *
+     * @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;
 private:
     class SortTrait : public OpenCLSort::SortTrait {
         int getDataSize() const {return 8;}
@@ -664,8 +673,9 @@ private:
     cl::Kernel pmeInterpolateForceKernel;
     std::map<std::string, std::string> pmeDefines;
     std::vector<std::pair<int, int> > exceptionAtoms;
-    double ewaldSelfEnergy, dispersionCoefficient, alpha;
+    double ewaldSelfEnergy, dispersionCoefficient, alpha, dispersionAlpha;
     int gridSizeX, gridSizeY, gridSizeZ;
+    int dispersionGridSizeX, dispersionGridSizeY, dispersionGridSizeZ;
     bool hasCoulomb, hasLJ, usePmeQueue;
     NonbondedMethod nonbondedMethod;
     static const int PmeOrder = 5;

platforms/opencl/include/OpenCLParallelKernels.h

@@ -431,13 +431,22 @@ public:
     void copyParametersToContext(ContextImpl& context, const NonbondedForce& force);
     /**
      * Get the parameters being used for PME.
-     * 
+     *
      * @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 getPMEParameters(double& alpha, int& nx, int& ny, int& nz) const;
+    /**
+     * 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
+     */
+    void getLJPMEParameters(double& dalpha, int& dnx, int& dny, int& dnz) const;
 private:
     class Task;
     OpenCLPlatform::PlatformData& data;

platforms/opencl/src/OpenCLKernels.cpp

@@ -1733,7 +1733,7 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
     else if (nonbondedMethod == PME) {
         // Compute the PME parameters.
 
-        NonbondedForceImpl::calcPMEParameters(system, force, alpha, gridSizeX, gridSizeY, gridSizeZ);
+        NonbondedForceImpl::calcPMEParameters(system, force, alpha, gridSizeX, gridSizeY, gridSizeZ, false);
         gridSizeX = OpenCLFFT3D::findLegalDimension(gridSizeX);
         gridSizeY = OpenCLFFT3D::findLegalDimension(gridSizeY);
         gridSizeZ = OpenCLFFT3D::findLegalDimension(gridSizeZ);
@@ -2205,6 +2205,19 @@ void OpenCLCalcNonbondedForceKernel::getPMEParameters(double& alpha, int& nx, in
     }
 }
 
+void OpenCLCalcNonbondedForceKernel::getLJPMEParameters(double& dalpha, int& dnx, int& dny, int& dnz) 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);
+    else {
+        dalpha = this->alpha;
+        dnx = gridSizeX;
+        dny = gridSizeY;
+        dnz = gridSizeZ;
+    }
+}
+
 class OpenCLCustomNonbondedForceInfo : public OpenCLForceInfo {
 public:
     OpenCLCustomNonbondedForceInfo(int requiredBuffers, const CustomNonbondedForce& force) : OpenCLForceInfo(requiredBuffers), force(force) {

platforms/opencl/src/OpenCLParallelKernels.cpp

@@ -583,6 +583,10 @@ 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);
+}
+
 class OpenCLParallelCalcCustomNonbondedForceKernel::Task : public OpenCLContext::WorkTask {
 public:
     Task(ContextImpl& context, OpenCLCalcCustomNonbondedForceKernel& kernel, bool includeForce,

platforms/reference/include/ReferenceKernels.h

@@ -604,12 +604,21 @@ public:
      * @param nz      the number of grid points along the Z axis
      */
     void getPMEParameters(double& alpha, int& nx, int& ny, int& nz) const;
+    /**
+     * 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
+     */
+    void getLJPMEParameters(double& dalpha, int& dnx, int& dny, int& dnz) const;
 private:
     int numParticles, num14;
     int **bonded14IndexArray;
     RealOpenMM **particleParamArray, **bonded14ParamArray;
-    RealOpenMM nonbondedCutoff, switchingDistance, rfDielectric, ewaldAlpha, dispersionCoefficient;
-    int kmax[3], gridSize[3];
+    RealOpenMM nonbondedCutoff, switchingDistance, rfDielectric, ewaldAlpha, ewaldDispersionAlpha, dispersionCoefficient;
+    int kmax[3], gridSize[3], dispersionGridSize[3];
     bool useSwitchingFunction;
     std::vector<std::set<int> > exclusions;
     NonbondedMethod nonbondedMethod;

platforms/reference/include/ReferenceLJCoulombIxn.h

@@ -38,14 +38,14 @@ class ReferenceLJCoulombIxn {
       bool useSwitch;
       bool periodic;
       bool ewald;
-      bool pme;
+      bool pme, ljpme;
       const OpenMM::NeighborList* neighborList;
       OpenMM::RealVec periodicBoxVectors[3];
       RealOpenMM cutoffDistance, switchingDistance;
       RealOpenMM krf, crf;
-      RealOpenMM alphaEwald;
+      RealOpenMM alphaEwald, alphaDispersionEwald;
       int numRx, numRy, numRz;
-      int meshDim[3];
+      int meshDim[3], dispersionMeshDim[3];
 
       // parameter indices
 
@@ -139,16 +139,28 @@ class 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 setUsePME(RealOpenMM alpha, int meshSize[3]);
-      
+
+
+      /**---------------------------------------------------------------------------------------
+
+         Set the force to use Particle-Mesh Ewald (PME) summation for dispersion.
+
+         @param dalpha    the dispersion Ewald separation parameter
+         @param dgridSize the dimensions of the dispersion mesh
+
+         --------------------------------------------------------------------------------------- */
+
+      void setUseLJPME(RealOpenMM dalpha, int dmeshSize[3]);
+
       /**---------------------------------------------------------------------------------------
       
          Calculate LJ Coulomb pair ixn

platforms/reference/include/ReferencePME.h

@@ -87,6 +87,28 @@ pme_exec(pme_t       pme,
          RealOpenMM *    energy);
 
 
+/*
+ * Evaluate reciprocal space PME dispersion energy and forces.
+ *
+ * Args:
+ *
+ * pme         Opaque pme_t object, must have been initialized with pme_init()
+ * x           Pointer to coordinate data array (nm)
+ * f           Pointer to force data array (will be written as kJ/mol/nm)
+ * c6s         Array of c6 coefficients (units of sqrt(kJ/mol).nm^3 )
+ * box         Simulation cell dimensions (nm)
+ * energy      Total energy (will be written in units of kJ/mol)
+ */
+int OPENMM_EXPORT
+pme_exec_dpme(pme_t       pme,
+              const std::vector<OpenMM::RealVec>& atomCoordinates,
+              std::vector<OpenMM::RealVec>& forces,
+              const std::vector<RealOpenMM>& c6s,
+              const OpenMM::RealVec  periodicBoxVectors[3],
+              RealOpenMM *    energy);
+
+
+
 
 /* Release all memory in pme structure */
 int OPENMM_EXPORT
@@ -94,4 +116,4 @@ pme_destroy(pme_t    pme);
 
 } // namespace OpenMM
 
-#endif // __ReferencePME_H__
+#endif // __ReferencePME_H__
\ No newline at end of file

platforms/reference/src/ReferenceKernels.cpp

@@ -969,9 +969,17 @@ void ReferenceCalcNonbondedForceKernel::initialize(const System& system, const N
     }
     else if (nonbondedMethod == PME) {
         double alpha;
-        NonbondedForceImpl::calcPMEParameters(system, force, alpha, gridSize[0], gridSize[1], gridSize[2]);
+        NonbondedForceImpl::calcPMEParameters(system, force, alpha, gridSize[0], gridSize[1], gridSize[2], false);
         ewaldAlpha = (RealOpenMM) alpha;
     }
+    else if (nonbondedMethod == LJPME) {
+        double alpha;
+        NonbondedForceImpl::calcPMEParameters(system, force, alpha, gridSize[0], gridSize[1], gridSize[2], false);
+        ewaldAlpha = (RealOpenMM) alpha;
+        NonbondedForceImpl::calcPMEParameters(system, force, alpha, dispersionGridSize[0], dispersionGridSize[1], dispersionGridSize[2], true);
+        ewaldDispersionAlpha = (RealOpenMM) alpha;
+        useSwitchingFunction = false;
+    }
     rfDielectric = (RealOpenMM)force.getReactionFieldDielectric();
     if (force.getUseDispersionCorrection())
         dispersionCoefficient = NonbondedForceImpl::calcDispersionCorrection(system, force);
@@ -987,11 +995,12 @@ double ReferenceCalcNonbondedForceKernel::execute(ContextImpl& context, bool inc
     bool periodic = (nonbondedMethod == CutoffPeriodic);
     bool ewald  = (nonbondedMethod == Ewald);
     bool pme  = (nonbondedMethod == PME);
+    bool ljpme = (nonbondedMethod == LJPME);
     if (nonbondedMethod != NoCutoff) {
-        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, extractBoxVectors(context), periodic || ewald || pme, nonbondedCutoff, 0.0);
+        computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, extractBoxVectors(context), periodic || ewald || pme || ljpme, nonbondedCutoff, 0.0);
         clj.setUseCutoff(nonbondedCutoff, *neighborList, rfDielectric);
     }
-    if (periodic || ewald || pme) {
+    if (periodic || ewald || pme || ljpme) {
         RealVec* boxVectors = extractBoxVectors(context);
         double minAllowedSize = 1.999999*nonbondedCutoff;
         if (boxVectors[0][0] < minAllowedSize || boxVectors[1][1] < minAllowedSize || boxVectors[2][2] < minAllowedSize)
@@ -1002,6 +1011,10 @@ double ReferenceCalcNonbondedForceKernel::execute(ContextImpl& context, bool inc
         clj.setUseEwald(ewaldAlpha, kmax[0], kmax[1], kmax[2]);
     if (pme)
         clj.setUsePME(ewaldAlpha, gridSize);
+    if (ljpme){
+        clj.setUsePME(ewaldAlpha, gridSize);
+        clj.setUseLJPME(ewaldDispersionAlpha, dispersionGridSize);
+    }
     if (useSwitchingFunction)
         clj.setUseSwitchingFunction(switchingDistance);
     clj.calculatePairIxn(numParticles, posData, particleParamArray, exclusions, 0, forceData, 0, includeEnergy ? &energy : NULL, includeDirect, includeReciprocal);
@@ -1059,14 +1072,23 @@ void ReferenceCalcNonbondedForceKernel::copyParametersToContext(ContextImpl& con
 }
 
 void ReferenceCalcNonbondedForceKernel::getPMEParameters(double& alpha, int& nx, int& ny, int& nz) const {
-    if (nonbondedMethod != PME)
-        throw OpenMMException("getPMEParametersInContext: This Context is not using PME");
+    if (nonbondedMethod != PME && nonbondedMethod != LJPME)
+        throw OpenMMException("getPMEParametersInContext: This Context is not using PME or LJPME");
     alpha = ewaldAlpha;
     nx = gridSize[0];
     ny = gridSize[1];
     nz = gridSize[2];
 }
 
+void ReferenceCalcNonbondedForceKernel::getLJPMEParameters(double& dalpha, int& dnx, int& dny, int& dnz) const {
+    if (nonbondedMethod != LJPME)
+        throw OpenMMException("getPMEParametersInContext: This Context is not using LJPME");
+    dalpha = ewaldDispersionAlpha;
+    dnx = dispersionGridSize[0];
+    dny = dispersionGridSize[1];
+    dnz = dispersionGridSize[2];
+}
+
 ReferenceCalcCustomNonbondedForceKernel::~ReferenceCalcCustomNonbondedForceKernel() {
     disposeRealArray(particleParamArray, numParticles);
     if (neighborList != NULL)

platforms/reference/src/SimTKReference/ReferencePME.cpp

@@ -513,6 +513,106 @@ pme_reciprocal_convolution(pme_t     pme,
 }
 
 
+static void
+dpme_reciprocal_convolution(pme_t     pme,
+                           const RealVec periodicBoxVectors[3],
+                           const RealVec recipBoxVectors[3],
+                           RealOpenMM *  energy)
+{
+    int kx,ky,kz;
+    int nx,ny,nz;
+    RealOpenMM mx,my,mz;
+    RealOpenMM mhx,mhy,mhz,m2;
+    RealOpenMM bx,by,bz;
+    RealOpenMM d1,d2;
+    RealOpenMM eterm,struct2,ets2;
+    RealOpenMM esum;
+    RealOpenMM denom;
+    RealOpenMM boxfactor;
+    RealOpenMM maxkx,maxky,maxkz;
+
+    t_complex *ptr;
+
+    nx = pme->ngrid[0];
+    ny = pme->ngrid[1];
+    nz = pme->ngrid[2];
+
+    boxfactor = (RealOpenMM) M_PI*sqrt(M_PI) / (6.0*periodicBoxVectors[0][0]*periodicBoxVectors[1][1]*periodicBoxVectors[2][2]);
+
+    esum = 0;
+
+    maxkx = (RealOpenMM) ((nx+1)/2);
+    maxky = (RealOpenMM) ((ny+1)/2);
+    maxkz = (RealOpenMM) ((nz+1)/2);
+
+    RealOpenMM bfac = M_PI / pme->ewaldcoeff;
+    RealOpenMM fac1 = 2.0*M_PI*M_PI*M_PI*sqrt(M_PI);
+    RealOpenMM fac2 = pme->ewaldcoeff*pme->ewaldcoeff*pme->ewaldcoeff;
+    RealOpenMM fac3 = -2.0*pme->ewaldcoeff*M_PI*M_PI;
+    RealOpenMM b, m, m3, expfac, expterm, erfcterm;
+
+    for (kx=0;kx<nx;kx++)
+    {
+        /* Calculate frequency. Grid indices in the upper half correspond to negative frequencies! */
+        mx  = (RealOpenMM) ((kx<maxkx) ? kx : (kx-nx));
+        mhx = mx*recipBoxVectors[0][0];
+        bx  = pme->bsplines_moduli[0][kx];
+
+        for (ky=0;ky<ny;ky++)
+        {
+            /* Calculate frequency. Grid indices in the upper half correspond to negative frequencies! */
+            my  = (RealOpenMM) ((ky<maxky) ? ky : (ky-ny));
+            mhy = mx*recipBoxVectors[1][0]+my*recipBoxVectors[1][1];
+            by  = pme->bsplines_moduli[1][ky];
+
+            for (kz=0;kz<nz;kz++)
+            {
+                /*
+                 * Unlike the Coulombic case, there's an m=0 term so all terms are considered here.
+                 */
+
+                /* Calculate frequency. Grid indices in the upper half correspond to negative frequencies! */
+                mz        = (RealOpenMM) ((kz<maxkz) ? kz : (kz-nz));
+                mhz       = mx*recipBoxVectors[2][0]+my*recipBoxVectors[2][1]+mz*recipBoxVectors[2][2];
+
+                /* Pointer to the grid cell in question */
+                ptr       = pme->grid + kx*ny*nz + ky*nz + kz;
+
+                /* Get grid data for this frequency */
+                d1        = ptr->re;
+                d2        = ptr->im;
+
+                /* Calculate the convolution - see the Essman/Darden paper for the equation! */
+                m2        = mhx*mhx+mhy*mhy+mhz*mhz;
+                bz        = pme->bsplines_moduli[2][kz];
+                denom     = boxfactor / (bx*by*bz);
+
+                m = sqrt(m2);
+                m3 = m*m2;
+                b = bfac*m;
+                expfac = -b*b;
+                erfcterm = erfc(b);
+                expterm = exp(expfac);
+
+                eterm     = (fac1*erfcterm*m3 + expterm*(fac2 + fac3*m2)) * denom;
+
+                /* write back convolution data to grid */
+                ptr->re   = d1*eterm;
+                ptr->im   = d2*eterm;
+
+                struct2   = (d1*d1+d2*d2);
+
+                /* Long-range PME contribution to the energy for this frequency */
+                ets2      = eterm*struct2;
+                esum     += ets2;
+            }
+        }
+    }
+    // Remember the C6 energy is attractive, hence the negative sign.
+    *energy = (RealOpenMM) (-esum);
+}
+
+
 static void
 pme_grid_interpolate_force(pme_t pme,
                            const RealVec recipBoxVectors[3],
@@ -704,6 +804,49 @@ int pme_exec(pme_t       pme,
 }
 
 
+int pme_exec_dpme(pme_t       pme,
+             const vector<RealVec>& atomCoordinates,
+             vector<RealVec>& forces,
+             const vector<RealOpenMM>& c6s,
+             const RealVec periodicBoxVectors[3],
+             RealOpenMM* energy)
+{
+    /* Routine is called with coordinates in x, a box, and charges in q */
+
+    RealVec recipBoxVectors[3];
+    invert_box_vectors(periodicBoxVectors, recipBoxVectors);
+
+    /* Before we can do the actual interpolation, we need to recalculate and update
+     * the indices for each particle in the charge grid (initialized in pme_init()),
+     * and what its fractional offset in this grid cell is.
+     */
+
+    /* Update charge grid indices and fractional offsets for each atom.
+     * The indices/fractions are stored internally in the pme datatype
+     */
+    pme_update_grid_index_and_fraction(pme,atomCoordinates,periodicBoxVectors,recipBoxVectors);
+
+    /* Calculate bsplines (and their differentials) from current fractional coordinates, store in pme structure */
+    pme_update_bsplines(pme);
+
+    /* Spread the charges on grid (using newly calculated bsplines in the pme structure) */
+    pme_grid_spread_charge(pme, c6s);
+
+    /* do 3d-fft */
+    fftpack_exec_3d(pme->fftplan,FFTPACK_FORWARD,pme->grid,pme->grid);
+
+    /* solve in k-space */
+    dpme_reciprocal_convolution(pme,periodicBoxVectors,recipBoxVectors,energy);
+
+    /* do 3d-invfft */
+    fftpack_exec_3d(pme->fftplan,FFTPACK_BACKWARD,pme->grid,pme->grid);
+
+    /* Get the particle forces from the grid and bsplines in the pme structure */
+    pme_grid_interpolate_force(pme,recipBoxVectors,c6s,forces);
+
+    return 0;
+}
+
 
 int
 pme_destroy(pme_t    pme)

plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.cpp

@@ -602,7 +602,7 @@ void ReferenceCalcAmoebaMultipoleForceKernel::initialize(const System& system, c
             nb.setEwaldErrorTolerance(force.getEwaldErrorTolerance());
             nb.setCutoffDistance(force.getCutoffDistance());
             int gridSizeX, gridSizeY, gridSizeZ;
-            NonbondedForceImpl::calcPMEParameters(system, nb, alphaEwald, gridSizeX, gridSizeY, gridSizeZ);
+            NonbondedForceImpl::calcPMEParameters(system, nb, alphaEwald, gridSizeX, gridSizeY, gridSizeZ, false);
             pmeGridDimension[0] = gridSizeX;
             pmeGridDimension[1] = gridSizeY;
             pmeGridDimension[2] = gridSizeZ;

plugins/cpupme/tests/TestCpuPme.cpp

@@ -108,7 +108,7 @@ void testPME(bool triclinic) {
     
     double alpha;
     int gridx, gridy, gridz;
-    NonbondedForceImpl::calcPMEParameters(system, *force, alpha, gridx, gridy, gridz);
+    NonbondedForceImpl::calcPMEParameters(system, *force, alpha, gridx, gridy, gridz, false);
     CpuCalcPmeReciprocalForceKernel pme(CalcPmeReciprocalForceKernel::Name(), platform);
     IO io;
     double sumSquaredCharges = 0;

tests/TestEwald.h

@@ -365,7 +365,7 @@ void testErrorTolerance(NonbondedForce::NonbondedMethod method) {
 
                 double expectedAlpha, actualAlpha;
                 int expectedSize[3], actualSize[3];
-                NonbondedForceImpl::calcPMEParameters(system, *force, expectedAlpha, expectedSize[0], expectedSize[1], expectedSize[2]);
+                NonbondedForceImpl::calcPMEParameters(system, *force, expectedAlpha, expectedSize[0], expectedSize[1], expectedSize[2], false);
                 force->getPMEParametersInContext(context, actualAlpha, actualSize[0], actualSize[1], actualSize[2]);
                 ASSERT_EQUAL_TOL(expectedAlpha, actualAlpha, 1e-5);
                 for (int i = 0; i < 3; i++) {