Fixed pressure calculation in MOnteCarloFlexibleBarostat (#5251)

openmmapi/src/MonteCarloFlexibleBarostatImpl.cpp

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

tests/TestMonteCarloFlexibleBarostat.h

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