Fixed errors in computeCurrentPressure() (#5114)

platforms/common/src/kernels/monteCarloBarostat.cc

@@ -79,7 +79,7 @@ KERNEL void computeMolecularKineticEnergy(int numMolecules, GLOBAL mixed4* RESTR
             molVel += mass*trimTo3(v);
             molMass += mass;
         }
-        molVel *= RECIP((mixed) molMass);
+        molVel *= (molMass == 0 ? 0 : RECIP((mixed) molMass));
 #if COMPONENTS == 1
         ke[0] += 0.5f*molMass*dot(molVel, molVel);
 #else

platforms/reference/src/SimTKReference/ReferenceMonteCarloBarostat.cpp

@@ -133,7 +133,8 @@ void ReferenceMonteCarloBarostat::computeMolecularKineticEnergy(const vector<Vec
             molVel += masses[atom]*velocities[atom];
             molMass += masses[atom];
         }
-        molVel /= molMass;
+        if (molMass != 0)
+            molVel /= molMass;
         if (components == 1)
             ke[0] += 0.5*molMass*molVel.dot(molVel);
         else {

tests/TestMonteCarloBarostat.h

@@ -164,9 +164,18 @@ void testMolecularGas() {
     OpenMM_SFMT::SFMT sfmt;
     init_gen_rand(0, sfmt);
     for (int i = 0; i < numMolecules; ++i) {
-        system.addParticle(1.0);
-        system.addParticle(2.0);
-        system.addParticle(3.0);
+        if (i == 0) {
+            // Make one molecule massless to make sure that doesn't cause any problems.
+
+            system.addParticle(0.0);
+            system.addParticle(0.0);
+            system.addParticle(0.0);
+        }
+        else {
+            system.addParticle(1.0);
+            system.addParticle(2.0);
+            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);

wrappers/python/src/swig_doxygen/swigInputConfig.py

@@ -493,6 +493,7 @@ UNITS = {
 ("MonteCarloMembraneBarostat", "MonteCarloMembraneBarostat") : (None, ("unit.bar", "unit.bar*unit.nanometer", "unit.kelvin", None, None, None)),
 ("MonteCarloMembraneBarostat", "getXYMode") : (None, ()),
 ("MonteCarloMembraneBarostat", "getZMode") : (None, ()),
+("*", "computeCurrentPressure") : ("unit.bar", ()),
 ("CustomIntegrator", "CustomIntegrator") : (None, ("unit.picosecond",)),
 ("BrownianIntegrator", "BrownianIntegrator") : (None, ("unit.kelvin", "unit.picosecond**-1", "unit.picosecond")),
 ("LangevinIntegrator", "LangevinIntegrator") : (None, ("unit.kelvin", "unit.picosecond**-1", "unit.picosecond")),

wrappers/python/tests/TestAPIUnits.py

@@ -1242,8 +1242,31 @@ class TestAPIUnits(unittest.TestCase):
         self.assertEqual(force.getDefaultTemperature(), 298.15*kelvin)
         self.assertAlmostEqualUnit(force.getDefaultCollisionFrequency(), 1/picosecond)
 
-    def testMonteCarloMembraneBarostat(self):
-        """ Tests the MonteCarloMembraneBarostat API features """
+    def testMonteCarloBarostat(self):
+        """ Tests the various Monte Carlo barostats API features """
+        def checkPressureUnits(force):
+            system = System()
+            system.addParticle(0.0)
+            system.addForce(force)
+            bonds = HarmonicBondForce()
+            bonds.setUsesPeriodicBoundaryConditions(True)
+            system.addForce(bonds)
+            context = Context(system, VerletIntegrator(0.001))
+            context.setPositions([Vec3(0, 0, 0)])
+            pressure = force.computeCurrentPressure(context)
+            self.assertTrue(is_quantity(pressure))
+            self.assertTrue(pressure.unit == bar)
+
+        force = MonteCarloBarostat(1.1*bar, 350*kelvin)
+        self.assertEqual(force.getDefaultPressure(), 1.1*bar)
+        self.assertEqual(force.getDefaultTemperature(), 350*kelvin)
+        checkPressureUnits(force)
+
+        force = MonteCarloFlexibleBarostat(1.1*bar, 350*kelvin)
+        self.assertEqual(force.getDefaultPressure(), 1.1*bar)
+        self.assertEqual(force.getDefaultTemperature(), 350*kelvin)
+        checkPressureUnits(force)
+
         force = MonteCarloMembraneBarostat(1.0, 1.5, 300, MonteCarloMembraneBarostat.XYAnisotropic, MonteCarloMembraneBarostat.ZFixed, 25)
         self.assertEqual(force.getDefaultPressure(), 1.0*bar)
         self.assertEqual(force.getDefaultSurfaceTension(), 1.5*bar*nanometer)
@@ -1251,6 +1274,7 @@ class TestAPIUnits(unittest.TestCase):
         self.assertEqual(force.getXYMode(), MonteCarloMembraneBarostat.XYAnisotropic)
         self.assertEqual(force.getZMode(), MonteCarloMembraneBarostat.ZFixed)
         self.assertEqual(force.getFrequency(), 25)
+        checkPressureUnits(force)
 
         force = MonteCarloMembraneBarostat(1.1*bar, 2.0*bar*nanometer, 350*kelvin, MonteCarloMembraneBarostat.XYAnisotropic, MonteCarloMembraneBarostat.ZFixed, 25)
         self.assertEqual(force.getDefaultPressure(), 1.1*bar)
@@ -1264,6 +1288,11 @@ class TestAPIUnits(unittest.TestCase):
         self.assertEqual(force.getDefaultSurfaceTension(), 2.5*bar*nanometer)
         self.assertEqual(force.getDefaultTemperature(), 298.15*kelvin)
 
+        force = MonteCarloAnisotropicBarostat(Vec3(1.1, 2.2, 3.3)*bar, 350*kelvin)
+        self.assertEqual(force.getDefaultPressure(), Vec3(1.1, 2.2, 3.3)*bar)
+        self.assertEqual(force.getDefaultTemperature(), 350*kelvin)
+        checkPressureUnits(force)
+
     def testDrudeSCFIntegrator(self):
         """ Tests the DrudeSCFIntegrator API features """
         integrator = DrudeSCFIntegrator(0.002)