Cleanup and bugfixes to MonteCarloAnisotropicBarostat

openmmapi/include/openmm/MonteCarloAnisotropicBarostat.h

@@ -9,8 +9,8 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2010 Stanford University and the Authors.           *
- * Authors: Peter Eastman                                                     *
+ * Portions copyright (c) 2010-2013 Stanford University and the Authors.      *
+ * Authors: Peter Eastman, Lee-Ping Wang                                      *
  * Contributors:                                                              *
  *                                                                            *
  * Permission is hereby granted, free of charge, to any person obtaining a    *
@@ -43,8 +43,11 @@ namespace OpenMM {
  * This class uses a Monte Carlo algorithm to adjust the size of the periodic box, simulating the
  * effect of constant pressure.
  *
- * Compared to MonteCarloBarostat, this class scales the three axes of the simulation cell independently.
- * The user supplies three doubles to specify the pressure along each axis.
+ * This class is similar to MonteCarloBarostat, but each Monte Carlo move is applied to only one axis
+ * of the periodic box (unlike MonteCarloBarostat, which scales the entire box isotropically).  This
+ * means that the box may change shape as well as size over the course of the simulation.  It also
+ * allows you to specify a different pressure for each axis of the box, or to keep the box size fixed
+ * along certain axes while still allowing it to change along others.
  * 
  * This class assumes the simulation is also being run at constant temperature, and requires you
  * to specify the system temperature (since it affects the acceptance probability for Monte Carlo
@@ -84,39 +87,33 @@ public:
      * @param defaultPressure   The default pressure acting on each axis (in bar)
      * @param temperature       the 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 scaleX            on/off switch for whether to scale the X axis
-     * @param scaleY            on/off switch for whether to scale the Y axis
-     * @param scaleZ            on/off switch for whether to scale the Z axis
+     * @param scaleX            whether to allow the X dimension of the periodic box to change size
+     * @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
      */
     MonteCarloAnisotropicBarostat(const Vec3& defaultPressure, double temperature, int frequency = 25, bool scaleX = 1, bool scaleY = 1, bool scaleZ = 1);
     /**
      * Get the default pressure (in bar).
      *
-     * @return the default pressure acting on the system, measured in bar.
+     * @return the default pressure acting along each axis, measured in bar.
      */
-    Vec3 getDefaultPressure() const {
+    const Vec3& getDefaultPressure() const {
         return defaultPressure;
     }
     /**
-     * Get the true/false flag for scaling the X-axis.
-     *
-     * @return the true/false flag for scaling the X-axis.
+     * Get whether to allow the X dimension of the periodic box to change size.
      */
     bool getScaleX() const {
       return scaleX;
     }
     /**
-     * Get the true/false flag for scaling the Y-axis.
-     *
-     * @return the true/false flag for scaling the Y-axis.
+     * Get whether to allow the Y dimension of the periodic box to change size.
      */
     bool getScaleY() const {
       return scaleY;
     }
     /**
-     * Get the true/false flag for scaling the Z-axis.
-     *
-     * @return the true/false flag for scaling the Z-axis.
+     * Get whether to allow the Z dimension of the periodic box to change size.
      */
     bool getScaleZ() const {
       return scaleZ;
@@ -172,15 +169,7 @@ private:
     double temperature;
     bool scaleX, scaleY, scaleZ;
     int frequency, randomNumberSeed;
-
-        double GetTemperature() const {
-            return temperature;
-        }
-
-        void SetTemperature(double temperature) {
-            this->temperature = temperature;
-        }
-    };
+};
 
 } // namespace OpenMM
 

openmmapi/include/openmm/MonteCarloBarostat.h

@@ -123,15 +123,7 @@ protected:
 private:
     double defaultPressure, temperature;
     int frequency, randomNumberSeed;
-
-        double GetTemperature() const {
-            return temperature;
-        }
-
-        void SetTemperature(double temperature) {
-            this->temperature = temperature;
-        }
-    };
+};
 
 } // namespace OpenMM
 

openmmapi/include/openmm/internal/AssertionUtilities.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-2012 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2013 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -54,6 +54,8 @@ void OPENMM_EXPORT throwException(const char* file, int line, const std::string&
 
 #define ASSERT_EQUAL_VEC(expected, found, tol) {double _norm_ = std::sqrt((expected).dot(expected)); double _scale_ = _norm_ > 1.0 ? _norm_ : 1.0; if ((std::abs(((expected)[0])-((found)[0]))/_scale_ > (tol)) || (std::abs(((expected)[1])-((found)[1]))/_scale_ > (tol)) || (std::abs(((expected)[2])-((found)[2]))/_scale_ > (tol))) {std::stringstream details; details << " Expected "<<(expected)<<", found "<<(found); throwException(__FILE__, __LINE__, details.str());}};
 
+#define ASSERT_USUALLY_TRUE(cond) {if (!(cond)) throwException(__FILE__, __LINE__, "(This test is stochastic and may occasionally fail)");};
+
 #define ASSERT_USUALLY_EQUAL_TOL(expected, found, tol) {double _scale_ = std::abs(expected) > 1.0 ? std::abs(expected) : 1.0; if (!(std::abs((expected)-(found))/_scale_ <= (tol))) {std::stringstream details; details << "Expected "<<(expected)<<", found "<<(found)<<" (This test is stochastic and may occasionally fail)"; throwException(__FILE__, __LINE__, details.str());}};
 
 #define ASSERT_VALID_INDEX(index, vector) {if (index < 0 || index >= (int) vector.size()) throwException(__FILE__, __LINE__, "Index out of range");};

openmmapi/src/MonteCarloAnisotropicBarostatImpl.cpp

@@ -6,8 +6,8 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2010 Stanford University and the Authors.           *
- * Authors: Peter Eastman                                                     *
+ * Portions copyright (c) 2010-2013 Stanford University and the Authors.      *
+ * Authors: Peter Eastman, Lee-Ping Wang                                      *
  * Contributors:                                                              *
  *                                                                            *
  * Permission is hereby granted, free of charge, to any person obtaining a    *
@@ -65,7 +65,7 @@ void MonteCarloAnisotropicBarostatImpl::initialize(ContextImpl& context) {
 void MonteCarloAnisotropicBarostatImpl::updateContextState(ContextImpl& context) {
     if (++step < owner.getFrequency() || owner.getFrequency() == 0)
         return;
-    if (owner.getScaleX() == 0 && owner.getScaleY() == 0 && owner.getScaleZ() == 0)
+    if (!owner.getScaleX() && !owner.getScaleY() && !owner.getScaleZ())
         return;
     step = 0;
     
@@ -75,23 +75,26 @@ void MonteCarloAnisotropicBarostatImpl::updateContextState(ContextImpl& context)
     double pressure;
     
     // Choose which axis to modify at random.
-    double rnd = genrand_real2(random)*3.0;
     int axis;
-    while (1) {
-        if (rnd < 1.0 && owner.getScaleX()) {
+    while (true) {
+        double rnd = genrand_real2(random)*3.0;
+        if (rnd < 1.0) {
+            if (owner.getScaleX()) {
                 axis = 0;
                 pressure = context.getParameter(MonteCarloAnisotropicBarostat::PressureX())*(AVOGADRO*1e-25);
                 break;
-        } else if (rnd < 2.0 && owner.getScaleY()) {
+            }
+        } else if (rnd < 2.0) {
+            if (owner.getScaleY()) {
                 axis = 1;
                 pressure = context.getParameter(MonteCarloAnisotropicBarostat::PressureY())*(AVOGADRO*1e-25);
                 break;
+            }
         } else if (owner.getScaleZ()) {
             axis = 2;
             pressure = context.getParameter(MonteCarloAnisotropicBarostat::PressureZ())*(AVOGADRO*1e-25);
             break;
         }
-	rnd = genrand_real2(random)*3.0;
     }
     
     // Modify the periodic box size.
@@ -101,9 +104,7 @@ void MonteCarloAnisotropicBarostatImpl::updateContextState(ContextImpl& context)
     double volume = box[0][0]*box[1][1]*box[2][2];
     double deltaVolume = volumeScale[axis]*2*(genrand_real2(random)-0.5);
     double newVolume = volume+deltaVolume;
-    Vec3 lengthScale;
-    for (int i=0; i<3; i++)
-        lengthScale[i] = 1.0;
+    Vec3 lengthScale(1.0, 1.0, 1.0);
     lengthScale[axis] = newVolume/volume;
     kernel.getAs<ApplyMonteCarloBarostatKernel>().scaleCoordinates(context, lengthScale[0], lengthScale[1], lengthScale[2]);
     context.getOwner().setPeriodicBoxVectors(box[0]*lengthScale[0], box[1]*lengthScale[1], box[2]*lengthScale[2]);

platforms/cuda/tests/TestCudaMonteCarloAnisotropicBarostat.cpp

@@ -6,8 +6,8 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2012 Stanford University and the Authors.      *
- * Authors: Peter Eastman                                                     *
+ * Portions copyright (c) 2008-2013 Stanford University and the Authors.      *
+ * Authors: Peter Eastman, Lee-Ping Wang                                      *
  * Contributors:                                                              *
  *                                                                            *
  * Permission is hereby granted, free of charge, to any person obtaining a    *
@@ -274,42 +274,35 @@ void testRandomSeed() {
     }
 }
 
-void testEinsteinCrystal() {
-  /* 
-     
-     Run a constant pressure simulation on an anisotropic Einstein crystal
-     using isotropic and anisotropic barostats.  There are a total of 15 simulations:
-
-     1) 3 pressures: 9.0, 10.0, 11.0 bar, for each of the following groups:
-     2) 3 groups of simulations that scale just one axis: x, y, z
-     3) 1 group of simulations that scales all three axes in the anisotropic barostat
-     4) 1 group of simulations that scales all three axes in the isotropic barostat
-
-     Results that we will check:
-
-     a) In each group of simulations, the volume should decrease with increasing pressure
-     b) In the three simulation groups that scale just one axis, the compressibility (i.e. incremental volume change
-     with increasing pressure) should go like kx > ky > kz (because the spring constant is largest in the z-direction)
-     c) The anisotropic barostat should produce the same result as the isotropic barostat when all three axes are scaled
-     
+/**
+ * Run a constant pressure simulation on an anisotropic Einstein crystal
+ * using isotropic and anisotropic barostats.  There are a total of 15 simulations:
+ *
+ * 1) 3 pressures: 9.0, 10.0, 11.0 bar, for each of the following groups:
+ * 2) 3 groups of simulations that scale just one axis: x, y, z
+ * 3) 1 group of simulations that scales all three axes in the anisotropic barostat
+ * 4) 1 group of simulations that scales all three axes in the isotropic barostat
+ *
+ * Results that we will check:
+ *
+ * a) In each group of simulations, the volume should decrease with increasing pressure
+ * b) In the three simulation groups that scale just one axis, the compressibility (i.e. incremental volume change
+ * with increasing pressure) should go like kx > ky > kz (because the spring constant is largest in the z-direction)
+ * c) The anisotropic barostat should produce the same result as the isotropic barostat when all three axes are scaled
  */
-
+void testEinsteinCrystal() {
     const int numParticles = 64;
-  const int frequency = 10;
+    const int frequency = 2;
     const int equil = 10000;
-  const int steps = 10000;
+    const int steps = 5000;
     const double pressure = 10.0;
     const double pressureInMD = pressure*(AVOGADRO*1e-25); // pressure in kJ/mol/nm^3
     const double temp = 300.0; // Only test one temperature since we're looking at three pressures.
-  const double pres3[] = {9.0, 10.0, 11.0};
+    const double pres3[] = {2.0, 8.0, 15.0};
     const double initialVolume = numParticles*BOLTZ*temp/pressureInMD;
     const double initialLength = std::pow(initialVolume, 1.0/3.0);
     vector<double> initialPositions(3);
-  // All results.
-  double results[] = {0.0, 0.0, 0.0, 0.0, 0.0, 
-		      0.0, 0.0, 0.0, 0.0, 0.0,
-		      0.0, 0.0, 0.0, 0.0, 0.0};
-  int simNum = 0;
+    vector<double> results;
     // Run four groups of anisotropic simulations; scaling just x, y, z, then all three.
     for (int a = 0; a < 4; a++) {
         // Test barostat for three different pressures.
@@ -325,6 +318,8 @@ void testEinsteinCrystal() {
             force->addPerParticleParameter("x0");
             force->addPerParticleParameter("y0");
             force->addPerParticleParameter("z0");
+            NonbondedForce* nb = new NonbondedForce();
+            nb->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
             for (int i = 0; i < numParticles; ++i) {
                 system.addParticle(1.0);
                 positions[i] = Vec3(((i/16)%4+0.5)*initialLength/4, ((i/4)%4+0.5)*initialLength/4, (i%4+0.5)*initialLength/4);
@@ -332,8 +327,10 @@ void testEinsteinCrystal() {
                 initialPositions[1] = positions[i][1];
                 initialPositions[2] = positions[i][2];
                 force->addParticle(i, initialPositions);
+                nb->addParticle(0, initialLength/6, 0.1);
             }
             system.addForce(force);
+            system.addForce(nb);
             // Create the barostat.
             MonteCarloAnisotropicBarostat* barostat = new MonteCarloAnisotropicBarostat(Vec3(pres3[p], pres3[p], pres3[p]), temp, frequency, (a==0||a==3), (a==1||a==3), (a==2||a==3));
             system.addForce(barostat);
@@ -352,8 +349,7 @@ void testEinsteinCrystal() {
                 integrator.step(frequency);
             }
             volume /= steps;
-      results[simNum] = volume;
-      simNum += 1;
+            results.push_back(volume);
         }
     }
     for (int p = 0; p < 3; p++) {
@@ -368,6 +364,8 @@ void testEinsteinCrystal() {
         force->addPerParticleParameter("x0");
         force->addPerParticleParameter("y0");
         force->addPerParticleParameter("z0");
+        NonbondedForce* nb = new NonbondedForce();
+        nb->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
         for (int i = 0; i < numParticles; ++i) {
             system.addParticle(1.0);
             positions[i] = Vec3(((i/16)%4+0.5)*initialLength/4, ((i/4)%4+0.5)*initialLength/4, (i%4+0.5)*initialLength/4);
@@ -375,13 +373,15 @@ void testEinsteinCrystal() {
             initialPositions[1] = positions[i][1];
             initialPositions[2] = positions[i][2];
             force->addParticle(i, initialPositions);
+            nb->addParticle(0, initialLength/6, 0.1);
         }
         system.addForce(force);
+        system.addForce(nb);
         // Create the barostat.
         MonteCarloBarostat* barostat = new MonteCarloBarostat(pres3[p], temp, frequency);
         system.addForce(barostat);
         barostat->setTemperature(temp);
-    LangevinIntegrator integrator(temp, 0.1, 0.01);
+        LangevinIntegrator integrator(temp, 0.1, 0.001);
         Context context(system, integrator, platform);
         context.setPositions(positions);
         // Let it equilibrate.
@@ -395,29 +395,22 @@ void testEinsteinCrystal() {
             integrator.step(frequency);
         }
         volume /= steps;
-    results[simNum] = volume;
-    simNum += 1;
+        results.push_back(volume);
     }
     
-  /*
-  for (int j = 0; j < 15; j++) {
-    printf("%.6f\n",results[j]);
-  }
-  */
-
     // Check to see if volumes decrease with increasing pressure.
-  ASSERT(results[0] > results[1]);
-  ASSERT(results[1] > results[2]);
-  ASSERT(results[3] > results[4]);
-  ASSERT(results[4] > results[5]);
-  ASSERT(results[6] > results[7]);
-  ASSERT(results[7] > results[8]);
+    ASSERT_USUALLY_TRUE(results[0] > results[1]);
+    ASSERT_USUALLY_TRUE(results[1] > results[2]);
+    ASSERT_USUALLY_TRUE(results[3] > results[4]);
+    ASSERT_USUALLY_TRUE(results[4] > results[5]);
+    ASSERT_USUALLY_TRUE(results[6] > results[7]);
+    ASSERT_USUALLY_TRUE(results[7] > results[8]);
 
     // Check to see if incremental volume changes with increasing pressure go like kx > ky > kz.
-  ASSERT((results[0] - results[1]) > (results[3] - results[4]));
-  ASSERT((results[1] - results[2]) > (results[4] - results[5]));
-  ASSERT((results[3] - results[4]) > (results[6] - results[7]));
-  ASSERT((results[4] - results[5]) > (results[7] - results[8]));
+    ASSERT_USUALLY_TRUE((results[0] - results[1]) > (results[3] - results[4]));
+    ASSERT_USUALLY_TRUE((results[1] - results[2]) > (results[4] - results[5]));
+    ASSERT_USUALLY_TRUE((results[3] - results[4]) > (results[6] - results[7]));
+    ASSERT_USUALLY_TRUE((results[4] - results[5]) > (results[7] - results[8]));
     
     // Check to see if the volumes are equal for isotropic and anisotropic (all axis).
     ASSERT_USUALLY_EQUAL_TOL(results[9], results[12], 3/std::sqrt((double) steps));

platforms/opencl/tests/TestOpenCLMonteCarloAnisotropicBarostat.cpp

@@ -6,8 +6,8 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2012 Stanford University and the Authors.      *
- * Authors: Peter Eastman                                                     *
+ * Portions copyright (c) 2008-2013 Stanford University and the Authors.      *
+ * Authors: Peter Eastman, Lee-Ping Wang                                      *
  * Contributors:                                                              *
  *                                                                            *
  * Permission is hereby granted, free of charge, to any person obtaining a    *
@@ -51,7 +51,7 @@
 using namespace OpenMM;
 using namespace std;
 
-static OpenCLPlatform platform;
+OpenCLPlatform platform;
 
 void testChangingBoxSize() {
     System system;
@@ -274,42 +274,35 @@ void testRandomSeed() {
     }
 }
 
-void testEinsteinCrystal() {
-  /* 
-     
-     Run a constant pressure simulation on an anisotropic Einstein crystal
-     using isotropic and anisotropic barostats.  There are a total of 15 simulations:
-
-     1) 3 pressures: 9.0, 10.0, 11.0 bar, for each of the following groups:
-     2) 3 groups of simulations that scale just one axis: x, y, z
-     3) 1 group of simulations that scales all three axes in the anisotropic barostat
-     4) 1 group of simulations that scales all three axes in the isotropic barostat
-
-     Results that we will check:
-
-     a) In each group of simulations, the volume should decrease with increasing pressure
-     b) In the three simulation groups that scale just one axis, the compressibility (i.e. incremental volume change
-     with increasing pressure) should go like kx > ky > kz (because the spring constant is largest in the z-direction)
-     c) The anisotropic barostat should produce the same result as the isotropic barostat when all three axes are scaled
-     
+/**
+ * Run a constant pressure simulation on an anisotropic Einstein crystal
+ * using isotropic and anisotropic barostats.  There are a total of 15 simulations:
+ *
+ * 1) 3 pressures: 9.0, 10.0, 11.0 bar, for each of the following groups:
+ * 2) 3 groups of simulations that scale just one axis: x, y, z
+ * 3) 1 group of simulations that scales all three axes in the anisotropic barostat
+ * 4) 1 group of simulations that scales all three axes in the isotropic barostat
+ *
+ * Results that we will check:
+ *
+ * a) In each group of simulations, the volume should decrease with increasing pressure
+ * b) In the three simulation groups that scale just one axis, the compressibility (i.e. incremental volume change
+ * with increasing pressure) should go like kx > ky > kz (because the spring constant is largest in the z-direction)
+ * c) The anisotropic barostat should produce the same result as the isotropic barostat when all three axes are scaled
  */
-
+void testEinsteinCrystal() {
     const int numParticles = 64;
-  const int frequency = 10;
+    const int frequency = 2;
     const int equil = 10000;
-  const int steps = 10000;
+    const int steps = 5000;
     const double pressure = 10.0;
     const double pressureInMD = pressure*(AVOGADRO*1e-25); // pressure in kJ/mol/nm^3
     const double temp = 300.0; // Only test one temperature since we're looking at three pressures.
-  const double pres3[] = {9.0, 10.0, 11.0};
+    const double pres3[] = {2.0, 8.0, 15.0};
     const double initialVolume = numParticles*BOLTZ*temp/pressureInMD;
     const double initialLength = std::pow(initialVolume, 1.0/3.0);
     vector<double> initialPositions(3);
-  // All results.
-  double results[] = {0.0, 0.0, 0.0, 0.0, 0.0, 
-		      0.0, 0.0, 0.0, 0.0, 0.0,
-		      0.0, 0.0, 0.0, 0.0, 0.0};
-  int simNum = 0;
+    vector<double> results;
     // Run four groups of anisotropic simulations; scaling just x, y, z, then all three.
     for (int a = 0; a < 4; a++) {
         // Test barostat for three different pressures.
@@ -325,6 +318,8 @@ void testEinsteinCrystal() {
             force->addPerParticleParameter("x0");
             force->addPerParticleParameter("y0");
             force->addPerParticleParameter("z0");
+            NonbondedForce* nb = new NonbondedForce();
+            nb->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
             for (int i = 0; i < numParticles; ++i) {
                 system.addParticle(1.0);
                 positions[i] = Vec3(((i/16)%4+0.5)*initialLength/4, ((i/4)%4+0.5)*initialLength/4, (i%4+0.5)*initialLength/4);
@@ -332,8 +327,10 @@ void testEinsteinCrystal() {
                 initialPositions[1] = positions[i][1];
                 initialPositions[2] = positions[i][2];
                 force->addParticle(i, initialPositions);
+                nb->addParticle(0, initialLength/6, 0.1);
             }
             system.addForce(force);
+            system.addForce(nb);
             // Create the barostat.
             MonteCarloAnisotropicBarostat* barostat = new MonteCarloAnisotropicBarostat(Vec3(pres3[p], pres3[p], pres3[p]), temp, frequency, (a==0||a==3), (a==1||a==3), (a==2||a==3));
             system.addForce(barostat);
@@ -352,8 +349,7 @@ void testEinsteinCrystal() {
                 integrator.step(frequency);
             }
             volume /= steps;
-      results[simNum] = volume;
-      simNum += 1;
+            results.push_back(volume);
         }
     }
     for (int p = 0; p < 3; p++) {
@@ -368,6 +364,8 @@ void testEinsteinCrystal() {
         force->addPerParticleParameter("x0");
         force->addPerParticleParameter("y0");
         force->addPerParticleParameter("z0");
+        NonbondedForce* nb = new NonbondedForce();
+        nb->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
         for (int i = 0; i < numParticles; ++i) {
             system.addParticle(1.0);
             positions[i] = Vec3(((i/16)%4+0.5)*initialLength/4, ((i/4)%4+0.5)*initialLength/4, (i%4+0.5)*initialLength/4);
@@ -375,13 +373,15 @@ void testEinsteinCrystal() {
             initialPositions[1] = positions[i][1];
             initialPositions[2] = positions[i][2];
             force->addParticle(i, initialPositions);
+            nb->addParticle(0, initialLength/6, 0.1);
         }
         system.addForce(force);
+        system.addForce(nb);
         // Create the barostat.
         MonteCarloBarostat* barostat = new MonteCarloBarostat(pres3[p], temp, frequency);
         system.addForce(barostat);
         barostat->setTemperature(temp);
-    LangevinIntegrator integrator(temp, 0.1, 0.01);
+        LangevinIntegrator integrator(temp, 0.1, 0.001);
         Context context(system, integrator, platform);
         context.setPositions(positions);
         // Let it equilibrate.
@@ -395,29 +395,22 @@ void testEinsteinCrystal() {
             integrator.step(frequency);
         }
         volume /= steps;
-    results[simNum] = volume;
-    simNum += 1;
+        results.push_back(volume);
     }
     
-  /*
-  for (int j = 0; j < 15; j++) {
-    printf("%.6f\n",results[j]);
-  }
-  */
-
     // Check to see if volumes decrease with increasing pressure.
-  ASSERT(results[0] > results[1]);
-  ASSERT(results[1] > results[2]);
-  ASSERT(results[3] > results[4]);
-  ASSERT(results[4] > results[5]);
-  ASSERT(results[6] > results[7]);
-  ASSERT(results[7] > results[8]);
+    ASSERT_USUALLY_TRUE(results[0] > results[1]);
+    ASSERT_USUALLY_TRUE(results[1] > results[2]);
+    ASSERT_USUALLY_TRUE(results[3] > results[4]);
+    ASSERT_USUALLY_TRUE(results[4] > results[5]);
+    ASSERT_USUALLY_TRUE(results[6] > results[7]);
+    ASSERT_USUALLY_TRUE(results[7] > results[8]);
 
     // Check to see if incremental volume changes with increasing pressure go like kx > ky > kz.
-  ASSERT((results[0] - results[1]) > (results[3] - results[4]));
-  ASSERT((results[1] - results[2]) > (results[4] - results[5]));
-  ASSERT((results[3] - results[4]) > (results[6] - results[7]));
-  ASSERT((results[4] - results[5]) > (results[7] - results[8]));
+    ASSERT_USUALLY_TRUE((results[0] - results[1]) > (results[3] - results[4]));
+    ASSERT_USUALLY_TRUE((results[1] - results[2]) > (results[4] - results[5]));
+    ASSERT_USUALLY_TRUE((results[3] - results[4]) > (results[6] - results[7]));
+    ASSERT_USUALLY_TRUE((results[4] - results[5]) > (results[7] - results[8]));
     
     // Check to see if the volumes are equal for isotropic and anisotropic (all axis).
     ASSERT_USUALLY_EQUAL_TOL(results[9], results[12], 3/std::sqrt((double) steps));

platforms/reference/tests/TestReferenceMonteCarloAnisotropicBarostat.cpp

@@ -6,8 +6,8 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2010 Stanford University and the Authors.      *
- * Authors: Peter Eastman                                                     *
+ * Portions copyright (c) 2008-2013 Stanford University and the Authors.      *
+ * Authors: Peter Eastman, Lee-Ping Wang                                      *
  * Contributors:                                                              *
  *                                                                            *
  * Permission is hereby granted, free of charge, to any person obtaining a    *
@@ -34,6 +34,8 @@
  */
 
 #include "openmm/internal/AssertionUtilities.h"
+#include "openmm/CustomExternalForce.h"
+#include "openmm/MonteCarloBarostat.h"
 #include "openmm/MonteCarloAnisotropicBarostat.h"
 #include "openmm/Context.h"
 #include "ReferencePlatform.h"
@@ -274,6 +276,151 @@ void testRandomSeed() {
     }
 }
 
+/**
+ * Run a constant pressure simulation on an anisotropic Einstein crystal
+ * using isotropic and anisotropic barostats.  There are a total of 15 simulations:
+ *
+ * 1) 3 pressures: 9.0, 10.0, 11.0 bar, for each of the following groups:
+ * 2) 3 groups of simulations that scale just one axis: x, y, z
+ * 3) 1 group of simulations that scales all three axes in the anisotropic barostat
+ * 4) 1 group of simulations that scales all three axes in the isotropic barostat
+ *
+ * Results that we will check:
+ *
+ * a) In each group of simulations, the volume should decrease with increasing pressure
+ * b) In the three simulation groups that scale just one axis, the compressibility (i.e. incremental volume change
+ * with increasing pressure) should go like kx > ky > kz (because the spring constant is largest in the z-direction)
+ * c) The anisotropic barostat should produce the same result as the isotropic barostat when all three axes are scaled
+ */
+void testEinsteinCrystal() {
+    const int numParticles = 64;
+    const int frequency = 2;
+    const int equil = 10000;
+    const int steps = 5000;
+    const double pressure = 10.0;
+    const double pressureInMD = pressure*(AVOGADRO*1e-25); // pressure in kJ/mol/nm^3
+    const double temp = 300.0; // Only test one temperature since we're looking at three pressures.
+    const double pres3[] = {2.0, 8.0, 15.0};
+    const double initialVolume = numParticles*BOLTZ*temp/pressureInMD;
+    const double initialLength = std::pow(initialVolume, 1.0/3.0);
+    ReferencePlatform platform;
+    vector<double> initialPositions(3);
+    vector<double> results;
+    // Run four groups of anisotropic simulations; scaling just x, y, z, then all three.
+    for (int a = 0; a < 4; a++) {
+        // Test barostat for three different pressures.
+        for (int p = 0; p < 3; p++) {
+            // Create a system of noninteracting particles attached by harmonic springs to their initial positions.
+            System system;
+            system.setDefaultPeriodicBoxVectors(Vec3(initialLength, 0, 0), Vec3(0, initialLength, 0), Vec3(0, 0, initialLength));
+            vector<Vec3> positions(numParticles);
+            OpenMM_SFMT::SFMT sfmt;
+            init_gen_rand(0, sfmt);
+            // Anisotropic force constants.
+            CustomExternalForce* force = new CustomExternalForce("0.005*(x-x0)^2 + 0.01*(y-y0)^2 + 0.02*(z-z0)^2");
+            force->addPerParticleParameter("x0");
+            force->addPerParticleParameter("y0");
+            force->addPerParticleParameter("z0");
+            NonbondedForce* nb = new NonbondedForce();
+            nb->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+            for (int i = 0; i < numParticles; ++i) {
+                system.addParticle(1.0);
+                positions[i] = Vec3(((i/16)%4+0.5)*initialLength/4, ((i/4)%4+0.5)*initialLength/4, (i%4+0.5)*initialLength/4);
+                initialPositions[0] = positions[i][0];
+                initialPositions[1] = positions[i][1];
+                initialPositions[2] = positions[i][2];
+                force->addParticle(i, initialPositions);
+                nb->addParticle(0, initialLength/6, 0.1);
+            }
+            system.addForce(force);
+            system.addForce(nb);
+            // Create the barostat.
+            MonteCarloAnisotropicBarostat* barostat = new MonteCarloAnisotropicBarostat(Vec3(pres3[p], pres3[p], pres3[p]), temp, frequency, (a==0||a==3), (a==1||a==3), (a==2||a==3));
+            system.addForce(barostat);
+            barostat->setTemperature(temp);
+            LangevinIntegrator integrator(temp, 0.1, 0.01);
+            Context context(system, integrator, platform);
+            context.setPositions(positions);
+            // Let it equilibrate.
+            integrator.step(equil);
+            // Now run it for a while and see if the volume is correct.
+            double volume = 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);
+            }
+            volume /= steps;
+            results.push_back(volume);
+        }
+    }
+    for (int p = 0; p < 3; p++) {
+        // Create a system of noninteracting particles attached by harmonic springs to their initial positions.
+        System system;
+        system.setDefaultPeriodicBoxVectors(Vec3(initialLength, 0, 0), Vec3(0, initialLength, 0), Vec3(0, 0, initialLength));
+        vector<Vec3> positions(numParticles);
+        OpenMM_SFMT::SFMT sfmt;
+        init_gen_rand(0, sfmt);
+        // Anisotropic force constants.
+        CustomExternalForce* force = new CustomExternalForce("0.005*(x-x0)^2 + 0.01*(y-y0)^2 + 0.02*(z-z0)^2");
+        force->addPerParticleParameter("x0");
+        force->addPerParticleParameter("y0");
+        force->addPerParticleParameter("z0");
+        NonbondedForce* nb = new NonbondedForce();
+        nb->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+        for (int i = 0; i < numParticles; ++i) {
+            system.addParticle(1.0);
+            positions[i] = Vec3(((i/16)%4+0.5)*initialLength/4, ((i/4)%4+0.5)*initialLength/4, (i%4+0.5)*initialLength/4);
+            initialPositions[0] = positions[i][0];
+            initialPositions[1] = positions[i][1];
+            initialPositions[2] = positions[i][2];
+            force->addParticle(i, initialPositions);
+            nb->addParticle(0, initialLength/6, 0.1);
+        }
+        system.addForce(force);
+        system.addForce(nb);
+        // Create the barostat.
+        MonteCarloBarostat* barostat = new MonteCarloBarostat(pres3[p], temp, frequency);
+        system.addForce(barostat);
+        barostat->setTemperature(temp);
+        LangevinIntegrator integrator(temp, 0.1, 0.001);
+        Context context(system, integrator, platform);
+        context.setPositions(positions);
+        // Let it equilibrate.
+        integrator.step(equil);
+        // Now run it for a while and see if the volume is correct.
+        double volume = 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);
+        }
+        volume /= steps;
+        results.push_back(volume);
+    }
+    
+    // Check to see if volumes decrease with increasing pressure.
+    ASSERT_USUALLY_TRUE(results[0] > results[1]);
+    ASSERT_USUALLY_TRUE(results[1] > results[2]);
+    ASSERT_USUALLY_TRUE(results[3] > results[4]);
+    ASSERT_USUALLY_TRUE(results[4] > results[5]);
+    ASSERT_USUALLY_TRUE(results[6] > results[7]);
+    ASSERT_USUALLY_TRUE(results[7] > results[8]);
+
+    // Check to see if incremental volume changes with increasing pressure go like kx > ky > kz.
+    ASSERT_USUALLY_TRUE((results[0] - results[1]) > (results[3] - results[4]));
+    ASSERT_USUALLY_TRUE((results[1] - results[2]) > (results[4] - results[5]));
+    ASSERT_USUALLY_TRUE((results[3] - results[4]) > (results[6] - results[7]));
+    ASSERT_USUALLY_TRUE((results[4] - results[5]) > (results[7] - results[8]));
+    
+    // Check to see if the volumes are equal for isotropic and anisotropic (all axis).
+    ASSERT_USUALLY_EQUAL_TOL(results[9], results[12], 3/std::sqrt((double) steps));
+    ASSERT_USUALLY_EQUAL_TOL(results[10], results[13], 3/std::sqrt((double) steps));
+    ASSERT_USUALLY_EQUAL_TOL(results[11], results[14], 3/std::sqrt((double) steps));
+}
+
 int main() {
     try {
         testIdealGas();
@@ -281,6 +428,7 @@ int main() {
         testIdealGasAxis(1);
         testIdealGasAxis(2);
         testRandomSeed();
+        testEinsteinCrystal();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;