Implemented cutoffs for 1-4 nonbonded interactions

platforms/reference/src/ReferenceKernels.cpp

@@ -197,6 +197,8 @@ void ReferenceCalcStandardMMForceFieldKernel::executeForces(const Stream& positi
         clj.setPeriodic(periodicBoxSize);
     clj.calculatePairIxn(numAtoms, posData, atomParamArray, exclusionArray, 0, forceData, 0, 0);
     ReferenceLJCoulomb14 nonbonded14;
+    if (nonbondedMethod != NoCutoff)
+        nonbonded14.setUseCutoff(nonbondedCutoff, 78.3);
     refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, 0, 0, 0, nonbonded14);
 }
 
@@ -233,6 +235,8 @@ double ReferenceCalcStandardMMForceFieldKernel::executeEnergy(const Stream& posi
         clj.setPeriodic(periodicBoxSize);
     clj.calculatePairIxn(numAtoms, posData, atomParamArray, exclusionArray, 0, forceData, 0, &energy);
     ReferenceLJCoulomb14 nonbonded14;
+    if (nonbondedMethod != NoCutoff)
+        nonbonded14.setUseCutoff(nonbondedCutoff, 78.3);
     for (int i = 0; i < arraySize; ++i)
         energyArray[i] = 0;
     refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, energyArray, 0, &energy, nonbonded14);

platforms/reference/src/SimTKReference/ReferenceLJCoulomb14.cpp

@@ -37,7 +37,7 @@
 
    --------------------------------------------------------------------------------------- */
 
-ReferenceLJCoulomb14::ReferenceLJCoulomb14( ){
+ReferenceLJCoulomb14::ReferenceLJCoulomb14( ) : cutoff(false) {
 
    // ---------------------------------------------------------------------------------------
 
@@ -63,6 +63,27 @@ ReferenceLJCoulomb14::~ReferenceLJCoulomb14( ){
 
 }
 
+  /**---------------------------------------------------------------------------------------
+
+     Set the force to use a cutoff.
+
+     @param distance            the cutoff distance
+     @param solventDielectric   the dielectric constant of the bulk solvent
+
+     @return ReferenceForce::DefaultReturn
+
+     --------------------------------------------------------------------------------------- */
+
+  int ReferenceLJCoulomb14::setUseCutoff( RealOpenMM distance, RealOpenMM solventDielectric ) {
+    
+    cutoff = true;
+    cutoffDistance = distance;
+    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);
+            
+    return ReferenceForce::DefaultReturn;
+  }
+  
 /**---------------------------------------------------------------------------------------
 
    Calculate parameters for LJ 1-4 ixn
@@ -171,12 +192,18 @@ int ReferenceLJCoulomb14::calculateBondIxn( int* atomIndices, RealOpenMM** atomC
    int atomBIndex = atomIndices[1];
    ReferenceForce::getDeltaR( atomCoordinates[atomBIndex], atomCoordinates[atomAIndex], deltaR[0] );  
 
+   if (cutoff && deltaR[0][ReferenceForce::RIndex] > cutoffDistance)
+       return ReferenceForce::DefaultReturn;
+   RealOpenMM r2        = deltaR[0][ReferenceForce::R2Index];
    RealOpenMM inverseR  = one/(deltaR[0][ReferenceForce::RIndex]);
    RealOpenMM sig2      = inverseR*parameters[0];
               sig2     *= sig2;
    RealOpenMM sig6      = sig2*sig2*sig2;
 
    RealOpenMM dEdR      = parameters[1]*( twelve*sig6 - six )*sig6;
+              if (cutoff)
+                  dEdR += parameters[2]*(inverseR-2.0*krf*r2);
+              else
                   dEdR += parameters[2]*inverseR;
               dEdR     *= inverseR*inverseR;
 
@@ -188,7 +215,11 @@ int ReferenceLJCoulomb14::calculateBondIxn( int* atomIndices, RealOpenMM** atomC
       forces[atomBIndex][ii] -= force;
    }
 
-   RealOpenMM energy = parameters[1]*( sig6 - one )*sig6 + parameters[2]*inverseR;
+   RealOpenMM energy = parameters[1]*( sig6 - one )*sig6;
+   if (cutoff)
+       energy += parameters[2]*(inverseR+krf*r2-crf);
+   else
+       energy += parameters[2]*inverseR;
 
    // accumulate energies
 

platforms/reference/src/SimTKReference/ReferenceLJCoulomb14.h

@@ -33,6 +33,10 @@ class ReferenceLJCoulomb14 : public ReferenceBondIxn {
 
    private:
 
+        bool cutoff;
+        RealOpenMM cutoffDistance;
+        RealOpenMM krf, crf;
+
    public:
 
       /**---------------------------------------------------------------------------------------
@@ -51,6 +55,19 @@ class ReferenceLJCoulomb14 : public ReferenceBondIxn {
 
        ~ReferenceLJCoulomb14( );
 
+      /**---------------------------------------------------------------------------------------
+      
+         Set the force to use a cutoff.
+      
+         @param distance            the cutoff distance
+         @param solventDielectric   the dielectric constant of the bulk solvent
+      
+         @return ReferenceForce::DefaultReturn
+      
+         --------------------------------------------------------------------------------------- */
+      
+      int setUseCutoff( RealOpenMM distance, RealOpenMM solventDielectric );
+       
       /**---------------------------------------------------------------------------------------
       
          Calculate parameters for LJ 1-4 ixn

platforms/reference/tests/TestReferenceStandardMMForceField.cpp

@@ -298,6 +298,83 @@ void testCutoff() {
     ASSERT_EQUAL_TOL(energy1+energy2, state.getPotentialEnergy(), TOL);
 }
 
+void testCutoff14() {
+    ReferencePlatform platform;
+    System system(5, 0);
+    VerletIntegrator integrator(0.01);
+    StandardMMForceField* forceField = new StandardMMForceField(5, 4, 0, 0, 0);
+    forceField->setBondParameters(0, 0, 1, 1, 0);
+    forceField->setBondParameters(1, 1, 2, 1, 0);
+    forceField->setBondParameters(2, 2, 3, 1, 0);
+    forceField->setBondParameters(3, 3, 4, 1, 0);
+    forceField->setNonbondedMethod(StandardMMForceField::CutoffNonPeriodic);
+    const double cutoff = 3.5;
+    forceField->setCutoffDistance(cutoff);
+    system.addForce(forceField);
+    OpenMMContext context(system, integrator, platform);
+    vector<Vec3> positions(5);
+    positions[0] = Vec3(0, 0, 0);
+    positions[1] = Vec3(1, 0, 0);
+    positions[2] = Vec3(2, 0, 0);
+    positions[3] = Vec3(3, 0, 0);
+    positions[4] = Vec3(4, 0, 0);
+    for (int i = 1; i < 5; ++i) {
+ 
+        // Test LJ forces
+        
+        forceField->setAtomParameters(0, 0, 1.5, 1);
+        for (int j = 1; j < 5; ++j)
+            forceField->setAtomParameters(j, 0, 1.5, 0);
+        forceField->setAtomParameters(i, 0, 1.5, 1);
+        context.reinitialize();
+        context.setPositions(positions);
+        State state = context.getState(State::Forces | State::Energy);
+        const vector<Vec3>& forces = state.getForces();
+        double r = positions[i][0];
+        double x = 1.5/r;
+        double e = 1.0;
+        double force = 4.0*e*(12*std::pow(x, 12.0)-6*std::pow(x, 6.0))/r;
+        double energy = 4.0*e*(std::pow(x, 12.0)-std::pow(x, 6.0));
+        if (i == 3) {
+            force *= 0.5;
+            energy *= 0.5;
+        }
+        if (i < 3 || r > cutoff) {
+            force = 0;
+            energy = 0;
+        }
+        ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces[0], TOL);
+        ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces[i], TOL);
+        ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
+
+        // Test Coulomb forces
+        
+        const double q = 0.7;
+        forceField->setAtomParameters(0, q, 1.5, 0);
+        forceField->setAtomParameters(i, q, 1.5, 0);
+        context.reinitialize();
+        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);
+        energy = 138.935485*q*q*(1.0/r+krf*r*r-crf);
+        if (i == 3) {
+            force /= 1.2;
+            energy /= 1.2;
+        }
+        if (i < 3 || r > cutoff) {
+            force = 0;
+            energy = 0;
+        }
+        ASSERT_EQUAL_VEC(Vec3(-force, 0, 0), forces2[0], TOL);
+        ASSERT_EQUAL_VEC(Vec3(force, 0, 0), forces2[i], TOL);
+        ASSERT_EQUAL_TOL(energy, state.getPotentialEnergy(), TOL);
+    }
+}
+
 void testPeriodic() {
     ReferencePlatform platform;
     System system(3, 0);
@@ -340,6 +417,7 @@ int main() {
         testLJ();
         testExclusionsAnd14();
         testCutoff();
+        testCutoff14();
         testPeriodic();
     }
     catch(const exception& e) {