Merge pull request #363 from peastman/pme

Added option to explicitly set PME parameters

docs/usersguide/theory.rst

@@ -417,7 +417,8 @@ and the number of nodes in the mesh along each dimension as
    {n}_{\text{mesh}}=\frac{2\alpha d}{{3d}^{1/5}}
 
 
-where *d* is the width of the periodic box along that dimension.  (Note that
+where *d* is the width of the periodic box along that dimension.  Alternatively,
+the user may choose to explicitly set values for these parameters.  (Note that
 some Platforms may choose to use a larger value of :math:`n_\text{mesh}` than that
 given by this equation.  For example, some FFT implementations require the mesh
 size to be a multiple of certain small prime numbers, so a Platform might round

openmmapi/include/openmm/NonbondedForce.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2013 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2014 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -182,15 +182,41 @@ public:
      * which is acceptable.  This value is used to select the reciprocal space cutoff and separation
      * parameter so that the average error level will be less than the tolerance.  There is not a
      * rigorous guarantee that all forces on all atoms will be less than the tolerance, however.
+     * 
+     * For PME calculations, if setPMEParameters() is used to set alpha to something other than 0,
+     * this value is ignored.
      */
     double getEwaldErrorTolerance() const;
     /**
-     * Get the error tolerance for Ewald summation.  This corresponds to the fractional error in the forces
+     * Set the error tolerance for Ewald summation.  This corresponds to the fractional error in the forces
      * which is acceptable.  This value is used to select the reciprocal space cutoff and separation
      * parameter so that the average error level will be less than the tolerance.  There is not a
      * rigorous guarantee that all forces on all atoms will be less than the tolerance, however.
+     * 
+     * For PME calculations, if setPMEParameters() is used to set alpha to something other than 0,
+     * this value is ignored.
      */
     void setEwaldErrorTolerance(double tol);
+    /**
+     * Get 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.
+     * 
+     * @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;
+    /**
+     * 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.
+     * 
+     * @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 setPMEParameters(double alpha, int nx, int ny, int nz);
     /**
      * 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.
@@ -332,9 +358,9 @@ private:
     class ParticleInfo;
     class ExceptionInfo;
     NonbondedMethod nonbondedMethod;
-    double cutoffDistance, switchingDistance, rfDielectric, ewaldErrorTol;
+    double cutoffDistance, switchingDistance, rfDielectric, ewaldErrorTol, alpha;
     bool useSwitchingFunction, useDispersionCorrection;
-    int recipForceGroup;
+    int recipForceGroup, nx, ny, nz;
     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/src/NonbondedForce.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-2013 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2014 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -48,7 +48,7 @@ using std::stringstream;
 using std::vector;
 
 NonbondedForce::NonbondedForce() : nonbondedMethod(NoCutoff), cutoffDistance(1.0), switchingDistance(-1.0), rfDielectric(78.3),
-        ewaldErrorTol(5e-4), useSwitchingFunction(false), useDispersionCorrection(true), recipForceGroup(-1) {
+        ewaldErrorTol(5e-4), alpha(0.0), useSwitchingFunction(false), useDispersionCorrection(true), recipForceGroup(-1), nx(0), ny(0), nz(0) {
 }
 
 NonbondedForce::NonbondedMethod NonbondedForce::getNonbondedMethod() const {
@@ -95,11 +95,24 @@ double NonbondedForce::getEwaldErrorTolerance() const {
     return ewaldErrorTol;
 }
 
-void NonbondedForce::setEwaldErrorTolerance(double tol)
-{
+void NonbondedForce::setEwaldErrorTolerance(double tol) {
     ewaldErrorTol = tol;
 }
 
+void NonbondedForce::getPMEParameters(double& alpha, int& nx, int& ny, int& nz) const {
+    alpha = this->alpha;
+    nx = this->nx;
+    ny = this->ny;
+    nz = this->nz;
+}
+
+void NonbondedForce::setPMEParameters(double alpha, int nx, int ny, int nz) {
+    this->alpha = alpha;
+    this->nx = nx;
+    this->ny = ny;
+    this->nz = nz;
+}
+
 int NonbondedForce::addParticle(double charge, double sigma, double epsilon) {
     particles.push_back(ParticleInfo(charge, sigma, epsilon));
     return particles.size()-1;

openmmapi/src/NonbondedForceImpl.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-2010 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2014 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -149,16 +149,19 @@ void NonbondedForceImpl::calcEwaldParameters(const System& system, const Nonbond
 }
 
 void NonbondedForceImpl::calcPMEParameters(const System& system, const NonbondedForce& force, double& alpha, int& xsize, int& ysize, int& zsize) {
-    Vec3 boxVectors[3];
-    system.getDefaultPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
-    double tol = force.getEwaldErrorTolerance();
-    alpha = (1.0/force.getCutoffDistance())*std::sqrt(-log(2.0*tol));
-    xsize = (int) ceil(2*alpha*boxVectors[0][0]/(3*pow(tol, 0.2)));
-    ysize = (int) ceil(2*alpha*boxVectors[1][1]/(3*pow(tol, 0.2)));
-    zsize = (int) ceil(2*alpha*boxVectors[2][2]/(3*pow(tol, 0.2)));
-    xsize = max(xsize, 5);
-    ysize = max(ysize, 5);
-    zsize = max(zsize, 5);
+    force.getPMEParameters(alpha, xsize, ysize, zsize);
+    if (alpha == 0.0) {
+        Vec3 boxVectors[3];
+        system.getDefaultPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
+        double tol = force.getEwaldErrorTolerance();
+        alpha = (1.0/force.getCutoffDistance())*std::sqrt(-log(2.0*tol));
+        xsize = (int) ceil(2*alpha*boxVectors[0][0]/(3*pow(tol, 0.2)));
+        ysize = (int) ceil(2*alpha*boxVectors[1][1]/(3*pow(tol, 0.2)));
+        zsize = (int) ceil(2*alpha*boxVectors[2][2]/(3*pow(tol, 0.2)));
+        xsize = max(xsize, 5);
+        ysize = max(ysize, 5);
+        zsize = max(zsize, 5);
+    }
 }
 
 int NonbondedForceImpl::findZero(const NonbondedForceImpl::ErrorFunction& f, int initialGuess) {

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-2012 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2014 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -252,6 +252,54 @@ void testErrorTolerance(NonbondedForce::NonbondedMethod method) {
     }
 }
 
+void testPMEParameters() {
+    // Create a cloud of random point charges.
+
+    const int numParticles = 51;
+    const double boxWidth = 4.7;
+    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(-1.0+i*2.0/(numParticles-1), 1.0, 0.0);
+        positions[i] = Vec3(boxWidth*genrand_real2(sfmt), boxWidth*genrand_real2(sfmt), boxWidth*genrand_real2(sfmt));
+    }
+    force->setNonbondedMethod(NonbondedForce::PME);
+    
+    // Compute the energy with an error tolerance of 1e-3.
+
+    force->setEwaldErrorTolerance(1e-3);
+    VerletIntegrator integrator1(0.01);
+    Context context1(system, integrator1, platform);
+    context1.setPositions(positions);
+    double energy1 = context1.getState(State::Energy).getPotentialEnergy();
+    
+    // Try again with an error tolerance of 1e-4.
+
+    force->setEwaldErrorTolerance(1e-4);
+    VerletIntegrator integrator2(0.01);
+    Context context2(system, integrator2, platform);
+    context2.setPositions(positions);
+    double energy2 = context2.getState(State::Energy).getPotentialEnergy();
+    
+    // Now explicitly set the parameters.  These should match the values that were
+    // used for tolerance 1e-3.
+
+    force->setPMEParameters(2.49291157051793, 32, 32, 32);
+    VerletIntegrator integrator3(0.01);
+    Context context3(system, integrator3, platform);
+    context3.setPositions(positions);
+    double energy3 = context3.getState(State::Energy).getPotentialEnergy();
+    ASSERT_EQUAL_TOL(energy1, energy3, 1e-6);
+    ASSERT(fabs((energy1-energy2)/energy1) > 1e-5);
+}
+
 int main(int argc, char* argv[]) {
     try {
         if (argc > 1)
@@ -261,6 +309,7 @@ int main(int argc, char* argv[]) {
 //        testEwald2Ions();
         testErrorTolerance(NonbondedForce::Ewald);
         testErrorTolerance(NonbondedForce::PME);
+        testPMEParameters();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;

platforms/opencl/tests/TestOpenCLEwald.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-2010 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2014 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -252,6 +252,54 @@ void testErrorTolerance(NonbondedForce::NonbondedMethod method) {
     }
 }
 
+void testPMEParameters() {
+    // Create a cloud of random point charges.
+
+    const int numParticles = 51;
+    const double boxWidth = 4.7;
+    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(-1.0+i*2.0/(numParticles-1), 1.0, 0.0);
+        positions[i] = Vec3(boxWidth*genrand_real2(sfmt), boxWidth*genrand_real2(sfmt), boxWidth*genrand_real2(sfmt));
+    }
+    force->setNonbondedMethod(NonbondedForce::PME);
+    
+    // Compute the energy with an error tolerance of 1e-3.
+
+    force->setEwaldErrorTolerance(1e-3);
+    VerletIntegrator integrator1(0.01);
+    Context context1(system, integrator1, platform);
+    context1.setPositions(positions);
+    double energy1 = context1.getState(State::Energy).getPotentialEnergy();
+    
+    // Try again with an error tolerance of 1e-4.
+
+    force->setEwaldErrorTolerance(1e-4);
+    VerletIntegrator integrator2(0.01);
+    Context context2(system, integrator2, platform);
+    context2.setPositions(positions);
+    double energy2 = context2.getState(State::Energy).getPotentialEnergy();
+    
+    // Now explicitly set the parameters.  These should match the values that were
+    // used for tolerance 1e-3.
+
+    force->setPMEParameters(2.49291157051793, 32, 32, 32);
+    VerletIntegrator integrator3(0.01);
+    Context context3(system, integrator3, platform);
+    context3.setPositions(positions);
+    double energy3 = context3.getState(State::Energy).getPotentialEnergy();
+    ASSERT_EQUAL_TOL(energy1, energy3, 1e-6);
+    ASSERT(fabs((energy1-energy2)/energy1) > 1e-5);
+}
+
 int main(int argc, char* argv[]) {
     try {
         if (argc > 1)
@@ -261,6 +309,7 @@ int main(int argc, char* argv[]) {
 //        testEwald2Ions();
         testErrorTolerance(NonbondedForce::Ewald);
         testErrorTolerance(NonbondedForce::PME);
+        testPMEParameters();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;

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-2010 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2014 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -354,6 +354,55 @@ void testErrorTolerance(NonbondedForce::NonbondedMethod method) {
     }
 }
 
+void testPMEParameters() {
+    // Create a cloud of random point charges.
+
+    const int numParticles = 51;
+    const double boxWidth = 4.7;
+    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(-1.0+i*2.0/(numParticles-1), 1.0, 0.0);
+        positions[i] = Vec3(boxWidth*genrand_real2(sfmt), boxWidth*genrand_real2(sfmt), boxWidth*genrand_real2(sfmt));
+    }
+    force->setNonbondedMethod(NonbondedForce::PME);
+    ReferencePlatform platform;
+    
+    // Compute the energy with an error tolerance of 1e-3.
+
+    force->setEwaldErrorTolerance(1e-3);
+    VerletIntegrator integrator1(0.01);
+    Context context1(system, integrator1, platform);
+    context1.setPositions(positions);
+    double energy1 = context1.getState(State::Energy).getPotentialEnergy();
+    
+    // Try again with an error tolerance of 1e-4.
+
+    force->setEwaldErrorTolerance(1e-4);
+    VerletIntegrator integrator2(0.01);
+    Context context2(system, integrator2, platform);
+    context2.setPositions(positions);
+    double energy2 = context2.getState(State::Energy).getPotentialEnergy();
+    
+    // Now explicitly set the parameters.  These should match the values that were
+    // used for tolerance 1e-3.
+
+    force->setPMEParameters(2.49291157051793, 32, 32, 32);
+    VerletIntegrator integrator3(0.01);
+    Context context3(system, integrator3, platform);
+    context3.setPositions(positions);
+    double energy3 = context3.getState(State::Energy).getPotentialEnergy();
+    ASSERT_EQUAL_TOL(energy1, energy3, 1e-6);
+    ASSERT(fabs((energy1-energy2)/energy1) > 1e-5);
+}
+
 int main() {
     try {
      testEwaldExact();
@@ -362,6 +411,7 @@ int main() {
 //     testWaterSystem();
      testErrorTolerance(NonbondedForce::Ewald);
      testErrorTolerance(NonbondedForce::PME);
+     testPMEParameters();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;