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";
@@ -596,6 +597,15 @@ public:
      * @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) {
+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

@@ -258,20 +258,29 @@ 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 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

@@ -614,6 +614,15 @@ 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 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

@@ -438,6 +438,15 @@ 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 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
 
@@ -149,6 +149,18 @@ class ReferenceLJCoulombIxn {
 
       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

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/ReferenceLJCoulombIxn.cpp

@@ -26,6 +26,7 @@
 #include <sstream>
 #include <complex>
 #include <algorithm>
+#include <iostream>
 
 #include "SimTKOpenMMUtilities.h"
 #include "ReferenceLJCoulombIxn.h"
@@ -47,7 +48,7 @@ using namespace OpenMM;
 
    --------------------------------------------------------------------------------------- */
 
-ReferenceLJCoulombIxn::ReferenceLJCoulombIxn() : cutoff(false), useSwitch(false), periodic(false), ewald(false), pme(false) {
+ReferenceLJCoulombIxn::ReferenceLJCoulombIxn() : cutoff(false), useSwitch(false), periodic(false), ewald(false), pme(false), ljpme(false) {
 
     // ---------------------------------------------------------------------------------------
 
@@ -73,7 +74,7 @@ ReferenceLJCoulombIxn::~ReferenceLJCoulombIxn() {
 
 }
 
-  /**---------------------------------------------------------------------------------------
+/**---------------------------------------------------------------------------------------
 
      Set the force to use a cutoff.
 
@@ -83,14 +84,14 @@ ReferenceLJCoulombIxn::~ReferenceLJCoulombIxn() {
 
      --------------------------------------------------------------------------------------- */
 
-  void ReferenceLJCoulombIxn::setUseCutoff(RealOpenMM distance, const OpenMM::NeighborList& neighbors, RealOpenMM solventDielectric) {
+void ReferenceLJCoulombIxn::setUseCutoff(RealOpenMM distance, const OpenMM::NeighborList& neighbors, RealOpenMM solventDielectric) {
 
     cutoff = true;
     cutoffDistance = distance;
     neighborList = &neighbors;
     krf = pow(cutoffDistance, -3.0)*(solventDielectric-1.0)/(2.0*solventDielectric+1.0);
     crf = (1.0/cutoffDistance)*(3.0*solventDielectric)/(2.0*solventDielectric+1.0);
-  }
+}
 
 /**---------------------------------------------------------------------------------------
 
@@ -105,7 +106,7 @@ void ReferenceLJCoulombIxn::setUseSwitchingFunction(RealOpenMM distance) {
     switchingDistance = distance;
 }
 
-  /**---------------------------------------------------------------------------------------
+/**---------------------------------------------------------------------------------------
 
      Set the force to use periodic boundary conditions.  This requires that a cutoff has
      also been set, and the smallest side of the periodic box is at least twice the cutoff
@@ -115,7 +116,7 @@ void ReferenceLJCoulombIxn::setUseSwitchingFunction(RealOpenMM distance) {
 
      --------------------------------------------------------------------------------------- */
 
-  void ReferenceLJCoulombIxn::setPeriodic(OpenMM::RealVec* vectors) {
+void ReferenceLJCoulombIxn::setPeriodic(OpenMM::RealVec* vectors) {
 
     assert(cutoff);
     assert(vectors[0][0] >= 2.0*cutoffDistance);
@@ -125,9 +126,9 @@ void ReferenceLJCoulombIxn::setUseSwitchingFunction(RealOpenMM distance) {
     periodicBoxVectors[0] = vectors[0];
     periodicBoxVectors[1] = vectors[1];
     periodicBoxVectors[2] = vectors[2];
-  }
+}
 
-  /**---------------------------------------------------------------------------------------
+/**---------------------------------------------------------------------------------------
 
      Set the force to use Ewald summation.
 
@@ -138,15 +139,15 @@ void ReferenceLJCoulombIxn::setUseSwitchingFunction(RealOpenMM distance) {
 
      --------------------------------------------------------------------------------------- */
 
-  void ReferenceLJCoulombIxn::setUseEwald(RealOpenMM alpha, int kmaxx, int kmaxy, int kmaxz) {
+void ReferenceLJCoulombIxn::setUseEwald(RealOpenMM alpha, int kmaxx, int kmaxy, int kmaxz) {
     alphaEwald = alpha;
     numRx = kmaxx;
     numRy = kmaxy;
     numRz = kmaxz;
     ewald = true;
-  }
+}
 
-  /**---------------------------------------------------------------------------------------
+/**---------------------------------------------------------------------------------------
 
      Set the force to use Particle-Mesh Ewald (PME) summation.
 
@@ -155,13 +156,30 @@ void ReferenceLJCoulombIxn::setUseSwitchingFunction(RealOpenMM distance) {
 
      --------------------------------------------------------------------------------------- */
 
-  void ReferenceLJCoulombIxn::setUsePME(RealOpenMM alpha, int meshSize[3]) {
+void ReferenceLJCoulombIxn::setUsePME(RealOpenMM alpha, int meshSize[3]) {
     alphaEwald = alpha;
     meshDim[0] = meshSize[0];
     meshDim[1] = meshSize[1];
     meshDim[2] = meshSize[2];
     pme = true;
-  }
+}
+
+/**---------------------------------------------------------------------------------------
+
+     Set the force to use Particle-Mesh Ewald (PME) summation for dispersion and terms.
+
+     @param dalpha  the dispersion Ewald separation parameter
+     @param dgridSize 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];
+    ljpme = true;
+}
 
 /**---------------------------------------------------------------------------------------
 
@@ -201,16 +219,27 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
     RealOpenMM totalSelfEwaldEnergy     = 0.0;
     RealOpenMM realSpaceEwaldEnergy     = 0.0;
     RealOpenMM recipEnergy              = 0.0;
+    RealOpenMM recipDispersionEnergy    = 0.0;
     RealOpenMM totalRecipEnergy         = 0.0;
     RealOpenMM vdwEnergy                = 0.0;
 
-// **************************************************************************************
-// SELF ENERGY
-// **************************************************************************************
+    // A couple of sanity checks for
+    if(ljpme && useSwitch)
+        throw OpenMMException("Switching cannot be used with LJPME");
+    if(ljpme && !pme)
+        throw OpenMMException("LJPME has been set, without PME being set");
+
+    // **************************************************************************************
+    // SELF ENERGY
+    // **************************************************************************************
 
     if (includeReciprocal) {
         for (int atomID = 0; atomID < numberOfAtoms; atomID++) {
             RealOpenMM selfEwaldEnergy       = (RealOpenMM) (ONE_4PI_EPS0*atomParameters[atomID][QIndex]*atomParameters[atomID][QIndex] * alphaEwald/SQRT_PI);
+            if(ljpme) {
+                // Dispersion self term
+                selfEwaldEnergy -= pow(alphaDispersionEwald, 6.0) * 64.0*pow(atomParameters[atomID][SigIndex], 6.0) * pow(atomParameters[atomID][EpsIndex], 2.0) / 12.0;
+            }
             totalSelfEwaldEnergy            -= selfEwaldEnergy;
             if (energyByAtom) {
                 energyByAtom[atomID]        -= selfEwaldEnergy;
@@ -222,9 +251,9 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
         *totalEnergy += totalSelfEwaldEnergy;
     }
 
-// **************************************************************************************
-// RECIPROCAL SPACE EWALD ENERGY AND FORCES
-// **************************************************************************************
+    // **************************************************************************************
+    // RECIPROCAL SPACE EWALD ENERGY AND FORCES
+    // **************************************************************************************
     // PME
 
     if (pme && includeReciprocal) {
@@ -245,8 +274,31 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
                 energyByAtom[n] += recipEnergy;
 
         pme_destroy(pmedata);
+
+        if (ljpme) {
+            // Dispersion reciprocal space terms
+            pme_init(&pmedata,alphaDispersionEwald,numberOfAtoms,dispersionMeshDim,5,1);
+
+            std::vector<RealVec> dpmeforces;
+            for (int i = 0; i < numberOfAtoms; i++){
+                charges[i] = 8.0*pow(atomParameters[i][SigIndex], 3.0) * atomParameters[i][EpsIndex];
+                dpmeforces.push_back(RealVec());
+            }
+            pme_exec_dpme(pmedata,atomCoordinates,dpmeforces,charges,periodicBoxVectors,&recipDispersionEnergy);
+            for (int i = 0; i < numberOfAtoms; i++){
+                forces[i][0] -= 2.0*dpmeforces[i][0];
+                forces[i][1] -= 2.0*dpmeforces[i][1];
+                forces[i][2] -= 2.0*dpmeforces[i][2];
             }
+            if (totalEnergy)
+                *totalEnergy += recipDispersionEnergy;
 
+            if (energyByAtom)
+                for (int n = 0; n < numberOfAtoms; n++)
+                    energyByAtom[n] += recipDispersionEnergy;
+            pme_destroy(pmedata);
+        }
+    }
     // Ewald method
 
     else if (ewald && includeReciprocal) {
@@ -258,7 +310,7 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
 
         // setup K-vectors
 
-  #define EIR(x, y, z) eir[(x)*numberOfAtoms*3+(y)*3+z]
+#define EIR(x, y, z) eir[(x)*numberOfAtoms*3+(y)*3+z]
         vector<d_complex> eir(kmax*numberOfAtoms*3);
         vector<d_complex> tab_xy(numberOfAtoms);
         vector<d_complex> tab_qxyz(numberOfAtoms);
@@ -350,15 +402,16 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
         }
     }
 
-// **************************************************************************************
-// SHORT-RANGE ENERGY AND FORCES
-// **************************************************************************************
+    // **************************************************************************************
+    // SHORT-RANGE ENERGY AND FORCES
+    // **************************************************************************************
 
     if (!includeDirect)
         return;
     RealOpenMM totalVdwEnergy            = 0.0f;
     RealOpenMM totalRealSpaceEwaldEnergy = 0.0f;
 
+
     for (int i = 0; i < (int) neighborList->size(); i++) {
         OpenMM::AtomPair pair = (*neighborList)[i];
         int ii = pair.first;
@@ -387,6 +440,36 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
         RealOpenMM eps       = atomParameters[ii][EpsIndex]*atomParameters[jj][EpsIndex];
         dEdR     += switchValue*eps*(twelve*sig6 - six)*sig6*inverseR*inverseR;
         vdwEnergy = eps*(sig6-one)*sig6;
+
+        if (ljpme) {
+            RealOpenMM dalphaR   = alphaDispersionEwald * r;
+            RealOpenMM dar2 = dalphaR*dalphaR;
+            RealOpenMM dar4 = dar2*dar2;
+            RealOpenMM dar6 = dar4*dar2;
+            RealOpenMM inverseR2 = inverseR*inverseR;
+            RealOpenMM c6i = 8.0*pow(atomParameters[ii][SigIndex], 3.0) * atomParameters[ii][EpsIndex];
+            RealOpenMM c6j = 8.0*pow(atomParameters[jj][SigIndex], 3.0) * atomParameters[jj][EpsIndex];
+            // For the energies and forces, we first add the regular Lorentz−Berthelot terms.  The C12 term is treated as usual
+            // but we then subtract out (remembering that the C6 term is negative) the multiplicative C6 term that has been
+            // computed in real space.  Finally, we add a potential shift term to account for the difference between the LB
+            // 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;
+            dalphaR   = alphaDispersionEwald * cutoffDistance;
+            dar2 = dalphaR*dalphaR;
+            dar4 = dar2*dar2;
+            // The multiplicative part of the potential shift
+            potentialshift -= c6i*c6j*inverseCut6*(1.0 - EXP(-dar2) * (1.0 + dar2 + 0.5*dar4));
+            vdwEnergy += emult + potentialshift;
+        }
+
         if (useSwitch) {
             dEdR -= vdwEnergy*switchDeriv*inverseR;
             vdwEnergy *= switchValue;
@@ -452,6 +535,24 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
                     realSpaceEwaldEnergy = (RealOpenMM) (alphaEwald*TWO_OVER_SQRT_PI*ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]);
                 }
 
+                if(ljpme){
+                    // Dispersion terms.  Here we just back out the reciprocal space terms, and don't add any extra real space terms.
+                    RealOpenMM dalphaR   = alphaDispersionEwald * r;
+                    RealOpenMM inverseR2 = inverseR*inverseR;
+                    RealOpenMM dar2 = dalphaR*dalphaR;
+                    RealOpenMM dar4 = dar2*dar2;
+                    RealOpenMM dar6 = dar4*dar2;
+                    RealOpenMM c6i = 8.0*pow(atomParameters[ii][SigIndex], 3.0) * atomParameters[ii][EpsIndex];
+                    RealOpenMM c6j = 8.0*pow(atomParameters[jj][SigIndex], 3.0) * atomParameters[jj][EpsIndex];
+                    realSpaceEwaldEnergy -= c6i*c6j*inverseR2*inverseR2*inverseR2*(1.0 - EXP(-dar2) * (1.0 + dar2 + 0.5*dar4));
+                    RealOpenMM dEdR = -6.0*c6i*c6j*inverseR2*inverseR2*inverseR2*inverseR2*(1.0 - EXP(-dar2) * (1.0 + dar2 + 0.5*dar4 + dar6/6.0));
+                    for (int kk = 0; kk < 3; kk++) {
+                        RealOpenMM force  = dEdR*deltaR[0][kk];
+                        forces[ii][kk]   -= force;
+                        forces[jj][kk]   += force;
+                    }
+                }
+
                 totalExclusionEnergy += realSpaceEwaldEnergy;
                 if (energyByAtom) {
                     energyByAtom[ii] -= realSpaceEwaldEnergy;
@@ -488,7 +589,7 @@ void ReferenceLJCoulombIxn::calculatePairIxn(int numberOfAtoms, vector<RealVec>&
                                              RealOpenMM* fixedParameters, vector<RealVec>& forces,
                                              RealOpenMM* energyByAtom, RealOpenMM* totalEnergy, bool includeDirect, bool includeReciprocal) const {
 
-   if (ewald || pme) {
+    if (ewald || pme || ljpme) {
         calculateEwaldIxn(numberOfAtoms, atomCoordinates, atomParameters, exclusions, fixedParameters, forces, energyByAtom,
                           totalEnergy, includeDirect, includeReciprocal);
         return;
@@ -512,7 +613,7 @@ void ReferenceLJCoulombIxn::calculatePairIxn(int numberOfAtoms, vector<RealVec>&
     }
 }
 
-  /**---------------------------------------------------------------------------------------
+/**---------------------------------------------------------------------------------------
 
      Calculate LJ Coulomb pair ixn between two atoms
 
@@ -608,5 +709,5 @@ void ReferenceLJCoulombIxn::calculateOneIxn(int ii, int jj, vector<RealVec>& ato
         energyByAtom[ii] += energy;
         energyByAtom[jj] += energy;
     }
-  }
+}
 

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++) {