Merge pull request #230 from peastman/master

Major simplification to constraint handling in reference platform

platforms/reference/src/SimTKReference/ReferenceSETTLEAlgorithm.cpp

@@ -35,19 +35,11 @@ using namespace OpenMM;
 using namespace std;
 
 ReferenceSETTLEAlgorithm::ReferenceSETTLEAlgorithm(const vector<int>& atom1, const vector<int>& atom2, const vector<int>& atom3,
-        const vector<RealOpenMM>& distance1, const vector<RealOpenMM>& distance2, vector<RealOpenMM>& masses, RealOpenMM tolerance) :
-        atom1(atom1), atom2(atom2), atom3(atom3), distance1(distance1), distance2(distance2), masses(masses), tolerance(tolerance) {
+        const vector<RealOpenMM>& distance1, const vector<RealOpenMM>& distance2, vector<RealOpenMM>& masses) :
+        atom1(atom1), atom2(atom2), atom3(atom3), distance1(distance1), distance2(distance2), masses(masses) {
 }
 
-RealOpenMM ReferenceSETTLEAlgorithm::getTolerance() const {
-    return tolerance;
-}
-
-void ReferenceSETTLEAlgorithm::setTolerance(RealOpenMM tolerance) {
-    this->tolerance = tolerance;
-}
-
-void ReferenceSETTLEAlgorithm::apply(vector<OpenMM::RealVec>& atomCoordinates, vector<OpenMM::RealVec>& atomCoordinatesP, vector<RealOpenMM>& inverseMasses) {
+void ReferenceSETTLEAlgorithm::apply(vector<OpenMM::RealVec>& atomCoordinates, vector<OpenMM::RealVec>& atomCoordinatesP, vector<RealOpenMM>& inverseMasses, RealOpenMM tolerance) {
     for (int index = 0; index < (int) atom1.size(); ++index) {
         RealVec apos0 = atomCoordinates[atom1[index]];
         RealVec xp0 = atomCoordinatesP[atom1[index]]-apos0;
@@ -190,7 +182,7 @@ void ReferenceSETTLEAlgorithm::apply(vector<OpenMM::RealVec>& atomCoordinates, v
     }
 }
 
-void ReferenceSETTLEAlgorithm::applyToVelocities(vector<OpenMM::RealVec>& atomCoordinates, vector<OpenMM::RealVec>& velocities, vector<RealOpenMM>& inverseMasses) {
+void ReferenceSETTLEAlgorithm::applyToVelocities(vector<OpenMM::RealVec>& atomCoordinates, vector<OpenMM::RealVec>& velocities, vector<RealOpenMM>& inverseMasses, RealOpenMM tolerance) {
     for (int index = 0; index < (int) atom1.size(); ++index) {
         RealVec apos0 = atomCoordinates[atom1[index]];
         RealVec apos1 = atomCoordinates[atom2[index]];

platforms/reference/src/SimTKReference/ReferenceStochasticDynamics.cpp

@@ -200,7 +200,7 @@ void ReferenceStochasticDynamics::updatePart2( int numberOfAtoms, vector<RealVec
    --------------------------------------------------------------------------------------- */
 
 void ReferenceStochasticDynamics::update(const OpenMM::System& system, vector<RealVec>& atomCoordinates,
-                                          vector<RealVec>& velocities, vector<RealVec>& forces, vector<RealOpenMM>& masses ){
+                                          vector<RealVec>& velocities, vector<RealVec>& forces, vector<RealOpenMM>& masses, RealOpenMM tolerance) {
 
    // ---------------------------------------------------------------------------------------
 
@@ -235,7 +235,7 @@ void ReferenceStochasticDynamics::update(const OpenMM::System& system, vector<Re
 
    ReferenceConstraintAlgorithm* referenceConstraintAlgorithm = getReferenceConstraintAlgorithm();
    if (referenceConstraintAlgorithm)
-      referenceConstraintAlgorithm->apply(atomCoordinates, xPrime, inverseMasses);
+      referenceConstraintAlgorithm->apply(atomCoordinates, xPrime, inverseMasses, tolerance);
 
    // copy xPrime -> atomCoordinates
 

platforms/reference/src/SimTKReference/ReferenceVariableStochasticDynamics.cpp

@@ -260,7 +260,7 @@ void ReferenceVariableStochasticDynamics::updatePart2( int numberOfAtoms, vector
 
 void ReferenceVariableStochasticDynamics::update(const OpenMM::System& system, vector<RealVec>& atomCoordinates,
                                           vector<RealVec>& velocities,
-                                          vector<RealVec>& forces, vector<RealOpenMM>& masses, RealOpenMM maxStepSize ){
+                                          vector<RealVec>& forces, vector<RealOpenMM>& masses, RealOpenMM maxStepSize, RealOpenMM tolerance) {
 
    // ---------------------------------------------------------------------------------------
 
@@ -279,7 +279,7 @@ void ReferenceVariableStochasticDynamics::update(const OpenMM::System& system, v
 
    ReferenceConstraintAlgorithm* referenceConstraintAlgorithm = getReferenceConstraintAlgorithm();
    if (referenceConstraintAlgorithm)
-      referenceConstraintAlgorithm->apply(atomCoordinates, xPrime, inverseMasses);
+      referenceConstraintAlgorithm->apply(atomCoordinates, xPrime, inverseMasses, tolerance);
 
    // copy xPrime -> atomCoordinates
 

platforms/reference/src/SimTKReference/ReferenceVariableVerletDynamics.cpp

@@ -104,7 +104,7 @@ void ReferenceVariableVerletDynamics::setAccuracy( RealOpenMM accuracy ) {
 
 void ReferenceVariableVerletDynamics::update(const OpenMM::System& system, vector<RealVec>& atomCoordinates,
                                           vector<RealVec>& velocities,
-                                          vector<RealVec>& forces, vector<RealOpenMM>& masses, RealOpenMM maxStepSize ){
+                                          vector<RealVec>& forces, vector<RealOpenMM>& masses, RealOpenMM maxStepSize, RealOpenMM tolerance) {
 
     // ---------------------------------------------------------------------------------------
 
@@ -155,7 +155,7 @@ void ReferenceVariableVerletDynamics::update(const OpenMM::System& system, vecto
     }
     ReferenceConstraintAlgorithm* referenceConstraintAlgorithm = getReferenceConstraintAlgorithm();
     if (referenceConstraintAlgorithm)
-        referenceConstraintAlgorithm->apply(atomCoordinates, xPrime, inverseMasses);
+        referenceConstraintAlgorithm->apply(atomCoordinates, xPrime, inverseMasses, tolerance);
 
    // Update the positions and velocities.
 

platforms/reference/src/SimTKReference/ReferenceVerletDynamics.cpp

@@ -95,7 +95,7 @@ ReferenceVerletDynamics::~ReferenceVerletDynamics( ){
 
 void ReferenceVerletDynamics::update(const OpenMM::System& system, vector<RealVec>& atomCoordinates,
                                           vector<RealVec>& velocities,
-                                          vector<RealVec>& forces, vector<RealOpenMM>& masses ){
+                                          vector<RealVec>& forces, vector<RealOpenMM>& masses, RealOpenMM tolerance) {
 
    // ---------------------------------------------------------------------------------------
 
@@ -131,7 +131,7 @@ void ReferenceVerletDynamics::update(const OpenMM::System& system, vector<RealVe
    }
    ReferenceConstraintAlgorithm* referenceConstraintAlgorithm = getReferenceConstraintAlgorithm();
    if (referenceConstraintAlgorithm)
-      referenceConstraintAlgorithm->apply(atomCoordinates, xPrime, inverseMasses);
+      referenceConstraintAlgorithm->apply(atomCoordinates, xPrime, inverseMasses, tolerance);
    
    // Update the positions and velocities.
    

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

@@ -56,7 +56,12 @@ static vector<RealVec>& extractForces(ContextImpl& context) {
     return *((vector<RealVec>*) data->forces);
 }
 
-static double computeShiftedKineticEnergy(ContextImpl& context, vector<double>& inverseMasses, double timeShift, ReferenceConstraintAlgorithm* constraints) {
+static ReferenceConstraints& extractConstraints(ContextImpl& context) {
+    ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
+    return *(ReferenceConstraints*) data->constraints;
+}
+
+static double computeShiftedKineticEnergy(ContextImpl& context, vector<double>& inverseMasses, double timeShift) {
     const System& system = context.getSystem();
     int numParticles = system.getNumParticles();
     vector<RealVec>& posData = extractPositions(context);
@@ -75,10 +80,7 @@ static double computeShiftedKineticEnergy(ContextImpl& context, vector<double>&
     
     // Apply constraints to them.
     
-    if (constraints != NULL) {
-        constraints->setTolerance(1e-4);
-        constraints->applyToVelocities(posData, shiftedVel, inverseMasses);
-    }
+    extractConstraints(context).applyToVelocities(posData, shiftedVel, inverseMasses, 1e-4);
     
     // Compute the kinetic energy.
     
@@ -224,8 +226,6 @@ void ReferenceCalcDrudeForceKernel::copyParametersToContext(ContextImpl& context
 }
 
 ReferenceIntegrateDrudeLangevinStepKernel::~ReferenceIntegrateDrudeLangevinStepKernel() {
-    if (constraints != NULL)
-        delete constraints;
 }
 
 void ReferenceIntegrateDrudeLangevinStepKernel::initialize(const System& system, const DrudeLangevinIntegrator& integrator, const DrudeForce& force) {
@@ -254,23 +254,6 @@ void ReferenceIntegrateDrudeLangevinStepKernel::initialize(const System& system,
         pairInvReducedMass.push_back((m1+m2)/(m1*m2));
     }
     normalParticles.insert(normalParticles.begin(), particles.begin(), particles.end());
-    
-    // Prepare constraints.
-    
-    if (system.getNumConstraints() > 0) {
-        vector<pair<int, int> > constraintIndices;
-        vector<RealOpenMM> constraintDistances;
-        for (int i = 0; i < system.getNumConstraints(); ++i) {
-            int particle1, particle2;
-            double distance;
-            system.getConstraintParameters(i, particle1, particle2, distance);
-            if (system.getParticleMass(particle1) != 0 || system.getParticleMass(particle2) != 0) {
-                constraintIndices.push_back(make_pair(particle1, particle2));
-                constraintDistances.push_back(distance);
-            }
-        }
-        constraints = new ReferenceConstraints(system, (RealOpenMM) integrator.getConstraintTolerance());
-    }
 }
 
 void ReferenceIntegrateDrudeLangevinStepKernel::execute(ContextImpl& context, const DrudeLangevinIntegrator& integrator) {
@@ -330,8 +313,7 @@ void ReferenceIntegrateDrudeLangevinStepKernel::execute(ContextImpl& context, co
     
     // Apply constraints.
     
-    if (constraints != NULL)
-        constraints->apply(pos, xPrime, particleInvMass);
+    extractConstraints(context).apply(pos, xPrime, particleInvMass, integrator.getConstraintTolerance());
     
     // Record the constrained positions and velocities.
     
@@ -348,12 +330,10 @@ void ReferenceIntegrateDrudeLangevinStepKernel::execute(ContextImpl& context, co
 }
 
 double ReferenceIntegrateDrudeLangevinStepKernel::computeKineticEnergy(ContextImpl& context, const DrudeLangevinIntegrator& integrator) {
-    return computeShiftedKineticEnergy(context, particleInvMass, 0.5*integrator.getStepSize(), constraints);
+    return computeShiftedKineticEnergy(context, particleInvMass, 0.5*integrator.getStepSize());
 }
 
 ReferenceIntegrateDrudeSCFStepKernel::~ReferenceIntegrateDrudeSCFStepKernel() {
-    if (constraints != NULL)
-        delete constraints;
     if (minimizerPos != NULL)
         lbfgs_free(minimizerPos);
 }
@@ -377,23 +357,6 @@ void ReferenceIntegrateDrudeSCFStepKernel::initialize(const System& system, cons
         particleInvMass.push_back(mass == 0.0 ? 0.0 : 1.0/mass);
     }
     
-    // Prepare constraints.
-    
-    if (system.getNumConstraints() > 0) {
-        vector<pair<int, int> > constraintIndices;
-        vector<RealOpenMM> constraintDistances;
-        for (int i = 0; i < system.getNumConstraints(); ++i) {
-            int particle1, particle2;
-            double distance;
-            system.getConstraintParameters(i, particle1, particle2, distance);
-            if (system.getParticleMass(particle1) != 0 || system.getParticleMass(particle2) != 0) {
-                constraintIndices.push_back(make_pair(particle1, particle2));
-                constraintDistances.push_back(distance);
-            }
-        }
-        constraints = new ReferenceConstraints(system, (RealOpenMM) integrator.getConstraintTolerance());
-    }
-    
     // Initialize the energy minimizer.
     
     minimizerPos = lbfgs_malloc(drudeParticles.size()*3);
@@ -424,8 +387,7 @@ void ReferenceIntegrateDrudeSCFStepKernel::execute(ContextImpl& context, const D
         
     // Apply constraints.
     
-    if (constraints != NULL)
-        constraints->apply(pos, xPrime, particleInvMass);
+    extractConstraints(context).apply(pos, xPrime, particleInvMass, integrator.getConstraintTolerance());
     
     // Record the constrained positions and velocities.
     
@@ -446,7 +408,7 @@ void ReferenceIntegrateDrudeSCFStepKernel::execute(ContextImpl& context, const D
 }
 
 double ReferenceIntegrateDrudeSCFStepKernel::computeKineticEnergy(ContextImpl& context, const DrudeSCFIntegrator& integrator) {
-    return computeShiftedKineticEnergy(context, particleInvMass, 0.5*integrator.getStepSize(), constraints);
+    return computeShiftedKineticEnergy(context, particleInvMass, 0.5*integrator.getStepSize());
 }
 
 struct MinimizerData {

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

@@ -86,7 +86,7 @@ private:
 class ReferenceIntegrateDrudeLangevinStepKernel : public IntegrateDrudeLangevinStepKernel {
 public:
     ReferenceIntegrateDrudeLangevinStepKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) :
-        IntegrateDrudeLangevinStepKernel(name, platform), data(data), constraints(NULL) {
+        IntegrateDrudeLangevinStepKernel(name, platform), data(data) {
     }
     ~ReferenceIntegrateDrudeLangevinStepKernel();
     /**
@@ -118,7 +118,6 @@ private:
     std::vector<double> particleInvMass;
     std::vector<double> pairInvTotalMass;
     std::vector<double> pairInvReducedMass;
-    ReferenceConstraintAlgorithm* constraints;
 };
 
 /**
@@ -127,7 +126,7 @@ private:
 class ReferenceIntegrateDrudeSCFStepKernel : public IntegrateDrudeSCFStepKernel {
 public:
     ReferenceIntegrateDrudeSCFStepKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) :
-        IntegrateDrudeSCFStepKernel(name, platform), data(data), constraints(NULL), minimizerPos(NULL) {
+        IntegrateDrudeSCFStepKernel(name, platform), data(data), minimizerPos(NULL) {
     }
     ~ReferenceIntegrateDrudeSCFStepKernel();
     /**
@@ -157,7 +156,6 @@ private:
     ReferencePlatform::PlatformData& data;
     std::vector<int> drudeParticles;
     std::vector<double> particleInvMass;
-    ReferenceConstraintAlgorithm* constraints;
     lbfgsfloatval_t *minimizerPos;
     lbfgs_parameter_t minimizerParams;
 };