Implemented new method for selecting Ewald and PME parameters

openmmapi/include/openmm/internal/NonbondedForceImpl.h

@@ -41,6 +41,8 @@
 
 namespace OpenMM {
 
+class System;
+
 /**
  * This is the internal implementation of NonbondedForce.
  */
@@ -62,7 +64,20 @@ public:
         return std::map<std::string, double>(); // This force field doesn't define any parameters.
     }
     std::vector<std::string> getKernelNames();
+    /**
+     * This is a utility routine that calculates the values to use for alpha and kmax when using
+     * Ewald summation.
+     */
+    static void calcEwaldParameters(const System& system, const NonbondedForce& force, double& alpha, int& kmaxx, int& kmaxy, int& kmaxz);
+    /**
+     * 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);
 private:
+    class ErrorFunction;
+    class EwaldErrorFunction;
+    static double findZero(const ErrorFunction& f, int initialGuess);
     NonbondedForce& owner;
     Kernel kernel;
 };

openmmapi/src/NonbondedForceImpl.cpp

@@ -33,13 +33,11 @@
 #include "openmm/internal/ContextImpl.h"
 #include "openmm/internal/NonbondedForceImpl.h"
 #include "openmm/kernels.h"
+#include <cmath>
 #include <sstream>
 
 using namespace OpenMM;
-using std::pair;
-using std::vector;
-using std::set;
-using std::stringstream;
+using namespace std;
 
 NonbondedForceImpl::NonbondedForceImpl(NonbondedForce& owner) : owner(owner) {
 }
@@ -99,3 +97,59 @@ std::vector<std::string> NonbondedForceImpl::getKernelNames() {
     names.push_back(CalcNonbondedForceKernel::Name());
     return names;
 }
+
+class NonbondedForceImpl::ErrorFunction {
+public:
+    virtual double getValue(int arg) const = 0;
+};
+
+class NonbondedForceImpl::EwaldErrorFunction : public ErrorFunction {
+public:
+    EwaldErrorFunction(double width, double alpha, double target) : width(width), alpha(alpha), target(target) {
+    }
+    double getValue(int arg) const {
+        double temp = arg*M_PI/(width*alpha);
+        return target-0.05*sqrt(width*alpha)*arg*exp(-temp*temp);
+    }
+private:
+    double width, alpha, target;
+};
+
+void NonbondedForceImpl::calcEwaldParameters(const System& system, const NonbondedForce& force, double& alpha, int& kmaxx, int& kmaxy, int& kmaxz) {
+    Vec3 boxVectors[3];
+    system.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
+    double tol = force.getEwaldErrorTolerance();
+    alpha = (1.0/force.getCutoffDistance())*std::sqrt(-log(0.5*tol));
+    kmaxx = findZero(EwaldErrorFunction(boxVectors[0][0], alpha, tol), 10);
+    kmaxy = findZero(EwaldErrorFunction(boxVectors[1][1], alpha, tol), 10);
+    kmaxz = findZero(EwaldErrorFunction(boxVectors[2][2], alpha, tol), 10);
+    if (kmaxx%2 == 0)
+        kmaxx++;
+    if (kmaxy%2 == 0)
+        kmaxy++;
+    if (kmaxz%2 == 0)
+        kmaxz++;
+}
+
+void NonbondedForceImpl::calcPMEParameters(const System& system, const NonbondedForce& force, double& alpha, int& xsize, int& ysize, int& zsize) {
+    Vec3 boxVectors[3];
+    system.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
+    double tol = force.getEwaldErrorTolerance();
+    alpha = (1.0/force.getCutoffDistance())*std::sqrt(-log(2.0*tol));
+    xsize = (int) ceil(alpha*boxVectors[0][0]/pow(0.5*tol, 0.2));
+    ysize = (int) ceil(alpha*boxVectors[1][1]/pow(0.5*tol, 0.2));
+    zsize = (int) ceil(alpha*boxVectors[2][2]/pow(0.5*tol, 0.2));
+}
+
+double NonbondedForceImpl::findZero(const NonbondedForceImpl::ErrorFunction& f, int initialGuess) {
+    int arg = initialGuess;
+    double value = f.getValue(arg);
+    if (value > 0.0) {
+        while (value > 0.0 && arg > 0)
+            value = f.getValue(--arg);
+        return arg+1;
+    }
+    while (value < 0.0)
+        value = f.getValue(++arg);
+    return arg;
+}

platforms/cuda/src/CudaKernels.cpp

@@ -28,6 +28,7 @@
 #include "openmm/LangevinIntegrator.h"
 #include "openmm/Context.h"
 #include "openmm/internal/ContextImpl.h"
+#include "openmm/internal/NonbondedForceImpl.h"
 #include "kernels/gputypes.h"
 #include "kernels/cudaKernels.h"
 #include <cmath>
@@ -343,29 +344,17 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
             method = PERIODIC;
         }
         if (force.getNonbondedMethod() == NonbondedForce::Ewald || force.getNonbondedMethod() == NonbondedForce::PME) {
-            double ewaldErrorTol = force.getEwaldErrorTolerance();
-            double alpha = (1.0/force.getCutoffDistance())*std::sqrt(-std::log(ewaldErrorTol));
-            double mx = boxVectors[0][0]/force.getCutoffDistance();
-            double my = boxVectors[1][1]/force.getCutoffDistance();
-            double mz = boxVectors[2][2]/force.getCutoffDistance();
-            double pi = 3.1415926535897932385;
-            int kmaxx = (int)std::ceil(-(mx/pi)*std::log(ewaldErrorTol));
-            int kmaxy = (int)std::ceil(-(my/pi)*std::log(ewaldErrorTol));
-            int kmaxz = (int)std::ceil(-(mz/pi)*std::log(ewaldErrorTol));
             if (force.getNonbondedMethod() == NonbondedForce::Ewald) {
-                if (kmaxx%2 == 0)
-                    kmaxx++;
-                if (kmaxy%2 == 0)
-                    kmaxy++;
-                if (kmaxz%2 == 0)
-                    kmaxz++;
+                double alpha;
+                int kmaxx, kmaxy, kmaxz;
+                NonbondedForceImpl::calcEwaldParameters(system, force, alpha, kmaxx, kmaxy, kmaxz);
                 gpuSetEwaldParameters(gpu, (float) alpha, kmaxx, kmaxy, kmaxz);
                 method = EWALD;
             }
             else {
-                int gridSizeX = -0.5*kmaxx*std::log(ewaldErrorTol);
-                int gridSizeY = -0.5*kmaxy*std::log(ewaldErrorTol);
-                int gridSizeZ = -0.5*kmaxz*std::log(ewaldErrorTol);
+                double alpha;
+                int gridSizeX, gridSizeY, gridSizeZ;
+                NonbondedForceImpl::calcPMEParameters(system, force, alpha, gridSizeX, gridSizeY, gridSizeZ);
                 gpuSetPMEParameters(gpu, (float) alpha, gridSizeX, gridSizeY, gridSizeZ);
                 method = PARTICLE_MESH_EWALD;
             }

platforms/cuda/tests/TestCudaEwald.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Portions copyright (c) 2008-2009 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -119,7 +119,7 @@ void testEwaldPME() {
     Context cudaContext2(system, integrator, cuda);
     cudaContext2.setPositions(positions);
     
-    tol = 1e-4;
+    tol = 1e-3;
     State cudaState2 = cudaContext2.getState(State::Energy);
     ASSERT_EQUAL_TOL(norm, (cudaState2.getPotentialEnergy()-cudaState.getPotentialEnergy())/delta, tol)
 
@@ -157,10 +157,9 @@ void testEwaldPME() {
     Context cudaContext3(system, integrator, cuda);
     cudaContext3.setPositions(positions);
      
-    tol = 1e-4;
+    tol = 1e-3;
     State cudaState3 = cudaContext3.getState(State::Energy);
     ASSERT_EQUAL_TOL(norm, (cudaState3.getPotentialEnergy()-cudaState.getPotentialEnergy())/delta, tol)
-
 }
 
 void testEwald2Ions() {
@@ -191,11 +190,62 @@ void testEwald2Ions() {
     ASSERT_EQUAL_TOL(-217.276, state.getPotentialEnergy(), 0.01/*10*TOL*/);
 }
 
+void testErrorTolerance(NonbondedForce::NonbondedMethod method) {
+    // Create a cloud of random point charges.
+
+    const int numParticles = 51;
+    const double boxWidth = 5.0;
+    System system;
+    system.setPeriodicBoxVectors(Vec3(boxWidth, 0, 0), Vec3(0, boxWidth, 0), Vec3(0, 0, boxWidth));
+    NonbondedForce* force = new NonbondedForce();
+    system.addForce(force);
+    vector<Vec3> positions(numParticles);
+    init_gen_rand(0);
+    for (int i = 0; i < numParticles; i++) {
+        system.addParticle(1.0);
+        force->addParticle(-1.0+i*2.0/(numParticles-1), 1.0, 0.0);
+        positions[i] = Vec3(boxWidth*genrand_real2(), boxWidth*genrand_real2(), boxWidth*genrand_real2());
+    }
+    force->setNonbondedMethod(method);
+    CudaPlatform platform;
+    VerletIntegrator integrator(0.01);
+
+    // For various values of the cutoff and error tolerance, see if the actual error is reasonable.
+
+    for (double cutoff = 1.0; cutoff < boxWidth/2; cutoff += 0.2) {
+        force->setCutoffDistance(cutoff);
+        vector<Vec3> refForces;
+        double norm = 0.0;
+        for (double tol = 5e-5; tol < 1e-3; tol *= 2.0) {
+            force->setEwaldErrorTolerance(tol);
+            Context context(system, integrator, platform);
+            context.setPositions(positions);
+            State state = context.getState(State::Forces);
+            if (refForces.size() == 0) {
+                refForces = state.getForces();
+                for (int i = 0; i < numParticles; i++)
+                    norm += refForces[i].dot(refForces[i]);
+                norm = sqrt(norm);
+            }
+            else {
+                double diff = 0.0;
+                for (int i = 0; i < numParticles; i++) {
+                    Vec3 delta = refForces[i]-state.getForces()[i];
+                    diff += delta.dot(delta);
+                }
+                diff = sqrt(diff)/norm;
+                ASSERT(diff < 5*tol);
+            }
+        }
+    }
+}
 
 int main() {
     try {
      testEwaldPME();
 //     testEwald2Ions();
+     testErrorTolerance(NonbondedForce::Ewald);
+     testErrorTolerance(NonbondedForce::PME);
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;

platforms/reference/src/ReferenceKernels.cpp

@@ -50,6 +50,7 @@
 #include "openmm/CMMotionRemover.h"
 #include "openmm/System.h"
 #include "openmm/internal/ContextImpl.h"
+#include "openmm/internal/NonbondedForceImpl.h"
 #include "openmm/Integrator.h"
 #include "SimTKUtilities/SimTKOpenMMUtilities.h"
 #include "lepton/CustomFunction.h"
@@ -470,29 +471,15 @@ void ReferenceCalcNonbondedForceKernel::initialize(const System& system, const N
         neighborList = NULL;
     else
         neighborList = new NeighborList();
-    if (nonbondedMethod == Ewald || nonbondedMethod == PME) {
-        RealOpenMM ewaldErrorTol = (RealOpenMM) force.getEwaldErrorTolerance();
-        ewaldAlpha = (RealOpenMM) (std::sqrt(-std::log(ewaldErrorTol))/nonbondedCutoff);
-        RealOpenMM mx = periodicBoxSize[0]/nonbondedCutoff;
-        RealOpenMM my = periodicBoxSize[1]/nonbondedCutoff;
-        RealOpenMM mz = periodicBoxSize[2]/nonbondedCutoff;
-        RealOpenMM pi = (RealOpenMM) 3.1415926535897932385;
-        kmax[0] = (int)std::ceil(-(mx/pi)*std::log(ewaldErrorTol));
-        kmax[1] = (int)std::ceil(-(my/pi)*std::log(ewaldErrorTol));
-        kmax[2] = (int)std::ceil(-(mz/pi)*std::log(ewaldErrorTol));
     if (nonbondedMethod == Ewald) {
-            if (kmax[0]%2 == 0)
-                kmax[0]++;
-            if (kmax[1]%2 == 0)
-                kmax[1]++;
-            if (kmax[2]%2 == 0)
-                kmax[2]++;
-        }
-        else {
-            gridSize[0] = -0.5*kmax[0]*std::log(ewaldErrorTol);
-            gridSize[1] = -0.5*kmax[1]*std::log(ewaldErrorTol);
-            gridSize[2] = -0.5*kmax[2]*std::log(ewaldErrorTol);
+        double alpha;
+        NonbondedForceImpl::calcEwaldParameters(system, force, alpha, kmax[0], kmax[1], kmax[2]);
+        ewaldAlpha = alpha;
     }
+    else if (nonbondedMethod == PME) {
+        double alpha;
+        NonbondedForceImpl::calcPMEParameters(system, force, alpha, gridSize[0], gridSize[1], gridSize[2]);
+        ewaldAlpha = alpha;
     }
     rfDielectric = (RealOpenMM)force.getReactionFieldDielectric();
 }

platforms/reference/tests/TestReferenceEwald.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008 Stanford University and the Authors.           *
+ * Portions copyright (c) 2008-2009 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -40,6 +40,7 @@
 #include "openmm/System.h"
 #include "openmm/VerletIntegrator.h"
 #include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
+#include "../src/sfmt/SFMT.h"
 #include "openmm/HarmonicBondForce.h"
 #include <iostream>
 #include <vector>
@@ -138,7 +139,7 @@ void testEwaldPME() {
 
 //   (1)   CHECK EXACT VALUE OF EWALD ENERGY (Against Gromacs output)
 
-    tol = 1e-5;
+    tol = 1e-4;
     ASSERT_EQUAL_TOL(-3.82047e+05, state1.getPotentialEnergy(), tol);
 
 //   (2)   CHECK WHETHER THE EWALD FORCES ARE THE SAME AS THE GROMACS OUTPUT
@@ -303,6 +304,55 @@ void testWaterSystem() {
 
 }
 
+void testErrorTolerance(NonbondedForce::NonbondedMethod method) {
+    // Create a cloud of random point charges.
+
+    const int numParticles = 51;
+    const double boxWidth = 5.0;
+    System system;
+    system.setPeriodicBoxVectors(Vec3(boxWidth, 0, 0), Vec3(0, boxWidth, 0), Vec3(0, 0, boxWidth));
+    NonbondedForce* force = new NonbondedForce();
+    system.addForce(force);
+    vector<Vec3> positions(numParticles);
+    init_gen_rand(0);
+    for (int i = 0; i < numParticles; i++) {
+        system.addParticle(1.0);
+        force->addParticle(-1.0+i*2.0/(numParticles-1), 1.0, 0.0);
+        positions[i] = Vec3(boxWidth*genrand_real2(), boxWidth*genrand_real2(), boxWidth*genrand_real2());
+    }
+    force->setNonbondedMethod(method);
+    ReferencePlatform platform;
+    VerletIntegrator integrator(0.01);
+
+    // For various values of the cutoff and error tolerance, see if the actual error is reasonable.
+
+    for (double cutoff = 1.0; cutoff < boxWidth/2; cutoff += 0.2) {
+        force->setCutoffDistance(cutoff);
+        vector<Vec3> refForces;
+        double norm = 0.0;
+        for (double tol = 5e-5; tol < 1e-3; tol *= 2.0) {
+            force->setEwaldErrorTolerance(tol);
+            Context context(system, integrator, platform);
+            context.setPositions(positions);
+            State state = context.getState(State::Forces);
+            if (refForces.size() == 0) {
+                refForces = state.getForces();
+                for (int i = 0; i < numParticles; i++)
+                    norm += refForces[i].dot(refForces[i]);
+                norm = sqrt(norm);
+            }
+            else {
+                double diff = 0.0;
+                for (int i = 0; i < numParticles; i++) {
+                    Vec3 delta = refForces[i]-state.getForces()[i];
+                    diff += delta.dot(delta);
+                }
+                diff = sqrt(diff)/norm;
+                ASSERT(diff < 5*tol);
+            }
+        }
+    }
+}
 
 int main() {
     try {
@@ -310,6 +360,8 @@ int main() {
      testEwaldPME();
 //     testEwald2Ions();
 //     testWaterSystem();
+     testErrorTolerance(NonbondedForce::Ewald);
+     testErrorTolerance(NonbondedForce::PME);
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;