Merge pull request #4 from peastman/ljpme

Revised calculation for dispersion PME grid size

docs-source/usersguide/theory.rst

@@ -458,7 +458,7 @@ The formula used to compute the number of nodes along each dimension of the mesh
 is slightly different from the one used for Coulomb interactions:
 
 .. math::
-   n_\mathit{mesh}=\frac{\alpha d}{{3\delta}^{1/20}}
+   n_\mathit{mesh}=\frac{\alpha d}{{3\delta}^{1/5}}
 
 As before, this is an empirical formula.  It will usually produce an average
 relative error in the forces less than or similar to :math:`\delta`\ , but that

openmmapi/src/NonbondedForceImpl.cpp

@@ -160,9 +160,9 @@ void NonbondedForceImpl::calcPMEParameters(const System& system, const Nonbonded
         double tol = force.getEwaldErrorTolerance();
         alpha = (1.0/force.getCutoffDistance())*std::sqrt(-log(2.0*tol));
         if (lj) {
-            xsize = (int) ceil(alpha*boxVectors[0][0]/(3*pow(tol, 0.05)));
-            ysize = (int) ceil(alpha*boxVectors[1][1]/(3*pow(tol, 0.05)));
-            zsize = (int) ceil(alpha*boxVectors[2][2]/(3*pow(tol, 0.05)));
+            xsize = (int) ceil(alpha*boxVectors[0][0]/(3*pow(tol, 0.2)));
+            ysize = (int) ceil(alpha*boxVectors[1][1]/(3*pow(tol, 0.2)));
+            zsize = (int) ceil(alpha*boxVectors[2][2]/(3*pow(tol, 0.2)));
         }
         else {
             xsize = (int) ceil(2*alpha*boxVectors[0][0]/(3*pow(tol, 0.2)));

platforms/cuda/src/CudaKernels.cpp

@@ -1623,12 +1623,18 @@ CudaCalcNonbondedForceKernel::~CudaCalcNonbondedForceKernel() {
         delete sort;
     if (fft != NULL)
         delete fft;
+    if (dispersionFft != NULL)
+        delete dispersionFft;
     if (pmeio != NULL)
         delete pmeio;
     if (hasInitializedFFT) {
         if (useCudaFFT) {
             cufftDestroy(fftForward);
             cufftDestroy(fftBackward);
+            if (doLJPME) {
+                cufftDestroy(dispersionFftForward);
+                cufftDestroy(dispersionFftBackward);                
+            }
         }
         if (usePmeStream) {
             cuStreamDestroy(pmeStream);
@@ -1911,6 +1917,8 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
                     if (useCudaFFT) {
                         cufftSetStream(fftForward, pmeStream);
                         cufftSetStream(fftBackward, pmeStream);
+                        cufftSetStream(dispersionFftForward, pmeStream);
+                        cufftSetStream(dispersionFftBackward, pmeStream);
                     }
                     CHECK_RESULT(cuEventCreate(&pmeSyncEvent, CU_EVENT_DISABLE_TIMING), "Error creating event for NonbondedForce");
                     int recipForceGroup = force.getReciprocalSpaceForceGroup();
@@ -2170,7 +2178,7 @@ double CudaCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeF
                     cufftExecR2C(dispersionFftForward, (float*) directPmeGrid->getDevicePointer(), (float2*) reciprocalPmeGrid->getDevicePointer());
             }
             else {
-                fft->execFFT(*directPmeGrid, *reciprocalPmeGrid, true);
+                dispersionFft->execFFT(*directPmeGrid, *reciprocalPmeGrid, true);
             }
 
             if (includeEnergy) {
@@ -2192,7 +2200,7 @@ double CudaCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeF
                     cufftExecC2R(dispersionFftBackward, (float2*) reciprocalPmeGrid->getDevicePointer(), (float*)  directPmeGrid->getDevicePointer());
             }
             else {
-                fft->execFFT(*reciprocalPmeGrid, *directPmeGrid, false);
+                dispersionFft->execFFT(*reciprocalPmeGrid, *directPmeGrid, false);
             }
 
             void* interpolateArgs[] = {&cu.getPosq().getDevicePointer(), &cu.getForce().getDevicePointer(), &directPmeGrid->getDevicePointer(), cu.getPeriodicBoxSizePointer(),

platforms/cuda/src/kernels/pme.cu

@@ -238,7 +238,7 @@ gridEvaluateEnergy(real2* __restrict__ halfcomplex_pmeGrid, mixed* __restrict__
 #endif
             energy += eterm*(grid.x*grid.x + grid.y*grid.y);
     }
-#ifdef USE_PME_STREAM
+#if defined(USE_PME_STREAM) && !defined(USE_LJPME)
     energyBuffer[blockIdx.x*blockDim.x+threadIdx.x] = 0.5f*energy;
 #else
     energyBuffer[blockIdx.x*blockDim.x+threadIdx.x] += 0.5f*energy;

tests/TestDispersionPME.h

@@ -43,8 +43,60 @@
 using namespace OpenMM;
 using namespace std;
 
+void testConvergence() {
+    // Create a cloud of random particles.
+
+    const int numParticles = 1000;
+    const double boxWidth = 6.0;
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(boxWidth, 0, 0), Vec3(0, boxWidth, 0), Vec3(0, 0, boxWidth));
+    NonbondedForce* force = new NonbondedForce();
+    system.addForce(force);
+    vector<Vec3> positions(numParticles);
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+
+    for (int i = 0; i < numParticles; i++) {
+        system.addParticle(1.0);
+        force->addParticle(0.0, 0.1+0.3*genrand_real2(sfmt), genrand_real2(sfmt));
+        while (true) {
+            Vec3 pos = Vec3(boxWidth*genrand_real2(sfmt), boxWidth*genrand_real2(sfmt), boxWidth*genrand_real2(sfmt));
+            double minDist = boxWidth;
+            for (int j = 0; j < i; j++) {
+                Vec3 delta = pos-positions[j];
+                minDist = min(minDist, sqrt(delta.dot(delta)));
+            }
+            if (minDist > 0.15) {
+                positions[i] = pos;
+                break;
+            }
+        }
+    }
+    
+    // Compute the energy with short and long cutoffs, and compare them to LJPME.
+    // The long cutoff should match the LJPME result better than the short cutoff.
+    
+    force->setNonbondedMethod(NonbondedForce::LJPME);
+    force->setCutoffDistance(1.0);
+    force->setEwaldErrorTolerance(1e-5);
+    force->setUseDispersionCorrection(false);
+    VerletIntegrator integrator(0.01);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    double pmeEnergy = context.getState(State::Energy, false, 1).getPotentialEnergy();
+    force->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    context.reinitialize();
+    context.setPositions(positions);
+    double shortEnergy = context.getState(State::Energy, false, 1).getPotentialEnergy();
+    force->setCutoffDistance(3.0);
+    context.reinitialize();
+    context.setPositions(positions);
+    double longEnergy = context.getState(State::Energy, false, 1).getPotentialEnergy();
+    ASSERT(fabs(longEnergy-pmeEnergy) < fabs(shortEnergy-pmeEnergy))
+}
+
 void testErrorTolerance() {
-    // Create a cloud of random point charges.
+    // Create a cloud of random particles.
 
     const int numParticles = 200;
     const double boxWidth = 5.0;
@@ -118,7 +170,7 @@ void testErrorTolerance() {
 }
 
 void testPMEParameters() {
-    // Create a cloud of random point charges.
+    // Create a cloud of random particles.
 
     const int numParticles = 51;
     const double boxWidth = 4.7;
@@ -136,10 +188,11 @@ void testPMEParameters() {
         positions[i] = Vec3(boxWidth*genrand_real2(sfmt), boxWidth*genrand_real2(sfmt), boxWidth*genrand_real2(sfmt));
     }
     force->setNonbondedMethod(NonbondedForce::LJPME);
+    force->setCutoffDistance(0.5);
     
-    // Compute the energy with an error tolerance of 1e-3.
+    // Compute the energy with an error tolerance of 0.1.
 
-    force->setEwaldErrorTolerance(1e-3);
+    force->setEwaldErrorTolerance(0.1);
     VerletIntegrator integrator1(0.01);
     Context context1(system, integrator1, platform);
     context1.setPositions(positions);
@@ -157,7 +210,7 @@ void testPMEParameters() {
     double energy2 = context2.getState(State::Energy).getPotentialEnergy();
     
     // Now explicitly set the parameters.  These should match the values that were
-    // used for tolerance 1e-3.
+    // used for tolerance 0.1.
 
     force->setLJPMEParameters(alpha, gridx, gridy, gridz);
     VerletIntegrator integrator3(0.01);
@@ -1806,6 +1859,7 @@ void runPlatformTests();
 int main(int argc, char* argv[]) {
     try {
         initializeTests(argc, argv);
+        testConvergence();
         testErrorTolerance();
         testPMEParameters();
         testWater2DpmeEnergiesForcesNoExclusions();