Merge fixes for 8.5.1 (#5253)

* Fix leaked variable causing spurious template constraint assignment (#5236)

* fix 5234 and add test

* clean up docstring and standardize test name

* Take line search energy difference on CPU before reducing precision (#5242)

* Avoid multiple forces running on the worker thread at once (#5243)

* Fixed pressure calculation in MOnteCarloFlexibleBarostat (#5251)

---------
Co-authored-by: Jeff Wagner <jwagnerjpl@gmail.com>
Co-authored-by: Evan Pretti <pretti@stanford.edu>

olla/include/openmm/kernels.h

@@ -1994,10 +1994,10 @@ public:
     /**
      * Initialize the kernel.
      *
-     * @param system     the System this kernel will be applied to
+     * @param context    the ContextImpl this kernel will be applied to
      * @param force      the CustomCPPForceImpl this kernel will be used for
      */
-    virtual void initialize(const System& system, CustomCPPForceImpl& force) = 0;
+    virtual void initialize(const ContextImpl& context, CustomCPPForceImpl& force) = 0;
     /**
      * Execute the kernel to calculate the forces and/or energy.
      *
@@ -2022,10 +2022,10 @@ public:
     /**
      * Initialize the kernel.
      *
-     * @param system     the System this kernel will be applied to
+     * @param context    the ContextImpl this kernel will be applied to
      * @param force      the PythonForce this kernel will be used for
      */
-    virtual void initialize(const System& system, const PythonForce& force) = 0;
+    virtual void initialize(const ContextImpl& context, const PythonForce& force) = 0;
     /**
      * Execute the kernel to calculate the forces and/or energy.
      *

openmmapi/src/CustomCPPForceImpl.cpp

@@ -4,7 +4,7 @@
  * This is part of the OpenMM molecular simulation toolkit.                   *
  * See https://openmm.org/development.                                        *
  *                                                                            *
- * Portions copyright (c) 2008-2021 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2026 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -40,7 +40,7 @@ CustomCPPForceImpl::CustomCPPForceImpl(const Force& owner) {
 
 void CustomCPPForceImpl::initialize(ContextImpl& context) {
     kernel = context.getPlatform().createKernel(CalcCustomCPPForceKernel::Name(), context);
-    kernel.getAs<CalcCustomCPPForceKernel>().initialize(context.getSystem(), *this);
+    kernel.getAs<CalcCustomCPPForceKernel>().initialize(context, *this);
 }
 
 double CustomCPPForceImpl::calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy, int groups) {

openmmapi/src/MonteCarloFlexibleBarostatImpl.cpp

@@ -4,7 +4,7 @@
  * This is part of the OpenMM molecular simulation toolkit.                   *
  * See https://openmm.org/development.                                        *
  *                                                                            *
- * Portions copyright (c) 2010-2025 Stanford University and the Authors.      *
+ * Portions copyright (c) 2010-2026 Stanford University and the Authors.      *
  * Authors: Peter Eastman, Sander Vandenhaute                                 *
  * Contributors:                                                              *
  *                                                                            *
@@ -159,7 +159,7 @@ void MonteCarloFlexibleBarostatImpl::computeCurrentPressure(ContextImpl& context
     // Compute each component of the pressure tensor.
 
     for (int component = 0; component < 6; component++)
-        pressure[component] = (2.0*ke[component] - computePressureComponent(context, delta, component))/(volume*AVOGADRO*1e-25);
+        pressure[component] = (2.0*ke[component]/volume - computePressureComponent(context, delta, component))/(AVOGADRO*1e-25);
 
     // Restore the context to its original state.
 
@@ -223,7 +223,7 @@ double MonteCarloFlexibleBarostatImpl::computePressureComponent(ContextImpl& con
     // Compute the potential energy contribution to this element of the pressure tensor.
 
     double volume = box[0][0]*box[1][1]*box[2][2];
-    return (energy2-energy1)/(volume*2*delta);
+    return (energy1-energy2)/(volume*2*delta);
 }
 
 void MonteCarloFlexibleBarostatImpl::setBoxVectors(ContextImpl& context, Vec3 a, Vec3 b, Vec3 c) {

openmmapi/src/PythonForceImpl.cpp

@@ -4,7 +4,7 @@
  * This is part of the OpenMM molecular simulation toolkit.                   *
  * See https://openmm.org/development.                                        *
  *                                                                            *
- * Portions copyright (c) 2025 Stanford University and the Authors.           *
+ * Portions copyright (c) 2025-2026 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -47,7 +47,7 @@ PythonForceImpl::~PythonForceImpl() {
 
 void PythonForceImpl::initialize(ContextImpl& context) {
     kernel = context.getPlatform().createKernel(CalcPythonForceKernel::Name(), context);
-    kernel.getAs<CalcPythonForceKernel>().initialize(context.getSystem(), owner);
+    kernel.getAs<CalcPythonForceKernel>().initialize(context, owner);
 }
 
 double PythonForceImpl::calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy, int groups) {

platforms/common/include/openmm/common/CommonKernels.h

@@ -7,7 +7,7 @@
  * This is part of the OpenMM molecular simulation toolkit.                   *
  * See https://openmm.org/development.                                        *
  *                                                                            *
- * Portions copyright (c) 2008-2025 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2026 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -1469,10 +1469,10 @@ public:
     /**
      * Initialize the kernel.
      *
-     * @param system     the System this kernel will be applied to
+     * @param context    the ContextImpl this kernel will be applied to
      * @param force      the CustomCPPForceImpl this kernel will be used for
      */
-    void initialize(const System& system, CustomCPPForceImpl& force);
+    void initialize(const ContextImpl& context, CustomCPPForceImpl& force);
     /**
      * Execute the kernel to calculate the forces and/or energy.
      *
@@ -1507,6 +1507,7 @@ private:
     std::vector<float> floatForces;
     int forceGroupFlag;
     double energy;
+    bool useWorkerThread;
 };
 
 /**
@@ -1520,10 +1521,10 @@ public:
     /**
      * Initialize the kernel.
      *
-     * @param system     the System this kernel will be applied to
+     * @param context    the ContextImpl this kernel will be applied to
      * @param force      the PythonForce this kernel will be used for
      */
-    void initialize(const System& system, const PythonForce& force);
+    void initialize(const ContextImpl& context, const PythonForce& force);
     /**
      * Execute the kernel to calculate the forces and/or energy.
      *
@@ -1558,7 +1559,7 @@ private:
     std::vector<double> forcesVec;
     int forceGroupFlag;
     double energy;
-    bool usePeriodic;
+    bool usePeriodic, useWorkerThread;
 };
 
 } // namespace OpenMM

platforms/common/include/openmm/common/CommonMinimizeKernel.h

@@ -98,7 +98,7 @@ private:
     int maxIterations;
     MinimizationReporter* reporter;
 
-    double kRestraint, energy;
+    double kRestraint, energy, energyStart;
     bool largeGrad;
 
     ComputeArray constraintIndices, constraintDistances;

platforms/common/src/CommonKernels.cpp

@@ -4,7 +4,7 @@
  * This is part of the OpenMM molecular simulation toolkit.                   *
  * See https://openmm.org/development.                                        *
  *                                                                            *
- * Portions copyright (c) 2008-2025 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2026 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -34,6 +34,7 @@
 #include "openmm/internal/CustomCompoundBondForceImpl.h"
 #include "openmm/internal/DPDIntegratorUtilities.h"
 #include "openmm/internal/OSRngSeed.h"
+#include "openmm/internal/PythonForceImpl.h"
 #include "openmm/internal/ThreadPool.h"
 #include "openmm/internal/timer.h"
 #include "CommonKernelSources.h"
@@ -4688,10 +4689,10 @@ public:
     CommonCalcCustomCPPForceKernel& owner;
 };
 
-void CommonCalcCustomCPPForceKernel::initialize(const System& system, CustomCPPForceImpl& force) {
+void CommonCalcCustomCPPForceKernel::initialize(const ContextImpl& context, CustomCPPForceImpl& force) {
     ContextSelector selector(cc);
     this->force = &force;
-    int numParticles = system.getNumParticles();
+    int numParticles = context.getSystem().getNumParticles();
     forcesVec.resize(numParticles);
     positionsVec.resize(numParticles);
     floatForces.resize(3*numParticles);
@@ -4706,14 +4707,18 @@ void CommonCalcCustomCPPForceKernel::initialize(const System& system, CustomCPPF
     addForcesKernel->addArg(cc.getLongForceBuffer());
     addForcesKernel->addArg(cc.getAtomIndexArray());
     forceGroupFlag = (1<<force.getOwner().getForceGroup());
-    if (cc.getNumContexts() == 1) {
+    useWorkerThread = (cc.getNumContexts() == 1);
+    for (const ForceImpl* impl : context.getForceImpls())
+        if (dynamic_cast<const CustomCPPForceImpl*>(impl) != NULL || dynamic_cast<const PythonForceImpl*>(impl) != NULL)
+            useWorkerThread = false;
+    if (useWorkerThread) {
         cc.addPreComputation(new StartCalculationPreComputation(*this));
         cc.addPostComputation(new AddForcesPostComputation(*this));
     }
 }
 
 double CommonCalcCustomCPPForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
-    if (cc.getNumContexts() == 1) {
+    if (useWorkerThread) {
         // This method does nothing.  The actual calculation is started by the pre-computation, continued on
         // the worker thread, and finished by the post-computation.
 
@@ -4765,7 +4770,7 @@ double CommonCalcCustomCPPForceKernel::addForces(bool includeForces, bool includ
 
     // Wait until executeOnWorkerThread() is finished.
 
-    if (cc.getNumContexts() == 1)
+    if (useWorkerThread)
         cc.getWorkThread().flush();
 
     // Add in the forces.
@@ -4811,11 +4816,11 @@ public:
     CommonCalcPythonForceKernel& owner;
 };
 
-void CommonCalcPythonForceKernel::initialize(const System& system, const PythonForce& force) {
+void CommonCalcPythonForceKernel::initialize(const ContextImpl& context, const PythonForce& force) {
     ContextSelector selector(cc);
     computation = &force.getComputation();
     usePeriodic = force.usesPeriodicBoundaryConditions();
-    int numParticles = system.getNumParticles();
+    int numParticles = context.getSystem().getNumParticles();
     positionsVec.resize(numParticles);
     forcesVec.resize(3*numParticles);
     int elementSize = (cc.getUseDoublePrecision() ? sizeof(double) : sizeof(float));
@@ -4829,14 +4834,18 @@ void CommonCalcPythonForceKernel::initialize(const System& system, const PythonF
     addForcesKernel->addArg(cc.getLongForceBuffer());
     addForcesKernel->addArg(cc.getAtomIndexArray());
     forceGroupFlag = (1<<force.getForceGroup());
-    if (cc.getNumContexts() == 1) {
+    useWorkerThread = (cc.getNumContexts() == 1);
+    for (const ForceImpl* impl : context.getForceImpls())
+        if (dynamic_cast<const CustomCPPForceImpl*>(impl) != NULL || dynamic_cast<const PythonForceImpl*>(impl) != NULL)
+            useWorkerThread = false;
+    if (useWorkerThread) {
         cc.addPreComputation(new StartCalculationPreComputation(*this));
         cc.addPostComputation(new AddForcesPostComputation(*this));
     }
 }
 
 double CommonCalcPythonForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
-    if (cc.getNumContexts() == 1) {
+    if (useWorkerThread) {
         // This method does nothing.  The actual calculation is started by the pre-computation, continued on
         // the worker thread, and finished by the post-computation.
 
@@ -4857,9 +4866,9 @@ double CommonCalcPythonForceKernel::execute(ContextImpl& context, bool includeFo
 void CommonCalcPythonForceKernel::beginComputation(bool includeForces, bool includeEnergy, int groups) {
     if ((groups&forceGroupFlag) == 0)
         return;
-    
+
     // The actual force computation will be done on a different thread.
-    
+
     cc.getWorkThread().addTask(new ExecuteTask(*this, includeForces));
 }
 
@@ -4886,8 +4895,8 @@ double CommonCalcPythonForceKernel::addForces(bool includeForces, bool includeEn
         return 0;
 
     // Wait until executeOnWorkerThread() is finished.
-    
-    if (cc.getNumContexts() == 1)
+
+    if (useWorkerThread)
         cc.getWorkThread().flush();
 
     // Add in the forces.
@@ -4896,7 +4905,7 @@ double CommonCalcPythonForceKernel::addForces(bool includeForces, bool includeEn
         ContextSelector selector(cc);
         addForcesKernel->execute(cc.getNumAtoms());
     }
-    
+
     // Return the energy.
     
     return energy;

platforms/common/src/CommonMinimizeKernel.cpp

@@ -173,8 +173,8 @@ void CommonMinimizeKernel::setup(ContextImpl& context) {
     returnFlag.initialize<int>(cc, 1, "returnFlag");
     returnValue.initialize(cc, 1, elementSize, "returnValue");
     gradNorm.initialize(cc, 1, elementSize, "gradNorm");
-    lineSearchData.initialize(cc, 4, elementSize, "lineSearchData");
-    lineSearchDataBackup.initialize(cc, 4, elementSize, "lineSearchDataBackup");
+    lineSearchData.initialize(cc, 3, elementSize, "lineSearchData");
+    lineSearchDataBackup.initialize(cc, 3, elementSize, "lineSearchDataBackup");
 
     // Compile kernels and set arguments.
 
@@ -341,7 +341,6 @@ void CommonMinimizeKernel::setup(ContextImpl& context) {
     lineSearchSetupKernel->addArg(gradNorm);
     lineSearchSetupKernel->addArg(lineSearchData);
     lineSearchSetupKernel->addArg(numVariables);
-    lineSearchSetupKernel->addArg(); // energyStart
 
     lineSearchStepKernel = program->createKernel("lineSearchStep");
     lineSearchStepKernel->addArg(x);
@@ -405,12 +404,7 @@ void CommonMinimizeKernel::lbfgs(ContextImpl& context) {
     for (int iteration = 1, end = 0;;) {
         // Prepare for a line search.
 
-        if (mixedIsDouble) {
-            lineSearchSetupKernel->setArg(9, energy);
-        }
-        else {
-            lineSearchSetupKernel->setArg(9, (float) energy);
-        }
+        energyStart = energy;
         lineSearchSetupKernel->execute(numVariables);
 
         // Take line search steps.
@@ -746,11 +740,10 @@ double CommonMinimizeKernel::downloadGradNormSync() {
 
 void CommonMinimizeKernel::runLineSearchKernels() {
     if (mixedIsDouble) {
-        lineSearchDotKernel->setArg(6, isfinite(energy) ? energy : (double) std::numeric_limits<float>::max());
+        lineSearchDotKernel->setArg(6, isfinite(energy) ? energy - energyStart : (double) std::numeric_limits<float>::max());
     }
     else {
-        float hostEnergy = (float) energy;
-        lineSearchDotKernel->setArg(6, isfinite(hostEnergy) ? hostEnergy : std::numeric_limits<float>::max());
+        lineSearchDotKernel->setArg(6, isfinite((float) energy) ? (float) (energy - energyStart) : std::numeric_limits<float>::max());
     }
     lineSearchDotKernel->execute(numVariables);
     lineSearchContinueKernel->execute(1);

platforms/common/src/kernels/minimize.cc

 #define LS_DOT_START 0
 #define LS_DOT 1
-#define LS_ENERGY 2
-#define LS_STEP 3
+#define LS_STEP 2
 
 #define LS_FAIL 0
 #define LS_SUCCEED 1
@@ -668,8 +667,7 @@ KERNEL void lineSearchSetup(
     GLOBAL int* RESTRICT returnFlag,
     GLOBAL mixed* RESTRICT gradNorm,
     GLOBAL mixed* RESTRICT lineSearchData,
-    const int numVariables,
-    const mixed energyStart
+    const int numVariables
 ) {
     LOCAL volatile mixed temp[TEMP_SIZE];
 
@@ -694,7 +692,6 @@ KERNEL void lineSearchSetup(
     if (GLOBAL_ID == 0) {
         *returnFlag = LS_CONTINUE;
         *gradNorm = 0;
-        lineSearchData[LS_ENERGY] = energyStart;
     }
 }
 
@@ -758,7 +755,6 @@ KERNEL void lineSearchStep(
 
         lineSearchDataBackup[LS_DOT_START] = lineSearchData[LS_DOT_START];
         lineSearchDataBackup[LS_DOT] = lineSearchData[LS_DOT] = 0;
-        lineSearchDataBackup[LS_ENERGY] = lineSearchData[LS_ENERGY];
         lineSearchDataBackup[LS_STEP] = lineSearchData[LS_STEP];
     }
 }
@@ -770,7 +766,7 @@ KERNEL void lineSearchDot(
     GLOBAL int* RESTRICT returnFlag,
     GLOBAL const mixed* RESTRICT returnValue,
     const int numVariables,
-    mixed energy
+    mixed deltaEnergy
 ) {
     LOCAL volatile mixed temp[TEMP_SIZE];
 
@@ -781,13 +777,13 @@ KERNEL void lineSearchDot(
     // Any restraint energy in returnValue hasn't been downloaded yet to be
     // passed back up in the energy parameter, so add it in here.
 
-    energy += *returnValue;
+    deltaEnergy += *returnValue;
 
     // The energy may be such that we don't need to do a dot product and can
     // immediately decide to scale the step, so mark this case with LS_SUCCEED.
     // This will be checked in the following kernel.
 
-    if (!(FABS_MIXED(energy) < FLT_MAX) || energy > lineSearchData[LS_ENERGY] + lineSearchData[LS_STEP] * LBFGS_FTOL * lineSearchData[LS_DOT_START]) {
+    if (!(FABS_MIXED(deltaEnergy) < FLT_MAX) || deltaEnergy > lineSearchData[LS_STEP] * LBFGS_FTOL * lineSearchData[LS_DOT_START]) {
         if (GLOBAL_ID == 0) {
             *returnFlag = LS_SUCCEED;
         }

platforms/reference/include/ReferenceKernels.h

@@ -1992,10 +1992,10 @@ public:
     /**
      * Initialize the kernel.
      *
-     * @param system     the System this kernel will be applied to
+     * @param context    the ContextImpl this kernel will be applied to
      * @param force      the CustomCPPForceImpl this kernel will be used for
      */
-    void initialize(const System& system, CustomCPPForceImpl& force);
+    void initialize(const ContextImpl& context, CustomCPPForceImpl& force);
     /**
      * Execute the kernel to calculate the forces and/or energy.
      *
@@ -2020,10 +2020,10 @@ public:
     /**
      * Initialize the kernel.
      *
-     * @param system     the System this kernel will be applied to
+     * @param context    the ContextImpl this kernel will be applied to
      * @param force      the PythonForce this kernel will be used for
      */
-    void initialize(const System& system, const PythonForce& force);
+    void initialize(const ContextImpl& context, const PythonForce& force);
     /**
      * Execute the kernel to calculate the forces and/or energy.
      *

platforms/reference/src/ReferenceKernels.cpp

@@ -3535,9 +3535,9 @@ void ReferenceCalcATMForceKernel::copyParametersToContext(ContextImpl& context,
     loadParams(numParticles, force);
 }
 
-void ReferenceCalcCustomCPPForceKernel::initialize(const System& system, CustomCPPForceImpl& force) {
+void ReferenceCalcCustomCPPForceKernel::initialize(const ContextImpl& context, CustomCPPForceImpl& force) {
     this->force = &force;
-    forces.resize(system.getNumParticles());
+    forces.resize(context.getSystem().getNumParticles());
 }
 
 double ReferenceCalcCustomCPPForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
@@ -3550,9 +3550,9 @@ double ReferenceCalcCustomCPPForceKernel::execute(ContextImpl& context, bool inc
     return energy;
 }
 
-void ReferenceCalcPythonForceKernel::initialize(const System& system, const PythonForce& force) {
+void ReferenceCalcPythonForceKernel::initialize(const ContextImpl& context, const PythonForce& force) {
     computation = &force.getComputation();
-    forces.resize(system.getNumParticles());
+    forces.resize(context.getSystem().getNumParticles());
     usePeriodic = force.usesPeriodicBoundaryConditions();
 }
 

tests/TestMonteCarloFlexibleBarostat.h

@@ -4,7 +4,7 @@
  * This is part of the OpenMM molecular simulation toolkit.                   *
  * See https://openmm.org/development.                                        *
  *                                                                            *
- * Portions copyright (c) 2008-2021 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2026 Stanford University and the Authors.      *
  * Authors: Peter Eastman, Lee-Ping Wang                                      *
  * Contributors:                                                              *
  *                                                                            *
@@ -180,12 +180,12 @@ void testMoleculeScaling(bool rigid) {
 }
 
 void testMolecularGas(bool rigid) {
-    const int numMolecules = 64;
+    const int numMolecules = 256;
     const int frequency = 5;
     const int steps = 5000;
     const double pressure = 3.0;
     const double pressureInMD = pressure*(AVOGADRO*1e-25); // pressure in kJ/mol/nm^3
-    const double temp =300.0;
+    const double temp = 300.0;
     const double initialVolume = numMolecules*BOLTZ*temp/pressureInMD;
     const double initialLength = std::pow(initialVolume, 1.0/3.0);
 
@@ -206,8 +206,8 @@ void testMolecularGas(bool rigid) {
         system.addParticle(1.0);
         system.addParticle(1.0);
         Vec3 pos(initialLength*genrand_real2(sfmt), 0.5*initialLength*genrand_real2(sfmt), 2*initialLength*genrand_real2(sfmt));
-        bonds->addBond(positions.size(), positions.size()+1, 0.1, 0.0);
-        bonds->addBond(positions.size(), positions.size()+2, 0.1, 0.0);
+        bonds->addBond(positions.size(), positions.size()+1, 0.1, 1.0);
+        bonds->addBond(positions.size(), positions.size()+2, 0.1, 1.0);
         positions.push_back(pos);
         positions.push_back(pos+Vec3(0.1, 0.0, 0.0));
         positions.push_back(pos+Vec3(0.0, 0.1, 0.0));

wrappers/python/openmm/app/forcefield.py

@@ -1423,7 +1423,7 @@ class ForceField(object):
                     continue
                 t1 = data.atomTemplateIndexes[atom1]
                 t2 = data.atomTemplateIndexes[atom2]
-                for constraint in template.constraints:
+                for constraint in matchedTemplateForResidue[atom1.residue].constraints:
                     if sorted((t1, t2)) == sorted((constraint[0], constraint[1])):
                         bond.isConstrained = True
 
@@ -4719,4 +4719,4 @@ class DrudeGenerator(object):
             sys.setParticleMass(particle, drudeMass)
             sys.setParticleMass(parent, sys.getParticleMass(parent)-transferMass)
 
-parsers["DrudeForce"] = DrudeGenerator.parseElement
\ No newline at end of file
+parsers["DrudeForce"] = DrudeGenerator.parseElement

wrappers/python/tests/TestForceField.py

@@ -262,6 +262,75 @@ class TestForceField(unittest.TestCase):
                 self.assertEqual(1.0, lengths[(0, 4)])
                 self.assertEqual(1.0, lengths[(0, 5)])
 
+    def testTemplateConstraintsMultipleMols(self):
+        """Test that constraints defined by a residue template don't leak into
+        other residues.
+        See https://github.com/openmm/openmm/issues/5234
+        """
+        MOL_A = """<ForceField>
+ <AtomTypes>
+  <Type name="A0" mass="12" class="A0" element="C"/>
+  <Type name="A1" mass="12" class="A1" element="C"/>
+ </AtomTypes>
+ <HarmonicBondForce>
+  <Bond class1="A0" class2="A1" length="0.15" k="200000"/>
+ </HarmonicBondForce>
+ <NonbondedForce coulomb14scale="1" lj14scale="1">
+  <Atom class="A0" sigma="0.3" epsilon="0.1" charge="0"/>
+  <Atom class="A1" sigma="0.3" epsilon="0.1" charge="0"/>
+ </NonbondedForce>
+ <Residues>
+  <Residue name="MOLA">
+   <Atom name="a0" type="A0"/><Atom name="a1" type="A1"/>
+   <Bond atomName1="a0" atomName2="a1"/>
+  </Residue>
+ </Residues>
+</ForceField>"""
+
+        MOL_B = """<ForceField>
+ <AtomTypes>
+  <Type name="B0" mass="12" class="B0" element="C"/>
+  <Type name="B1" mass="1" class="B1" element="H"/>
+ </AtomTypes>
+ <HarmonicBondForce>
+  <Bond class1="B0" class2="B1" length="0.11" k="300000"/>
+ </HarmonicBondForce>
+ <NonbondedForce coulomb14scale="1" lj14scale="1">
+  <Atom class="B0" sigma="0.3" epsilon="0.1" charge="0"/>
+  <Atom class="B1" sigma="0.1" epsilon="0.01" charge="0"/>
+ </NonbondedForce>
+ <Residues>
+  <Residue name="MOLB">
+   <Atom name="b0" type="B0"/><Atom name="b1" type="B1"/>
+   <Bond atomName1="b0" atomName2="b1"/>
+   <Constraint atomName1="b0" atomName2="b1" distance="0.11"/>
+  </Residue>
+ </Residues>
+</ForceField>"""
+
+        top = Topology()
+        c = top.addChain()
+        C, H = Element.getBySymbol('C'), Element.getBySymbol('H')
+        r1 = top.addResidue('MOLA', c)
+        a0, a1 = top.addAtom('a0', C, r1), top.addAtom('a1', C, r1)
+        top.addBond(a0, a1)
+        r2 = top.addResidue('MOLB', c)
+        b0, b1 = top.addAtom('b0', C, r2), top.addAtom('b1', H, r2)
+        top.addBond(b0, b1)
+
+        ff = ForceField()
+        ff.loadFile(StringIO(MOL_A))
+        ff.loadFile(StringIO(MOL_B))
+        sys = ff.createSystem(top,
+                              nonbondedMethod=NoCutoff,
+                              constraints=None)
+
+        self.assertEqual(sys.getNumConstraints(), 1)
+        constraintParameters = sys.getConstraintParameters(0)
+        self.assertEqual(constraintParameters[0], 2)
+        self.assertEqual(constraintParameters[1], 3)
+        self.assertEqual(constraintParameters[2], 0.11 * nanometer)
+
     def test_ImplicitSolvent(self):
         """Test the four types of implicit solvents using the implicitSolvent
         parameter.