More barostats support scaling particles independently (#5278)

* More barostats support scaling particles independently

* Documentation update

* MonteCarloMembraneBarostat can scale particles independently

docs-source/usersguide/theory/02_standard_forces.rst

@@ -772,7 +772,11 @@ accordingly.
 Each Monte Carlo step modifies particle positions by scaling the centroid of
 each molecule, then applying the resulting displacement to each particle in the
 molecule.  This ensures that each molecule is translated as a unit, so bond
-lengths and constrained distances are unaffected.
+lengths and constrained distances are unaffected.  Alternatively, there is an
+option to scale the position of each particle independently.  This option is
+intended for unusual situations, such as when the entire system is a single
+molecule.  In most other cases it is less efficient than scaling molecules,
+and it should never be used with a system that contains constraints.
 
 MonteCarloBarostat assumes the simulation is being run at constant temperature
 as well as pressure, and the simulation temperature affects the step acceptance

openmmapi/include/openmm/MonteCarloAnisotropicBarostat.h

@@ -7,7 +7,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, Lee-Ping Wang                                      *
  * Contributors:                                                              *
  *                                                                            *
@@ -96,8 +96,10 @@ public:
      * @param scaleY              whether to allow the Y dimension of the periodic box to change size
      * @param scaleZ              whether to allow the Z dimension of the periodic box to change size
      * @param frequency           the frequency at which Monte Carlo pressure changes should be attempted (in time steps)
+     * @param scaleMoleculesAsRigid   if true, coordinate scaling keeps molecules rigid, scaling only the center of mass
+     *                                of each one.  If false, every atom is scaled independently.
      */
-    MonteCarloAnisotropicBarostat(const Vec3& defaultPressure, double defaultTemperature, bool scaleX = true, bool scaleY = true, bool scaleZ = true, int frequency = 25);
+    MonteCarloAnisotropicBarostat(const Vec3& defaultPressure, double defaultTemperature, bool scaleX = true, bool scaleY = true, bool scaleZ = true, int frequency = 25, bool scaleMoleculesAsRigid = true);
     /**
      * Get the default pressure (in bar).
      *
@@ -185,11 +187,28 @@ public:
     bool usesPeriodicBoundaryConditions() const {
         return false;
     }
+    /**
+     * Get whether scaling is applied to the centroid of each molecule while keeping
+     * the molecules rigid, or to each atom independently.
+     *
+     * @returns true if scaling is applied to molecule centroids, false if it is applied to each atom independently.
+     */
+    bool getScaleMoleculesAsRigid() const {
+        return scaleMoleculesAsRigid;
+    }
+    /**
+     * Set whether scaling is applied to the centroid of each molecule while keeping
+     * the molecules rigid, or to each atom independently.
+     */
+    void setScaleMoleculesAsRigid(bool rigid) {
+        scaleMoleculesAsRigid = rigid;
+    }
     /**
      * Compute the instantaneous pressure along each axis of a system to which this barostat
      * is applied.
      * 
-     * The pressure is computed from the molecular virial, using a finite difference to
+     * The pressure is computed from the molecular virial if getScaleMoleculesAsRigid()
+     * is true, or the atomic virial if it is false.  It uses a finite difference to
      * calculate the derivative of potential energy with respect to volume.  For most systems
      * in equilibrium, the time average of the instantaneous pressure should equal the
      * pressure applied by the barostat.  Fluctuations around the average value can be
@@ -203,6 +222,7 @@ public:
 protected:
     ForceImpl* createImpl() const;
 private:
+    bool scaleMoleculesAsRigid;
     Vec3 defaultPressure;
     double defaultTemperature;
     bool scaleX, scaleY, scaleZ;

openmmapi/include/openmm/MonteCarloBarostat.h

@@ -7,7 +7,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                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -71,8 +71,10 @@ public:
      * @param defaultPressure     the default pressure acting on the system (in bar)
      * @param defaultTemperature  the default temperature at which the system is being maintained (in Kelvin)
      * @param frequency           the frequency at which Monte Carlo pressure changes should be attempted (in time steps)
+     * @param scaleMoleculesAsRigid   if true, coordinate scaling keeps molecules rigid, scaling only the center of mass
+     *                                of each one.  If false, every atom is scaled independently.
      */
-    MonteCarloBarostat(double defaultPressure, double defaultTemperature, int frequency = 25);
+    MonteCarloBarostat(double defaultPressure, double defaultTemperature, int frequency = 25, bool scaleMoleculesAsRigid = true);
     /**
      * Get the default pressure acting on the system (in bar).
      *
@@ -142,10 +144,27 @@ public:
     bool usesPeriodicBoundaryConditions() const {
         return false;
     }
+    /**
+     * Get whether scaling is applied to the centroid of each molecule while keeping
+     * the molecules rigid, or to each atom independently.
+     *
+     * @returns true if scaling is applied to molecule centroids, false if it is applied to each atom independently.
+     */
+    bool getScaleMoleculesAsRigid() const {
+        return scaleMoleculesAsRigid;
+    }
+    /**
+     * Set whether scaling is applied to the centroid of each molecule while keeping
+     * the molecules rigid, or to each atom independently.
+     */
+    void setScaleMoleculesAsRigid(bool rigid) {
+        scaleMoleculesAsRigid = rigid;
+    }
     /**
      * Compute the instantaneous pressure of a system to which this barostat is applied.
      * 
-     * The pressure is computed from the molecular virial, using a finite difference to
+     * The pressure is computed from the molecular virial if getScaleMoleculesAsRigid()
+     * is true, or the atomic virial if it is false.  It uses a finite difference to
      * calculate the derivative of potential energy with respect to volume.  For most systems
      * in equilibrium, the time average of the instantaneous pressure should equal the
      * pressure applied by the barostat.  Fluctuations around the average value can be
@@ -159,6 +178,7 @@ public:
 protected:
     ForceImpl* createImpl() const;
 private:
+    bool scaleMoleculesAsRigid;
     double defaultPressure, defaultTemperature;
     int frequency, randomNumberSeed;
 };

openmmapi/include/openmm/MonteCarloMembraneBarostat.h

@@ -7,7 +7,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                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -127,8 +127,11 @@ public:
      * @param xymode                 the mode specifying the behavior of the X and Y axes
      * @param zmode                  the mode specifying the behavior of the Z axis
      * @param frequency              the frequency at which Monte Carlo volume changes should be attempted (in time steps)
+     * @param scaleMoleculesAsRigid  if true, coordinate scaling keeps molecules rigid, scaling only the center of mass
+     *                               of each one.  If false, every atom is scaled independently.
      */
-    MonteCarloMembraneBarostat(double defaultPressure, double defaultSurfaceTension, double defaultTemperature, XYMode xymode, ZMode zmode, int frequency = 25);
+    MonteCarloMembraneBarostat(double defaultPressure, double defaultSurfaceTension, double defaultTemperature, XYMode xymode, ZMode zmode,
+                               int frequency = 25, bool scaleMoleculesAsRigid = true);
     /**
      * Get the default pressure acting on the system (in bar).
      *
@@ -237,11 +240,28 @@ public:
     bool usesPeriodicBoundaryConditions() const {
         return false;
     }
+    /**
+     * Get whether scaling is applied to the centroid of each molecule while keeping
+     * the molecules rigid, or to each atom independently.
+     *
+     * @returns true if scaling is applied to molecule centroids, false if it is applied to each atom independently.
+     */
+    bool getScaleMoleculesAsRigid() const {
+        return scaleMoleculesAsRigid;
+    }
+    /**
+     * Set whether scaling is applied to the centroid of each molecule while keeping
+     * the molecules rigid, or to each atom independently.
+     */
+    void setScaleMoleculesAsRigid(bool rigid) {
+        scaleMoleculesAsRigid = rigid;
+    }
     /**
      * Compute the instantaneous pressure along each axis of a system to which this barostat
      * is applied.
      * 
-     * The pressure is computed from the molecular virial, using a finite difference to
+     * The pressure is computed from the molecular virial if getScaleMoleculesAsRigid()
+     * is true, or the atomic virial if it is false.  It uses a finite difference to
      * calculate the derivative of potential energy with respect to volume.  For most systems
      * in equilibrium, the time average of the instantaneous pressure should equal the
      * pressure applied by the barostat.  Fluctuations around the average value can be
@@ -259,6 +279,7 @@ private:
     XYMode xymode;
     ZMode zmode;
     int frequency, randomNumberSeed;
+    bool scaleMoleculesAsRigid;
 };
 
 } // namespace OpenMM

openmmapi/src/MonteCarloAnisotropicBarostat.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                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -32,8 +32,8 @@
 
 using namespace OpenMM;
 
-MonteCarloAnisotropicBarostat::MonteCarloAnisotropicBarostat(const Vec3& defaultPressure, double defaultTemperature, bool scaleX, bool scaleY, bool scaleZ, int frequency) :
-        scaleX(scaleX), scaleY(scaleY), scaleZ(scaleZ) {
+MonteCarloAnisotropicBarostat::MonteCarloAnisotropicBarostat(const Vec3& defaultPressure, double defaultTemperature, bool scaleX, bool scaleY, bool scaleZ, int frequency, bool scaleMoleculesAsRigid) :
+        scaleX(scaleX), scaleY(scaleY), scaleZ(scaleZ), scaleMoleculesAsRigid(scaleMoleculesAsRigid) {
     setDefaultPressure(defaultPressure);
     setDefaultTemperature(defaultTemperature);
     setFrequency(frequency);

openmmapi/src/MonteCarloAnisotropicBarostatImpl.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, Lee-Ping Wang                                      *
  * Contributors:                                                              *
  *                                                                            *
@@ -48,7 +48,7 @@ void MonteCarloAnisotropicBarostatImpl::initialize(ContextImpl& context) {
     if (!context.getSystem().usesPeriodicBoundaryConditions())
         throw OpenMMException("A barostat cannot be used with a non-periodic system");
     kernel = context.getPlatform().createKernel(ApplyMonteCarloBarostatKernel::Name(), context);
-    kernel.getAs<ApplyMonteCarloBarostatKernel>().initialize(context.getSystem(), owner, 3);
+    kernel.getAs<ApplyMonteCarloBarostatKernel>().initialize(context.getSystem(), owner, 3, owner.getScaleMoleculesAsRigid());
     Vec3 box[3];
     context.getPeriodicBoxVectors(box[0], box[1], box[2]);
     double volume = box[0][0]*box[1][1]*box[2][2];
@@ -113,9 +113,14 @@ void MonteCarloAnisotropicBarostatImpl::updateContextState(ContextImpl& context,
 
     // Compute the energy of the modified system.
 
+    double numberOfScaledParticles;
+    if (owner.getScaleMoleculesAsRigid())
+        numberOfScaledParticles = context.getMolecules().size();
+    else
+        numberOfScaledParticles = context.getSystem().getNumParticles();
     double finalEnergy = context.getOwner().getState(State::Energy, false, groups).getPotentialEnergy();
     double kT = BOLTZ*context.getParameter(MonteCarloAnisotropicBarostat::Temperature());
-    double w = finalEnergy-initialEnergy + pressure*deltaVolume - context.getMolecules().size()*kT*log(newVolume/volume);
+    double w = finalEnergy-initialEnergy + pressure*deltaVolume - numberOfScaledParticles*kT*log(newVolume/volume);
     if (w > 0 && SimTKOpenMMUtilities::getUniformlyDistributedRandomNumber() > exp(-w/kT)) {
         // Reject the step.
         

openmmapi/src/MonteCarloBarostat.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                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -34,7 +34,8 @@
 
 using namespace OpenMM;
 
-MonteCarloBarostat::MonteCarloBarostat(double defaultPressure, double defaultTemperature, int frequency) {
+MonteCarloBarostat::MonteCarloBarostat(double defaultPressure, double defaultTemperature, int frequency, bool scaleMoleculesAsRigid) :
+            scaleMoleculesAsRigid(scaleMoleculesAsRigid) {
     setDefaultPressure(defaultPressure);
     setDefaultTemperature(defaultTemperature);
     setFrequency(frequency);

openmmapi/src/MonteCarloBarostatImpl.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                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -48,7 +48,7 @@ void MonteCarloBarostatImpl::initialize(ContextImpl& context) {
     if (!context.getSystem().usesPeriodicBoundaryConditions())
         throw OpenMMException("A barostat cannot be used with a non-periodic system");
     kernel = context.getPlatform().createKernel(ApplyMonteCarloBarostatKernel::Name(), context);
-    kernel.getAs<ApplyMonteCarloBarostatKernel>().initialize(context.getSystem(), owner, 1);
+    kernel.getAs<ApplyMonteCarloBarostatKernel>().initialize(context.getSystem(), owner, 1, owner.getScaleMoleculesAsRigid());
     Vec3 box[3];
     context.getPeriodicBoxVectors(box[0], box[1], box[2]);
     double volume = box[0][0]*box[1][1]*box[2][2];
@@ -82,10 +82,15 @@ void MonteCarloBarostatImpl::updateContextState(ContextImpl& context, bool& forc
 
     // Compute the energy of the modified system.
 
+    double numberOfScaledParticles;
+    if (owner.getScaleMoleculesAsRigid())
+        numberOfScaledParticles = context.getMolecules().size();
+    else
+        numberOfScaledParticles = context.getSystem().getNumParticles();
     double finalEnergy = context.getOwner().getState(State::Energy, false, groups).getPotentialEnergy();
     double pressure = context.getParameter(MonteCarloBarostat::Pressure())*(AVOGADRO*1e-25);
     double kT = BOLTZ*context.getParameter(MonteCarloBarostat::Temperature());
-    double w = finalEnergy-initialEnergy + pressure*deltaVolume - context.getMolecules().size()*kT*log(newVolume/volume);
+    double w = finalEnergy-initialEnergy + pressure*deltaVolume - numberOfScaledParticles*kT*log(newVolume/volume);
     if (w > 0 && SimTKOpenMMUtilities::getUniformlyDistributedRandomNumber() > exp(-w/kT)) {
         // Reject the step.
 

openmmapi/src/MonteCarloMembraneBarostat.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                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -32,8 +32,8 @@
 
 using namespace OpenMM;
 
-MonteCarloMembraneBarostat::MonteCarloMembraneBarostat(double defaultPressure, double defaultSurfaceTension, double defaultTemperature, XYMode xymode, ZMode zmode, int frequency) :
-        xymode(xymode), zmode(zmode) {
+MonteCarloMembraneBarostat::MonteCarloMembraneBarostat(double defaultPressure, double defaultSurfaceTension, double defaultTemperature, XYMode xymode, ZMode zmode,
+        int frequency, bool scaleMoleculesAsRigid) : xymode(xymode), zmode(zmode), scaleMoleculesAsRigid(scaleMoleculesAsRigid) {
     setDefaultPressure(defaultPressure);
     setDefaultSurfaceTension(defaultSurfaceTension);
     setDefaultTemperature(defaultTemperature);

openmmapi/src/MonteCarloMembraneBarostatImpl.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, Lee-Ping Wang                                      *
  * Contributors:                                                              *
  *                                                                            *
@@ -48,7 +48,7 @@ void MonteCarloMembraneBarostatImpl::initialize(ContextImpl& context) {
     if (!context.getSystem().usesPeriodicBoundaryConditions())
         throw OpenMMException("A barostat cannot be used with a non-periodic system");
     kernel = context.getPlatform().createKernel(ApplyMonteCarloBarostatKernel::Name(), context);
-    kernel.getAs<ApplyMonteCarloBarostatKernel>().initialize(context.getSystem(), owner, 3);
+    kernel.getAs<ApplyMonteCarloBarostatKernel>().initialize(context.getSystem(), owner, 3, owner.getScaleMoleculesAsRigid());
     Vec3 box[3];
     context.getPeriodicBoxVectors(box[0], box[1], box[2]);
     double volume = box[0][0]*box[1][1]*box[2][2];
@@ -114,9 +114,14 @@ void MonteCarloMembraneBarostatImpl::updateContextState(ContextImpl& context, bo
 
     // Compute the energy of the modified system.
 
+    double numberOfScaledParticles;
+    if (owner.getScaleMoleculesAsRigid())
+        numberOfScaledParticles = context.getMolecules().size();
+    else
+        numberOfScaledParticles = context.getSystem().getNumParticles();
     double finalEnergy = context.getOwner().getState(State::Energy, false, groups).getPotentialEnergy();
     double kT = BOLTZ*context.getParameter(MonteCarloMembraneBarostat::Temperature());
-    double w = finalEnergy-initialEnergy + pressure*deltaVolume - tension*deltaArea - context.getMolecules().size()*kT*log(newVolume/volume);
+    double w = finalEnergy-initialEnergy + pressure*deltaVolume - tension*deltaArea - numberOfScaledParticles*kT*log(newVolume/volume);
     if (w > 0 && SimTKOpenMMUtilities::getUniformlyDistributedRandomNumber() > exp(-w/kT)) {
         // Reject the step.
         

platforms/reference/tests/TestReferenceMonteCarloMembraneBarostat.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,13 +34,14 @@
 #include "openmm/internal/AssertionUtilities.h"
 #include "openmm/MonteCarloMembraneBarostat.h"
 #include "openmm/Context.h"
-#include "ReferencePlatform.h"
+#include "openmm/CustomIntegrator.h"
 #include "openmm/HarmonicAngleForce.h"
 #include "openmm/HarmonicBondForce.h"
 #include "openmm/NonbondedForce.h"
 #include "openmm/System.h"
 #include "openmm/LangevinIntegrator.h"
 #include "openmm/VerletIntegrator.h"
+#include "ReferencePlatform.h"
 #include "sfmt/SFMT.h"
 #include "SimTKOpenMMRealType.h"
 #include <iostream>
@@ -118,6 +119,73 @@ void testIdealGas(MonteCarloMembraneBarostat::XYMode xymode, MonteCarloMembraneB
     }
 }
 
+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);
+    MonteCarloMembraneBarostat* barostat = new MonteCarloMembraneBarostat(1.0, 0.0, 300.0, MonteCarloMembraneBarostat::XYAnisotropic, MonteCarloMembraneBarostat::ZFree, 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 testRandomSeed() {
     const int numParticles = 8;
     const double temp = 100.0;
@@ -256,6 +324,8 @@ int main() {
         testIdealGas(MonteCarloMembraneBarostat::XYAnisotropic, MonteCarloMembraneBarostat::ZFree);
         testIdealGas(MonteCarloMembraneBarostat::XYAnisotropic, MonteCarloMembraneBarostat::ZFixed);
         testIdealGas(MonteCarloMembraneBarostat::XYAnisotropic, MonteCarloMembraneBarostat::ConstantVolume);
+        testMoleculeScaling(true);
+        testMoleculeScaling(false);
         testRandomSeed();
         testCrystal();
     }

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/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().