Ignore constraints that connect two massless particles (feature request 1915)

openmmapi/src/ContextImpl.cpp

@@ -69,7 +69,9 @@ ContextImpl::ContextImpl(Context& owner, const System& system, Integrator& integ
         int particle1, particle2;
         double distance;
         system.getConstraintParameters(i, particle1, particle2, distance);
-        if (system.getParticleMass(particle1) == 0.0 || system.getParticleMass(particle2) == 0.0)
+        double mass1 = system.getParticleMass(particle1);
+        double mass2 = system.getParticleMass(particle2);
+        if ((mass1 == 0.0 && mass2 != 0.0) || (mass2 == 0.0 && mass1 != 0.0))
             throw OpenMMException("A constraint cannot involve a massless particle");
     }
     

platforms/cuda/src/CudaIntegrationUtilities.cpp

@@ -123,15 +123,21 @@ CudaIntegrationUtilities::CudaIntegrationUtilities(CudaContext& context, const S
 
     // Record the set of constraints and how many constraints each atom is involved in.
 
-    int numConstraints = system.getNumConstraints();
-    vector<int> atom1(numConstraints);
-    vector<int> atom2(numConstraints);
-    vector<double> distance(numConstraints);
+    vector<int> atom1;
+    vector<int> atom2;
+    vector<double> distance;
     vector<int> constraintCount(context.getNumAtoms(), 0);
-    for (int i = 0; i < numConstraints; i++) {
-        system.getConstraintParameters(i, atom1[i], atom2[i], distance[i]);
-        constraintCount[atom1[i]]++;
-        constraintCount[atom2[i]]++;
+    for (int i = 0; i < system.getNumConstraints(); i++) {
+        int p1, p2;
+        double d;
+        system.getConstraintParameters(i, p1, p2, d);
+        if (system.getParticleMass(p1) != 0 || system.getParticleMass(p2) != 0) {
+            atom1.push_back(p1);
+            atom2.push_back(p2);
+            distance.push_back(d);
+            constraintCount[p1]++;
+            constraintCount[p2]++;
+        }
     }
 
     // Identify clusters of three atoms that can be treated with SETTLE.  First, for every

platforms/cuda/tests/TestCudaBrownianIntegrator.cpp

@@ -172,6 +172,37 @@ void testConstraints() {
     }
 }
 
+void testConstrainedMasslessParticles() {
+    System system;
+    system.addParticle(0.0);
+    system.addParticle(1.0);
+    system.addConstraint(0, 1, 1.5);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(-1, 0, 0);
+    positions[1] = Vec3(1, 0, 0);
+    BrownianIntegrator integrator(300.0, 2.0, 0.01);
+    bool failed = false;
+    try {
+        // This should throw an exception.
+        
+        Context context(system, integrator, platform);
+    }
+    catch (exception& ex) {
+        failed = true;
+    }
+    ASSERT(failed);
+    
+    // Now make both particles massless, which should work.
+    
+    system.setParticleMass(1, 0.0);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    context.setVelocitiesToTemperature(300.0);
+    integrator.step(1);
+    State state = context.getState(State::Velocities);
+    ASSERT_EQUAL(0.0, state.getVelocities()[0][0]);
+}
+
 void testRandomSeed() {
     const int numParticles = 8;
     const double temp = 100.0;
@@ -237,6 +268,7 @@ int main(int argc, char* argv[]) {
         testSingleBond();
         testTemperature();
         testConstraints();
+        testConstrainedMasslessParticles();
         testRandomSeed();
     }
     catch(const exception& e) {

platforms/cuda/tests/TestCudaCustomIntegrator.cpp

@@ -233,6 +233,43 @@ void testVelocityConstraints() {
     }
 }
 
+void testConstrainedMasslessParticles() {
+    System system;
+    system.addParticle(0.0);
+    system.addParticle(1.0);
+    system.addConstraint(0, 1, 1.5);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(-1, 0, 0);
+    positions[1] = Vec3(1, 0, 0);
+    CustomIntegrator integrator(0.002);
+    integrator.addPerDofVariable("oldx", 0);
+    integrator.addComputePerDof("v", "v+dt*f/m");
+    integrator.addComputePerDof("oldx", "x");
+    integrator.addComputePerDof("x", "x+dt*v");
+    integrator.addConstrainPositions();
+    integrator.addComputePerDof("v", "(x-oldx)/dt");
+    bool failed = false;
+    try {
+        // This should throw an exception.
+        
+        Context context(system, integrator, platform);
+    }
+    catch (exception& ex) {
+        failed = true;
+    }
+    ASSERT(failed);
+    
+    // Now make both particles massless, which should work.
+    
+    system.setParticleMass(1, 0.0);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    context.setVelocitiesToTemperature(300.0);
+    integrator.step(1);
+    State state = context.getState(State::Velocities | State::Positions);
+    ASSERT_EQUAL(0.0, state.getVelocities()[0][0]);
+}
+
 /**
  * Test an integrator with an AndersenThermostat to see if updateContextState()
  * is being handled correctly.
@@ -724,6 +761,7 @@ int main(int argc, char* argv[]) {
         testSingleBond();
         testConstraints();
         testVelocityConstraints();
+        testConstrainedMasslessParticles();
         testWithThermostat();
         testMonteCarlo();
         testSum();

platforms/cuda/tests/TestCudaLangevinIntegrator.cpp

@@ -177,6 +177,37 @@ void testConstraints() {
     }
 }
 
+void testConstrainedMasslessParticles() {
+    System system;
+    system.addParticle(0.0);
+    system.addParticle(1.0);
+    system.addConstraint(0, 1, 1.5);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(-1, 0, 0);
+    positions[1] = Vec3(1, 0, 0);
+    LangevinIntegrator integrator(300.0, 2.0, 0.01);
+    bool failed = false;
+    try {
+        // This should throw an exception.
+        
+        Context context(system, integrator, platform);
+    }
+    catch (exception& ex) {
+        failed = true;
+    }
+    ASSERT(failed);
+    
+    // Now make both particles massless, which should work.
+    
+    system.setParticleMass(1, 0.0);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    context.setVelocitiesToTemperature(300.0);
+    integrator.step(1);
+    State state = context.getState(State::Velocities);
+    ASSERT_EQUAL(0.0, state.getVelocities()[0][0]);
+}
+
 void testRandomSeed() {
     const int numParticles = 8;
     const double temp = 100.0;
@@ -241,6 +272,7 @@ int main(int argc, char* argv[]) {
         testSingleBond();
         testTemperature();
         testConstraints();
+        testConstrainedMasslessParticles();
         testRandomSeed();
     }
     catch(const exception& e) {

platforms/cuda/tests/TestCudaVariableLangevinIntegrator.cpp

@@ -172,6 +172,37 @@ void testConstraints() {
     }
 }
 
+void testConstrainedMasslessParticles() {
+    System system;
+    system.addParticle(0.0);
+    system.addParticle(1.0);
+    system.addConstraint(0, 1, 1.5);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(-1, 0, 0);
+    positions[1] = Vec3(1, 0, 0);
+    VariableLangevinIntegrator integrator(300.0, 2.0, 0.01);
+    bool failed = false;
+    try {
+        // This should throw an exception.
+        
+        Context context(system, integrator, platform);
+    }
+    catch (exception& ex) {
+        failed = true;
+    }
+    ASSERT(failed);
+    
+    // Now make both particles massless, which should work.
+    
+    system.setParticleMass(1, 0.0);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    context.setVelocitiesToTemperature(300.0);
+    integrator.step(1);
+    State state = context.getState(State::Velocities);
+    ASSERT_EQUAL(0.0, state.getVelocities()[0][0]);
+}
+
 void testRandomSeed() {
     const int numParticles = 8;
     const double temp = 100.0;
@@ -295,6 +326,7 @@ int main(int argc, char* argv[]) {
         testSingleBond();
         testTemperature();
         testConstraints();
+        testConstrainedMasslessParticles();
         testRandomSeed();
         testArgonBox();
     }

platforms/cuda/tests/TestCudaVariableVerletIntegrator.cpp

@@ -211,6 +211,37 @@ void testConstrainedClusters() {
     ASSERT(context.getState(State::Positions).getTime() > 0.1);
 }
 
+void testConstrainedMasslessParticles() {
+    System system;
+    system.addParticle(0.0);
+    system.addParticle(1.0);
+    system.addConstraint(0, 1, 1.5);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(-1, 0, 0);
+    positions[1] = Vec3(1, 0, 0);
+    VariableVerletIntegrator integrator(0.01);
+    bool failed = false;
+    try {
+        // This should throw an exception.
+        
+        Context context(system, integrator, platform);
+    }
+    catch (exception& ex) {
+        failed = true;
+    }
+    ASSERT(failed);
+    
+    // Now make both particles massless, which should work.
+    
+    system.setParticleMass(1, 0.0);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    context.setVelocitiesToTemperature(300.0);
+    integrator.step(1);
+    State state = context.getState(State::Velocities);
+    ASSERT_EQUAL(0.0, state.getVelocities()[0][0]);
+}
+
 void testArgonBox() {
     const int gridSize = 8;
     const double mass = 40.0;            // Ar atomic mass
@@ -274,6 +305,7 @@ int main(int argc, char* argv[]) {
         testSingleBond();
         testConstraints();
         testConstrainedClusters();
+        testConstrainedMasslessParticles();
         testArgonBox();
     }
     catch(const exception& e) {

platforms/cuda/tests/TestCudaVerletIntegrator.cpp

@@ -203,6 +203,37 @@ void testConstrainedClusters() {
     }
 }
 
+void testConstrainedMasslessParticles() {
+    System system;
+    system.addParticle(0.0);
+    system.addParticle(1.0);
+    system.addConstraint(0, 1, 1.5);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(-1, 0, 0);
+    positions[1] = Vec3(1, 0, 0);
+    VerletIntegrator integrator(0.01);
+    bool failed = false;
+    try {
+        // This should throw an exception.
+        
+        Context context(system, integrator, platform);
+    }
+    catch (exception& ex) {
+        failed = true;
+    }
+    ASSERT(failed);
+    
+    // Now make both particles massless, which should work.
+    
+    system.setParticleMass(1, 0.0);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    context.setVelocitiesToTemperature(300.0);
+    integrator.step(1);
+    State state = context.getState(State::Velocities);
+    ASSERT_EQUAL(0.0, state.getVelocities()[0][0]);
+}
+
 int main(int argc, char* argv[]) {
     try {
         if (argc > 1)
@@ -210,6 +241,7 @@ int main(int argc, char* argv[]) {
         testSingleBond();
         testConstraints();
         testConstrainedClusters();
+        testConstrainedMasslessParticles();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;

platforms/opencl/src/OpenCLIntegrationUtilities.cpp

@@ -141,15 +141,21 @@ OpenCLIntegrationUtilities::OpenCLIntegrationUtilities(OpenCLContext& context, c
 
     // Record the set of constraints and how many constraints each atom is involved in.
 
-    int numConstraints = system.getNumConstraints();
-    vector<int> atom1(numConstraints);
-    vector<int> atom2(numConstraints);
-    vector<double> distance(numConstraints);
+    vector<int> atom1;
+    vector<int> atom2;
+    vector<double> distance;
     vector<int> constraintCount(context.getNumAtoms(), 0);
-    for (int i = 0; i < numConstraints; i++) {
-        system.getConstraintParameters(i, atom1[i], atom2[i], distance[i]);
-        constraintCount[atom1[i]]++;
-        constraintCount[atom2[i]]++;
+    for (int i = 0; i < system.getNumConstraints(); i++) {
+        int p1, p2;
+        double d;
+        system.getConstraintParameters(i, p1, p2, d);
+        if (system.getParticleMass(p1) != 0 || system.getParticleMass(p2) != 0) {
+            atom1.push_back(p1);
+            atom2.push_back(p2);
+            distance.push_back(d);
+            constraintCount[p1]++;
+            constraintCount[p2]++;
+        }
     }
 
     // Identify clusters of three atoms that can be treated with SETTLE.  First, for every

platforms/opencl/tests/TestOpenCLBrownianIntegrator.cpp

@@ -172,6 +172,37 @@ void testConstraints() {
     }
 }
 
+void testConstrainedMasslessParticles() {
+    System system;
+    system.addParticle(0.0);
+    system.addParticle(1.0);
+    system.addConstraint(0, 1, 1.5);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(-1, 0, 0);
+    positions[1] = Vec3(1, 0, 0);
+    BrownianIntegrator integrator(300.0, 2.0, 0.01);
+    bool failed = false;
+    try {
+        // This should throw an exception.
+        
+        Context context(system, integrator, platform);
+    }
+    catch (exception& ex) {
+        failed = true;
+    }
+    ASSERT(failed);
+    
+    // Now make both particles massless, which should work.
+    
+    system.setParticleMass(1, 0.0);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    context.setVelocitiesToTemperature(300.0);
+    integrator.step(1);
+    State state = context.getState(State::Velocities);
+    ASSERT_EQUAL(0.0, state.getVelocities()[0][0]);
+}
+
 void testRandomSeed() {
     const int numParticles = 8;
     const double temp = 100.0;
@@ -237,6 +268,7 @@ int main(int argc, char* argv[]) {
         testSingleBond();
         testTemperature();
         testConstraints();
+        testConstrainedMasslessParticles();
         testRandomSeed();
     }
     catch(const exception& e) {

platforms/opencl/tests/TestOpenCLCustomIntegrator.cpp

@@ -233,6 +233,43 @@ void testVelocityConstraints() {
     }
 }
 
+void testConstrainedMasslessParticles() {
+    System system;
+    system.addParticle(0.0);
+    system.addParticle(1.0);
+    system.addConstraint(0, 1, 1.5);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(-1, 0, 0);
+    positions[1] = Vec3(1, 0, 0);
+    CustomIntegrator integrator(0.002);
+    integrator.addPerDofVariable("oldx", 0);
+    integrator.addComputePerDof("v", "v+dt*f/m");
+    integrator.addComputePerDof("oldx", "x");
+    integrator.addComputePerDof("x", "x+dt*v");
+    integrator.addConstrainPositions();
+    integrator.addComputePerDof("v", "(x-oldx)/dt");
+    bool failed = false;
+    try {
+        // This should throw an exception.
+        
+        Context context(system, integrator, platform);
+    }
+    catch (exception& ex) {
+        failed = true;
+    }
+    ASSERT(failed);
+    
+    // Now make both particles massless, which should work.
+    
+    system.setParticleMass(1, 0.0);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    context.setVelocitiesToTemperature(300.0);
+    integrator.step(1);
+    State state = context.getState(State::Velocities | State::Positions);
+    ASSERT_EQUAL(0.0, state.getVelocities()[0][0]);
+}
+
 /**
  * Test an integrator with an AndersenThermostat to see if updateContextState()
  * is being handled correctly.
@@ -724,6 +761,7 @@ int main(int argc, char* argv[]) {
         testSingleBond();
         testConstraints();
         testVelocityConstraints();
+        testConstrainedMasslessParticles();
         testWithThermostat();
         testMonteCarlo();
         testSum();

platforms/opencl/tests/TestOpenCLLangevinIntegrator.cpp

@@ -177,6 +177,37 @@ void testConstraints() {
     }
 }
 
+void testConstrainedMasslessParticles() {
+    System system;
+    system.addParticle(0.0);
+    system.addParticle(1.0);
+    system.addConstraint(0, 1, 1.5);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(-1, 0, 0);
+    positions[1] = Vec3(1, 0, 0);
+    LangevinIntegrator integrator(300.0, 2.0, 0.01);
+    bool failed = false;
+    try {
+        // This should throw an exception.
+        
+        Context context(system, integrator, platform);
+    }
+    catch (exception& ex) {
+        failed = true;
+    }
+    ASSERT(failed);
+    
+    // Now make both particles massless, which should work.
+    
+    system.setParticleMass(1, 0.0);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    context.setVelocitiesToTemperature(300.0);
+    integrator.step(1);
+    State state = context.getState(State::Velocities);
+    ASSERT_EQUAL(0.0, state.getVelocities()[0][0]);
+}
+
 void testRandomSeed() {
     const int numParticles = 8;
     const double temp = 100.0;
@@ -241,6 +272,7 @@ int main(int argc, char* argv[]) {
         testSingleBond();
         testTemperature();
         testConstraints();
+        testConstrainedMasslessParticles();
         testRandomSeed();
     }
     catch(const exception& e) {

platforms/opencl/tests/TestOpenCLVariableLangevinIntegrator.cpp

@@ -172,6 +172,37 @@ void testConstraints() {
     }
 }
 
+void testConstrainedMasslessParticles() {
+    System system;
+    system.addParticle(0.0);
+    system.addParticle(1.0);
+    system.addConstraint(0, 1, 1.5);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(-1, 0, 0);
+    positions[1] = Vec3(1, 0, 0);
+    VariableLangevinIntegrator integrator(300.0, 2.0, 0.01);
+    bool failed = false;
+    try {
+        // This should throw an exception.
+        
+        Context context(system, integrator, platform);
+    }
+    catch (exception& ex) {
+        failed = true;
+    }
+    ASSERT(failed);
+    
+    // Now make both particles massless, which should work.
+    
+    system.setParticleMass(1, 0.0);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    context.setVelocitiesToTemperature(300.0);
+    integrator.step(1);
+    State state = context.getState(State::Velocities);
+    ASSERT_EQUAL(0.0, state.getVelocities()[0][0]);
+}
+
 void testRandomSeed() {
     const int numParticles = 8;
     const double temp = 100.0;
@@ -295,6 +326,7 @@ int main(int argc, char* argv[]) {
         testSingleBond();
         testTemperature();
         testConstraints();
+        testConstrainedMasslessParticles();
         testRandomSeed();
         testArgonBox();
     }

platforms/opencl/tests/TestOpenCLVariableVerletIntegrator.cpp

@@ -210,6 +210,38 @@ void testConstrainedClusters() {
     ASSERT(context.getState(State::Positions).getTime() > 0.1);
 }
 
+void testConstrainedMasslessParticles() {
+    System system;
+    system.addParticle(0.0);
+    system.addParticle(1.0);
+    system.addConstraint(0, 1, 1.5);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(-1, 0, 0);
+    positions[1] = Vec3(1, 0, 0);
+    VariableVerletIntegrator integrator(0.01);
+    bool failed = false;
+    try {
+        // This should throw an exception.
+        
+        Context context(system, integrator, platform);
+    }
+    catch (exception& ex) {
+        failed = true;
+    }
+    ASSERT(failed);
+    
+    // Now make both particles massless, which should work.
+    
+    system.setParticleMass(1, 0.0);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    context.setVelocitiesToTemperature(300.0);
+    integrator.step(1);
+    State state = context.getState(State::Velocities);
+    ASSERT_EQUAL(0.0, state.getVelocities()[0][0]);
+}
+
+
 void testArgonBox() {
     const int gridSize = 8;
     const double mass = 40.0;            // Ar atomic mass
@@ -273,6 +305,7 @@ int main(int argc, char* argv[]) {
         testSingleBond();
         testConstraints();
         testConstrainedClusters();
+        testConstrainedMasslessParticles();
         testArgonBox();
     }
     catch(const exception& e) {

platforms/opencl/tests/TestOpenCLVerletIntegrator.cpp

@@ -201,6 +201,37 @@ void testConstrainedClusters() {
     }
 }
 
+void testConstrainedMasslessParticles() {
+    System system;
+    system.addParticle(0.0);
+    system.addParticle(1.0);
+    system.addConstraint(0, 1, 1.5);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(-1, 0, 0);
+    positions[1] = Vec3(1, 0, 0);
+    VerletIntegrator integrator(0.01);
+    bool failed = false;
+    try {
+        // This should throw an exception.
+        
+        Context context(system, integrator, platform);
+    }
+    catch (exception& ex) {
+        failed = true;
+    }
+    ASSERT(failed);
+    
+    // Now make both particles massless, which should work.
+    
+    system.setParticleMass(1, 0.0);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    context.setVelocitiesToTemperature(300.0);
+    integrator.step(1);
+    State state = context.getState(State::Velocities);
+    ASSERT_EQUAL(0.0, state.getVelocities()[0][0]);
+}
+
 int main(int argc, char* argv[]) {
     try {
         if (argc > 1)
@@ -208,6 +239,7 @@ int main(int argc, char* argv[]) {
         testSingleBond();
         testConstraints();
         testConstrainedClusters();
+        testConstrainedMasslessParticles();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;

platforms/reference/src/ReferenceKernels.cpp

@@ -315,17 +315,18 @@ void ReferenceApplyConstraintsKernel::initialize(const System& system) {
         masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
         inverseMasses[i] = 1.0/masses[i];
     }
-    numConstraints = system.getNumConstraints();
-    constraintIndices.resize(numConstraints);
-    constraintDistances.resize(numConstraints);
-    for (int i = 0; i < numConstraints; ++i) {
+    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);
-        constraintIndices[i].first = particle1;
-        constraintIndices[i].second = particle2;
-        constraintDistances[i] = static_cast<RealOpenMM>(distance);
+        if (system.getParticleMass(particle1) != 0 || system.getParticleMass(particle2) != 0) {
+            constraintIndices.push_back(make_pair(particle1, particle2));
+            constraintDistances.push_back(distance);
         }
+    }
+    numConstraints = constraintIndices.size();
 }
 
 ReferenceApplyConstraintsKernel::~ReferenceApplyConstraintsKernel() {
@@ -1715,17 +1716,18 @@ void ReferenceIntegrateVerletStepKernel::initialize(const System& system, const
     masses.resize(numParticles);
     for (int i = 0; i < numParticles; ++i)
         masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
-    numConstraints = system.getNumConstraints();
-    vector<pair<int, int> > constraintIndices(numConstraints);
-    vector<RealOpenMM> constraintDistances(numConstraints);
-    for (int i = 0; i < numConstraints; ++i) {
+    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);
-        constraintIndices[i].first = particle1;
-        constraintIndices[i].second = particle2;
-        constraintDistances[i] = static_cast<RealOpenMM>(distance);
+        if (system.getParticleMass(particle1) != 0 || system.getParticleMass(particle2) != 0) {
+            constraintIndices.push_back(make_pair(particle1, particle2));
+            constraintDistances.push_back(distance);
+        }
     }
+    numConstraints = constraintIndices.size();
     vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
     findAnglesForCCMA(system, angles);
     constraints = new ReferenceCCMAAlgorithm(system.getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
@@ -1767,17 +1769,18 @@ void ReferenceIntegrateLangevinStepKernel::initialize(const System& system, cons
     masses.resize(numParticles);
     for (int i = 0; i < numParticles; ++i)
         masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
-    numConstraints = system.getNumConstraints();
-    vector<pair<int, int> > constraintIndices(numConstraints);
-    vector<RealOpenMM> constraintDistances(numConstraints);
-    for (int i = 0; i < numConstraints; ++i) {
+    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);
-        constraintIndices[i].first = particle1;
-        constraintIndices[i].second = particle2;
-        constraintDistances[i] = static_cast<RealOpenMM>(distance);
+        if (system.getParticleMass(particle1) != 0 || system.getParticleMass(particle2) != 0) {
+            constraintIndices.push_back(make_pair(particle1, particle2));
+            constraintDistances.push_back(distance);
         }
+    }
+    numConstraints = constraintIndices.size();
     SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
     vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
     findAnglesForCCMA(system, angles);
@@ -1829,17 +1832,18 @@ void ReferenceIntegrateBrownianStepKernel::initialize(const System& system, cons
     masses.resize(numParticles);
     for (int i = 0; i < numParticles; ++i)
         masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
-    numConstraints = system.getNumConstraints();
-    vector<pair<int, int> > constraintIndices(numConstraints);
-    vector<RealOpenMM> constraintDistances(numConstraints);
-    for (int i = 0; i < numConstraints; ++i) {
+    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);
-        constraintIndices[i].first = particle1;
-        constraintIndices[i].second = particle2;
-        constraintDistances[i] = static_cast<RealOpenMM>(distance);
+        if (system.getParticleMass(particle1) != 0 || system.getParticleMass(particle2) != 0) {
+            constraintIndices.push_back(make_pair(particle1, particle2));
+            constraintDistances.push_back(distance);
+        }
     }
+    numConstraints = constraintIndices.size();
     SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
     vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
     findAnglesForCCMA(system, angles);
@@ -1890,17 +1894,18 @@ void ReferenceIntegrateVariableLangevinStepKernel::initialize(const System& syst
     masses.resize(numParticles);
     for (int i = 0; i < numParticles; ++i)
         masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
-    numConstraints = system.getNumConstraints();
-    vector<pair<int, int> > constraintIndices(numConstraints);
-    vector<RealOpenMM> constraintDistances(numConstraints);
-    for (int i = 0; i < numConstraints; ++i) {
+    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);
-        constraintIndices[i].first = particle1;
-        constraintIndices[i].second = particle2;
-        constraintDistances[i] = static_cast<RealOpenMM>(distance);
+        if (system.getParticleMass(particle1) != 0 || system.getParticleMass(particle2) != 0) {
+            constraintIndices.push_back(make_pair(particle1, particle2));
+            constraintDistances.push_back(distance);
         }
+    }
+    numConstraints = constraintIndices.size();
     SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
     vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
     findAnglesForCCMA(system, angles);
@@ -1952,17 +1957,18 @@ void ReferenceIntegrateVariableVerletStepKernel::initialize(const System& system
     masses.resize(numParticles);
     for (int i = 0; i < numParticles; ++i)
         masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
-    numConstraints = system.getNumConstraints();
-    vector<pair<int, int> > constraintIndices(numConstraints);
-    vector<RealOpenMM> constraintDistances(numConstraints);
-    for (int i = 0; i < numConstraints; ++i) {
+    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);
-        constraintIndices[i].first = particle1;
-        constraintIndices[i].second = particle2;
-        constraintDistances[i] = static_cast<RealOpenMM>(distance);
+        if (system.getParticleMass(particle1) != 0 || system.getParticleMass(particle2) != 0) {
+            constraintIndices.push_back(make_pair(particle1, particle2));
+            constraintDistances.push_back(distance);
+        }
     }
+    numConstraints = constraintIndices.size();
     vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
     findAnglesForCCMA(system, angles);
     constraints = new ReferenceCCMAAlgorithm(system.getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
@@ -2008,17 +2014,18 @@ void ReferenceIntegrateCustomStepKernel::initialize(const System& system, const
     masses.resize(numParticles);
     for (int i = 0; i < numParticles; ++i)
         masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
-    numConstraints = system.getNumConstraints();
-    vector<pair<int, int> > constraintIndices(numConstraints);
-    vector<RealOpenMM> constraintDistances(numConstraints);
-    for (int i = 0; i < numConstraints; ++i) {
+    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);
-        constraintIndices[i].first = particle1;
-        constraintIndices[i].second = particle2;
-        constraintDistances[i] = static_cast<RealOpenMM>(distance);
+        if (system.getParticleMass(particle1) != 0 || system.getParticleMass(particle2) != 0) {
+            constraintIndices.push_back(make_pair(particle1, particle2));
+            constraintDistances.push_back(distance);
+        }
     }
+    numConstraints = constraintIndices.size();
     perDofValues.resize(integrator.getNumPerDofVariables());
     for (int i = 0; i < (int) perDofValues.size(); i++)
         perDofValues[i].resize(numParticles);

platforms/reference/tests/TestReferenceBrownianIntegrator.cpp

@@ -164,6 +164,38 @@ void testConstraints() {
     }
 }
 
+void testConstrainedMasslessParticles() {
+    ReferencePlatform platform;
+    System system;
+    system.addParticle(0.0);
+    system.addParticle(1.0);
+    system.addConstraint(0, 1, 1.5);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(-1, 0, 0);
+    positions[1] = Vec3(1, 0, 0);
+    BrownianIntegrator integrator(300.0, 2.0, 0.01);
+    bool failed = false;
+    try {
+        // This should throw an exception.
+        
+        Context context(system, integrator, platform);
+    }
+    catch (exception& ex) {
+        failed = true;
+    }
+    ASSERT(failed);
+    
+    // Now make both particles massless, which should work.
+    
+    system.setParticleMass(1, 0.0);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    context.setVelocitiesToTemperature(300.0);
+    integrator.step(1);
+    State state = context.getState(State::Velocities | State::Positions);
+    ASSERT_EQUAL(0.0, state.getVelocities()[0][0]);
+}
+
 void testRandomSeed() {
     const int numParticles = 8;
     const double temp = 100.0;
@@ -228,6 +260,7 @@ int main() {
         testSingleBond();
         testTemperature();
         testConstraints();
+        testConstrainedMasslessParticles();
         testRandomSeed();
     }
     catch(const exception& e) {

platforms/reference/tests/TestReferenceCustomIntegrator.cpp

@@ -221,6 +221,44 @@ void testVelocityConstraints() {
     }
 }
 
+void testConstrainedMasslessParticles() {
+    ReferencePlatform platform;
+    System system;
+    system.addParticle(0.0);
+    system.addParticle(1.0);
+    system.addConstraint(0, 1, 1.5);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(-1, 0, 0);
+    positions[1] = Vec3(1, 0, 0);
+    CustomIntegrator integrator(0.002);
+    integrator.addPerDofVariable("oldx", 0);
+    integrator.addComputePerDof("v", "v+dt*f/m");
+    integrator.addComputePerDof("oldx", "x");
+    integrator.addComputePerDof("x", "x+dt*v");
+    integrator.addConstrainPositions();
+    integrator.addComputePerDof("v", "(x-oldx)/dt");
+    bool failed = false;
+    try {
+        // This should throw an exception.
+        
+        Context context(system, integrator, platform);
+    }
+    catch (exception& ex) {
+        failed = true;
+    }
+    ASSERT(failed);
+    
+    // Now make both particles massless, which should work.
+    
+    system.setParticleMass(1, 0.0);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    context.setVelocitiesToTemperature(300.0);
+    integrator.step(1);
+    State state = context.getState(State::Velocities | State::Positions);
+    ASSERT_EQUAL(0.0, state.getVelocities()[0][0]);
+}
+
 /**
  * Test an integrator with an AndersenThermostat to see if updateContextState()
  * is being handled correctly.
@@ -651,6 +689,7 @@ int main() {
         testSingleBond();
         testConstraints();
         testVelocityConstraints();
+        testConstrainedMasslessParticles();
         testWithThermostat();
         testMonteCarlo();
         testSum();

platforms/reference/tests/TestReferenceLangevinIntegrator.cpp

@@ -171,6 +171,38 @@ void testConstraints() {
     }
 }
 
+void testConstrainedMasslessParticles() {
+    ReferencePlatform platform;
+    System system;
+    system.addParticle(0.0);
+    system.addParticle(1.0);
+    system.addConstraint(0, 1, 1.5);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(-1, 0, 0);
+    positions[1] = Vec3(1, 0, 0);
+    LangevinIntegrator integrator(300.0, 2.0, 0.01);
+    bool failed = false;
+    try {
+        // This should throw an exception.
+        
+        Context context(system, integrator, platform);
+    }
+    catch (exception& ex) {
+        failed = true;
+    }
+    ASSERT(failed);
+    
+    // Now make both particles massless, which should work.
+    
+    system.setParticleMass(1, 0.0);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    context.setVelocitiesToTemperature(300.0);
+    integrator.step(1);
+    State state = context.getState(State::Velocities | State::Positions);
+    ASSERT_EQUAL(0.0, state.getVelocities()[0][0]);
+}
+
 void testRandomSeed() {
     const int numParticles = 8;
     const double temp = 100.0;
@@ -235,6 +267,7 @@ int main() {
         testSingleBond();
         testTemperature();
         testConstraints();
+        testConstrainedMasslessParticles();
         testRandomSeed();
     }
     catch(const exception& e) {

platforms/reference/tests/TestReferenceVariableLangevinIntegrator.cpp

@@ -173,6 +173,38 @@ void testConstraints() {
     }
 }
 
+void testConstrainedMasslessParticles() {
+    ReferencePlatform platform;
+    System system;
+    system.addParticle(0.0);
+    system.addParticle(1.0);
+    system.addConstraint(0, 1, 1.5);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(-1, 0, 0);
+    positions[1] = Vec3(1, 0, 0);
+    VariableLangevinIntegrator integrator(300.0, 2.0, 0.01);
+    bool failed = false;
+    try {
+        // This should throw an exception.
+        
+        Context context(system, integrator, platform);
+    }
+    catch (exception& ex) {
+        failed = true;
+    }
+    ASSERT(failed);
+    
+    // Now make both particles massless, which should work.
+    
+    system.setParticleMass(1, 0.0);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    context.setVelocitiesToTemperature(300.0);
+    integrator.step(1);
+    State state = context.getState(State::Velocities | State::Positions);
+    ASSERT_EQUAL(0.0, state.getVelocities()[0][0]);
+}
+
 void testRandomSeed() {
     const int numParticles = 8;
     const double temp = 100.0;
@@ -296,6 +328,7 @@ int main() {
         testSingleBond();
         testTemperature();
         testConstraints();
+        testConstrainedMasslessParticles();
         testRandomSeed();
         testArgonBox();
     }