Merge OpenMM master into dtk26.04

plugins/drude/platforms/common/src/CommonDrudeKernels.cpp

@@ -502,7 +502,10 @@ void CommonIntegrateDrudeSCFStepKernel::minimize(ContextImpl& context, double to
     int numDrude = drudeParams.getSize();
     int paddedNumAtoms = cc.getPaddedNumAtoms();
     for (int iteration = 0; iteration < 50; iteration++) {
-        context.calcForcesAndEnergy(true, false, context.getIntegrator().getIntegrationForceGroups());
+        {
+            ContextDeselector deselector(cc);
+            context.calcForcesAndEnergy(true, false, context.getIntegrator().getIntegrationForceGroups());
+        }
         minimizeKernel->execute(drudeParams.getSize());
         cc.getLongForceBuffer().download(forces);
         double totalForce = 0;

plugins/rpmd/platforms/common/src/CommonRpmdKernels.cpp

@@ -297,7 +297,10 @@ void CommonIntegrateRPMDStepKernel::computeForces(ContextImpl& context) {
         context.getPeriodicBoxVectors(finalBox[0], finalBox[1], finalBox[2]);
         if (initialBox[0] != finalBox[0] || initialBox[1] != finalBox[1] || initialBox[2] != finalBox[2])
             throw OpenMMException("Standard barostats cannot be used with RPMDIntegrator.  Use RPMDMonteCarloBarostat instead.");
-        context.calcForcesAndEnergy(true, false, groupsNotContracted);
+        {
+            ContextDeselector deselector(cc);
+            context.calcForcesAndEnergy(true, false, groupsNotContracted);
+        }
         copyFromContextKernel->setArg(7, i);
         copyFromContextKernel->execute(cc.getNumAtoms());
     }
@@ -322,7 +325,10 @@ void CommonIntegrateRPMDStepKernel::computeForces(ContextImpl& context) {
                 copyToContextKernel->setArg(5, i);
                 copyToContextKernel->execute(cc.getNumAtoms());
                 context.computeVirtualSites();
-                context.calcForcesAndEnergy(true, false, groupFlags);
+                {
+                    ContextDeselector deselector(cc);
+                    context.calcForcesAndEnergy(true, false, groupFlags);
+                }
                 copyFromContextKernel->setArg(7, i);
                 copyFromContextKernel->execute(cc.getNumAtoms());
             }

serialization/src/MonteCarloAnisotropicBarostatProxy.cpp

@@ -4,7 +4,7 @@
  * This is part of the OpenMM molecular simulation toolkit.                   *
  * See https://openmm.org/development.                                        *
  *                                                                            *
- * Portions copyright (c) 2010-2016 Stanford University and the Authors.      *
+ * Portions copyright (c) 2010-2026 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -54,6 +54,7 @@ void MonteCarloAnisotropicBarostatProxy::serialize(const void* object, Serializa
     node.setDoubleProperty("temperature", force.getDefaultTemperature());
     node.setIntProperty("frequency", force.getFrequency());
     node.setIntProperty("randomSeed", force.getRandomNumberSeed());
+    node.setBoolProperty("rigidScaling", force.getScaleMoleculesAsRigid());
 }
 
 void* MonteCarloAnisotropicBarostatProxy::deserialize(const SerializationNode& node) const {
@@ -63,7 +64,7 @@ void* MonteCarloAnisotropicBarostatProxy::deserialize(const SerializationNode& n
     try {
         Vec3 pressure(node.getDoubleProperty("pressurex"), node.getDoubleProperty("pressurey"), node.getDoubleProperty("pressurez"));
         force = new MonteCarloAnisotropicBarostat(pressure, node.getDoubleProperty("temperature"), node.getBoolProperty("scalex"),
-                node.getBoolProperty("scaley"), node.getBoolProperty("scalez"), node.getIntProperty("frequency"));
+                node.getBoolProperty("scaley"), node.getBoolProperty("scalez"), node.getIntProperty("frequency"), node.getBoolProperty("rigidScaling", true));
         force->setForceGroup(node.getIntProperty("forceGroup", 0));
         force->setName(node.getStringProperty("name", force->getName()));
         force->setRandomNumberSeed(node.getIntProperty("randomSeed"));

serialization/src/MonteCarloBarostatProxy.cpp

@@ -4,7 +4,7 @@
  * This is part of the OpenMM molecular simulation toolkit.                   *
  * See https://openmm.org/development.                                        *
  *                                                                            *
- * Portions copyright (c) 2010-2016 Stanford University and the Authors.      *
+ * Portions copyright (c) 2010-2026 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -48,6 +48,7 @@ void MonteCarloBarostatProxy::serialize(const void* object, SerializationNode& n
     node.setDoubleProperty("temperature", force.getDefaultTemperature());
     node.setIntProperty("frequency", force.getFrequency());
     node.setIntProperty("randomSeed", force.getRandomNumberSeed());
+    node.setBoolProperty("rigidScaling", force.getScaleMoleculesAsRigid());
 }
 
 void* MonteCarloBarostatProxy::deserialize(const SerializationNode& node) const {
@@ -55,7 +56,8 @@ void* MonteCarloBarostatProxy::deserialize(const SerializationNode& node) const
         throw OpenMMException("Unsupported version number");
     MonteCarloBarostat* force = NULL;
     try {
-        force = new MonteCarloBarostat(node.getDoubleProperty("pressure"), node.getDoubleProperty("temperature"), node.getIntProperty("frequency"));
+        force = new MonteCarloBarostat(node.getDoubleProperty("pressure"), node.getDoubleProperty("temperature"), node.getIntProperty("frequency"),
+                                       node.getBoolProperty("rigidScaling", true));
         force->setForceGroup(node.getIntProperty("forceGroup", 0));
         force->setName(node.getStringProperty("name", force->getName()));
         force->setRandomNumberSeed(node.getIntProperty("randomSeed"));

serialization/src/MonteCarloMembraneBarostatProxy.cpp

@@ -4,7 +4,7 @@
  * This is part of the OpenMM molecular simulation toolkit.                   *
  * See https://openmm.org/development.                                        *
  *                                                                            *
- * Portions copyright (c) 2010-2016 Stanford University and the Authors.      *
+ * Portions copyright (c) 2010-2026 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -51,6 +51,7 @@ void MonteCarloMembraneBarostatProxy::serialize(const void* object, Serializatio
     node.setIntProperty("zmode", force.getZMode());
     node.setIntProperty("frequency", force.getFrequency());
     node.setIntProperty("randomSeed", force.getRandomNumberSeed());
+    node.setBoolProperty("rigidScaling", force.getScaleMoleculesAsRigid());
 }
 
 void* MonteCarloMembraneBarostatProxy::deserialize(const SerializationNode& node) const {
@@ -61,7 +62,7 @@ void* MonteCarloMembraneBarostatProxy::deserialize(const SerializationNode& node
         MonteCarloMembraneBarostat::XYMode xymode = (MonteCarloMembraneBarostat::XYMode) node.getIntProperty("xymode");
         MonteCarloMembraneBarostat::ZMode zmode = (MonteCarloMembraneBarostat::ZMode) node.getIntProperty("zmode");
         force = new MonteCarloMembraneBarostat(node.getDoubleProperty("pressure"), node.getDoubleProperty("surfaceTension"),
-                node.getDoubleProperty("temperature"), xymode, zmode, node.getIntProperty("frequency"));
+                node.getDoubleProperty("temperature"), xymode, zmode, node.getIntProperty("frequency"), node.getBoolProperty("rigidScaling", true));
         force->setForceGroup(node.getIntProperty("forceGroup", 0));
         force->setName(node.getStringProperty("name", force->getName()));
         force->setRandomNumberSeed(node.getIntProperty("randomSeed"));

serialization/src/g_fmt.cpp

@@ -33,10 +33,10 @@ extern "C" {
 #endif
 
  char *
-g_fmt(register char *b, double x)
+g_fmt(char *b, double x)
 {
-	register int i, k;
-	register char *s;
+	int i, k;
+	char *s;
 	int decpt, j, sign;
 	char *b0, *s0, *se;
 

serialization/tests/TestSerializeMonteCarloAnisotropicBarostat.cpp

@@ -4,7 +4,7 @@
  * This is part of the OpenMM molecular simulation toolkit.                   *
  * See https://openmm.org/development.                                        *
  *                                                                            *
- * Portions copyright (c) 2010-2016 Stanford University and the Authors.      *
+ * Portions copyright (c) 2010-2026 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -39,7 +39,7 @@ using namespace std;
 void testSerialization() {
     // Create a Force.
 
-    MonteCarloAnisotropicBarostat force(Vec3(15.1, 18.2, 19.3), 250.0, true, false, true, 14);
+    MonteCarloAnisotropicBarostat force(Vec3(15.1, 18.2, 19.3), 250.0, true, false, true, 14, false);
     force.setForceGroup(3);
     force.setName("custom name");
     force.setRandomNumberSeed(3);
@@ -62,6 +62,7 @@ void testSerialization() {
     ASSERT_EQUAL(force.getScaleZ(), force2.getScaleZ());
     ASSERT_EQUAL(force.getFrequency(), force2.getFrequency());
     ASSERT_EQUAL(force.getRandomNumberSeed(), force2.getRandomNumberSeed());
+    ASSERT_EQUAL(force.getScaleMoleculesAsRigid(), force2.getScaleMoleculesAsRigid());
 }
 
 int main() {

serialization/tests/TestSerializeMonteCarloBarostat.cpp

@@ -4,7 +4,7 @@
  * This is part of the OpenMM molecular simulation toolkit.                   *
  * See https://openmm.org/development.                                        *
  *                                                                            *
- * Portions copyright (c) 2010-2016 Stanford University and the Authors.      *
+ * Portions copyright (c) 2010-2026 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -39,7 +39,7 @@ using namespace std;
 void testSerialization() {
     // Create a Force.
 
-    MonteCarloBarostat force(25.5, 250.0, 14);
+    MonteCarloBarostat force(25.5, 250.0, 14, false);
     force.setForceGroup(3);
     force.setName("custom name");
     force.setRandomNumberSeed(3);
@@ -59,6 +59,7 @@ void testSerialization() {
     ASSERT_EQUAL(force.getDefaultTemperature(), force2.getDefaultTemperature());
     ASSERT_EQUAL(force.getFrequency(), force2.getFrequency());
     ASSERT_EQUAL(force.getRandomNumberSeed(), force2.getRandomNumberSeed());
+    ASSERT_EQUAL(force.getScaleMoleculesAsRigid(), force2.getScaleMoleculesAsRigid());
 }
 
 int main() {

serialization/tests/TestSerializeMonteCarloMembraneBarostat.cpp

@@ -4,7 +4,7 @@
  * This is part of the OpenMM molecular simulation toolkit.                   *
  * See https://openmm.org/development.                                        *
  *                                                                            *
- * Portions copyright (c) 2010-2016 Stanford University and the Authors.      *
+ * Portions copyright (c) 2010-2026 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -39,7 +39,7 @@ using namespace std;
 void testSerialization() {
     // Create a Force.
 
-    MonteCarloMembraneBarostat force(25.5, 11.2, 250.0, MonteCarloMembraneBarostat::XYAnisotropic, MonteCarloMembraneBarostat::ZFixed, 14);
+    MonteCarloMembraneBarostat force(25.5, 11.2, 250.0, MonteCarloMembraneBarostat::XYAnisotropic, MonteCarloMembraneBarostat::ZFixed, 14, false);
     force.setForceGroup(3);
     force.setName("custom name");
     force.setRandomNumberSeed(3);
@@ -62,6 +62,7 @@ void testSerialization() {
     ASSERT_EQUAL(force.getZMode(), force2.getZMode());
     ASSERT_EQUAL(force.getFrequency(), force2.getFrequency());
     ASSERT_EQUAL(force.getRandomNumberSeed(), force2.getRandomNumberSeed());
+    ASSERT_EQUAL(force.getScaleMoleculesAsRigid(), force2.getScaleMoleculesAsRigid());
 }
 
 int main() {

tests/TestMonteCarloAnisotropicBarostat.h

@@ -4,7 +4,7 @@
  * This is part of the OpenMM molecular simulation toolkit.                   *
  * See https://openmm.org/development.                                        *
  *                                                                            *
- * Portions copyright (c) 2008-2016 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2026 Stanford University and the Authors.      *
  * Authors: Peter Eastman, Lee-Ping Wang                                      *
  * Contributors:                                                              *
  *                                                                            *
@@ -29,6 +29,7 @@
 
 #include "openmm/internal/AssertionUtilities.h"
 #include "openmm/CustomExternalForce.h"
+#include "openmm/CustomIntegrator.h"
 #include "openmm/HarmonicBondForce.h"
 #include "openmm/MonteCarloBarostat.h"
 #include "openmm/MonteCarloAnisotropicBarostat.h"
@@ -184,8 +185,75 @@ void testIdealGasAxis(int axis) {
     }
 }
 
-void testMolecularGas() {
-    const int numMolecules = 64;
+void testMoleculeScaling(bool rigid) {
+    int numMolecules = 10;
+    double initialWidth = 3.0;
+
+    // Create a system of diatomic molecules.
+
+    System system;
+    Vec3 initialBox[] = {Vec3(initialWidth, 0, 0), Vec3(0, initialWidth, 0), Vec3(0, 0, initialWidth)};
+    system.setDefaultPeriodicBoxVectors(initialBox[0], initialBox[1], initialBox[2]);
+    HarmonicBondForce* bonds = new HarmonicBondForce();
+    system.addForce(bonds);
+    NonbondedForce* nonbonded = new NonbondedForce();
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    system.addForce(nonbonded);
+    MonteCarloAnisotropicBarostat* barostat = new MonteCarloAnisotropicBarostat(Vec3(1.0, 1.0, 1.0), 300.0, true, true, true, 1, rigid);
+    system.addForce(barostat);
+    vector<Vec3> positions;
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    for (int i = 0; i < numMolecules; i++) {
+        system.addParticle(1.0);
+        system.addParticle(1.0);
+        bonds->addBond(2*i, 2*i+1, 0.2, 1000.0);
+        nonbonded->addParticle(0.0, 0.1, 1.0);
+        nonbonded->addParticle(0.0, 0.1, 1.0);
+        Vec3 pos1(initialWidth*genrand_real2(sfmt), initialWidth*genrand_real2(sfmt), initialWidth*genrand_real2(sfmt));
+        Vec3 delta(genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5);
+        delta /= sqrt(delta.dot(delta));
+        positions.push_back(pos1);
+        positions.push_back(pos1+delta);
+    }
+
+    // Use an integrator that applies the barostat but nothing else.
+
+    CustomIntegrator integrator(1.0);
+    integrator.addUpdateContextState();
+
+    // Let the barostat make some moves.
+
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    integrator.step(100);
+    State state = context.getState(State::Positions);
+
+    // The box size should have changed.
+
+    Vec3 finalBox[3];
+    state.getPeriodicBoxVectors(finalBox[0], finalBox[1], finalBox[2]);
+    for (int i = 0; i < 3; i++)
+        ASSERT(finalBox[i][i] != initialBox[i][i]);
+
+    // See if the molecules were scaled correctly.
+
+    Vec3 boxScale(finalBox[0][0]/initialBox[0][0], finalBox[1][1]/initialBox[1][1], finalBox[2][2]/initialBox[2][2]);
+    for (int i = 0; i < numMolecules; i++) {
+        Vec3 delta1 = positions[2*i+1]-positions[2*i];
+        Vec3 delta2 = state.getPositions()[2*i+1]-state.getPositions()[2*i];
+        if (rigid) {
+            ASSERT_EQUAL_VEC(delta1, delta2, 1e-5);
+        }
+        else {
+            Vec3 expected(delta1[0]*boxScale[0], delta1[1]*boxScale[1], delta1[2]*boxScale[2]);
+            ASSERT_EQUAL_VEC(expected, delta2, 1e-5);
+        }
+    }
+}
+
+void testMolecularGas(bool rigid) {
+    const int numMolecules = 256;
     const int frequency = 5;
     const int steps = 5000;
     const double pressure = 3.0;
@@ -198,7 +266,7 @@ void testMolecularGas() {
 
     System system;
     system.setDefaultPeriodicBoxVectors(Vec3(initialLength, 0, 0), Vec3(0, 0.5*initialLength, 0), Vec3(0, 0, 2*initialLength));
-    MonteCarloAnisotropicBarostat* barostat = new MonteCarloAnisotropicBarostat(Vec3(pressure, pressure, pressure), temp, true, true, true, frequency);
+    MonteCarloAnisotropicBarostat* barostat = new MonteCarloAnisotropicBarostat(Vec3(pressure, pressure, pressure), temp, true, true, true, frequency, rigid);
     system.addForce(barostat);
     HarmonicBondForce* bonds = new HarmonicBondForce();
     bonds->setUsesPeriodicBoundaryConditions(true);
@@ -211,8 +279,8 @@ void testMolecularGas() {
         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));
-        system.addConstraint(positions.size(), positions.size()+1, 0.1);
-        system.addConstraint(positions.size(), positions.size()+2, 0.1);
+        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));
@@ -565,7 +633,10 @@ int main(int argc, char* argv[]) {
         testIdealGasAxis(0);
         testIdealGasAxis(1);
         testIdealGasAxis(2);
-        testMolecularGas();
+        testMoleculeScaling(true);
+        testMoleculeScaling(false);
+        testMolecularGas(true);
+        testMolecularGas(false);
         testRandomSeed();
         testTriclinic();
         //testEinsteinCrystal();

tests/TestMonteCarloBarostat.h

@@ -4,7 +4,7 @@
  * This is part of the OpenMM molecular simulation toolkit.                   *
  * See https://openmm.org/development.                                        *
  *                                                                            *
- * Portions copyright (c) 2008-2016 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2026 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -30,6 +30,7 @@
 #include "openmm/internal/AssertionUtilities.h"
 #include "openmm/MonteCarloBarostat.h"
 #include "openmm/Context.h"
+#include "openmm/CustomIntegrator.h"
 #include "openmm/HarmonicBondForce.h"
 #include "openmm/NonbondedForce.h"
 #include "openmm/System.h"
@@ -141,9 +142,75 @@ void testIdealGas() {
     }
 }
 
-void testMolecularGas() {
-    const int numMolecules = 64;
-    const int frequency = 10;
+void testMoleculeScaling(bool rigid) {
+    int numMolecules = 10;
+    double initialWidth = 3.0;
+
+    // Create a system of diatomic molecules.
+
+    System system;
+    Vec3 initialBox[] = {Vec3(initialWidth, 0, 0), Vec3(0, initialWidth, 0), Vec3(0, 0, initialWidth)};
+    system.setDefaultPeriodicBoxVectors(initialBox[0], initialBox[1], initialBox[2]);
+    HarmonicBondForce* bonds = new HarmonicBondForce();
+    system.addForce(bonds);
+    NonbondedForce* nonbonded = new NonbondedForce();
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    system.addForce(nonbonded);
+    MonteCarloBarostat* barostat = new MonteCarloBarostat(1.0, 300.0, 1, rigid);
+    system.addForce(barostat);
+    vector<Vec3> positions;
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    for (int i = 0; i < numMolecules; i++) {
+        system.addParticle(1.0);
+        system.addParticle(1.0);
+        bonds->addBond(2*i, 2*i+1, 0.2, 1000.0);
+        nonbonded->addParticle(0.0, 0.1, 1.0);
+        nonbonded->addParticle(0.0, 0.1, 1.0);
+        Vec3 pos1(initialWidth*genrand_real2(sfmt), initialWidth*genrand_real2(sfmt), initialWidth*genrand_real2(sfmt));
+        Vec3 delta(genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5, genrand_real2(sfmt)-0.5);
+        delta /= sqrt(delta.dot(delta));
+        positions.push_back(pos1);
+        positions.push_back(pos1+delta);
+    }
+
+    // Use an integrator that applies the barostat but nothing else.
+
+    CustomIntegrator integrator(1.0);
+    integrator.addUpdateContextState();
+
+    // Let the barostat make some moves.
+
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    integrator.step(100);
+    State state = context.getState(State::Positions);
+
+    // The box size should have changed.
+
+    Vec3 finalBox[3];
+    state.getPeriodicBoxVectors(finalBox[0], finalBox[1], finalBox[2]);
+    ASSERT(finalBox[0][0] != initialBox[0][0]);
+
+    // See if the molecules were scaled correctly.
+
+    Vec3 boxScale(finalBox[0][0]/initialBox[0][0], finalBox[1][1]/initialBox[1][1], finalBox[2][2]/initialBox[2][2]);
+    for (int i = 0; i < numMolecules; i++) {
+        Vec3 delta1 = positions[2*i+1]-positions[2*i];
+        Vec3 delta2 = state.getPositions()[2*i+1]-state.getPositions()[2*i];
+        if (rigid) {
+            ASSERT_EQUAL_VEC(delta1, delta2, 1e-5);
+        }
+        else {
+            Vec3 expected(delta1[0]*boxScale[0], delta1[1]*boxScale[1], delta1[2]*boxScale[2]);
+            ASSERT_EQUAL_VEC(expected, delta2, 1e-5);
+        }
+    }
+}
+
+void testMolecularGas(bool rigid) {
+    const int numMolecules = 256;
+    const int frequency = 5;
     const int steps = 1000;
     const double pressure = 1.5;
     const double pressureInMD = pressure*(AVOGADRO*1e-25); // pressure in kJ/mol/nm^3
@@ -155,7 +222,7 @@ void testMolecularGas() {
 
     System system;
     system.setDefaultPeriodicBoxVectors(Vec3(initialLength, 0, 0), Vec3(0, 0.5*initialLength, 0), Vec3(0, 0, 2*initialLength));
-    MonteCarloBarostat* barostat = new MonteCarloBarostat(pressure, temp, frequency);
+    MonteCarloBarostat* barostat = new MonteCarloBarostat(pressure, temp, frequency, rigid);
     system.addForce(barostat);
     HarmonicBondForce* bonds = new HarmonicBondForce();
     bonds->setUsesPeriodicBoundaryConditions(true);
@@ -177,8 +244,8 @@ void testMolecularGas() {
             system.addParticle(3.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, 10.0);
-        system.addConstraint(positions.size(), positions.size()+2, 0.1);
+        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));
@@ -381,7 +448,10 @@ int main(int argc, char* argv[]) {
         initializeTests(argc, argv);
         testChangingBoxSize();
         testIdealGas();
-        testMolecularGas();
+        testMoleculeScaling(true);
+        testMoleculeScaling(false);
+        testMolecularGas(true);
+        testMolecularGas(false);
         testRandomSeed();
         // Don't run testWater() or testLJPressure() here, because they're very slow on
         // Reference platform.  Individual platforms can run them from runPlatformTests().

wrappers/python/openmm/app/__init__.py

@@ -38,6 +38,8 @@ from .charmmpsffile import CharmmPsfFile, CharmmPSFWarning
 from .simulatedtempering import SimulatedTempering
 from .metadynamics import Metadynamics, BiasVariable
 from .replicaexchangesampler import ReplicaExchangeSampler
+from .replicaexchangereporter import ReplicaExchangeReporter
+from .expandedensemblesampler import ExpandedEnsembleSampler
 
 # Enumerated values
 

wrappers/python/openmm/app/expandedensemblesampler.py

+from __future__ import print_function
+
+"""
+expandedensemblesampler.py: Performs multistate sampling with the expanded ensemble method
+
+This is part of the OpenMM molecular simulation toolkit.
+See https://openmm.org/development.
+
+Portions copyright (c) 2015-2026 Stanford University and the Authors.
+Authors: Peter Eastman
+Contributors:
+
+Permission is hereby granted, free of charge, to any person obtaining a 
+copy of this software and associated documentation files (the "Software"),
+to deal in the Software without restriction, including without limitation
+the rights to use, copy, modify, merge, publish, distribute, sublicense,
+and/or sell copies of the Software, and to permit persons to whom the
+Software is furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
+DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
+OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
+USE OR OTHER DEALINGS IN THE SOFTWARE.
+"""
+__author__ = "Peter Eastman"
+__version__ = "1.0"
+
+import openmm
+import openmm.unit as unit
+from openmm.app.internal import safesave
+from openmm.app.internal.multistatesampler import MultistateSampler
+import math
+import pickle
+import random
+
+class ExpandedEnsembleSampler(object):
+    """
+    ExpandedEnsembleSampler uses the expanded ensemble method to simulate a system in a collection of thermodynamic
+    states.  It supports both the temperature version, in which the states correspond to different temperatures, and the
+    Hamiltonian version, in which they correspond to different potential functions.  It also can combine them to vary
+    temperature and potential function at the same time.
+
+    You provide a Simulation describing the system to simulate and a list of states.  It simulates the system while
+    periodically moving between states in a way that ensures correct sampling.  The method is based on Gibbs sampling,
+    which alternates sampling of conformations x and thermodynamic states s.  It first performs molecular dynamics to
+    sample P(x|s), the distribution of conformations for a fixed state.  It then chooses a new state from the
+    distribution P(s|x) for the current conformation.  For more information on this method, see
+    https://doi.org/10.33011/livecoms.4.1.1583.
+
+    Expanded ensemble sampling is used for a variety of purposes.  Here are some of the more common ones.
+
+    - The temperature version is most often used to accelerate sampling of rare events.  A particular transition might
+      only happen rarely at physiological temperature.  You therefore simulate both the temperature you care about and
+      also higher temperatures at which the transition happens more easily.  Allowing it to move between temperatures
+      produces accurate sampling at low temperature of the entire phase space.
+
+    - The Hamiltonian version may be used to sample an alchemical transition between two endpoints.  It produces
+      efficient sampling at every point along the transition pathway from which a free energy difference can be computed.
+
+    Each thermodynamic state (not to be confused with a State object) is represented by a dict containing property
+    values.  All states must specify values for the same properties.  They describe the ways in which the states differ
+    from each other.  The following properties are supported.
+
+    - 'temperature': the simulation temperature
+    - 'groups': a set containing the force groups to include when computing the energy, for example {0, 2}.  It also may
+      be a weighted sum specified as a dict.  For example, {0:1.0, 2:0.5} means to compute the energy of group 0 plus
+      half the energy of group 2.
+    - Context parameters, specified by the parameter name
+    - Global variables defined by a CustomIntegrator, specified by the variable name
+
+    For example, a typical use of temperature sampling might define the states as
+
+    >>> states = [{'temperature':t*kelvin} for t in np.geomspace(300.0, 500.0, 10)]
+
+    A typical use of Hamiltonian sampling might include forces that depend on a parameter `lambda`.  The states could
+    be defined as
+
+    >>> states = [{'lambda':x} for x in np.linspace(0.0, 1.0, 11)]
+
+    States may also specify multiple properties.  This example includes all combinations of two different parameters.
+
+    >>> states = [{'lambda1':x, 'lambda2':y} for x in np.linspace(0.0, 1.0, 6) for y in np.linspace(0.0, 1.0, 6)]
+
+    To run an expanded ensemble simulation, first create a Simulation object in the normal way.  Then create a
+    ExpandedEnsembleSampler:
+
+    >>> sampler = ExpandedEnsembleSampler(states, simulation, stepsPerIteration)
+
+    The third argument is the number of time steps to integrate between attempted state changes.  Then perform a
+    simulation by calling simulation.step() in the usual way.  The ExpandedEnsembleSampler adds a reporter to the
+    Simulation that manages state changes as it runs.
+
+    The probability distribution P(s|x) is often very nonuniform, which can lead to the simulation spending much more
+    time in some states than others.  This is addressed by including a weight factor for each state when computing the
+    probability.  Ideally the weights should be chosen so it spends equal time in every state.  You can specify the
+    weights yourself with a constructor argument, but usually it is easier to let the ExpandedEnsembleSampler choose
+    them automatically.  This is done with a combination of the Wang-Landau (https://doi.org/10.1103/PhysRevLett.86.2050)
+    and Self-Adjusted Mixture Sampling (http://dx.doi.org/10.1080/10618600.2015.1113975) algorithms.  You should monitor
+    the weights and discard the initial part of the simulation where they are changing rapidly.  Until the weights have
+    converged, the simulation will not sample the correct distribution.
+
+    An ExpandedEnsembleSampler can act as a reporter, generating output to files at regular intervals.  The reporting
+    interval should generally be the same as the one used for writing a trajectory so they will be synchronized with
+    each other.  The following files can optionally be written.
+
+    - A log file that records the current thermodynamic state and weights.  To analyze the results of a simulation, it
+    is essential to know what state it was in at every point in time.
+    - A file recording the current reduced energy (potential energy divided by kT) of every thermodynamic state.  This
+    information is useful for calculating free energies.
+    - A checkpoint file containing all information necessary to resume the simulation.  This includes internal fields
+    of the ExpandedEnsembleSampler itself, as well as a State object recording the current positions, velocities,
+    parameters, etc.
+
+    To resume a simulation from the saved checkpoint, pass `resume=True` to the constructor.  It will load all necessary
+    information and configure the ExpandedEnsembleSampler correctly.
+
+    Attributes
+    ----------
+    states: list[dict]
+        The states to sample
+    simulation: openmm.app.Simulation
+        The Simulation defining the System, Context, and Integrator to use
+    stepsPerIteration: int
+        The number of time steps to integrate between attempted state changes
+    reinitializeVelocities: bool
+        If true, the simulation's velocities are reinitialized from a Boltzmann distribution every time its state
+        changes.  This may sometimes improve stability, but also decreases efficiency.
+    currentIteration: int
+        The number of iterations that have been completed, each consisting of stepsPerIteration time steps followed by
+        an attempted state change
+    currentStateIndex: int
+        The current state of the simulation, specified as an index into states.
+    """
+
+    def __init__(self, states: list[dict], simulation: "openmm.app.Simulation", stepsPerIteration: int,
+                 reinitializeVelocities: bool = False, weights: list[float] | None = None, reportInterval: int = 1000,
+                 logFile: str | object | None = None, energyFile: str | object | None = None,
+                 checkpointFile: str | None = None, resume: bool = False):
+        """Create a new ExpandedEnsembleSampler.
+
+        Parameters
+        ----------
+        states: list[dict]
+            The states to sample
+        simulation: openmm.app.Simulation
+            The Simulation defining the System, Context, and Integrator to use
+        stepsPerIteration: int
+            The number of time steps to integrate between attempted state changes
+        reinitializeVelocities: bool
+            If true, the simulation's velocities are reinitialized from a Boltzmann distribution every time its state
+            changes.  This may sometimes improve stability, but also decreases efficiency.
+        weights: list[float] | None
+            The weights to use for each state.  If None, weights are chosen automatically as the simulation runs.
+        reportInterval: int
+            The frequency at which to write output, measured in time steps.  This must be a multiple of
+            stepsPerIteration.
+        logFile: str | object | None
+            An optional file to write a log to.  This may be either a file-like object or a string containing the path
+            to the file.
+        energyFile: str | object | None
+            An optional file to write reduced energies to.  This may be either a file-like object or a string containing
+            the path to the file.
+        checkpointFile: str | None
+            The path to an optional file for saving checkpointing information
+        resume: bool
+            Specifies whether to resume an earlier simulation.  If True, the checkpoint will be loaded and future output
+            will be appended to the existing files.
+        """
+        self.states = states
+        self.simulation = simulation
+        self.stepsPerIteration = stepsPerIteration
+        self.reinitializeVelocities = reinitializeVelocities
+        self.currentIteration = 0
+        self._stage = 0
+        self.currentStateIndex = 0
+        self._sampler = MultistateSampler(states, simulation.context)
+        if 'temperature' in states[0]:
+            temperature = [s['temperature'] for s in states]
+            for i, t in enumerate(temperature):
+                if not unit.is_quantity(t):
+                    temperature[i] = t*unit.kelvin
+            self._kT = [unit.MOLAR_GAS_CONSTANT_R*t for t in temperature]
+        else:
+            self._kT = None
+            if not hasattr(simulation.integrator, 'getTemperature'):
+                raise ValueError('Cannot determine temperature because the integrator does not have a getTemperature() method. '
+                                 'Specify the temperature in each state dict.')
+        self.reportInterval = reportInterval
+
+        # Initialize the weights.
+
+        if weights is None:
+            self._weights = [0.0]*len(states)
+            self._updateWeights = True
+            self._weightUpdateFactor = 1.0
+            self._histogram = [0]*len(states)
+            self._hasMadeTransition = False
+        else:
+            self._weights = weights
+            self._updateWeights = False
+
+        # Restore from the checkpoint file.
+
+        if resume:
+            if checkpointFile is None:
+                raise ValueError('resume=True, but no checkpoint file was specified to restore from.')
+            with open(checkpointFile, 'rb') as input:
+                checkpoint = pickle.load(input)
+                self._applyCheckpoint(checkpoint)
+
+        # Apply the current state to the simulation.
+
+        self._sampler.applyState(self.currentStateIndex)
+
+        # Open output files.
+
+        if (logFile is not None or energyFile is not None) and reportInterval%stepsPerIteration != 0:
+            raise ValueError('The reporting interval must be a multiple iteration length.')
+        mode = 'a' if resume else 'w'
+        self._openedLogFile = isinstance(logFile, str)
+        if self._openedLogFile:
+            self._log = open(logFile, mode, 1)
+        else:
+            self._log = logFile
+        self._openedEnergyFile = isinstance(energyFile, str)
+        if self._openedEnergyFile:
+            self._energy = open(energyFile, mode, 1)
+        else:
+            self._energy = energyFile
+        self._checkpointFile = checkpointFile
+
+        # Add a reporter to the simulation which will handle the updates and reports.
+
+        class EEReporter(object):
+            def __init__(self, sampler):
+                self.sampler = sampler
+
+            def describeNextReport(self, simulation):
+                sampler = self.sampler
+                steps1 = sampler.stepsPerIteration - simulation.currentStep%sampler.stepsPerIteration
+                steps2 = sampler.reportInterval - simulation.currentStep%sampler.reportInterval
+                steps = min(steps1, steps2)
+                isUpdateAttempt = (steps1 == steps)
+                if isUpdateAttempt and not sampler.reinitializeVelocities:
+                    return {'steps': steps, 'periodic':None, 'include':['velocities']}
+                else:
+                    return {'steps': steps, 'periodic':None, 'include':[]}
+
+            def report(self, simulation, state):
+                sampler = self.sampler
+                if simulation.currentStep%sampler.stepsPerIteration == 0:
+                    reducedEnergy = sampler.attemptStateChange(state)
+                    if simulation.currentStep%sampler.reportInterval == 0:
+                        sampler._writeReport(reducedEnergy)
+
+        simulation.reporters.append(EEReporter(self))
+
+        # Write header lines to the output files.
+
+        if not resume:
+            if self._log is not None:
+                headers = ['"Steps"', '"Iteration"', '"State"'] + [f'"Weight {i}"' for i in range(len(self.states))]
+                print('%s' % (',').join(headers), file=self._log)
+            if self._energy is not None:
+                headers = ['"Steps"'] + [f'"u{i}"' for i in range(len(self.states))]
+                print('%s' % (',').join(headers), file=self._energy)
+
+    def __del__(self):
+        if self._openedLogFile:
+            self._log.close()
+        if self._openedEnergyFile:
+            self._energy.close()
+
+    @property
+    def weights(self):
+        return [x-self._weights[0] for x in self._weights]
+
+    def step(self, steps: int):
+        """Advance the simulation by integrating a specified number of time steps."""
+        self.simulation.step(steps)
+
+    def attemptStateChange(self, state: openmm.State) -> list[float]:
+        """Attempt to move to a different state.  This is called automatically during the simulation, and there is not
+        normally a reason to call it directly.
+
+        You can create subclasses that override this method to select states in different ways."""
+
+        # Compute the probability for each state.  This is done in log space to avoid overflow.
+
+        self.currentIteration += 1
+        energies = self._sampler.computeAllEnergies()
+        if self._kT is None:
+            kT = [unit.MOLAR_GAS_CONSTANT_R*self.simulation.integrator.getTemperature()]*len(self.states)
+        else:
+            kT = self._kT
+        reducedEnergy = [energies[i]/kT[i] for i in range(len(self._weights))]
+        logProbability = [(self._weights[i]-reducedEnergy[i]) for i in range(len(self._weights))]
+        maxLogProb = max(logProbability)
+        offset = maxLogProb + math.log(sum(math.exp(x-maxLogProb) for x in logProbability))
+        probability = [math.exp(x-offset) for x in logProbability]
+
+        # Select a new state.
+
+        prevState = self.currentStateIndex
+        r = random.random()
+        for j in range(len(probability)):
+            if r < probability[j]:
+                if j != self.currentStateIndex:
+                    # Select this state.
+
+                    self._hasMadeTransition = True
+                    self.currentStateIndex = j
+                if self._updateWeights:
+                    # Update the weights.
+
+                    updateFactor = 1/self.currentIteration
+                    if self._stage == 0:
+                        if self._weightUpdateFactor >= updateFactor:
+                            updateFactor = self._weightUpdateFactor
+                        else:
+                            self._stage = 1
+                    for k in range(len(probability)):
+                        self._weights[k] -= updateFactor*probability[k]
+                    if self._stage == 0:
+                        # Early in the simulation we reduce the update factor when the histogram is sufficiently flat.
+                        # Later we will switch over to just using 1/iteration.
+
+                        self._histogram[j] += 1
+                        minCounts = min(self._histogram)
+                        if minCounts > 20 and minCounts >= 0.2*sum(self._histogram)/len(self._histogram):
+                            # Reduce the weight update factor and reset the histogram.
+
+                            self._weightUpdateFactor *= 0.5
+                            self._histogram = [0]*len(self.states)
+                            self._weights = [x-self._weights[0] for x in self._weights]
+                        elif not self._hasMadeTransition and probability[self.currentStateIndex] > 0.99 and sum(self._histogram) >= 10:
+                            # Rapidly increase the weight update factor at the start of the simulation to find
+                            # a reasonable starting value.
+
+                            self._weightUpdateFactor *= 2.0
+                            self._histogram = [0]*len(self.states)
+                break
+            r -= probability[j]
+
+        # Apply the state.
+
+        self._sampler.applyState(self.currentStateIndex)
+
+        # Reinitialize or rescale the velocities.
+
+        if self.reinitializeVelocities:
+            if 'temperature' in self.states[self.currentStateIndex]:
+                temperature = self.states[self.currentStateIndex]['temperature']
+            else:
+                temperature = self.simulation.integrator.getTemperature()
+            self.simulation.context.setVelocitiesToTemperature(temperature)
+        elif kT[prevState] != kT[self.currentStateIndex]:
+            scale = math.sqrt(kT[self.currentStateIndex]/kT[prevState])
+            self.simulation.context.setVelocities(scale*state.getVelocities(asNumpy=True))
+        return reducedEnergy
+
+    def _writeReport(self, reducedEnergy: list[float]):
+        """Write output to the files."""
+        if self._log is not None:
+            print(f'{self.simulation.currentStep},{self.currentIteration},{self.currentStateIndex},' + ','.join('%g' % v for v in self.weights), file=self._log)
+        if self._energy is not None:
+            print(f'{self.simulation.currentStep},' + ','.join('%g' % v for v in reducedEnergy), file=self._energy)
+        if self._checkpointFile is not None:
+            checkpoint = self._createCheckpoint()
+            safesave.save(pickle.dumps(checkpoint), self._checkpointFile)
+
+    def _createCheckpoint(self) -> dict:
+        """Create a dict containing the information to save to a checkpoint file."""
+        checkpoint = {}
+        checkpoint['currentIteration'] = self.currentIteration
+        checkpoint['stage'] = self._stage
+        checkpoint['currentStateIndex'] = self.currentStateIndex
+        checkpoint['weights'] = self._weights
+        checkpoint['weightUpdateFactor'] = self._weightUpdateFactor
+        checkpoint['histogram'] = self._histogram
+        checkpoint['hasMadeTransition'] = self._hasMadeTransition
+        checkpoint['state'] = self.simulation.context.getState(positions=True, velocities=True, parameters=True, integratorParameters=True)
+        return checkpoint
+
+    def _applyCheckpoint(self, checkpoint: dict):
+        """Record the information that was loaded from a checkpoint file."""
+        self.currentIteration = checkpoint['currentIteration']
+        self._stage = checkpoint['stage']
+        self.currentStateIndex = checkpoint['currentStateIndex']
+        self._weights = checkpoint['weights']
+        self._weightUpdateFactor = checkpoint['weightUpdateFactor']
+        self._histogram = checkpoint['histogram']
+        self._hasMadeTransition = checkpoint['hasMadeTransition']
+        self.simulation.context.setState(checkpoint['state'])

wrappers/python/openmm/app/forcefield.py

@@ -1286,6 +1286,12 @@ class ForceField(object):
         -------
         system
             the newly created System
+
+        Notes
+        -----
+        The constructed System contains one particle for each atom in the Topology,
+        and the particles are in the same order as the atoms.  An atom's index in the
+        Topology is therefore also the index of its corresponding particle in the System.
         """
         args['switchDistance'] = switchDistance
         args['flexibleConstraints'] = flexibleConstraints

wrappers/python/openmm/app/internal/multistatesampler.py

@@ -209,7 +209,7 @@ class MultistateSampler(object):
         -------
         an array containing the potential energies of all states in the order they appear in self.states
         """
-        energies = [0*kilojoules_per_mole for _ in self.states]
+        energies = [0 for _ in self.states]*kilojoules_per_mole
         for i, subset in enumerate(self.subsets):
             if self.groups is None:
                 # States don't depend on force groups, so we can just evaluate the energy of one state.
@@ -241,6 +241,6 @@ class MultistateSampler(object):
         an array containing the shifted potential energies of all states in the order they appear in self.states
         """
         if len(self.subsets) == 1 and self.groups is None:
-            return [0*kilojoules_per_mole]*len(self.states)
+            return [0]*len(self.states)*kilojoules_per_mole
         energies = self.computeAllEnergies()
         return [e-energies[0] for e in energies]
\ No newline at end of file

wrappers/python/openmm/app/internal/safesave.py

@@ -4,9 +4,9 @@ safesave.py: Helper module to ensure atomic overwrite/backup of existing files.
 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: Evan Pretti
-Contributors:
+Contributors: Peter Eastman
 
 Permission is hereby granted, free of charge, to any person obtaining a
 copy of this software and associated documentation files (the "Software"),
@@ -51,11 +51,8 @@ def _getTempFilename(prefix):
 
     for index in itertools.count():
         name = f'{prefix}.{index}.tmp'
-        try:
-            with open(name, 'x'):
-                return name
-        except FileExistsError:
-            pass
+        if not os.path.exists(name):
+            return name
 
 def save(data, filename):
     """

wrappers/python/openmm/app/replicaexchangereporter.py

+"""
+replicaexchangereporter.py: Writes output for replica exchange simulations
+
+This is part of the OpenMM molecular simulation toolkit.
+See https://openmm.org/development.
+
+Portions copyright (c) 2026 Stanford University and the Authors.
+Authors: Peter Eastman
+Contributors:
+
+Permission is hereby granted, free of charge, to any person obtaining a
+copy of this software and associated documentation files (the "Software"),
+to deal in the Software without restriction, including without limitation
+the rights to use, copy, modify, merge, publish, distribute, sublicense,
+and/or sell copies of the Software, and to permit persons to whom the
+Software is furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
+DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
+OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
+USE OR OTHER DEALINGS IN THE SOFTWARE.
+"""
+
+import openmm as mm
+import openmm.app as app
+import openmm.unit as unit
+from openmm.app.internal import safesave
+import os
+
+class ReplicaExchangeReporter(object):
+    """
+    This is a reporter for use with ReplicaExchangeSampler.  It outputs the most important information required for
+    nearly all replica exchange simulations.  It must be added to the sampler's list of reporters:
+
+    >>> sampler.reporters.append(ReplicaExchangeReporter(directory, interval, sampler))
+
+    As the simulation runs, it creates the following files in a user specified directory.
+
+    - A CSV file containing a log of what state each replica was in at each iteration.
+    - A set of XML files containing serialized State objects for each replica. These serve as checkpoints, allowing a
+      simulation to be resumed later.
+    - (optional) A trajectory file for each thermodynamic state.  The coordinates in these files jump discontinuously
+      whenever an exchange happens.
+    - (optional) A trajectory file for each replica.  The coordinates in these files are continuous, but the
+      thermodynamic states they explore change with time.
+    - (optional) A CSV file containing the reduced energy (potential energy divided by kT) of every replica in every
+      state at each iteration.  This information is useful for calculating free energies.
+
+    The reduced energies are written as a matrix in flattened order for each iteration.  You can load the file with
+    NumPy using the command
+
+    >>> energy = np.loadtxt('energy.csv', delimiter=',').reshape(-1, n_states, n_states)
+
+    `energy[i][j][k]` is the reduced energy of state k for replica j in iteration i.
+
+    To resume a simulation from the saved files, pass `resume=True` to the constructor.  It will load all necessary
+    information and configure the ReplicaExchangeSampler correctly.
+    """
+    def __init__(self, directory: str, reportInterval: int, sampler: "app.ReplicaExchangeSampler",
+                 trajectoryPerState: bool = True, trajectoryPerReplica: bool = False, trajectoryFormat: str = 'xtc',
+                 enforcePeriodicBox: bool | None = None, energy: bool = False, resume: bool = False):
+        """
+        Create a ReplicaExchangeReporter.
+
+        Parameters
+        ----------
+        directory: str
+            The path to the directory where the files will be written.  The directory will be created if it does not
+            already exist.  The directory must be empty unless resume is True, in which case it must contain the files
+            from an earlier simulation.
+        reportInterval: int
+            The frequency at which to write output, measured in iterations
+        sampler: openmm.app.ReplicaExchangeSampler
+            The sampler this reporter will be added to
+        trajectoryPerState: bool
+            If True, a trajectory file will be saved for each thermodynamic state.
+        trajectoryPerReplica: bool
+            If True, a trajectory file will be saved for each replica.
+        trajectoryFormat: str
+            The format in which to save trajectories.  Supported options are 'xtc' and 'dcd'.
+        enforcePeriodicBox: bool
+            Specifies whether particle positions should be translated so the center of every molecule
+            lies in the same periodic box.  If None (the default), it will automatically decide whether
+            to translate molecules based on whether the system being simulated uses periodic boundary
+            conditions.
+        energy: bool
+            If True, a CSV file will be saved containing reduced energies.
+        resume: bool
+            Specifies whether to resume an earlier simulation.  If True, the checkpoint and log data will be loaded into
+            the ReplicaExchangeSampler, and future output will be appended to the existing files.
+        """
+        trajectoryFormat = trajectoryFormat.lower()
+        self.directory = directory
+        self.reportInterval = reportInterval
+        self.trajectoryPerState = trajectoryPerState
+        self.trajectoryPerReplica = trajectoryPerReplica
+        self.format = format
+        self._log = None
+        self._energy = None
+        numStates = len(sampler.states)
+
+        # Validate the inputs and create the output directory if necessary.
+
+        if trajectoryFormat not in ('xtc', 'dcd'):
+            raise ValueError(f'Unsupported trajectory format: {trajectoryFormat}.  Allowed values are "xtc" and "dcd".')
+        if resume:
+            if not os.path.isdir(directory):
+                raise ValueError(f'Cannot resume because the directory does not exist: {directory}')
+
+            def checkExists(filename):
+                if not os.path.isfile(os.path.join(directory, filename)):
+                    raise ValueError(f'Cannot resume because the file {filename} does not exist.')
+
+            checkExists('log.csv')
+            if energy:
+                checkExists('energy.csv')
+            for i in range(numStates):
+                checkExists(f'checkpoint_{i}.xml')
+                if trajectoryPerState:
+                    checkExists(f'state_{i}.{trajectoryFormat}')
+                if trajectoryPerReplica:
+                    checkExists(f'replica_{i}.{trajectoryFormat}')
+        else:
+            if os.path.isdir(directory):
+                if len(os.listdir(directory)) > 0:
+                    raise ValueError(f'The directory {directory} is not empty.')
+            else:
+                os.makedirs(directory)
+
+        # Load state assignments and checkpoints if resuming another simulation.
+
+        if resume:
+            for i in range(numStates):
+                with open(os.path.join(directory, f'checkpoint_{i}.xml')) as input:
+                    sampler.replicaConformation[i] = mm.XmlSerializer.deserialize(input.read())
+            log = None
+            with open(os.path.join(directory, 'log.csv')) as input:
+                for line in input:
+                    log = line
+            if log is None:
+                raise ValueError('Cannot resume because log file is empty.')
+            fields = log.split(',')
+            sampler.replicaStateIndex = [int(x) for x in fields[2:(numStates+2)]]
+            sampler._previousReplicaStateIndex = sampler.replicaStateIndex[:]
+            sampler.currentIteration = int(fields[0])
+
+        # Create reporters and open files for writing output.
+
+        def createReporter(filename):
+            path = os.path.join(directory, filename)
+            interval = reportInterval*sampler.stepsPerIteration
+            if trajectoryFormat == 'xtc':
+                return app.XTCReporter(path, reportInterval*sampler.stepsPerIteration, append=resume, enforcePeriodicBox=enforcePeriodicBox)
+            return app.DCDReporter(path, interval, append=resume, enforcePeriodicBox=enforcePeriodicBox)
+
+        self._stateReporters = []
+        self._replicaReporters = []
+        for i in range(numStates):
+            if trajectoryPerState:
+                self._stateReporters.append(createReporter(f'state_{i}.{trajectoryFormat}'))
+            if trajectoryPerReplica:
+                self._replicaReporters.append(createReporter(f'replica_{i}.{trajectoryFormat}'))
+        self._log = open(os.path.join(directory, 'log.csv'), 'a' if resume else 'w')
+        if not resume:
+            print(','.join(['Iteration', 'Step']+[f'Replica_{i}_State' for i in range(numStates)]), file=self._log)
+            self._log.flush()
+        if energy:
+            self._energy = open(os.path.join(directory, 'energy.csv'), 'a' if resume else 'w')
+
+    def __del__(self):
+        if self._log is not None:
+            self._log.close()
+        if self._energy is not None:
+            self._energy.close()
+
+    def __call__(self, sampler: "app.ReplicaExchangeSampler"):
+        """This is invoked by the ReplicaExchangeSampler at the end of every iteration to generate output."""
+        if sampler.currentIteration % self.reportInterval != 0:
+            return
+        logData = [sampler.currentIteration, sampler.simulation.currentStep]+sampler.replicaStateIndex
+        print(','.join([str(x) for x in logData]), file=self._log)
+        self._log.flush()
+        if self._energy is not None:
+            import numpy as np
+            energy = np.array(sampler.replicaStateEnergy)
+            if sampler._kT is None:
+                energy /= unit.MOLAR_GAS_CONSTANT_R*sampler.simulation.integrator.getTemperature()
+            else:
+                energy /= sampler._kT
+            print(','.join([str(x) for x in energy.flatten()]), file=self._energy)
+            self._energy.flush()
+        for i, conf in enumerate(sampler.replicaConformation):
+            if len(self._stateReporters) > 0:
+                self._stateReporters[sampler.replicaStateIndex[i]].report(sampler.simulation, conf)
+            if len(self._replicaReporters) > 0:
+                self._replicaReporters[i].report(sampler.simulation, conf)
+            safesave.save(mm.XmlSerializer.serialize(conf), os.path.join(self.directory, f'checkpoint_{i}.xml'))

wrappers/python/openmm/app/replicaexchangesampler.py

@@ -35,7 +35,7 @@ import random
 
 class ReplicaExchangeSampler(object):
     """
-    ReplicaExchangeSampler uses replica exchange to simulate a system in a collection of different states.  It supports
+    ReplicaExchangeSampler uses replica exchange to simulate a system in a collection of thermodynamic states.  It supports
     both temperature replica exchange, in which the states correspond to different temperatures, and Hamiltonian replica
     exchange, in which they correspond to different potential functions.  It also can combine them to vary temperature
     and potential function at the same time.
@@ -54,9 +54,9 @@ class ReplicaExchangeSampler(object):
     - Hamiltonian replica exchange may be used to sample an alchemical transition between two endpoints.  It produces
       efficient sampling at every point along the transition pathway from which a free energy difference can be computed.
 
-    Each state (not to be confused with a State object) is represented by a dict containing property values.  All states
-    must specify values for the same properties.  They describe the ways in which the states differ from each other.  The
-    following properties are supported.
+    Each thermodynamic state (not to be confused with a State object) is represented by a dict containing property
+    values.  All states must specify values for the same properties.  They describe the ways in which the states differ
+    from each other.  The following properties are supported.
 
     - 'temperature': the simulation temperature
     - 'groups': a set containing the force groups to include when computing the energy, for example {0, 2}.  It also may
@@ -90,7 +90,16 @@ class ReplicaExchangeSampler(object):
 
     Because replica exchange involves simulating many distinct replicas, the standard reporters used to generate output
     from a simulation are often not appropriate.  ReplicaExchangeSampler instead provides its own reporting mechanism.
-    Define a function that takes a ReplicaExchangeSampler as its only argument:
+    For most purposes, you can use the ReplicaExchangeReporter class.  Create a reporter and add it to the sampler:
+
+    >>> sampler.reporters.append(ReplicaExchangeReporter(directory, interval, sampler))
+
+    It can write out a variety of types of information including a log of state assignments, trajectories for each
+    replica and/or state, reduced energies, and checkpoint files for resuming a simulation.  See the documentation on
+    ReplicaExchangeReporter for details.
+
+    If you want to output other information, you can create your own reporters.  Define a function that takes a
+    ReplicaExchangeSampler as its only argument:
 
     >>> def report(sampler):
     >>>     # Generate whatever output you want here

wrappers/python/openmm/app/simulatedtempering.py

@@ -6,7 +6,7 @@ simulatedtempering.py: Implements simulated tempering
 This is part of the OpenMM molecular simulation toolkit.
 See https://openmm.org/development.
 
-Portions copyright (c) 2015 Stanford University and the Authors.
+Portions copyright (c) 2015-2026 Stanford University and the Authors.
 Authors: Peter Eastman
 Contributors:
 
@@ -32,28 +32,14 @@ __author__ = "Peter Eastman"
 __version__ = "1.0"
 
 import openmm.unit as unit
-import math
-import random
 from sys import stdout
 
-try:
-    import bz2
-    have_bz2 = True
-except: have_bz2 = False
-
-try:
-    import gzip
-    have_gzip = True
-except: have_gzip = False
-
-try:
-    import numpy
-    have_numpy = True
-except: have_numpy = False
 
 class SimulatedTempering(object):
     """SimulatedTempering implements the simulated tempering algorithm for accelerated sampling.
-    
+
+    @deprecated This class exists mainly for historical reasons.  It now is just a thin wrapper around ExpandedEnsembleSampler.  It is still supported for backward compatibility, but using ExpandedEnsembleSampler directly is recommended since it is much more flexible.
+
     It runs a simulation while allowing the temperature to vary.  At high temperatures, it can more easily cross
     energy barriers to explore a wider area of conformation space.  At low temperatures, it can thoroughly
     explore each local region.  For details, see Marinari, E. and Parisi, G., Europhys. Lett. 19(6). pp. 451-458 (1992).
@@ -106,160 +92,27 @@ class SimulatedTempering(object):
             The interval (in time steps) at which to write information to the report file
         reportFile: string or file
             The file to write reporting information to, specified as a file name or file object
-        """        
-        self.simulation = simulation
+        """
         if temperatures is None:
             if unit.is_quantity(minTemperature):
                 minTemperature = minTemperature.value_in_unit(unit.kelvin)
             if unit.is_quantity(maxTemperature):
                 maxTemperature = maxTemperature.value_in_unit(unit.kelvin)
-            self.temperatures = [minTemperature*((float(maxTemperature)/minTemperature)**(i/float(numTemperatures-1))) for i in range(numTemperatures)]*unit.kelvin
-        else:
-            numTemperatures = len(temperatures)
-            self.temperatures = [(t.value_in_unit(unit.kelvin) if unit.is_quantity(t) else t)*unit.kelvin for t in temperatures]
-            if any(self.temperatures[i] >= self.temperatures[i+1] for i in range(numTemperatures-1)):
-                raise ValueError('The temperatures must be in strictly increasing order')
-        self.tempChangeInterval = tempChangeInterval
-        self.reportInterval = reportInterval
-        self.inverseTemperatures = [1.0/(unit.MOLAR_GAS_CONSTANT_R*t) for t in self.temperatures]
-
-        # If necessary, open the file we will write reports to.
-
-        self._openedFile = isinstance(reportFile, str)
-        if self._openedFile:
-            # Detect the desired compression scheme from the filename extension
-            # and open all files unbuffered
-            if reportFile.endswith('.gz'):
-                if not have_gzip:
-                    raise RuntimeError("Cannot write .gz file because Python could not import gzip library")
-                self._out = gzip.GzipFile(fileobj=open(reportFile, 'wb', 0))
-            elif reportFile.endswith('.bz2'):
-                if not have_bz2:
-                    raise RuntimeError("Cannot write .bz2 file because Python could not import bz2 library")
-                self._out = bz2.BZ2File(reportFile, 'w', 0)
-            else:
-                self._out = open(reportFile, 'w', 1)
-        else:
-            self._out = reportFile
-        
-        # Initialize the weights.
-        
-        if weights is None:
-            self._weights = [0.0]*numTemperatures
-            self._updateWeights = True
-            self._weightUpdateFactor = 1.0
-            self._histogram = [0]*numTemperatures
-            self._hasMadeTransition = False
-        else:
-            self._weights = weights
-            self._updateWeights = False
-
-        # Select the initial temperature.
-        
-        self.currentTemperature = 0
-        self.simulation.integrator.setTemperature(self.temperatures[self.currentTemperature])
-        for param in self.simulation.context.getParameters():
-            if 'MonteCarloTemperature' in param:
-                self.simulation.context.setParameter(param, self.temperatures[self.currentTemperature])
-        
-        # Add a reporter to the simulation which will handle the updates and reports.
-        
-        class STReporter(object):
-            def __init__(self, st):
-                self.st = st
+            temperatures = [minTemperature*((float(maxTemperature)/minTemperature)**(i/float(numTemperatures-1))) for i in range(numTemperatures)]*unit.kelvin
+        self.temperatures = [(t.value_in_unit(unit.kelvin) if unit.is_quantity(t) else t)*unit.kelvin for t in temperatures]
+        states = [{'temperature':t} for t in temperatures]
+        from . import ExpandedEnsembleSampler
+        self._sampler = ExpandedEnsembleSampler(states, simulation, tempChangeInterval, weights=weights,
+                                                reportInterval=reportInterval, logFile=reportFile)
 
-            def describeNextReport(self, simulation):
-                st = self.st
-                steps1 = st.tempChangeInterval - simulation.currentStep%st.tempChangeInterval
-                steps2 = st.reportInterval - simulation.currentStep%st.reportInterval
-                steps = min(steps1, steps2)
-                isUpdateAttempt = (steps1 == steps)
-                if isUpdateAttempt:
-                    return {'steps': steps, 'periodic':None, 'include':['velocities', 'energy']}
-                else:
-                    return {'steps': steps, 'periodic':None, 'include':[]}
-
-            def report(self, simulation, state):
-                st = self.st
-                if simulation.currentStep%st.tempChangeInterval == 0:
-                    st._attemptTemperatureChange(state)
-                if simulation.currentStep%st.reportInterval == 0:
-                    st._writeReport()
-        
-        simulation.reporters.append(STReporter(self))
-        
-        # Write out the header line.
-        
-        headers = ['Steps', 'Temperature (K)']
-        for t in self.temperatures:
-            headers.append('%gK Weight' % t.value_in_unit(unit.kelvin))
-        print('#"%s"' % ('"\t"').join(headers), file=self._out)
-
-    def __del__(self):
-        if self._openedFile:
-            self._out.close()
-    
     @property
     def weights(self):
-        return [x-self._weights[0] for x in self._weights]
+        return self._sampler.weights
+
+    @property
+    def currentTemperature(self):
+        return self._sampler.currentStateIndex
 
     def step(self, steps):
         """Advance the simulation by integrating a specified number of time steps."""
-        self.simulation.step(steps)
-    
-    def _attemptTemperatureChange(self, state):
-        """Attempt to move to a different temperature."""
-        
-        # Compute the probability for each temperature.  This is done in log space to avoid overflow.
-
-        logProbability = [(self._weights[i]-self.inverseTemperatures[i]*state.getPotentialEnergy()) for i in range(len(self._weights))]
-        maxLogProb = max(logProbability)
-        offset = maxLogProb + math.log(sum(math.exp(x-maxLogProb) for x in logProbability))
-        probability = [math.exp(x-offset) for x in logProbability]
-        r = random.random()
-        for j in range(len(probability)):
-            if r < probability[j]:
-                if j != self.currentTemperature:
-                    # Rescale the velocities.
-                    
-                    scale = math.sqrt(self.temperatures[j]/self.temperatures[self.currentTemperature])
-                    if have_numpy:
-                        velocities = scale*state.getVelocities(asNumpy=True).value_in_unit(unit.nanometers/unit.picoseconds)
-                    else:
-                        velocities = [v*scale for v in state.getVelocities().value_in_unit(unit.nanometers/unit.picoseconds)]
-                    self.simulation.context.setVelocities(velocities)
-
-                    # Select this temperature.
-                    
-                    self._hasMadeTransition = True
-                    self.currentTemperature = j
-                    self.simulation.integrator.setTemperature(self.temperatures[j])
-                    for param in self.simulation.context.getParameters():
-                        if 'MonteCarloTemperature' in param:
-                            self.simulation.context.setParameter(param, self.temperatures[self.currentTemperature])
-                if self._updateWeights:
-                    # Update the weight factors.
-                    
-                    self._weights[j] -= self._weightUpdateFactor
-                    self._histogram[j] += 1
-                    minCounts = min(self._histogram)
-                    if minCounts > 20 and minCounts >= 0.2*sum(self._histogram)/len(self._histogram):
-                        # Reduce the weight update factor and reset the histogram.
-                        
-                        self._weightUpdateFactor *= 0.5
-                        self._histogram = [0]*len(self.temperatures)
-                        self._weights = [x-self._weights[0] for x in self._weights]
-                    elif not self._hasMadeTransition and probability[self.currentTemperature] > 0.99:
-                        # Rapidly increase the weight update factor at the start of the simulation to find
-                        # a reasonable starting value.
-                        
-                        self._weightUpdateFactor *= 2.0
-                        self._histogram = [0]*len(self.temperatures)
-                return
-            r -= probability[j]
-
-    def _writeReport(self):
-        """Write out a line to the report."""
-        temperature = self.temperatures[self.currentTemperature].value_in_unit(unit.kelvin)
-        values = [temperature]+self.weights
-        print(('%d\t' % self.simulation.currentStep) + '\t'.join('%g' % v for v in values), file=self._out)
+        self._sampler.step(steps)

wrappers/python/openmm/app/topology.py

@@ -74,6 +74,10 @@ class Topology(object):
     pairs are bonded to each other, and the dimensions of the crystallographic unit cell.
 
     Atom and residue names should follow the PDB 3.0 nomenclature for all molecules for which one exists.
+
+    When a System is constructed from a Topology (for example with ForceField.createSystem()),
+    each atom is converted to a corresponding particle in the System.  The particles are in the
+    same order as the atoms, so an atom's index in the Topology is also the index of its particle.
     """
 
     _standardBonds = {}