Added computeCurrentPressure() to MonteCarloBarostat (#4881)

* Added computeCurrentPressure() to MonteCarloBarostat

* Use instantaneous temperature to compute pressure

* Added computeCurrentPressure() to MonteCarloAnisotropicBarostat

* Added computeCurrentPressure() to MonteCarloMembraneBarostat

* Fixed compilation error

* Fixed error in typemap

* Added documentation on computing pressure

* Fixed CUDA compilation errors

* Made test case more robust

* Made a test case more robust

* Added computeCurrentPressure() to MonteCarloFlexibleBarostat

* Fixed compilation error

* More documentation on computing pressure

platforms/common/src/CommonKernels.cpp

@@ -3731,7 +3731,8 @@ void CommonApplyAndersenThermostatKernel::execute(ContextImpl& context) {
     kernel->execute(cc.getNumAtoms());
 }
 
-void CommonApplyMonteCarloBarostatKernel::initialize(const System& system, const Force& thermostat, bool rigidMolecules) {
+void CommonApplyMonteCarloBarostatKernel::initialize(const System& system, const Force& thermostat, int components, bool rigidMolecules) {
+    this->components = components;
     this->rigidMolecules = rigidMolecules;
     ContextSelector selector(cc);
     savedPositions.initialize(cc, cc.getPaddedNumAtoms(), cc.getUseDoublePrecision() ? sizeof(mm_double4) : sizeof(mm_float4), "savedPositions");
@@ -3744,8 +3745,18 @@ void CommonApplyMonteCarloBarostatKernel::initialize(const System& system, const
     catch (...) {
         // The CUDA platform doesn't have a floating point force buffer, so we don't need to copy it.
     }
-    ComputeProgram program = cc.compileProgram(CommonKernelSources::monteCarloBarostat);
+    energyBuffers.resize(components);
+    for (int i = 0; i < components; i++)
+        energyBuffers[i].initialize(cc, cc.getNumThreadBlocks(), cc.getUseDoublePrecision() || cc.getUseMixedPrecision() ? sizeof(double) : sizeof(float), "energyBuffer");
+    vector<float> zeros(energyBuffers[0].getSize(), 0.0f);
+    for (int i = 0; i < components; i++)
+        energyBuffers[i].upload(zeros, true);
+    map<string, string> defines;
+    defines["WORK_GROUP_SIZE"] = cc.intToString(cc.ThreadBlockSize);
+    defines["COMPONENTS"] = cc.intToString(components);
+    ComputeProgram program = cc.compileProgram(CommonKernelSources::monteCarloBarostat, defines);
     kernel = program->createKernel("scalePositions");
+    kineticEnergyKernel = program->createKernel("computeMolecularKineticEnergy");
 }
 
 void CommonApplyMonteCarloBarostatKernel::saveCoordinates(ContextImpl& context) {
@@ -3804,12 +3815,18 @@ void CommonApplyMonteCarloBarostatKernel::scaleCoordinates(ContextImpl& context,
         kernel->addArg(cc.getPosq());
         kernel->addArg(moleculeAtoms);
         kernel->addArg(moleculeStartIndex);
+        kineticEnergyKernel->addArg(numMolecules);
+        kineticEnergyKernel->addArg(cc.getVelm());
+        kineticEnergyKernel->addArg(moleculeAtoms);
+        kineticEnergyKernel->addArg(moleculeStartIndex);
+        for (int i = 0; i < components; i++)
+            kineticEnergyKernel->addArg(energyBuffers[i]);
     }
     kernel->setArg(0, (float) scaleX);
     kernel->setArg(1, (float) scaleY);
     kernel->setArg(2, (float) scaleZ);
     setPeriodicBoxArgs(cc, kernel, 4);
-    kernel->execute(cc.getNumAtoms());
+    kernel->execute(numMolecules);
 }
 
 void CommonApplyMonteCarloBarostatKernel::restoreCoordinates(ContextImpl& context) {
@@ -3826,6 +3843,27 @@ void CommonApplyMonteCarloBarostatKernel::restoreCoordinates(ContextImpl& contex
         cc.setAtomIndex(lastAtomOrder);
 }
 
+void CommonApplyMonteCarloBarostatKernel::computeKineticEnergy(ContextImpl& context, vector<double>& ke) {
+    ContextSelector selector(cc);
+    ke.resize(components);
+    kineticEnergyKernel->execute(numMolecules);
+    for (int j = 0; j < components; j++) {
+        ke[j] = 0.0;
+        if (energyBuffers[j].getElementSize() == sizeof(float)) {
+            vector<float> buffer;
+            energyBuffers[j].download(buffer);
+            for (int i = 0; i < buffer.size(); i++)
+                ke[j] += buffer[i];
+        }
+        else {
+            vector<double> buffer;
+            energyBuffers[j].download(buffer);
+            for (int i = 0; i < buffer.size(); i++)
+                ke[j] += buffer[i];
+        }
+    }
+}
+
 class CommonCalcATMForceKernel::ReorderListener : public ComputeContext::ReorderListener {
 public:
     ReorderListener(ComputeContext& cc, ArrayInterface& invAtomOrder) : cc(cc), invAtomOrder(invAtomOrder) {

platforms/common/src/kernels/monteCarloBarostat.cc

@@ -39,3 +39,72 @@ KERNEL void scalePositions(float scaleX, float scaleY, float scaleZ, int numMole
         }
     }
 }
+
+/**
+ * Compute the kinetic energy of molecular centers of mass.  Depending on the value
+ * of the COMPONENTS macro, this may compute either 1) the total kinetic energy,
+ * 2) three components of kinetic energy corresponding to the three axes, or
+ * 3) six components corresponding to all elements of the box tensor.
+ */
+KERNEL void computeMolecularKineticEnergy(int numMolecules, GLOBAL mixed4* RESTRICT velm, GLOBAL const int* RESTRICT moleculeAtoms,
+        GLOBAL const int* RESTRICT moleculeStartIndex,
+#if COMPONENTS == 1
+        GLOBAL mixed* RESTRICT energyBuffer
+#else
+        GLOBAL mixed* RESTRICT energyBuffer1, GLOBAL mixed* RESTRICT energyBuffer2, GLOBAL mixed* RESTRICT energyBuffer3
+#if COMPONENTS == 6
+        , GLOBAL mixed* RESTRICT energyBuffer4, GLOBAL mixed* RESTRICT energyBuffer5, GLOBAL mixed* RESTRICT energyBuffer6
+#endif
+#endif
+    ) {
+#if COMPONENTS == 1
+    GLOBAL mixed* buffers[] = {energyBuffer};
+#elif COMPONENTS == 3
+    GLOBAL mixed* buffers[] = {energyBuffer1, energyBuffer2, energyBuffer3};
+#else
+    GLOBAL mixed* buffers[] = {energyBuffer1, energyBuffer2, energyBuffer3, energyBuffer4, energyBuffer5, energyBuffer6};
+#endif
+    mixed ke[COMPONENTS];
+    for (int i = 0; i < COMPONENTS; i++)
+        ke[i] = 0;
+    for (int index = GLOBAL_ID; index < numMolecules; index += GLOBAL_SIZE) {
+        int first = moleculeStartIndex[index];
+        int last = moleculeStartIndex[index+1];
+        int numAtoms = last-first;
+        mixed3 molVel = make_mixed3(0);
+        float molMass = 0;
+        for (int atom = first; atom < last; atom++) {
+            mixed4 v = velm[moleculeAtoms[atom]];
+            float mass = (v.w == 0 ? 0 : 1/v.w);
+            molVel += mass*trimTo3(v);
+            molMass += mass;
+        }
+        molVel *= RECIP((mixed) numAtoms);
+#if COMPONENTS == 1
+        ke[0] += 0.5f*molMass*dot(molVel, molVel);
+#else
+        ke[0] += 0.5f*molMass*molVel.x*molVel.x;
+        ke[1] += 0.5f*molMass*molVel.y*molVel.y;
+        ke[2] += 0.5f*molMass*molVel.z*molVel.z;
+#if COMPONENTS == 6
+        ke[3] += 0.5f*molMass*molVel.x*molVel.y;
+        ke[4] += 0.5f*molMass*molVel.x*molVel.z;
+        ke[5] += 0.5f*molMass*molVel.y*molVel.z;
+#endif
+#endif
+    }
+
+    // Sum the contributions from all the threads in this block and write the result.
+
+    LOCAL mixed tempBuffer[WORK_GROUP_SIZE];
+    for (int j = 0; j < COMPONENTS; j++) {
+        tempBuffer[LOCAL_ID] = ke[j];
+        for (int i = 1; i < WORK_GROUP_SIZE; i *= 2) {
+            SYNC_THREADS;
+            if (LOCAL_ID%(i*2) == 0 && LOCAL_ID+i < WORK_GROUP_SIZE)
+                tempBuffer[LOCAL_ID] += tempBuffer[LOCAL_ID+i];
+        }
+        if (LOCAL_ID == 0)
+            buffers[j][GROUP_ID] = tempBuffer[0];
+    }
+}

platforms/cuda/tests/TestCudaMonteCarloAnisotropicBarostat.cpp

@@ -33,4 +33,5 @@
 #include "TestMonteCarloAnisotropicBarostat.h"
 
 void runPlatformTests() {
+    testLJPressure();
 }

platforms/cuda/tests/TestCudaMonteCarloBarostat.cpp

@@ -34,4 +34,5 @@
 
 void runPlatformTests() {
     testWater();
+    testLJPressure();
 }

platforms/hip/tests/TestHipMonteCarloAnisotropicBarostat.cpp

@@ -34,4 +34,5 @@
 #include "TestMonteCarloAnisotropicBarostat.h"
 
 void runPlatformTests() {
+    testLJPressure();
 }

platforms/hip/tests/TestHipMonteCarloBarostat.cpp

@@ -35,4 +35,5 @@
 
 void runPlatformTests() {
     testWater();
+    testLJPressure();
 }

platforms/opencl/tests/TestOpenCLMonteCarloAnisotropicBarostat.cpp

@@ -33,4 +33,5 @@
 #include "TestMonteCarloAnisotropicBarostat.h"
 
 void runPlatformTests() {
+    testLJPressure();
 }

platforms/opencl/tests/TestOpenCLMonteCarloBarostat.cpp

@@ -34,4 +34,5 @@
 
 void runPlatformTests() {
     testWater();
+    testLJPressure();
 }

platforms/reference/include/ReferenceKernels.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2024 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2025 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -1587,8 +1587,10 @@ public:
      * @param system          the System this kernel will be applied to
      * @param barostat        the MonteCarloBarostat this kernel will be used for
      * @param rigidMolecules  whether molecules should be kept rigid while scaling coordinates
+     * @param components      the number of box components the barostat operates one (1 for isotropic scaling,
+     *                        3 for anisotropic, 6 for both lengths and angles)
      */
-    void initialize(const System& system, const Force& barostat, bool rigidMolecules=true);
+    void initialize(const System& system, const Force& barostat, int components, bool rigidMolecules=true);
     /**
      * Save the coordinates before attempting a Monte Carlo step.  This allows us to restore them
      * if the step is rejected.
@@ -1616,8 +1618,19 @@ public:
      * @param context    the context in which to execute this kernel
      */
     void restoreCoordinates(ContextImpl& context);
+    /**
+     * Compute the kinetic energy of the system.  If initialize() was called with rigidMolecules=true, this
+     * should include only the translational center of mass motion of molecules.  Otherwise it should include
+     * the total kinetic energy of all particles.  This is used when computing instantaneous pressure.
+     * 
+     * @param context    the context in which to execute this kernel
+     * @param ke         a vector to store the kinetic energy components into.  On output, its length will
+     *                   equal the number of components passed to initialize().
+     */
+    void computeKineticEnergy(ContextImpl& context, std::vector<double>& ke);
 private:
     bool rigidMolecules;
+    int components;
     ReferenceMonteCarloBarostat* barostat;
 };
 

platforms/reference/include/ReferenceMonteCarloBarostat.h

@@ -37,6 +37,7 @@ class ReferenceMonteCarloBarostat {
 
        std::vector<double> savedAtomPositions[3];
        std::vector<std::vector<int> > molecules;
+       std::vector<double> masses;
 
    public:
 
@@ -46,7 +47,7 @@ class ReferenceMonteCarloBarostat {
 
          --------------------------------------------------------------------------------------- */
 
-       ReferenceMonteCarloBarostat(int numAtoms, const std::vector<std::vector<int> >& molecules);
+       ReferenceMonteCarloBarostat(int numAtoms, const std::vector<std::vector<int> >& molecules, const std::vector<double>& masses);
 
       /**---------------------------------------------------------------------------------------
 
@@ -90,6 +91,17 @@ class ReferenceMonteCarloBarostat {
 
       void restorePositions(std::vector<OpenMM::Vec3>& atomPositions);
 
+      /**---------------------------------------------------------------------------------------
+
+         Compute the kinetic energy of molecular center of mass translations.
+
+         @param velocities      atom velocities
+         @param ke              a vector to store the kinetic energy components into
+         @param components      the number of components to calculate
+
+         --------------------------------------------------------------------------------------- */
+
+      void computeMolecularKineticEnergy(const std::vector<OpenMM::Vec3>& velocities, std::vector<double>& ke, int components);
 };
 
 } // namespace OpenMM

platforms/reference/src/ReferenceKernels.cpp

@@ -2921,19 +2921,24 @@ ReferenceApplyMonteCarloBarostatKernel::~ReferenceApplyMonteCarloBarostatKernel(
         delete barostat;
 }
 
-void ReferenceApplyMonteCarloBarostatKernel::initialize(const System& system, const Force& barostat, bool rigidMolecules) {
+void ReferenceApplyMonteCarloBarostatKernel::initialize(const System& system, const Force& barostat, int components, bool rigidMolecules) {
+    this->components = components;
     this->rigidMolecules = rigidMolecules;
 }
 
 void ReferenceApplyMonteCarloBarostatKernel::saveCoordinates(ContextImpl& context) {
     if (barostat == NULL) {
+        const System& system = context.getSystem();
+        vector<double> masses;
+        for (int i = 0; i < system.getNumParticles(); i++)
+            masses.push_back(system.getParticleMass(i));
         if (rigidMolecules)
-            barostat = new ReferenceMonteCarloBarostat(context.getSystem().getNumParticles(), context.getMolecules());
+            barostat = new ReferenceMonteCarloBarostat(system.getNumParticles(), context.getMolecules(), masses);
         else {
-            vector<vector<int> > molecules(context.getSystem().getNumParticles());
+            vector<vector<int> > molecules(system.getNumParticles());
             for (int i = 0; i < molecules.size(); i++)
                 molecules[i].push_back(i);
-            barostat = new ReferenceMonteCarloBarostat(context.getSystem().getNumParticles(), molecules);
+            barostat = new ReferenceMonteCarloBarostat(system.getNumParticles(), molecules, masses);
         }
     }
     vector<Vec3>& posData = extractPositions(context);
@@ -2951,6 +2956,10 @@ void ReferenceApplyMonteCarloBarostatKernel::restoreCoordinates(ContextImpl& con
     barostat->restorePositions(posData);
 }
 
+void ReferenceApplyMonteCarloBarostatKernel::computeKineticEnergy(ContextImpl& context, vector<double>& ke) {
+    barostat->computeMolecularKineticEnergy(extractVelocities(context), ke, components);
+}
+
 void ReferenceRemoveCMMotionKernel::initialize(const System& system, const CMMotionRemover& force) {
     frequency = force.getFrequency();
     masses.resize(system.getNumParticles());

platforms/reference/src/SimTKReference/ReferenceMonteCarloBarostat.cpp

@@ -37,7 +37,8 @@ using namespace OpenMM;
 
   --------------------------------------------------------------------------------------- */
 
-ReferenceMonteCarloBarostat::ReferenceMonteCarloBarostat(int numAtoms, const vector<vector<int> >& molecules) : molecules(molecules) {
+ReferenceMonteCarloBarostat::ReferenceMonteCarloBarostat(int numAtoms, const vector<vector<int> >& molecules, const std::vector<double>& masses) :
+            molecules(molecules), masses(masses) {
     savedAtomPositions[0].resize(numAtoms);
     savedAtomPositions[1].resize(numAtoms);
     savedAtomPositions[2].resize(numAtoms);
@@ -120,3 +121,30 @@ void ReferenceMonteCarloBarostat::restorePositions(vector<Vec3>& atomPositions)
         for (int j = 0; j < 3; j++)
             atomPositions[i][j] = savedAtomPositions[j][i];
 }
+
+void ReferenceMonteCarloBarostat::computeMolecularKineticEnergy(const vector<Vec3>& velocities, vector<double>& ke, int components) {
+    ke.resize(components);
+    for (int i = 0; i < components; i++)
+        ke[i] = 0.0;
+    for (auto& molecule : molecules) {
+        Vec3 molVel;
+        double molMass = 0.0;
+        for (int atom : molecule) {
+            molVel += masses[atom]*velocities[atom];
+            molMass += masses[atom];
+        }
+        molVel /= molecule.size();
+        if (components == 1)
+            ke[0] += 0.5*molMass*molVel.dot(molVel);
+        else {
+            ke[0] += 0.5*molMass*molVel[0]*molVel[0];
+            ke[1] += 0.5*molMass*molVel[1]*molVel[1];
+            ke[2] += 0.5*molMass*molVel[2]*molVel[2];
+            if (components == 6) {
+                ke[3] += 0.5*molMass*molVel[1]*molVel[0];
+                ke[4] += 0.5*molMass*molVel[2]*molVel[0];
+                ke[5] += 0.5*molMass*molVel[2]*molVel[1];
+            }
+        }
+    }
+}

platforms/reference/tests/TestReferenceMonteCarloMembraneBarostat.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2016 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2025 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -37,6 +37,7 @@
 #include "openmm/MonteCarloMembraneBarostat.h"
 #include "openmm/Context.h"
 #include "ReferencePlatform.h"
+#include "openmm/HarmonicAngleForce.h"
 #include "openmm/HarmonicBondForce.h"
 #include "openmm/NonbondedForce.h"
 #include "openmm/System.h"
@@ -182,6 +183,73 @@ void testRandomSeed() {
     }
 }
 
+void testCrystal() {
+    const int gridSize = 5;
+    const int numParticles = gridSize*gridSize*gridSize;
+    const int frequency = 5;
+    const int steps = 10000;
+    const double pressure = 15;
+    const double tension = 15;
+    const double temp = 300.0;
+    const double spacing = 1.0;
+    const double bondK = 20.0;
+    const double angleK = 20.0;
+
+    // Create a periodic crystal.
+
+    System system;
+    vector<vector<vector<int> > > index(gridSize, vector<vector<int> >(gridSize, vector<int>(gridSize)));
+    system.setDefaultPeriodicBoxVectors(Vec3(gridSize*spacing, 0, 0), Vec3(0, gridSize*spacing, 0), Vec3(0, 0, gridSize*spacing));
+    vector<Vec3> positions;
+    for (int i = 0; i < gridSize; i++)
+        for (int j = 0; j < gridSize; j++)
+            for (int k = 0; k < gridSize; k++) {
+                index[i][j][k] = system.getNumParticles();
+                system.addParticle(1.0);
+                positions.push_back(Vec3(i, j, k)*spacing);
+            }
+    HarmonicBondForce* bonds = new HarmonicBondForce();
+    HarmonicAngleForce* angles = new HarmonicAngleForce();
+    bonds->setUsesPeriodicBoundaryConditions(true);
+    angles->setUsesPeriodicBoundaryConditions(true);
+    system.addForce(bonds);
+    system.addForce(angles);
+    for (int i = 0; i < gridSize; i++)
+        for (int j = 0; j < gridSize; j++)
+            for (int k = 0; k < gridSize; k++) {
+                bonds->addBond(index[i][j][k], index[(i+1)%gridSize][j][k], spacing, bondK);
+                bonds->addBond(index[i][j][k], index[i][(j+1)%gridSize][k], spacing, bondK);
+                bonds->addBond(index[i][j][k], index[i][j][(k+1)%gridSize], spacing, bondK);
+                angles->addAngle(index[(i+1)%gridSize][j][k], index[i][j][k], index[i][(j+1)%gridSize][k], M_PI/2, angleK);
+                angles->addAngle(index[(i+1)%gridSize][j][k], index[i][j][k], index[i][j][(k+1)%gridSize], M_PI/2, angleK);
+                angles->addAngle(index[i][(j+1)%gridSize][k], index[i][j][k], index[i][j][(k+1)%gridSize], M_PI/2, angleK);
+            }
+    MonteCarloMembraneBarostat* barostat = new MonteCarloMembraneBarostat(pressure, tension, temp, MonteCarloMembraneBarostat::XYAnisotropic, MonteCarloMembraneBarostat::ZFree, 1);
+    system.addForce(barostat);
+
+    // Simulate it, seeing if the pressure is correct.
+
+    LangevinIntegrator integrator(temp, 10.0, 0.002);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    context.setVelocitiesToTemperature(temp);
+    integrator.step(1000);
+    Vec3 size;
+    Vec3 avgPressure;
+    for (int j = 0; j < steps; ++j) {
+        Vec3 box[3];
+        context.getState(0).getPeriodicBoxVectors(box[0], box[1], box[2]);
+        size += Vec3(box[0][0], box[1][1], box[2][2]);
+        integrator.step(frequency);
+        avgPressure += barostat->computeCurrentPressure(context);
+    }
+    size /= steps;
+    avgPressure /= steps;
+    ASSERT_USUALLY_EQUAL_TOL(pressure-tension/size[0], avgPressure[0], 0.15);
+    ASSERT_USUALLY_EQUAL_TOL(pressure-tension/size[1], avgPressure[1], 0.15);
+    ASSERT_USUALLY_EQUAL_TOL(pressure, avgPressure[2], 0.15);
+}
+
 int main() {
     try {
         testIdealGas(MonteCarloMembraneBarostat::XYIsotropic, MonteCarloMembraneBarostat::ZFree);
@@ -191,6 +259,7 @@ int main() {
         testIdealGas(MonteCarloMembraneBarostat::XYAnisotropic, MonteCarloMembraneBarostat::ZFixed);
         testIdealGas(MonteCarloMembraneBarostat::XYAnisotropic, MonteCarloMembraneBarostat::ConstantVolume);
         testRandomSeed();
+        testCrystal();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;

plugins/rpmd/openmmapi/src/RPMDMonteCarloBarostatImpl.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2010-2023 Stanford University and the Authors.      *
+ * Portions copyright (c) 2010-2025 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -55,7 +55,7 @@ void RPMDMonteCarloBarostatImpl::initialize(ContextImpl& context) {
     if (!integrator->getApplyThermostat())
         throw OpenMMException("RPMDMonteCarloBarostat requires the integrator's thermostat to be enabled");;
     kernel = context.getPlatform().createKernel(ApplyMonteCarloBarostatKernel::Name(), context);
-    kernel.getAs<ApplyMonteCarloBarostatKernel>().initialize(context.getSystem(), owner);
+    kernel.getAs<ApplyMonteCarloBarostatKernel>().initialize(context.getSystem(), owner, 1);
     savedPositions.resize(integrator->getNumCopies());
     Vec3 box[3];
     context.getPeriodicBoxVectors(box[0], box[1], box[2]);

tests/TestMonteCarloAnisotropicBarostat.h

@@ -51,7 +51,7 @@ void testIdealGas() {
     const int numParticles = 64;
     const int frequency = 10;
     const int steps = 1000;
-    const double pressure = 1.5;
+    const double pressure = 3.0;
     const double pressureInMD = pressure*(AVOGADRO*1e-25); // pressure in kJ/mol/nm^3
     const double temp[] = {300.0, 600.0, 1000.0};
     const double initialVolume = numParticles*BOLTZ*temp[1]/pressureInMD;
@@ -84,20 +84,26 @@ void testIdealGas() {
         
         // Let it equilibrate.
         
-        integrator.step(10000);
+        integrator.step(5000);
         
         // Now run it for a while and see if the volume is correct.
         
         double volume = 0.0;
+        Vec3 avgPressure;
         for (int j = 0; j < steps; ++j) {
             Vec3 box[3];
             context.getState(0).getPeriodicBoxVectors(box[0], box[1], box[2]);
             volume += box[0][0]*box[1][1]*box[2][2];
             integrator.step(frequency);
+            avgPressure += barostat->computeCurrentPressure(context);
         }
         volume /= steps;
+        avgPressure /= steps;
         double expected = (numParticles+1)*BOLTZ*temp[i]/pressureInMD;
         ASSERT_USUALLY_EQUAL_TOL(expected, volume, 3/std::sqrt((double) steps));
+        ASSERT_USUALLY_EQUAL_TOL(pressure, avgPressure[0], 0.2);
+        ASSERT_USUALLY_EQUAL_TOL(pressure, avgPressure[1], 0.2);
+        ASSERT_USUALLY_EQUAL_TOL(pressure, avgPressure[2], 0.2);
     }
 }
 
@@ -106,7 +112,7 @@ void testIdealGasAxis(int axis) {
     const int numParticles = 64;
     const int frequency = 10;
     const int steps = 1000;
-    const double pressure = 1.5;
+    const double pressure = 3.0;
     const double pressureInMD = pressure*(AVOGADRO*1e-25); // pressure in kJ/mol/nm^3
     const double temp[] = {300.0, 600.0, 1000.0};
     const double initialVolume = numParticles*BOLTZ*temp[1]/pressureInMD;
@@ -145,11 +151,12 @@ void testIdealGasAxis(int axis) {
         
         // Let it equilibrate.
         
-        integrator.step(10000);
+        integrator.step(5000);
         
         // Now run it for a while and see if the volume is correct.
         
         double volume = 0.0;
+        Vec3 avgPressure;
         for (int j = 0; j < steps; ++j) {
             Vec3 box[3];
             context.getState(0).getPeriodicBoxVectors(box[0], box[1], box[2]);
@@ -158,10 +165,15 @@ void testIdealGasAxis(int axis) {
             boxZ = box[2][2];
             volume += box[0][0]*box[1][1]*box[2][2];
             integrator.step(frequency);
+            avgPressure += barostat->computeCurrentPressure(context);
         }
         volume /= steps;
+        avgPressure /= steps;
         double expected = (numParticles+1)*BOLTZ*temp[i]/pressureInMD;
         ASSERT_USUALLY_EQUAL_TOL(expected, volume, 3/std::sqrt((double) steps));
+        ASSERT_USUALLY_EQUAL_TOL(pressure, avgPressure[0], 0.2);
+        ASSERT_USUALLY_EQUAL_TOL(pressure, avgPressure[1], 0.2);
+        ASSERT_USUALLY_EQUAL_TOL(pressure, avgPressure[2], 0.2);
         if (!scaleX) {
             ASSERT(boxX == initialLength);
         }
@@ -174,6 +186,57 @@ void testIdealGasAxis(int axis) {
     }
 }
 
+void testMolecularGas() {
+    const int numMolecules = 64;
+    const int frequency = 5;
+    const int steps = 5000;
+    const double pressure = 3.0;
+    const double pressureInMD = pressure*(AVOGADRO*1e-25); // pressure in kJ/mol/nm^3
+    const double temp = 300.0;
+    const double initialVolume = numMolecules*BOLTZ*temp/pressureInMD;
+    const double initialLength = std::pow(initialVolume, 1.0/3.0);
+
+    // Create a gas of noninteracting molecules.
+
+    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);
+    system.addForce(barostat);
+    HarmonicBondForce* bonds = new HarmonicBondForce();
+    bonds->setUsesPeriodicBoundaryConditions(true);
+    system.addForce(bonds);
+    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);
+        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);
+        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));
+    }
+
+    // Simulate it and see if the pressure is correct.
+
+    LangevinIntegrator integrator(temp, 0.1, 0.005);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    integrator.step(5000);
+    Vec3 avgPressure;
+    for (int j = 0; j < steps; ++j) {
+        integrator.step(frequency);
+        avgPressure += barostat->computeCurrentPressure(context);
+    }
+    avgPressure /= steps;
+    ASSERT_USUALLY_EQUAL_TOL(pressure, avgPressure[0], 0.2);
+    ASSERT_USUALLY_EQUAL_TOL(pressure, avgPressure[1], 0.2);
+    ASSERT_USUALLY_EQUAL_TOL(pressure, avgPressure[2], 0.2);
+}
+
 void testRandomSeed() {
     const int numParticles = 8;
     const double temp = 100.0;
@@ -275,7 +338,7 @@ void testTriclinic() {
 
     // Let it equilibrate.
 
-    integrator.step(10000);
+    integrator.step(5000);
 
     // Now run it for a while and see if the volume is correct.
 
@@ -447,6 +510,54 @@ void testEinsteinCrystal() {
     ASSERT_USUALLY_EQUAL_TOL(results[11], results[14], 3/std::sqrt((double) steps));
 }
 
+void testLJPressure() {
+    const int gridSize = 8;
+    const int frequency = 10;
+    const int steps = 1000;
+    const double spacing = 1.2;
+    const double temp = 150.0;
+    const double pressure = 10.0;
+
+    // Create a gas of LJ particles.  This is much more compressible than water, which means
+    // the fluctuations in pressure are much smaller and it's easier to compute the average.
+
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(gridSize*spacing, 0, 0), Vec3(0, gridSize*spacing, 0), Vec3(0, 0, gridSize*spacing));
+    NonbondedForce* nonbonded = new NonbondedForce();
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    vector<Vec3> positions;
+    for (int i = 0; i < gridSize; ++i) {
+        for (int j = 0; j < gridSize; ++j) {
+            for (int k = 0; k < gridSize; ++k) {
+                system.addParticle(1.0);
+                nonbonded->addParticle(0.0, 1.0, 1.0);
+                positions.push_back(Vec3(spacing*i, spacing*j, spacing*k));
+            }
+        }
+    }
+    system.addForce(nonbonded);
+    MonteCarloAnisotropicBarostat* barostat = new MonteCarloAnisotropicBarostat(Vec3(pressure, pressure, pressure), temp, true, true, true, frequency);
+    system.addForce(barostat);
+
+    // Simulate it and see if the average instantaneous pressure matches the heat bath pressure.
+
+    LangevinIntegrator integrator(temp, 1.0, 0.01);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    integrator.step(1000);
+    Vec3 avgPressure;
+    for (int j = 0; j < steps; ++j) {
+        Vec3 box[3];
+        context.getState(0).getPeriodicBoxVectors(box[0], box[1], box[2]);
+        avgPressure += barostat->computeCurrentPressure(context);
+        integrator.step(frequency);
+    }
+    avgPressure /= steps;
+    ASSERT_USUALLY_EQUAL_TOL(pressure, avgPressure[0], 0.1);
+    ASSERT_USUALLY_EQUAL_TOL(pressure, avgPressure[1], 0.1);
+    ASSERT_USUALLY_EQUAL_TOL(pressure, avgPressure[2], 0.1);
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
@@ -456,6 +567,7 @@ int main(int argc, char* argv[]) {
         testIdealGasAxis(0);
         testIdealGasAxis(1);
         testIdealGasAxis(2);
+        testMolecularGas();
         testRandomSeed();
         testTriclinic();
         //testEinsteinCrystal();

tests/TestMonteCarloBarostat.h

@@ -125,6 +125,7 @@ void testIdealGas() {
         // Now run it for a while and see if the volume is correct.
 
         double volume = 0.0;
+        double avgPressure = 0.0;
         for (int j = 0; j < steps; ++j) {
             Vec3 box[3];
             context.getState(0).getPeriodicBoxVectors(box[0], box[1], box[2]);
@@ -132,13 +133,65 @@ void testIdealGas() {
             ASSERT_EQUAL_TOL(0.5*box[0][0], box[1][1], 1e-5);
             ASSERT_EQUAL_TOL(2*box[0][0], box[2][2], 1e-5);
             integrator.step(frequency);
+            avgPressure += barostat->computeCurrentPressure(context);
         }
         volume /= steps;
+        avgPressure /= steps;
         double expected = (numParticles+1)*BOLTZ*temp[i]/pressureInMD;
         ASSERT_USUALLY_EQUAL_TOL(expected, volume, 3/std::sqrt((double) steps));
+        ASSERT_USUALLY_EQUAL_TOL(pressure, avgPressure, 0.1);
     }
 }
 
+void testMolecularGas() {
+    const int numMolecules = 64;
+    const int frequency = 10;
+    const int steps = 1000;
+    const double pressure = 1.5;
+    const double pressureInMD = pressure*(AVOGADRO*1e-25); // pressure in kJ/mol/nm^3
+    const double temp = 300.0;
+    const double initialVolume = numMolecules*BOLTZ*temp/pressureInMD;
+    const double initialLength = std::pow(initialVolume, 1.0/3.0);
+
+    // Create a gas of noninteracting molecules.
+
+    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);
+    system.addForce(barostat);
+    HarmonicBondForce* bonds = new HarmonicBondForce();
+    bonds->setUsesPeriodicBoundaryConditions(true);
+    system.addForce(bonds);
+    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);
+        system.addParticle(1.0);
+        Vec3 pos(initialLength*genrand_real2(sfmt), 0.5*initialLength*genrand_real2(sfmt), 2*initialLength*genrand_real2(sfmt));
+        bonds->addBond(positions.size(), positions.size()+1, 0.1, 10.0);
+        system.addConstraint(positions.size(), positions.size()+2, 0.1);
+        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));
+    }
+
+    // Simulate it and see if the pressure is correct..
+
+    LangevinIntegrator integrator(temp, 0.1, 0.01);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    integrator.step(10000);
+    double avgPressure = 0.0;
+    for (int j = 0; j < steps; ++j) {
+        integrator.step(frequency);
+        avgPressure += barostat->computeCurrentPressure(context);
+    }
+    avgPressure /= steps;
+    ASSERT_USUALLY_EQUAL_TOL(pressure, avgPressure, 0.1);
+}
+
 void testRandomSeed() {
     const int numParticles = 8;
     const double temp = 100.0;
@@ -268,6 +321,52 @@ void testWater() {
     ASSERT_USUALLY_EQUAL_TOL(1.0, density, 0.02);
 }
 
+void testLJPressure() {
+    const int gridSize = 8;
+    const int frequency = 10;
+    const int steps = 1000;
+    const double spacing = 1.2;
+    const double temp = 150.0;
+    const double pressure = 10.0;
+
+    // Create a gas of LJ particles.  This is much more compressible than water, which means
+    // the fluctuations in pressure are much smaller and it's easier to compute the average.
+
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(gridSize*spacing, 0, 0), Vec3(0, gridSize*spacing, 0), Vec3(0, 0, gridSize*spacing));
+    NonbondedForce* nonbonded = new NonbondedForce();
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    vector<Vec3> positions;
+    for (int i = 0; i < gridSize; ++i) {
+        for (int j = 0; j < gridSize; ++j) {
+            for (int k = 0; k < gridSize; ++k) {
+                system.addParticle(1.0);
+                nonbonded->addParticle(0.0, 1.0, 1.0);
+                positions.push_back(Vec3(spacing*i, spacing*j, spacing*k));
+            }
+        }
+    }
+    system.addForce(nonbonded);
+    MonteCarloBarostat* barostat = new MonteCarloBarostat(pressure, temp, frequency);
+    system.addForce(barostat);
+
+    // Simulate it and see if the average instantaneous pressure matches the heat bath pressure.
+
+    LangevinIntegrator integrator(temp, 1.0, 0.01);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    integrator.step(1000);
+    double avgPressure = 0.0;
+    for (int j = 0; j < steps; ++j) {
+        Vec3 box[3];
+        context.getState(0).getPeriodicBoxVectors(box[0], box[1], box[2]);
+        avgPressure += barostat->computeCurrentPressure(context);
+        integrator.step(frequency);
+    }
+    avgPressure /= steps;
+    ASSERT_USUALLY_EQUAL_TOL(pressure, avgPressure, 0.1);
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
@@ -275,9 +374,10 @@ int main(int argc, char* argv[]) {
         initializeTests(argc, argv);
         testChangingBoxSize();
         testIdealGas();
+        testMolecularGas();
         testRandomSeed();
-        // Don't run testWater() here, because it's very slow on Reference platform.
-        // Individual platforms can run it from runPlatformTests().
+        // Don't run testWater() or testLJPressure() here, because they're very slow on
+        // Reference platform.  Individual platforms can run them from runPlatformTests().
         runPlatformTests();
     }
     catch(const exception& e) {

tests/TestMonteCarloFlexibleBarostat.h

@@ -50,7 +50,7 @@ void testIdealGas() {
     const int numParticles = 64;
     const int frequency = 10;
     const int steps = 1000;
-    const double pressure = 1.5;
+    const double pressure = 3.0;
     const double pressureInMD = pressure*(AVOGADRO*1e-25); // pressure in kJ/mol/nm^3
     const double temp[] = {300.0, 600.0, 1000.0};
     const double initialVolume = numParticles*BOLTZ*temp[1]/pressureInMD;
@@ -83,20 +83,33 @@ void testIdealGas() {
 
         // Let it equilibrate.
 
-        integrator.step(10000);
+        integrator.step(5000);
 
         // Now run it for a while and see if the volume is correct.
 
         double volume = 0.0;
+        vector<double> avgPressure(6, 0.0);
         for (int j = 0; j < steps; ++j) {
             Vec3 box[3];
             context.getState(0).getPeriodicBoxVectors(box[0], box[1], box[2]);
             volume += box[0][0]*box[1][1]*box[2][2];
             integrator.step(frequency);
+            vector<double> pressure;
+            barostat->computeCurrentPressure(context, pressure);
+            for (int j = 0; j < 6; j++)
+                avgPressure[j] += pressure[j];
         }
         volume /= steps;
         double expected = (numParticles+1)*BOLTZ*temp[i]/pressureInMD;
         ASSERT_USUALLY_EQUAL_TOL(expected, volume, 3/std::sqrt((double) steps));
+        for (int i = 0; i < 3; i++) {
+            avgPressure[i] /= steps;
+            ASSERT_USUALLY_EQUAL_TOL(pressure, avgPressure[i], 0.2);
+        }
+        for (int i = 3; i < 6; i++) {
+            avgPressure[i] /= steps;
+            ASSERT_USUALLY_EQUAL_TOL(0.0, avgPressure[i], 0.4);
+        }
     }
 }
 
@@ -168,6 +181,65 @@ void testMoleculeScaling(bool rigid) {
     }
 }
 
+void testMolecularGas(bool rigid) {
+    const int numMolecules = 64;
+    const int frequency = 5;
+    const int steps = 5000;
+    const double pressure = 3.0;
+    const double pressureInMD = pressure*(AVOGADRO*1e-25); // pressure in kJ/mol/nm^3
+    const double temp =300.0;
+    const double initialVolume = numMolecules*BOLTZ*temp/pressureInMD;
+    const double initialLength = std::pow(initialVolume, 1.0/3.0);
+
+    // Create a gas of noninteracting molecules.
+
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(initialLength, 0, 0), Vec3(0, 0.5*initialLength, 0), Vec3(0, 0, 2*initialLength));
+    MonteCarloFlexibleBarostat* barostat = new MonteCarloFlexibleBarostat(pressure, temp, frequency, rigid);
+    system.addForce(barostat);
+    HarmonicBondForce* bonds = new HarmonicBondForce();
+    bonds->setUsesPeriodicBoundaryConditions(true);
+    system.addForce(bonds);
+    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);
+        system.addParticle(1.0);
+        Vec3 pos(initialLength*genrand_real2(sfmt), 0.5*initialLength*genrand_real2(sfmt), 2*initialLength*genrand_real2(sfmt));
+        bonds->addBond(positions.size(), positions.size()+1, 0.1, 0.0);
+        bonds->addBond(positions.size(), positions.size()+2, 0.1, 0.0);
+        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));
+    }
+
+    // Simulate it and see if the pressure is correct.
+
+    LangevinIntegrator integrator(temp, 0.1, 0.005);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    context.setVelocitiesToTemperature(temp);
+    integrator.step(5000);
+    vector<double> avgPressure(6, 0.0);
+    for (int j = 0; j < steps; ++j) {
+        integrator.step(frequency);
+        vector<double> pressure;
+        barostat->computeCurrentPressure(context, pressure);
+        for (int j = 0; j < 6; j++)
+            avgPressure[j] += pressure[j];
+    }
+    for (int i = 0; i < 3; i++) {
+        avgPressure[i] /= steps;
+        ASSERT_USUALLY_EQUAL_TOL(pressure, avgPressure[i], 0.2);
+    }
+    for (int i = 3; i < 6; i++) {
+        avgPressure[i] /= steps;
+        ASSERT_USUALLY_EQUAL_TOL(0.0, avgPressure[i], 0.4);
+    }
+}
+
 void testRandomSeed() {
     const int numParticles = 8;
     const double temp = 100.0;
@@ -238,6 +310,8 @@ int main(int argc, char* argv[]) {
         testIdealGas();
         testMoleculeScaling(true);
         testMoleculeScaling(false);
+        testMolecularGas(true);
+        testMolecularGas(false);
         testRandomSeed();
         runPlatformTests();
     }

wrappers/generateWrappers.py

@@ -79,6 +79,8 @@ class WrapperGenerator:
                             'const std::vector<int>& OpenMM::NoseHooverIntegrator::getAllThermostatedIndividualParticles',
                             'const std::vector<std::tuple<int, int, double> >& OpenMM::NoseHooverIntegrator::getAllThermostatedPairs',
                             'virtual void OpenMM::NoseHooverIntegrator::stateChanged',
+                            'Vec3 OpenMM::MonteCarloAnisotropicBarostat::computeCurrentPressure',
+                            'Vec3 OpenMM::MonteCarloMembraneBarostat::computeCurrentPressure',
                             'virtual bool OpenMM::TabulatedFunction::operator==',
                             'bool OpenMM::Continuous1DFunction::operator==',
                             'bool OpenMM::Continuous2DFunction::operator==',

wrappers/python/src/swig_doxygen/swigInputConfig.py

@@ -95,9 +95,6 @@ SKIP_METHODS = [('State', 'getPositions'),
                 ('IntegrateDrudeSCFStepKernel',),
                 ('XmlSerializer',  'serialize'),
                 ('XmlSerializer',  'deserialize'),
-                ('LocalCoordinatesSite',  'getOriginWeights', 0),
-                ('LocalCoordinatesSite',  'getXWeights', 0),
-                ('LocalCoordinatesSite',  'getYWeights', 0),
                 ("NoseHooverIntegrator", "getAllThermostatedIndividualParticles"),
                 ("NoseHooverIntegrator", "getAllThermostatedPairs"),
 ]

wrappers/python/src/swig_doxygen/swig_lib/python/typemaps.i

@@ -545,7 +545,7 @@ int Py_SequenceToVecVecVecDouble(PyObject* obj, std::vector<std::vector<std::vec
 
 /* The following typemaps handle the ways a Vec3 can be returned from a function. */
 %typemap(out, fragment="Vec3_to_PyVec3") Vec3 {
-    $result = Vec3_to_PyVec3(*$1);
+    $result = Vec3_to_PyVec3($1);
 }
 
 %typemap(out, fragment="Vec3_to_PyVec3") const Vec3& {