Constraint tolerance is now adjustable

openmmapi/include/Integrator.h

@@ -66,6 +66,18 @@ public:
     void setStepSize(double size) {
         stepSize = size;
     }
+    /**
+     * Get the distance tolerance within which constraints are maintained, measured in nm.
+     */
+    double getConstraintTolerance() const {
+        return constraintTol;
+    }
+    /**
+     * Set the distance tolerance within which constraints are maintained, measured in nm.
+     */
+    void setConstraintTolerance(double tol) {
+        constraintTol = tol;
+    }
     /**
      * Advance a simulation through time by taking a series of time steps.
      * 
@@ -85,7 +97,7 @@ protected:
      */
     virtual std::vector<std::string> getKernelNames() = 0;
 private:
-    double stepSize;
+    double stepSize, constraintTol;
 };
 
 } // namespace OpenMM

openmmapi/include/System.h

@@ -46,7 +46,7 @@ class OPENMM_EXPORT Force;
  * <ol>
  * <li>The set of particles in the system</li>
  * <li>The forces acting on them</li>
- * <li>Pairs of particles whose separation should be connstrained to a fixed value</li>
+ * <li>Pairs of particles whose separation should be constrained to a fixed value</li>
  * </ol>
  * 
  * The particles and constraints are defined directly by the System object.

openmmapi/src/BrownianIntegrator.cpp

@@ -43,26 +43,13 @@ BrownianIntegrator::BrownianIntegrator(double temperature, double frictionCoeff,
     setTemperature(temperature);
     setFriction(frictionCoeff);
     setStepSize(stepSize);
+    setConstraintTolerance(1e-4);
 }
 
 void BrownianIntegrator::initialize(OpenMMContextImpl& contextRef) {
     context = &contextRef;
     kernel = context->getPlatform().createKernel(IntegrateBrownianStepKernel::Name(), contextRef);
-    const System& system = context->getSystem();
-    vector<double> masses(system.getNumParticles());
-    vector<std::vector<int> > constraintIndices(system.getNumConstraints());
-    vector<double> constraintLengths(system.getNumConstraints());
-    for (int i = 0; i < system.getNumParticles(); ++i)
-        masses[i] = system.getParticleMass(i);
-    for (int i = 0; i < system.getNumConstraints(); ++i) {
-        int particle1, particle2;
-        double distance;
-        system.getConstraintParameters(i, particle1, particle2, distance);
-        constraintIndices[i].push_back(particle1);
-        constraintIndices[i].push_back(particle2);
-        constraintLengths[i] = distance;
-    }
-    dynamic_cast<IntegrateBrownianStepKernel&>(kernel.getImpl()).initialize(system, *this);
+    dynamic_cast<IntegrateBrownianStepKernel&>(kernel.getImpl()).initialize(contextRef.getSystem(), *this);
 }
 
 vector<string> BrownianIntegrator::getKernelNames() {

openmmapi/src/LangevinIntegrator.cpp

@@ -43,6 +43,7 @@ LangevinIntegrator::LangevinIntegrator(double temperature, double frictionCoeff,
     setTemperature(temperature);
     setFriction(frictionCoeff);
     setStepSize(stepSize);
+    setConstraintTolerance(1e-4);
 }
 
 void LangevinIntegrator::initialize(OpenMMContextImpl& contextRef) {

openmmapi/src/VerletIntegrator.cpp

@@ -41,6 +41,7 @@ using std::vector;
 
 VerletIntegrator::VerletIntegrator(double stepSize) {
     setStepSize(stepSize);
+    setConstraintTolerance(1e-4);
 }
 
 void VerletIntegrator::initialize(OpenMMContextImpl& contextRef) {

platforms/cuda/src/CudaKernels.cpp

@@ -311,7 +311,7 @@ void CudaCalcGBSAOBCForceKernel::initialize(const System& system, const GBSAOBCF
 void CudaCalcGBSAOBCForceKernel::executeForces(OpenMMContextImpl& context) {
 }
 
-static void initializeIntegration(const System& system, CudaPlatform::PlatformData& data) {
+static void initializeIntegration(const System& system, CudaPlatform::PlatformData& data, const Integrator& integrator) {
     
     // Set masses.
     
@@ -340,7 +340,7 @@ static void initializeIntegration(const System& system, CudaPlatform::PlatformDa
         invMass1[i] = 1.0f/mass[particle1Index];
         invMass2[i] = 1.0f/mass[particle2Index];
     }
-    gpuSetShakeParameters(gpu, particle1, particle2, distance, invMass1, invMass2);
+    gpuSetShakeParameters(gpu, particle1, particle2, distance, invMass1, invMass2, integrator.getConstraintTolerance());
     
     // Initialize any terms that haven't already been handled by a Force.
     
@@ -377,7 +377,7 @@ CudaIntegrateVerletStepKernel::~CudaIntegrateVerletStepKernel() {
 }
 
 void CudaIntegrateVerletStepKernel::initialize(const System& system, const VerletIntegrator& integrator) {
-    initializeIntegration(system, data);
+    initializeIntegration(system, data, integrator);
     prevStepSize = -1.0;
 }
 
@@ -406,7 +406,7 @@ CudaIntegrateLangevinStepKernel::~CudaIntegrateLangevinStepKernel() {
 }
 
 void CudaIntegrateLangevinStepKernel::initialize(const System& system, const LangevinIntegrator& integrator) {
-    initializeIntegration(system, data);
+    initializeIntegration(system, data, integrator);
     prevStepSize = -1.0;
 }
 
@@ -441,7 +441,7 @@ CudaIntegrateBrownianStepKernel::~CudaIntegrateBrownianStepKernel() {
 }
 
 void CudaIntegrateBrownianStepKernel::initialize(const System& system, const BrownianIntegrator& integrator) {
-    initializeIntegration(system, data);
+    initializeIntegration(system, data, integrator);
     prevStepSize = -1.0;
 }
 

platforms/reference/src/ReferenceKernels.cpp

@@ -477,6 +477,7 @@ void ReferenceIntegrateVerletStepKernel::execute(OpenMMContextImpl& context, con
         dynamics->setReferenceShakeAlgorithm(shake);
         prevStepSize = stepSize;
     }
+    shake->setTolerance(integrator.getConstraintTolerance());
     dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses);
 }
 
@@ -537,6 +538,7 @@ void ReferenceIntegrateLangevinStepKernel::execute(OpenMMContextImpl& context, c
         prevFriction = friction;
         prevStepSize = stepSize;
     }
+    shake->setTolerance(integrator.getConstraintTolerance());
     dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses);
 }
 
@@ -596,6 +598,7 @@ void ReferenceIntegrateBrownianStepKernel::execute(OpenMMContextImpl& context, c
         prevFriction = friction;
         prevStepSize = stepSize;
     }
+    shake->setTolerance(integrator.getConstraintTolerance());
     dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses);
 }