Began implementing switching function for Lennard-Jones

openmmapi/include/openmm/CustomNonbondedForce.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-2012 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2013 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -179,11 +179,13 @@ public:
      */
     void setCutoffDistance(double distance);
     /**
-     * Get whether a switching is applied to the interaction.
+     * Get whether a switching function is applied to the interaction.  If the nonbonded method is set
+     * to NoCutoff, this option is ignored.
      */
     bool getUseSwitchingFunction() const;
     /**
-     * Set whether a switching is applied to the interaction.
+     * Set whether a switching function is applied to the interaction.  If the nonbonded method is set
+     * to NoCutoff, this option is ignored.
      */
     void setUseSwitchingFunction(bool use);
     /**

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-2012 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2013 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -65,6 +65,12 @@ namespace OpenMM {
  * This class provides a convenience method for this case called createExceptionsFromBonds().  You pass to it
  * a list of bonds and the scale factors to use for 1-4 interactions.  It identifies all pairs of particles which
  * are separated by 1, 2, or 3 bonds, then automatically creates exceptions for them.
+ * 
+ * When using a cutoff, by default Lennard-Jones interactions are sharply truncated at the cutoff distance.
+ * Optionally you can instead use a switching function to make the interaction smoothly go to zero over a finite
+ * distance range.  To enable this, call setUseSwitchingFunction().  You must also call setSwitchingDistance()
+ * to specify the distance at which the interaction should begin to decrease.  The switching distance must be
+ * less than the cutoff distance.
  */
 
 class OPENMM_EXPORT NonbondedForce : public Force {
@@ -138,6 +144,26 @@ public:
      * @param distance    the cutoff distance, measured in nm
      */
     void setCutoffDistance(double distance);
+    /**
+     * Get whether a switching function is applied to the Lennard-Jones interaction.  If the nonbonded method is set
+     * to NoCutoff, this option is ignored.
+     */
+    bool getUseSwitchingFunction() const;
+    /**
+     * Set whether a switching function is applied to the Lennard-Jones interaction.  If the nonbonded method is set
+     * to NoCutoff, this option is ignored.
+     */
+    void setUseSwitchingFunction(bool use);
+    /**
+     * Get the distance at which the switching function begins to reduce the Lennard-Jones interaction.  This must be
+     * less than the cutoff distance.
+     */
+    double getSwitchingDistance() const;
+    /**
+     * Set the distance at which the switching function begins to reduce the Lennard-Jones interaction.  This must be
+     * less than the cutoff distance.
+     */
+    void setSwitchingDistance(double distance);
     /**
      * Get the dielectric constant to use for the solvent in the reaction field approximation.
      */
@@ -301,8 +327,8 @@ private:
     class ParticleInfo;
     class ExceptionInfo;
     NonbondedMethod nonbondedMethod;
-    double cutoffDistance, rfDielectric, ewaldErrorTol;
-    bool useDispersionCorrection;
+    double cutoffDistance, switchingDistance, rfDielectric, ewaldErrorTol;
+    bool useSwitchingFunction, useDispersionCorrection;
     int recipForceGroup;
     void addExclusionsToSet(const std::vector<std::set<int> >& bonded12, std::set<int>& exclusions, int baseParticle, int fromParticle, int currentLevel) const;
     std::vector<ParticleInfo> particles;

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-2010 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2013 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -47,7 +47,8 @@ using std::string;
 using std::stringstream;
 using std::vector;
 
-NonbondedForce::NonbondedForce() : nonbondedMethod(NoCutoff), cutoffDistance(1.0), rfDielectric(78.3), ewaldErrorTol(5e-4), useDispersionCorrection(true), recipForceGroup(-1) {
+NonbondedForce::NonbondedForce() : nonbondedMethod(NoCutoff), cutoffDistance(1.0), switchingDistance(-1.0), rfDielectric(78.3),
+        ewaldErrorTol(5e-4), useSwitchingFunction(false), useDispersionCorrection(true), recipForceGroup(-1) {
 }
 
 NonbondedForce::NonbondedMethod NonbondedForce::getNonbondedMethod() const {
@@ -66,6 +67,22 @@ void NonbondedForce::setCutoffDistance(double distance) {
     cutoffDistance = distance;
 }
 
+bool NonbondedForce::getUseSwitchingFunction() const {
+    return useSwitchingFunction;
+}
+
+void NonbondedForce::setUseSwitchingFunction(bool use) {
+    useSwitchingFunction = use;
+}
+
+double NonbondedForce::getSwitchingDistance() const {
+    return switchingDistance;
+}
+
+void NonbondedForce::setSwitchingDistance(double distance) {
+    switchingDistance = distance;
+}
+
 double NonbondedForce::getReactionFieldDielectric() const {
     return rfDielectric;
 }

openmmapi/src/NonbondedForceImpl.cpp

@@ -57,6 +57,10 @@ void NonbondedForceImpl::initialize(ContextImpl& context) {
     const System& system = context.getSystem();
     if (owner.getNumParticles() != system.getNumParticles())
         throw OpenMMException("NonbondedForce must have exactly as many particles as the System it belongs to.");
+    if (owner.getUseSwitchingFunction()) {
+        if (owner.getSwitchingDistance() < 0 || owner.getSwitchingDistance() >= owner.getCutoffDistance())
+            throw OpenMMException("NonbondedForce: Switching distance must satisfy 0 <= r_switch < r_cutoff");
+    }
     vector<set<int> > exceptions(owner.getNumParticles());
     for (int i = 0; i < owner.getNumExceptions(); i++) {
         int particle1, particle2;

platforms/reference/src/ReferenceKernels.cpp

@@ -867,10 +867,15 @@ void ReferenceCalcNonbondedForceKernel::initialize(const System& system, const N
     }
     nonbondedMethod = CalcNonbondedForceKernel::NonbondedMethod(force.getNonbondedMethod());
     nonbondedCutoff = (RealOpenMM) force.getCutoffDistance();
-    if (nonbondedMethod == NoCutoff)
+    if (nonbondedMethod == NoCutoff) {
         neighborList = NULL;
-    else
+        useSwitchingFunction = false;
+    }
+    else {
         neighborList = new NeighborList();
+        useSwitchingFunction = force.getUseSwitchingFunction();
+        switchingDistance = force.getSwitchingDistance();
+    }
     if (nonbondedMethod == Ewald) {
         double alpha;
         NonbondedForceImpl::calcEwaldParameters(system, force, alpha, kmax[0], kmax[1], kmax[2]);
@@ -911,6 +916,8 @@ double ReferenceCalcNonbondedForceKernel::execute(ContextImpl& context, bool inc
         clj.setUseEwald(ewaldAlpha, kmax[0], kmax[1], kmax[2]);
     if (pme)
         clj.setUsePME(ewaldAlpha, gridSize);
+    if (useSwitchingFunction)
+        clj.setUseSwitchingFunction(switchingDistance);
     clj.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, forceData, 0, includeEnergy ? &energy : NULL, includeDirect, includeReciprocal);
     if (includeDirect) {
         ReferenceBondForce refBondForce;

platforms/reference/src/ReferenceKernels.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-2009 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2013 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -574,8 +574,9 @@ private:
     int numParticles, num14;
     int **exclusionArray, **bonded14IndexArray;
     RealOpenMM **particleParamArray, **bonded14ParamArray;
-    RealOpenMM nonbondedCutoff, rfDielectric, ewaldAlpha, dispersionCoefficient;
+    RealOpenMM nonbondedCutoff, switchingDistance, rfDielectric, ewaldAlpha, dispersionCoefficient;
     int kmax[3], gridSize[3];
+    bool useSwitchingFunction;
     std::vector<std::set<int> > exclusions;
     NonbondedMethod nonbondedMethod;
     NeighborList* neighborList;

platforms/reference/src/SimTKReference/ReferenceLJCoulombIxn.cpp

 
-/* Portions copyright (c) 2006 Stanford University and Simbios.
+/* Portions copyright (c) 2006-2013 Stanford University and Simbios.
  * Contributors: Pande Group
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -46,7 +46,7 @@ using OpenMM::RealVec;
 
    --------------------------------------------------------------------------------------- */
 
-ReferenceLJCoulombIxn::ReferenceLJCoulombIxn( ) : cutoff(false), periodic(false), ewald(false), pme(false) {
+ReferenceLJCoulombIxn::ReferenceLJCoulombIxn( ) : cutoff(false), useSwitch(false), periodic(false), ewald(false), pme(false) {
 
    // ---------------------------------------------------------------------------------------
 
@@ -91,6 +91,19 @@ ReferenceLJCoulombIxn::~ReferenceLJCoulombIxn( ){
     crf = (1.0/cutoffDistance)*(3.0*solventDielectric)/(2.0*solventDielectric+1.0);
   }
 
+/**---------------------------------------------------------------------------------------
+
+   Set the force to use a switching function on the Lennard-Jones interaction.
+
+   @param distance            the switching distance
+
+   --------------------------------------------------------------------------------------- */
+
+void ReferenceLJCoulombIxn::setUseSwitchingFunction( RealOpenMM distance ) {
+    useSwitch = true;
+    switchingDistance = distance;
+}
+
   /**---------------------------------------------------------------------------------------
 
      Set the force to use periodic boundary conditions.  This requires that a cutoff has
@@ -357,6 +370,12 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
        ReferenceForce::getDeltaRPeriodic( atomCoordinates[jj], atomCoordinates[ii], periodicBoxSize, deltaR[0] );
        RealOpenMM r         = deltaR[0][ReferenceForce::RIndex];
        RealOpenMM inverseR  = one/(deltaR[0][ReferenceForce::RIndex]);
+       RealOpenMM switchValue = 1, switchDeriv = 0;
+       if (useSwitch && r > switchingDistance) {
+           RealOpenMM t = (r-switchingDistance)/(cutoffDistance-switchingDistance);
+           switchValue = 1+t*t*t*(-10+t*(15-t*6));
+           switchDeriv = t*t*(-30+t*(60-t*30))/(cutoffDistance-switchingDistance);
+       }
        RealOpenMM alphaR    = alphaEwald * r;
 
 
@@ -368,7 +387,12 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
                   sig2     *= sig2;
        RealOpenMM sig6      = sig2*sig2*sig2;
        RealOpenMM eps       = atomParameters[ii][EpsIndex]*atomParameters[jj][EpsIndex];
-                  dEdR     += eps*( twelve*sig6 - six )*sig6*inverseR*inverseR;
+                  dEdR     += switchValue*eps*( twelve*sig6 - six )*sig6*inverseR*inverseR;
+       vdwEnergy = eps*(sig6-one)*sig6;
+       if (useSwitch) {
+           dEdR -= vdwEnergy*switchDeriv*inverseR;
+           vdwEnergy *= switchValue;
+       }
 
        // accumulate forces
 
@@ -381,7 +405,6 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<RealVec>
        // accumulate energies
 
        realSpaceEwaldEnergy        = (RealOpenMM) (ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR*erfc(alphaR));
-       vdwEnergy                   = eps*(sig6-one)*sig6;
 
        totalVdwEnergy             += vdwEnergy;
        totalRealSpaceEwaldEnergy  += realSpaceEwaldEnergy;
@@ -554,18 +577,36 @@ void ReferenceLJCoulombIxn::calculateOneIxn( int ii, int jj, vector<RealVec>& at
 
     RealOpenMM r2        = deltaR[0][ReferenceForce::R2Index];
     RealOpenMM inverseR  = one/(deltaR[0][ReferenceForce::RIndex]);
+    RealOpenMM switchValue = 1, switchDeriv = 0;
+    if (useSwitch) {
+        RealOpenMM r = deltaR[0][ReferenceForce::RIndex];
+        if (r > switchingDistance) {
+            RealOpenMM t = (r-switchingDistance)/(cutoffDistance-switchingDistance);
+            switchValue = 1+t*t*t*(-10+t*(15-t*6));
+            switchDeriv = t*t*(-30+t*(60-t*30))/(cutoffDistance-switchingDistance);
+        }
+    }
     RealOpenMM sig       = atomParameters[ii][SigIndex] +  atomParameters[jj][SigIndex];
     RealOpenMM sig2      = inverseR*sig;
                sig2     *= sig2;
     RealOpenMM sig6      = sig2*sig2*sig2;
 
     RealOpenMM eps       = atomParameters[ii][EpsIndex]*atomParameters[jj][EpsIndex];
-    RealOpenMM dEdR      = eps*( twelve*sig6 - six )*sig6;
+    RealOpenMM dEdR      = switchValue*eps*( twelve*sig6 - six )*sig6;
     if (cutoff)
         dEdR += (RealOpenMM) (ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*(inverseR-2.0f*krf*r2));
     else
         dEdR += (RealOpenMM) (ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR);
     dEdR     *= inverseR*inverseR;
+    RealOpenMM energy = eps*(sig6-one)*sig6;
+    if (useSwitch) {
+        dEdR -= energy*switchDeriv*inverseR;
+        energy *= switchValue;
+    }
+    if (cutoff)
+        energy += (RealOpenMM) (ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*(inverseR+krf*r2-crf));
+    else
+        energy += (RealOpenMM) (ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR);
 
     // accumulate forces
 
@@ -575,16 +616,8 @@ void ReferenceLJCoulombIxn::calculateOneIxn( int ii, int jj, vector<RealVec>& at
        forces[jj][kk]   -= force;
     }
 
-    RealOpenMM energy = 0.0;
-
     // accumulate energies
 
-    if( totalEnergy || energyByAtom ) {
-        if (cutoff)
-            energy = (RealOpenMM) (ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*(inverseR+krf*r2-crf));
-        else
-            energy = (RealOpenMM) (ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR);
-        energy += eps*(sig6-one)*sig6;
     if( totalEnergy )
        *totalEnergy += energy;
     if( energyByAtom ){
@@ -592,5 +625,4 @@ void ReferenceLJCoulombIxn::calculateOneIxn( int ii, int jj, vector<RealVec>& at
        energyByAtom[jj] += energy;
     }
   }
-  }
 

platforms/reference/src/SimTKReference/ReferenceLJCoulombIxn.h

 
-/* Portions copyright (c) 2006 Stanford University and Simbios.
+/* Portions copyright (c) 2006-2013 Stanford University and Simbios.
  * Contributors: Pande Group
  *
  * Permission is hereby granted, free of charge, to any person obtaining
@@ -35,12 +35,13 @@ class ReferenceLJCoulombIxn {
    private:
        
       bool cutoff;
+      bool useSwitch;
       bool periodic;
       bool ewald;
       bool pme;
       const OpenMM::NeighborList* neighborList;
       RealOpenMM periodicBoxSize[3];
-      RealOpenMM cutoffDistance;
+      RealOpenMM cutoffDistance, switchingDistance;
       RealOpenMM krf, crf;
       RealOpenMM alphaEwald;
       int numRx, numRy, numRz;
@@ -101,6 +102,16 @@ class ReferenceLJCoulombIxn {
       
       void setUseCutoff( RealOpenMM distance, const OpenMM::NeighborList& neighbors, RealOpenMM solventDielectric );
 
+      /**---------------------------------------------------------------------------------------
+      
+         Set the force to use a switching function on the Lennard-Jones interaction.
+      
+         @param distance            the switching distance
+      
+         --------------------------------------------------------------------------------------- */
+      
+      void setUseSwitchingFunction( RealOpenMM distance );
+      
       /**---------------------------------------------------------------------------------------
       
          Set the force to use periodic boundary conditions.  This requires that a cutoff has

platforms/reference/tests/TestReferenceNonbondedForce.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-2013 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -419,6 +419,62 @@ void testDispersionCorrection() {
     ASSERT_EQUAL_TOL(expected, energy1-energy2, 1e-4);
 }
 
+void testSwitchingFunction(NonbondedForce::NonbondedMethod method) {
+    ReferencePlatform platform;
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(6, 0, 0), Vec3(0, 6, 0), Vec3(0, 0, 6));
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    VerletIntegrator integrator(0.01);
+    NonbondedForce* nonbonded = new NonbondedForce();
+    nonbonded->addParticle(0, 1.2, 1);
+    nonbonded->addParticle(0, 1.4, 2);
+    nonbonded->setNonbondedMethod(method);
+    nonbonded->setCutoffDistance(2.0);
+    nonbonded->setUseSwitchingFunction(true);
+    nonbonded->setSwitchingDistance(1.5);
+    nonbonded->setUseDispersionCorrection(false);
+    system.addForce(nonbonded);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(0, 0, 0);
+    double eps = SQRT_TWO;
+    
+    // Compute the interaction at various distances.
+    
+    for (double r = 1.0; r < 2.5; r += 0.1) {
+        positions[1] = Vec3(r, 0, 0);
+        context.setPositions(positions);
+        State state = context.getState(State::Forces | State::Energy);
+        
+        // See if the energy is correct.
+        
+        double x = 1.3/r;
+        double expectedEnergy = 4.0*eps*(std::pow(x, 12.0)-std::pow(x, 6.0));
+        double switchValue;
+        if (r <= 1.5)
+            switchValue = 1;
+        else if (r >= 2.0)
+            switchValue = 0;
+        else {
+            double t = (r-1.5)/0.5;
+            switchValue = 1+t*t*t*(-10+t*(15-t*6));
+        }
+        ASSERT_EQUAL_TOL(switchValue*expectedEnergy, state.getPotentialEnergy(), TOL);
+        
+        // See if the force is the gradient of the energy.
+        
+        double delta = 1e-3;
+        positions[1] = Vec3(r-delta, 0, 0);
+        context.setPositions(positions);
+        double e1 = context.getState(State::Energy).getPotentialEnergy();
+        positions[1] = Vec3(r+delta, 0, 0);
+        context.setPositions(positions);
+        double e2 = context.getState(State::Energy).getPotentialEnergy();
+        ASSERT_EQUAL_TOL((e2-e1)/(2*delta), state.getForces()[0][0], 1e-3);
+    }
+}
+
 int main() {
     try {
         testCoulomb();
@@ -428,6 +484,8 @@ int main() {
         testCutoff14();
         testPeriodic();
         testDispersionCorrection();
+        testSwitchingFunction(NonbondedForce::CutoffNonPeriodic);
+        testSwitchingFunction(NonbondedForce::PME);
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;