Allow the user to specify the reaction field dielectric constant

openmmapi/include/openmm/NonbondedForce.h

@@ -133,6 +133,14 @@ public:
      * is NoCutoff, this value will have no effect.
      */
     void setCutoffDistance(double distance);
+    /**
+     * Get the dielectric constant to use for the solvent in the reaction field approximation.
+     */
+    double getReactionFieldDielectric() const;
+    /**
+     * Set the dielectric constant to use for the solvent in the reaction field approximation.
+     */
+    void setReactionFieldDielectric(double dielectric);
     /**
      * Get 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
@@ -264,7 +272,7 @@ private:
     class ParticleInfo;
     class ExceptionInfo;
     NonbondedMethod nonbondedMethod;
-    double cutoffDistance, ewaldErrorTol;
+    double cutoffDistance, rfDielectric, ewaldErrorTol;
     Vec3 periodicBoxVectors[3];
     void addExclusionsToSet(const std::vector<std::set<int> >& bonded12, std::set<int>& exclusions, int baseParticle, int fromParticle, int currentLevel) const;
 

openmmapi/src/NonbondedForce.cpp

@@ -46,7 +46,7 @@ using std::string;
 using std::stringstream;
 using std::vector;
 
-NonbondedForce::NonbondedForce() : nonbondedMethod(NoCutoff), cutoffDistance(1.0), ewaldErrorTol(1e-4) {
+NonbondedForce::NonbondedForce() : nonbondedMethod(NoCutoff), cutoffDistance(1.0), rfDielectric(78.3), ewaldErrorTol(1e-4) {
     periodicBoxVectors[0] = Vec3(2, 0, 0);
     periodicBoxVectors[1] = Vec3(0, 2, 0);
     periodicBoxVectors[2] = Vec3(0, 0, 2);
@@ -68,6 +68,14 @@ void NonbondedForce::setCutoffDistance(double distance) {
     cutoffDistance = distance;
 }
 
+double NonbondedForce::getReactionFieldDielectric() const {
+    return rfDielectric;
+}
+
+void NonbondedForce::setReactionFieldDielectric(double dielectric) {
+    rfDielectric = dielectric;
+}
+
 double NonbondedForce::getEwaldErrorTolerance() const {
     return ewaldErrorTol;
 }

platforms/cuda/src/CudaKernels.cpp

@@ -290,7 +290,7 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
         }
         CudaNonbondedMethod method = NO_CUTOFF;
         if (force.getNonbondedMethod() != NonbondedForce::NoCutoff) {
-            gpuSetNonbondedCutoff(gpu, (float)force.getCutoffDistance(), 78.3f);
+            gpuSetNonbondedCutoff(gpu, (float)force.getCutoffDistance(), force.getReactionFieldDielectric());
             method = CUTOFF;
         }
         if (force.getNonbondedMethod() == NonbondedForce::CutoffPeriodic) {

platforms/cuda/tests/TestCudaNonbondedForce.cpp

@@ -203,6 +203,8 @@ void testCutoff() {
     forceField->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
     const double cutoff = 2.9;
     forceField->setCutoffDistance(cutoff);
+    const double eps = 50.0;
+    forceField->setReactionFieldDielectric(eps);
     system.addForce(forceField);
     OpenMMContext context(system, integrator, platform);
     vector<Vec3> positions(3);
@@ -212,7 +214,6 @@ void testCutoff() {
     context.setPositions(positions);
     State state = context.getState(State::Forces | State::Energy);
     const vector<Vec3>& forces = state.getForces();
-    const double eps = 78.3;
     const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0);
     const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0);
     const double force1 = 138.935485*(1.0)*(0.25-2.0*krf*2.0);
@@ -237,6 +238,8 @@ void testCutoff14() {
     }
     const double cutoff = 3.5;
     nonbonded->setCutoffDistance(cutoff);
+    const double eps = 30.0;
+    nonbonded->setReactionFieldDielectric(eps);
     vector<pair<int, int> > bonds;
     bonds.push_back(pair<int, int>(0, 1));
     bonds.push_back(pair<int, int>(1, 2));
@@ -303,7 +306,6 @@ void testCutoff14() {
         context.setPositions(positions);
         state = context.getState(State::Forces | State::Energy);
         const vector<Vec3>& forces2 = state.getForces();
-        const double eps = 78.3;
         const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0);
         const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0);
         force = 138.935485*q*q*(1.0/(r*r)-2.0*krf*r);

platforms/reference/src/ReferenceKernels.cpp

@@ -410,6 +410,7 @@ void ReferenceCalcNonbondedForceKernel::initialize(const System& system, const N
         if (kmax[2]%2 == 0)
             kmax[2]++;
     }
+    rfDielectric = force.getReactionFieldDielectric();
 }
 
 void ReferenceCalcNonbondedForceKernel::executeForces(OpenMMContextImpl& context) {
@@ -421,7 +422,7 @@ void ReferenceCalcNonbondedForceKernel::executeForces(OpenMMContextImpl& context
     bool pme  = (nonbondedMethod == PME);
     if (nonbondedMethod != NoCutoff) {
         computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, (periodic || ewald || pme) ? periodicBoxSize : NULL, nonbondedCutoff, 0.0);
-        clj.setUseCutoff(nonbondedCutoff, *neighborList, 78.3f);
+        clj.setUseCutoff(nonbondedCutoff, *neighborList, rfDielectric);
     }
     if (periodic||ewald||pme)
         clj.setPeriodic(periodicBoxSize);
@@ -433,7 +434,7 @@ void ReferenceCalcNonbondedForceKernel::executeForces(OpenMMContextImpl& context
     ReferenceBondForce refBondForce;
     ReferenceLJCoulomb14 nonbonded14;
     if (nonbondedMethod != NoCutoff)
-        nonbonded14.setUseCutoff(nonbondedCutoff, 78.3f);
+        nonbonded14.setUseCutoff(nonbondedCutoff, rfDielectric);
     refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, 0, 0, 0, nonbonded14);
 }
 
@@ -447,7 +448,7 @@ double ReferenceCalcNonbondedForceKernel::executeEnergy(OpenMMContextImpl& conte
     bool pme  = (nonbondedMethod == PME);
     if (nonbondedMethod != NoCutoff) {
         computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, (periodic || ewald || pme) ? periodicBoxSize : NULL, nonbondedCutoff, 0.0);
-        clj.setUseCutoff(nonbondedCutoff, *neighborList, 78.3f);
+        clj.setUseCutoff(nonbondedCutoff, *neighborList, rfDielectric);
     }
     if (periodic || ewald || pme)
         clj.setPeriodic(periodicBoxSize);
@@ -459,7 +460,7 @@ double ReferenceCalcNonbondedForceKernel::executeEnergy(OpenMMContextImpl& conte
     ReferenceBondForce refBondForce;
     ReferenceLJCoulomb14 nonbonded14;
     if (nonbondedMethod != NoCutoff)
-        nonbonded14.setUseCutoff(nonbondedCutoff, 78.3f);
+        nonbonded14.setUseCutoff(nonbondedCutoff, rfDielectric);
     RealOpenMM* energyArray = new RealOpenMM[num14];
     for (int i = 0; i < num14; ++i)
         energyArray[i] = 0;

platforms/reference/src/ReferenceKernels.h

@@ -267,7 +267,7 @@ private:
     int numParticles, num14;
     int **exclusionArray, **bonded14IndexArray;
     RealOpenMM **particleParamArray, **bonded14ParamArray;
-    RealOpenMM nonbondedCutoff, periodicBoxSize[3], ewaldAlpha;
+    RealOpenMM nonbondedCutoff, periodicBoxSize[3], rfDielectric, ewaldAlpha;
     int kmax[3];
     std::vector<std::set<int> > exclusions;
     NonbondedMethod nonbondedMethod;

platforms/reference/tests/TestReferenceNonbondedForce.cpp

@@ -199,6 +199,8 @@ void testCutoff() {
     forceField->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
     const double cutoff = 2.9;
     forceField->setCutoffDistance(cutoff);
+    const double eps = 50.0;
+    forceField->setReactionFieldDielectric(eps);
     system.addForce(forceField);
     OpenMMContext context(system, integrator, platform);
     vector<Vec3> positions(3);
@@ -208,7 +210,6 @@ void testCutoff() {
     context.setPositions(positions);
     State state = context.getState(State::Forces | State::Energy);
     const vector<Vec3>& forces = state.getForces();
-    const double eps = 78.3;
     const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0);
     const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0);
     const double force1 = 138.935485*(1.0)*(0.25-2.0*krf*2.0);
@@ -233,6 +234,8 @@ void testCutoff14() {
     nonbonded->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
     const double cutoff = 3.5;
     nonbonded->setCutoffDistance(cutoff);
+    const double eps = 30.0;
+    nonbonded->setReactionFieldDielectric(eps);
     vector<pair<int, int> > bonds;
     bonds.push_back(pair<int, int>(0, 1));
     bonds.push_back(pair<int, int>(1, 2));
@@ -299,7 +302,6 @@ void testCutoff14() {
         context.setPositions(positions);
         state = context.getState(State::Forces | State::Energy);
         const vector<Vec3>& forces2 = state.getForces();
-        const double eps = 78.3;
         const double krf = (1.0/(cutoff*cutoff*cutoff))*(eps-1.0)/(2.0*eps+1.0);
         const double crf = (1.0/cutoff)*(3.0*eps)/(2.0*eps+1.0);
         force = 138.935485*q*q*(1.0/(r*r)-2.0*krf*r);