Reference implementation of hard wall constraint for DrudeLangevinIntegrator

plugins/drude/openmmapi/include/openmm/DrudeLangevinIntegrator.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-2015 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -45,6 +45,10 @@ namespace OpenMM {
  * A second thermostat, typically with a much lower temperature, is applied to the relative internal
  * displacement of each pair.
  *
+ * This integrator can optionally set an upper limit on how far any Drude particle is ever allowed to
+ * get from its parent particle.  This can sometimes help to improve stability.  The limit is enforced
+ * with a hard wall constraint.
+ * 
  * This Integrator requires the System to include a DrudeForce, which it uses to identify the Drude
  * particles.
  */
@@ -129,6 +133,16 @@ public:
     void setDrudeFriction(double coeff) {
         drudeFriction = coeff;
     }
+    /**
+     * Get the maximum distance a Drude particle can ever move from its parent particle, measured in nm.  This is implemented
+     * with a hard wall constraint.  If this distance is set to 0 (the default), the hard wall constraint is omitted.
+     */
+    double getMaxDrudeDistance() const;
+    /**
+     * Set the maximum distance a Drude particle can ever move from its parent particle, measured in nm.  This is implemented
+     * with a hard wall constraint.  If this distance is set to 0 (the default), the hard wall constraint is omitted.
+     */
+    void setMaxDrudeDistance(double distance);
     /**
      * Get the random number seed.  See setRandomNumberSeed() for details.
      */
@@ -181,7 +195,7 @@ protected:
      */
     double computeKineticEnergy();
 private:
-    double temperature, friction, drudeTemperature, drudeFriction;
+    double temperature, friction, drudeTemperature, drudeFriction, maxDrudeDistance;
     int randomNumberSeed;
     Kernel kernel;
 };

plugins/drude/openmmapi/src/DrudeLangevinIntegrator.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-2015 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -46,11 +46,22 @@ DrudeLangevinIntegrator::DrudeLangevinIntegrator(double temperature, double fric
     setFriction(frictionCoeff);
     setDrudeTemperature(drudeTemperature);
     setDrudeFriction(drudeFrictionCoeff);
+    setMaxDrudeDistance(0);
     setStepSize(stepSize);
     setConstraintTolerance(1e-5);
     setRandomNumberSeed(0);
 }
 
+double DrudeLangevinIntegrator::getMaxDrudeDistance() const {
+    return maxDrudeDistance;
+}
+
+void DrudeLangevinIntegrator::setMaxDrudeDistance(double distance) {
+    if (distance < 0)
+        throw OpenMMException("setMaxDrudeDistance: Distance cannot be negative");
+    maxDrudeDistance = distance;
+}
+
 void DrudeLangevinIntegrator::initialize(ContextImpl& contextRef) {
     if (owner != NULL && &contextRef.getOwner() != owner)
         throw OpenMMException("This Integrator is already bound to a context");

plugins/drude/platforms/reference/src/ReferenceDrudeKernels.cpp

@@ -235,7 +235,6 @@ void ReferenceIntegrateDrudeLangevinStepKernel::initialize(const System& system,
     // Identify particle pairs and ordinary particles.
     
     set<int> particles;
-    vector<RealOpenMM> particleMass;
     for (int i = 0; i < system.getNumParticles(); i++) {
         particles.insert(i);
         double mass = system.getParticleMass(i);
@@ -325,6 +324,75 @@ void ReferenceIntegrateDrudeLangevinStepKernel::execute(ContextImpl& context, co
             pos[i] = xPrime[i];
         }
     }
+
+    // Apply hard wall constraints.
+
+    const RealOpenMM maxDrudeDistance = integrator.getMaxDrudeDistance();
+    if (maxDrudeDistance > 0) {
+        const RealOpenMM hardwallscaleDrude = sqrt(kTDrude);
+        for (int i = 0; i < (int) pairParticles.size(); i++) {
+            int p1 = pairParticles[i].first;
+            int p2 = pairParticles[i].second;
+            RealVec delta = pos[p1]-pos[p2];
+            RealOpenMM r = sqrt(delta.dot(delta));
+            RealOpenMM rInv = 1/r;
+            if (rInv*maxDrudeDistance < 1.0) {
+                // The constraint has been violated, so make the inter-particle distance "bounce"
+                // off the hard wall.
+                
+                if (rInv*maxDrudeDistance < 0.5)
+                    throw OpenMMException("Drude particle moved too far beyond hard wall constraint");
+                RealVec bondDir = delta*rInv;
+                RealVec vel1 = vel[p1];
+                RealVec vel2 = vel[p2];
+                RealOpenMM mass1 = particleMass[p1];
+                RealOpenMM mass2 = particleMass[p2];
+                RealOpenMM deltaR = r-maxDrudeDistance;
+                RealOpenMM deltaT = dt;
+                RealOpenMM dotvr1 = vel1.dot(bondDir);
+                RealVec vb1 = bondDir*dotvr1;
+                RealVec vp1 = vel1-vb1;
+                if (mass2 == 0) {
+                    // The parent particle is massless, so move only the Drude particle.
+
+                    if (dotvr1 != 0.0)
+                        deltaT = deltaR/abs(dotvr1);
+                    if (deltaT > dt)
+                        deltaT = dt;
+                    dotvr1 = -dotvr1*hardwallscaleDrude/(abs(dotvr1)*sqrt(mass1));
+                    RealOpenMM dr = -deltaR + deltaT*dotvr1;
+                    pos[p1] += bondDir*dr;
+                    vel[p1] = vp1 + bondDir*dotvr1;
+                }
+                else {
+                    // Move both particles.
+
+                    RealOpenMM invTotalMass = pairInvTotalMass[i];
+                    RealOpenMM dotvr2 = vel2.dot(bondDir);
+                    RealVec vb2 = bondDir*dotvr2;
+                    RealVec vp2 = vel2-vb2;
+                    RealOpenMM vbCMass = (mass1*dotvr1 + mass2*dotvr2)*invTotalMass;
+                    dotvr1 -= vbCMass;
+                    dotvr2 -= vbCMass;
+                    if (dotvr1 != dotvr2)
+                        deltaT = deltaR/abs(dotvr1-dotvr2);
+                    if (deltaT > dt)
+                        deltaT = dt;
+                    RealOpenMM vBond = hardwallscaleDrude/sqrt(mass1);
+                    dotvr1 = -dotvr1*vBond*mass2*invTotalMass/abs(dotvr1);
+                    dotvr2 = -dotvr2*vBond*mass1*invTotalMass/abs(dotvr2);
+                    RealOpenMM dr1 = -deltaR*mass2*invTotalMass + deltaT*dotvr1;
+                    RealOpenMM dr2 = deltaR*mass1*invTotalMass + deltaT*dotvr2;
+                    dotvr1 += vbCMass;
+                    dotvr2 += vbCMass;
+                    pos[p1] += bondDir*dr1;
+                    pos[p2] += bondDir*dr2;
+                    vel[p1] = vp1 + bondDir*dotvr1;
+                    vel[p2] = vp2 + bondDir*dotvr2;
+                }
+            }
+        }
+    }
     ReferenceVirtualSites::computePositions(context.getSystem(), pos);
     data.time += integrator.getStepSize();
     data.stepCount++;

plugins/drude/platforms/reference/src/ReferenceDrudeKernels.h

@@ -115,6 +115,7 @@ private:
     ReferencePlatform::PlatformData& data;
     std::vector<int> normalParticles;
     std::vector<std::pair<int, int> > pairParticles;
+    std::vector<double> particleMass;
     std::vector<double> particleInvMass;
     std::vector<double> pairInvTotalMass;
     std::vector<double> pairInvReducedMass;
@@ -158,6 +159,7 @@ private:
     std::vector<double> particleInvMass;
     lbfgsfloatval_t *minimizerPos;
     lbfgs_parameter_t minimizerParams;
+    double maxDrudeDistance;
 };
 
 } // namespace OpenMM

plugins/drude/platforms/reference/tests/TestReferenceDrudeLangevinIntegrator.cpp

@@ -60,6 +60,7 @@ void testSinglePair() {
     const double mass2 = 0.1;
     const double totalMass = mass1+mass2;
     const double reducedMass = (mass1*mass2)/(mass1+mass2);
+    const double maxDistance = 0.05;
     System system;
     system.addParticle(mass1);
     system.addParticle(mass2);
@@ -70,6 +71,7 @@ void testSinglePair() {
     positions[0] = Vec3(0, 0, 0);
     positions[1] = Vec3(0, 0, 0);
     DrudeLangevinIntegrator integ(temperature, 20.0, temperatureDrude, 20.0, 0.003);
+    integ.setMaxDrudeDistance(maxDistance);
     Platform& platform = Platform::getPlatformByName("Reference");
     Context context(system, integ, platform);
     context.setPositions(positions);
@@ -84,12 +86,15 @@ void testSinglePair() {
     int numSteps = 10000;
     for (int i = 0; i < numSteps; i++) {
         integ.step(10);
-        State state = context.getState(State::Velocities);
+        State state = context.getState(State::Velocities | State::Positions);
         const vector<Vec3>& vel = state.getVelocities();
         Vec3 velCM = vel[0]*(mass1/totalMass) + vel[1]*(mass2/totalMass);
         keCM += 0.5*totalMass*velCM.dot(velCM);
         Vec3 velInternal = vel[0]-vel[1];
         keInternal += 0.5*reducedMass*velInternal.dot(velInternal);
+        Vec3 delta = state.getPositions()[0]-state.getPositions()[1];
+        double distance = sqrt(delta.dot(delta));
+        ASSERT(distance <= maxDistance*(1+1e-6));
     }
     ASSERT_USUALLY_EQUAL_TOL(3*0.5*BOLTZ*temperature, keCM/numSteps, 0.1);
     ASSERT_USUALLY_EQUAL_TOL(3*0.5*BOLTZ*temperatureDrude, keInternal/numSteps, 0.01);