Merge pull request #12 from leeping/master

Created MonteCarloAnisotropicBarostat

olla/include/openmm/kernels.h

@@ -1107,17 +1107,20 @@ public:
      * @param system     the System this kernel will be applied to
      * @param barostat   the MonteCarloBarostat this kernel will be used for
      */
-    virtual void initialize(const System& system, const MonteCarloBarostat& barostat) = 0;
+    virtual void initialize(const System& system, const Force& barostat) = 0;
     /**
      * Attempt a Monte Carlo step, scaling particle positions (or cluster centers) by a specified value.
+     * This version scales the x, y, and z positions independently.
      * This is called BEFORE the periodic box size is modified.  It should begin by translating each particle
      * or cluster into the first periodic box, so that coordinates will still be correct after the box size
      * is changed.
      *
      * @param context    the context in which to execute this kernel
-     * @param scale      the scale factor by which to multiply particle positions
+     * @param scaleX     the scale factor by which to multiply particle x-coordinate
+     * @param scaleY     the scale factor by which to multiply particle y-coordinate
+     * @param scaleZ     the scale factor by which to multiply particle z-coordinate
      */
-    virtual void scaleCoordinates(ContextImpl& context, double scale) = 0;
+    virtual void scaleCoordinates(ContextImpl& context, double scaleX, double scaleY, double scaleZ) = 0;
     /**
      * Reject the most recent Monte Carlo step, restoring the particle positions to where they were before
      * scaleCoordinates() was last called.

openmmapi/include/OpenMM.h

@@ -53,6 +53,7 @@
 #include "openmm/Integrator.h"
 #include "openmm/LangevinIntegrator.h"
 #include "openmm/LocalEnergyMinimizer.h"
+#include "openmm/MonteCarloAnisotropicBarostat.h"
 #include "openmm/MonteCarloBarostat.h"
 #include "openmm/NonbondedForce.h"
 #include "openmm/Context.h"

openmmapi/include/openmm/MonteCarloAnisotropicBarostat.h

+#ifndef OPENMM_MONTECARLOANISOTROPICBAROSTAT_H_
+#define OPENMM_MONTECARLOANISOTROPICBAROSTAT_H_
+
+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * 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    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "Force.h"
+#include "Vec3.h"
+#include <string>
+#include "internal/windowsExport.h"
+
+namespace OpenMM {
+
+/**
+ * This class uses a Monte Carlo algorithm to adjust the size of the periodic box, simulating the
+ * effect of constant pressure.
+ *
+ * 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
+ * moves).  It does not actually perform temperature regulation, however.  You must use another
+ * mechanism along with it to maintain the temperature, such as LangevinIntegrator or AndersenThermostat.
+ */
+
+class OPENMM_EXPORT MonteCarloAnisotropicBarostat : public Force {
+public:
+    /**
+     * This is the name of the parameter which stores the current pressure acting on
+     * the X-axis (in bar).
+     */
+    static const std::string& PressureX() {
+        static const std::string key = "MonteCarloPressureX";
+        return key;
+    }
+    /**
+     * This is the name of the parameter which stores the current pressure acting on
+     * the Y-axis (in bar).
+     */
+    static const std::string& PressureY() {
+        static const std::string key = "MonteCarloPressureY";
+        return key;
+    }
+    /**
+     * This is the name of the parameter which stores the current pressure acting on
+     * the Z-axis (in bar).
+     */
+    static const std::string& PressureZ() {
+        static const std::string key = "MonteCarloPressureZ";
+        return key;
+    }
+    /**
+     * Create a MonteCarloAnisotropicBarostat.
+     *
+     * @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            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 along each axis, measured in bar.
+     */
+    const Vec3& getDefaultPressure() const {
+        return defaultPressure;
+    }
+    /**
+     * Get whether to allow the X dimension of the periodic box to change size.
+     */
+    bool getScaleX() const {
+      return scaleX;
+    }
+    /**
+     * Get whether to allow the Y dimension of the periodic box to change size.
+     */
+    bool getScaleY() const {
+      return scaleY;
+    }
+    /**
+     * Get whether to allow the Z dimension of the periodic box to change size.
+     */
+    bool getScaleZ() const {
+      return scaleZ;
+    }
+    /**
+     * Get the frequency (in time steps) at which Monte Carlo pressure changes should be attempted.  If this is set to
+     * 0, the barostat is disabled.
+     */
+    int getFrequency() const {
+        return frequency;
+    }
+    /**
+     * Set the frequency (in time steps) at which Monte Carlo pressure changes should be attempted.  If this is set to
+     * 0, the barostat is disabled.
+     */
+    void setFrequency(int freq) {
+        frequency = freq;
+    }
+    /**
+     * Get the temperature at which the system is being maintained, measured in Kelvin.
+     */
+    double getTemperature() const {
+        return temperature;
+    }
+    /**
+     * Set the temperature at which the system is being maintained.
+     *
+     * @param temp     the system temperature, measured in Kelvin.
+     */
+    void setTemperature(double temp) {
+        temperature = temp;
+    }
+    /**
+     * Get the random number seed.  See setRandomNumberSeed() for details.
+     */
+    int getRandomNumberSeed() const {
+        return randomNumberSeed;
+    }
+    /**
+     * Set the random number seed.  It is guaranteed that if two simulations are run
+     * with different random number seeds, the sequence of Monte Carlo steps will be different.  On
+     * the other hand, no guarantees are made about the behavior of simulations that use the same seed.
+     * In particular, Platforms are permitted to use non-deterministic algorithms which produce different
+     * results on successive runs, even if those runs were initialized identically.
+     */
+    void setRandomNumberSeed(int seed) {
+        randomNumberSeed = seed;
+    }
+protected:
+    ForceImpl* createImpl() const;
+private:
+    Vec3 defaultPressure;
+    double temperature;
+    bool scaleX, scaleY, scaleZ;
+    int frequency, randomNumberSeed;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_MONTECARLOANISOTROPICBAROSTAT_H_*/

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/include/openmm/internal/MonteCarloAnisotropicBarostatImpl.h

+#ifndef OPENMM_MONTECARLOANISOTROPICBAROSTATIMPL_H_
+#define OPENMM_MONTECARLOANISOTROPICBAROSTATIMPL_H_
+
+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * 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                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "ForceImpl.h"
+#include "openmm/MonteCarloAnisotropicBarostat.h"
+#include "openmm/Kernel.h"
+#include "sfmt/SFMT.h"
+#include <string>
+
+namespace OpenMM {
+
+/**
+ * This is the internal implementation of MonteCarloAnisotropicBarostat.
+ */
+
+class MonteCarloAnisotropicBarostatImpl : public ForceImpl {
+public:
+    MonteCarloAnisotropicBarostatImpl(const MonteCarloAnisotropicBarostat& owner);
+    void initialize(ContextImpl& context);
+    const MonteCarloAnisotropicBarostat& getOwner() const {
+        return owner;
+    }
+    void updateContextState(ContextImpl& context);
+    double calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy, int groups) {
+        // This force doesn't apply forces to particles.
+        return 0.0;
+    }
+    std::map<std::string, double> getDefaultParameters();
+    std::vector<std::string> getKernelNames();
+private:
+    const MonteCarloAnisotropicBarostat& owner;
+    int step, numAttempted[3], numAccepted[3];
+    double volumeScale[3];
+    OpenMM_SFMT::SFMT random;
+    Kernel kernel;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_MONTECARLOANISOTROPICBAROSTATIMPL_H_*/

openmmapi/src/MonteCarloAnisotropicBarostat.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * 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                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "openmm/MonteCarloAnisotropicBarostat.h"
+#include "openmm/internal/MonteCarloAnisotropicBarostatImpl.h"
+#include <ctime>
+
+using namespace OpenMM;
+
+MonteCarloAnisotropicBarostat::MonteCarloAnisotropicBarostat(const Vec3& defaultPressure, double temperature, int frequency, bool scaleX, bool scaleY, bool scaleZ) :
+        defaultPressure(defaultPressure), temperature(temperature), frequency(frequency), scaleX(scaleX), scaleY(scaleY), scaleZ(scaleZ) {
+    setRandomNumberSeed((int) time(NULL));
+}
+
+ForceImpl* MonteCarloAnisotropicBarostat::createImpl() const {
+    return new MonteCarloAnisotropicBarostatImpl(*this);
+}

openmmapi/src/MonteCarloAnisotropicBarostatImpl.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * 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    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "openmm/internal/MonteCarloAnisotropicBarostatImpl.h"
+#include "openmm/internal/ContextImpl.h"
+#include "openmm/Context.h"
+#include "openmm/kernels.h"
+#include <cmath>
+#include <vector>
+
+using namespace OpenMM;
+using namespace OpenMM_SFMT;
+using std::vector;
+
+const float BOLTZMANN = 1.380658e-23f; // (J/K)
+const float AVOGADRO = 6.0221367e23f;
+const float RGAS = BOLTZMANN*AVOGADRO; // (J/(mol K))
+const float BOLTZ = RGAS/1000;         // (kJ/(mol K))
+
+MonteCarloAnisotropicBarostatImpl::MonteCarloAnisotropicBarostatImpl(const MonteCarloAnisotropicBarostat& owner) : owner(owner), step(0) {
+}
+
+void MonteCarloAnisotropicBarostatImpl::initialize(ContextImpl& context) {
+    kernel = context.getPlatform().createKernel(ApplyMonteCarloBarostatKernel::Name(), context);
+    kernel.getAs<ApplyMonteCarloBarostatKernel>().initialize(context.getSystem(), owner);
+    Vec3 box[3];
+    context.getPeriodicBoxVectors(box[0], box[1], box[2]);
+    double volume = box[0][0]*box[1][1]*box[2][2];
+    for (int i=0; i<3; i++) {
+        volumeScale[i] = 0.01*volume;
+        numAttempted[i] = 0;
+        numAccepted[i] = 0;
+    }
+    init_gen_rand(owner.getRandomNumberSeed(), random);
+}
+
+void MonteCarloAnisotropicBarostatImpl::updateContextState(ContextImpl& context) {
+    if (++step < owner.getFrequency() || owner.getFrequency() == 0)
+        return;
+    if (!owner.getScaleX() && !owner.getScaleY() && !owner.getScaleZ())
+        return;
+    step = 0;
+    
+    // Compute the current potential energy.
+    
+    double initialEnergy = context.getOwner().getState(State::Energy).getPotentialEnergy();
+    double pressure;
+    
+    // Choose which axis to modify at random.
+    int axis;
+    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) {
+            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;
+        }
+    }
+    
+    // Modify the periodic box size.
+    
+    Vec3 box[3];
+    context.getPeriodicBoxVectors(box[0], box[1], box[2]);
+    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(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]);
+    
+    // Compute the energy of the modified system.
+    
+    double finalEnergy = context.getOwner().getState(State::Energy).getPotentialEnergy();
+    double kT = BOLTZ*owner.getTemperature();
+    double w = finalEnergy-initialEnergy + pressure*deltaVolume - context.getMolecules().size()*kT*std::log(newVolume/volume);
+    if (w > 0 && genrand_real2(random) > std::exp(-w/kT)) {
+        // Reject the step.
+        
+        kernel.getAs<ApplyMonteCarloBarostatKernel>().restoreCoordinates(context);
+        context.getOwner().setPeriodicBoxVectors(box[0], box[1], box[2]);
+        volume = newVolume;
+    }
+    else
+        numAccepted[axis]++;
+    numAttempted[axis]++;
+    if (numAttempted[axis] >= 10) {
+        if (numAccepted[axis] < 0.25*numAttempted[axis]) {
+            volumeScale[axis] /= 1.1;
+            numAttempted[axis] = 0;
+            numAccepted[axis] = 0;
+        }
+        else if (numAccepted[axis] > 0.75*numAttempted[axis]) {
+            volumeScale[axis] = std::min(volumeScale[axis]*1.1, volume*0.3);
+            numAttempted[axis] = 0;
+            numAccepted[axis] = 0;
+        }
+    }
+}
+
+std::map<std::string, double> MonteCarloAnisotropicBarostatImpl::getDefaultParameters() {
+    std::map<std::string, double> parameters;
+    parameters[MonteCarloAnisotropicBarostat::PressureX()] = getOwner().getDefaultPressure()[0];
+    parameters[MonteCarloAnisotropicBarostat::PressureY()] = getOwner().getDefaultPressure()[1];
+    parameters[MonteCarloAnisotropicBarostat::PressureZ()] = getOwner().getDefaultPressure()[2];
+    return parameters;
+}
+
+std::vector<std::string> MonteCarloAnisotropicBarostatImpl::getKernelNames() {
+    std::vector<std::string> names;
+    names.push_back(ApplyMonteCarloBarostatKernel::Name());
+    return names;
+}
+

openmmapi/src/MonteCarloBarostatImpl.cpp

@@ -77,7 +77,7 @@ void MonteCarloBarostatImpl::updateContextState(ContextImpl& context) {
     double deltaVolume = volumeScale*2*(genrand_real2(random)-0.5);
     double newVolume = volume+deltaVolume;
     double lengthScale = std::pow(newVolume/volume, 1.0/3.0);
-    kernel.getAs<ApplyMonteCarloBarostatKernel>().scaleCoordinates(context, lengthScale);
+    kernel.getAs<ApplyMonteCarloBarostatKernel>().scaleCoordinates(context, lengthScale, lengthScale, lengthScale);
     context.getOwner().setPeriodicBoxVectors(box[0]*lengthScale, box[1]*lengthScale, box[2]*lengthScale);
 
     // Compute the energy of the modified system.

platforms/cuda/include/CudaKernels.h

@@ -1259,17 +1259,20 @@ public:
      * @param system     the System this kernel will be applied to
      * @param barostat   the MonteCarloBarostat this kernel will be used for
      */
-    void initialize(const System& system, const MonteCarloBarostat& barostat);
+    void initialize(const System& system, const Force& barostat);
     /**
      * Attempt a Monte Carlo step, scaling particle positions (or cluster centers) by a specified value.
+     * This version scales the x, y, and z positions independently.
      * This is called BEFORE the periodic box size is modified.  It should begin by translating each particle
      * or cluster into the first periodic box, so that coordinates will still be correct after the box size
      * is changed.
      *
      * @param context    the context in which to execute this kernel
-     * @param scale      the scale factor by which to multiply particle positions
+     * @param scaleX     the scale factor by which to multiply particle x-coordinate
+     * @param scaleY     the scale factor by which to multiply particle y-coordinate
+     * @param scaleZ     the scale factor by which to multiply particle z-coordinate
      */
-    void scaleCoordinates(ContextImpl& context, double scale);
+    void scaleCoordinates(ContextImpl& context, double scaleX, double scaleY, double scaleZ);
     /**
      * Reject the most recent Monte Carlo step, restoring the particle positions to where they were before
      * scaleCoordinates() was last called.

platforms/cuda/src/CudaKernels.cpp

@@ -5402,14 +5402,14 @@ CudaApplyMonteCarloBarostatKernel::~CudaApplyMonteCarloBarostatKernel() {
         delete moleculeStartIndex;
 }
 
-void CudaApplyMonteCarloBarostatKernel::initialize(const System& system, const MonteCarloBarostat& thermostat) {
+void CudaApplyMonteCarloBarostatKernel::initialize(const System& system, const Force& thermostat) {
     cu.setAsCurrent();
     savedPositions = new CudaArray(cu, cu.getPaddedNumAtoms(), cu.getUseDoublePrecision() ? sizeof(double4) : sizeof(float4), "savedPositions");
     CUmodule module = cu.createModule(CudaKernelSources::monteCarloBarostat);
     kernel = cu.getKernel(module, "scalePositions");
 }
 
-void CudaApplyMonteCarloBarostatKernel::scaleCoordinates(ContextImpl& context, double scale) {
+void CudaApplyMonteCarloBarostatKernel::scaleCoordinates(ContextImpl& context, double scaleX, double scaleY, double scaleZ) {
     cu.setAsCurrent();
     if (!hasInitializedKernels) {
         hasInitializedKernels = true;
@@ -5442,8 +5442,10 @@ void CudaApplyMonteCarloBarostatKernel::scaleCoordinates(ContextImpl& context, d
         m<<"Error saving positions for MC barostat: "<<cu.getErrorString(result)<<" ("<<result<<")";
         throw OpenMMException(m.str());
     }        
-    float scalef = (float) scale;
-    void* args[] = {&scalef, &numMolecules, cu.getPeriodicBoxSizePointer(), cu.getInvPeriodicBoxSizePointer(),
+    float scalefX = (float) scaleX;
+    float scalefY = (float) scaleY;
+    float scalefZ = (float) scaleZ;
+    void* args[] = {&scalefX, &scalefY, &scalefZ, &numMolecules, cu.getPeriodicBoxSizePointer(), cu.getInvPeriodicBoxSizePointer(), 
 		    &cu.getPosq().getDevicePointer(), &moleculeAtoms->getDevicePointer(), &moleculeStartIndex->getDevicePointer()};
     cu.executeKernel(kernel, args, cu.getNumAtoms());
     for (int i = 0; i < (int) cu.getPosCellOffsets().size(); i++)

platforms/cuda/src/kernels/monteCarloBarostat.cu

 /**
- * Scale the particle positions.
+ * Scale the particle positions with each axis independent
  */
 
-extern "C" __global__ void scalePositions(float scale, int numMolecules, real4 periodicBoxSize, real4 invPeriodicBoxSize, real4* __restrict__ posq,
+extern "C" __global__ void scalePositions(float scaleX, float scaleY, float scaleZ, int numMolecules, real4 periodicBoxSize, real4 invPeriodicBoxSize, real4* __restrict__ posq,
         const int* __restrict__ moleculeAtoms, const int* __restrict__ moleculeStartIndex) {
     for (int index = blockIdx.x*blockDim.x+threadIdx.x; index < numMolecules; index += blockDim.x*gridDim.x) {
         int first = moleculeStartIndex[index];
@@ -35,9 +35,9 @@ extern "C" __global__ void scalePositions(float scale, int numMolecules, real4 p
 
         // Now scale the position of the molecule center.
 
-        delta.x = center.x*(scale-1)-delta.x;
-        delta.y = center.y*(scale-1)-delta.y;
-        delta.z = center.z*(scale-1)-delta.z;
+        delta.x = center.x*(scaleX-1)-delta.x;
+        delta.y = center.y*(scaleY-1)-delta.y;
+        delta.z = center.z*(scaleZ-1)-delta.z;
         for (int atom = first; atom < last; atom++) {
             real4 pos = posq[moleculeAtoms[atom]];
             pos.x += delta.x;

platforms/cuda/tests/TestCudaMonteCarloAnisotropicBarostat.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * 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    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests the CUDA implementation of MonteCarloAnisotropicBarostat.
+ */
+
+#include "openmm/internal/AssertionUtilities.h"
+#include "openmm/CustomExternalForce.h"
+#include "openmm/MonteCarloBarostat.h"
+#include "openmm/MonteCarloAnisotropicBarostat.h"
+#include "openmm/Context.h"
+#include "CudaPlatform.h"
+#include "openmm/NonbondedForce.h"
+#include "openmm/System.h"
+#include "openmm/LangevinIntegrator.h"
+#include "openmm/VerletIntegrator.h"
+#include "sfmt/SFMT.h"
+#include "SimTKOpenMMRealType.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+CudaPlatform platform;
+
+void testChangingBoxSize() {
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(4, 0, 0), Vec3(0, 5, 0), Vec3(0, 0, 6));
+    system.addParticle(1.0);
+    NonbondedForce* nb = new NonbondedForce();
+    nb->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    nb->setCutoffDistance(2.0);
+    nb->addParticle(1, 0.5, 0.5);
+    system.addForce(nb);
+    LangevinIntegrator integrator(300.0, 1.0, 0.01);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions;
+    positions.push_back(Vec3());
+    context.setPositions(positions);
+    Vec3 x, y, z;
+    context.getState(State::Forces).getPeriodicBoxVectors(x, y, z);
+    ASSERT_EQUAL_VEC(Vec3(4, 0, 0), x, 0);
+    ASSERT_EQUAL_VEC(Vec3(0, 5, 0), y, 0);
+    ASSERT_EQUAL_VEC(Vec3(0, 0, 6), z, 0);
+    context.setPeriodicBoxVectors(Vec3(7, 0, 0), Vec3(0, 8, 0), Vec3(0, 0, 9));
+    context.getState(State::Forces).getPeriodicBoxVectors(x, y, z);
+    ASSERT_EQUAL_VEC(Vec3(7, 0, 0), x, 0);
+    ASSERT_EQUAL_VEC(Vec3(0, 8, 0), y, 0);
+    ASSERT_EQUAL_VEC(Vec3(0, 0, 9), z, 0);
+    
+    // Shrinking the box too small should produce an exception.
+    
+    context.setPeriodicBoxVectors(Vec3(7, 0, 0), Vec3(0, 3.9, 0), Vec3(0, 0, 9));
+    bool ok = true;
+    try {
+        context.getState(State::Forces).getPeriodicBoxVectors(x, y, z);
+        ok = false;
+    }
+    catch (exception& ex) {
+    }
+    ASSERT(ok);
+}
+
+void testIdealGas() {
+    const int numParticles = 64;
+    const int frequency = 10;
+    const int steps = 1000;
+    const double pressure = 1.5;
+    const double pressureInMD = pressure*(AVOGADRO*1e-25);
+    const double temp[] = {300.0, 600.0, 1000.0};
+    const double initialVolume = numParticles*BOLTZ*temp[1]/pressureInMD;
+    const double initialLength = std::pow(initialVolume, 1.0/3.0);
+    
+    // Create a gas of noninteracting particles.
+    
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(initialLength, 0, 0), Vec3(0, 0.5*initialLength, 0), Vec3(0, 0, 2*initialLength));
+    vector<Vec3> positions(numParticles);
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    for (int i = 0; i < numParticles; ++i) {
+        system.addParticle(1.0);
+        positions[i] = Vec3(initialLength*genrand_real2(sfmt), 0.5*initialLength*genrand_real2(sfmt), 2*initialLength*genrand_real2(sfmt));
+    }
+    MonteCarloAnisotropicBarostat* barostat = new MonteCarloAnisotropicBarostat(Vec3(pressure, pressure, pressure), temp[0], frequency);
+    system.addForce(barostat);
+    
+    // Test it for three different temperatures.
+    
+    for (int i = 0; i < 3; i++) {
+        barostat->setTemperature(temp[i]);
+        LangevinIntegrator integrator(temp[i], 0.1, 0.01);
+        Context context(system, integrator, platform);
+        context.setPositions(positions);
+        
+        // Let it equilibrate.
+        
+        integrator.step(10000);
+        
+        // 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;
+        double expected = (numParticles+1)*BOLTZ*temp[i]/pressureInMD;
+        ASSERT_USUALLY_EQUAL_TOL(expected, volume, 3/std::sqrt((double) steps));
+    }
+}
+
+void testIdealGasAxis(int axis) {
+    // Test scaling just one axis.
+    const int numParticles = 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, 600.0, 1000.0};
+    const double initialVolume = numParticles*BOLTZ*temp[1]/pressureInMD;
+    const double initialLength = std::pow(initialVolume, 1.0/3.0);
+    const bool scaleX = (axis == 0);
+    const bool scaleY = (axis == 1);
+    const bool scaleZ = (axis == 2);
+    double boxX;
+    double boxY;
+    double boxZ;
+    
+    // Create a gas of noninteracting particles.
+    
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(initialLength, 0, 0), Vec3(0, 0.5*initialLength, 0), Vec3(0, 0, 2*initialLength));
+    vector<Vec3> positions(numParticles);
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    for (int i = 0; i < numParticles; ++i) {
+        system.addParticle(1.0);
+        positions[i] = Vec3(initialLength*genrand_real2(sfmt), 0.5*initialLength*genrand_real2(sfmt), 2*initialLength*genrand_real2(sfmt));
+    }
+    MonteCarloAnisotropicBarostat* barostat = new MonteCarloAnisotropicBarostat(Vec3(pressure, pressure, pressure), temp[0], frequency, scaleX, scaleY, scaleZ);
+    system.addForce(barostat);
+    
+    // Test it for three different temperatures.
+    
+    for (int i = 0; i < 3; i++) {
+        barostat->setTemperature(temp[i]);
+        LangevinIntegrator integrator(temp[i], 0.1, 0.01);
+        Context context(system, integrator, platform);
+        context.setPositions(positions);
+        
+        // Let it equilibrate.
+        
+        integrator.step(10000);
+        
+        // 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]);
+            boxX = box[0][0];
+            boxY = box[1][1];
+            boxZ = box[2][2];
+            volume += box[0][0]*box[1][1]*box[2][2];
+            integrator.step(frequency);
+        }
+        volume /= steps;
+        double expected = (numParticles+1)*BOLTZ*temp[i]/pressureInMD;
+        ASSERT_USUALLY_EQUAL_TOL(expected, volume, 3/std::sqrt((double) steps));
+        if (!scaleX) {
+            ASSERT(boxX == initialLength);
+        }
+        if (!scaleY) {
+            ASSERT(boxY == 0.5*initialLength);
+        }
+        if (!scaleZ) {
+            ASSERT(boxZ == 2*initialLength);
+        }
+    }
+}
+
+void testRandomSeed() {
+    const int numParticles = 8;
+    const double temp = 100.0;
+    const double pressure = 1.5;
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(8, 0, 0), Vec3(0, 8, 0), Vec3(0, 0, 8));
+    VerletIntegrator integrator(0.01);
+    NonbondedForce* forceField = new NonbondedForce();
+    forceField->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    for (int i = 0; i < numParticles; ++i) {
+        system.addParticle(2.0);
+        forceField->addParticle((i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
+    }
+    system.addForce(forceField);
+    MonteCarloAnisotropicBarostat* barostat = new MonteCarloAnisotropicBarostat(Vec3(pressure, pressure, pressure), temp, 1);
+    system.addForce(barostat);
+    vector<Vec3> positions(numParticles);
+    vector<Vec3> velocities(numParticles);
+    for (int i = 0; i < numParticles; ++i) {
+        positions[i] = Vec3((i%2 == 0 ? 2 : -2), (i%4 < 2 ? 2 : -2), (i < 4 ? 2 : -2));
+        velocities[i] = Vec3(0, 0, 0);
+    }
+    
+    // Try twice with the same random seed.
+    
+    barostat->setRandomNumberSeed(5);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+    integrator.step(10);
+    State state1 = context.getState(State::Positions);
+    context.reinitialize();
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+    integrator.step(10);
+    State state2 = context.getState(State::Positions);
+    
+    // Try twice with a different random seed.
+    
+    barostat->setRandomNumberSeed(10);
+    context.reinitialize();
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+    integrator.step(10);
+    State state3 = context.getState(State::Positions);
+    context.reinitialize();
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+    integrator.step(10);
+    State state4 = context.getState(State::Positions);
+    
+    // Compare the results.
+    
+    for (int i = 0; i < numParticles; i++) {
+        for (int j = 0; j < 3; j++) {
+            ASSERT(state1.getPositions()[i][j] == state2.getPositions()[i][j]);
+            ASSERT(state3.getPositions()[i][j] == state4.getPositions()[i][j]);
+            ASSERT(state1.getPositions()[i][j] != state3.getPositions()[i][j]);
+        }
+    }
+}
+
+/**
+ * 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);
+    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(int argc, char* argv[]) {
+    try {
+        if (argc > 1)
+            platform.setPropertyDefaultValue("CudaPrecision", string(argv[1]));
+        testChangingBoxSize();
+        testIdealGas();
+        testIdealGasAxis(0);
+        testIdealGasAxis(1);
+        testIdealGasAxis(2);
+        testRandomSeed();
+        testEinsteinCrystal();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}

platforms/opencl/include/OpenCLKernels.h

@@ -1273,17 +1273,20 @@ public:
      * @param system     the System this kernel will be applied to
      * @param barostat   the MonteCarloBarostat this kernel will be used for
      */
-    void initialize(const System& system, const MonteCarloBarostat& barostat);
+    void initialize(const System& system, const Force& barostat);
     /**
      * Attempt a Monte Carlo step, scaling particle positions (or cluster centers) by a specified value.
+     * This version scales the x, y, and z positions independently.
      * This is called BEFORE the periodic box size is modified.  It should begin by translating each particle
      * or cluster into the first periodic box, so that coordinates will still be correct after the box size
      * is changed.
      *
      * @param context    the context in which to execute this kernel
-     * @param scale      the scale factor by which to multiply particle positions
+     * @param scaleX     the scale factor by which to multiply particle x-coordinate
+     * @param scaleY     the scale factor by which to multiply particle y-coordinate
+     * @param scaleZ     the scale factor by which to multiply particle z-coordinate
      */
-    void scaleCoordinates(ContextImpl& context, double scale);
+    void scaleCoordinates(ContextImpl& context, double scaleX, double scaleY, double scaleZ);
     /**
      * Reject the most recent Monte Carlo step, restoring the particle positions to where they were before
      * scaleCoordinates() was last called.

platforms/opencl/src/OpenCLKernels.cpp

@@ -5626,13 +5626,13 @@ OpenCLApplyMonteCarloBarostatKernel::~OpenCLApplyMonteCarloBarostatKernel() {
         delete moleculeStartIndex;
 }
 
-void OpenCLApplyMonteCarloBarostatKernel::initialize(const System& system, const MonteCarloBarostat& thermostat) {
+void OpenCLApplyMonteCarloBarostatKernel::initialize(const System& system, const Force& thermostat) {
     savedPositions = OpenCLArray::create<mm_float4>(cl, cl.getPaddedNumAtoms(), "savedPositions");
     cl::Program program = cl.createProgram(OpenCLKernelSources::monteCarloBarostat);
     kernel = cl::Kernel(program, "scalePositions");
 }
 
-void OpenCLApplyMonteCarloBarostatKernel::scaleCoordinates(ContextImpl& context, double scale) {
+void OpenCLApplyMonteCarloBarostatKernel::scaleCoordinates(ContextImpl& context, double scaleX, double scaleY, double scaleZ) {
     if (!hasInitializedKernels) {
         hasInitializedKernels = true;
 
@@ -5656,15 +5656,17 @@ void OpenCLApplyMonteCarloBarostatKernel::scaleCoordinates(ContextImpl& context,
 
         // Initialize the kernel arguments.
         
-        kernel.setArg<cl_int>(1, numMolecules);
-        kernel.setArg<cl::Buffer>(4, cl.getPosq().getDeviceBuffer());
-        kernel.setArg<cl::Buffer>(5, moleculeAtoms->getDeviceBuffer());
-        kernel.setArg<cl::Buffer>(6, moleculeStartIndex->getDeviceBuffer());
+        kernel.setArg<cl_int>(3, numMolecules);
+        kernel.setArg<cl::Buffer>(6, cl.getPosq().getDeviceBuffer());
+        kernel.setArg<cl::Buffer>(7, moleculeAtoms->getDeviceBuffer());
+        kernel.setArg<cl::Buffer>(8, moleculeStartIndex->getDeviceBuffer());
     }
     cl.getQueue().enqueueCopyBuffer(cl.getPosq().getDeviceBuffer(), savedPositions->getDeviceBuffer(), 0, 0, cl.getPosq().getSize()*sizeof(mm_float4));
-    kernel.setArg<cl_float>(0, (cl_float) scale);
-    setPeriodicBoxSizeArg(cl, kernel, 2);
-    setInvPeriodicBoxSizeArg(cl, kernel, 3);
+    kernel.setArg<cl_float>(0, (cl_float) scaleX);
+    kernel.setArg<cl_float>(1, (cl_float) scaleY);
+    kernel.setArg<cl_float>(2, (cl_float) scaleZ);
+    setPeriodicBoxSizeArg(cl, kernel, 4);
+    setInvPeriodicBoxSizeArg(cl, kernel, 5);
     cl.executeKernel(kernel, cl.getNumAtoms());
     for (int i = 0; i < (int) cl.getPosCellOffsets().size(); i++)
         cl.getPosCellOffsets()[i] = mm_int4(0, 0, 0, 0);

platforms/opencl/src/kernels/monteCarloBarostat.cl

 /**
- * Scale the particle positions.
+ * Scale the particle positions with each axis independent.
  */
 
-__kernel void scalePositions(float scale, int numMolecules, real4 periodicBoxSize, real4 invPeriodicBoxSize, __global real4* restrict posq,
+__kernel void scalePositions(float scaleX, float scaleY, float scaleZ, int numMolecules, real4 periodicBoxSize, real4 invPeriodicBoxSize, __global real4* restrict posq,
         __global const int* restrict moleculeAtoms, __global const int* restrict moleculeStartIndex) {
     for (int index = get_global_id(0); index < numMolecules; index += get_global_size(0)) {
         int first = moleculeStartIndex[index];
@@ -22,11 +22,12 @@ __kernel void scalePositions(float scale, int numMolecules, real4 periodicBoxSiz
         int ycell = (int) floor(center.y*invPeriodicBoxSize.y);
         int zcell = (int) floor(center.z*invPeriodicBoxSize.z);
         real4 delta = (real4) (xcell*periodicBoxSize.x, ycell*periodicBoxSize.y, zcell*periodicBoxSize.z, 0);
+        real4 scaleXYZ = (real4) (scaleX, scaleY, scaleZ, 1);
         center -= delta;
 
         // Now scale the position of the molecule center.
 
-        delta = center*(scale-1)-delta;
+        delta = center*(scaleXYZ-1)-delta;
         for (int atom = first; atom < last; atom++) {
             real4 pos = posq[moleculeAtoms[atom]];
             pos.xyz += delta.xyz;

platforms/opencl/tests/TestOpenCLMonteCarloAnisotropicBarostat.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * 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    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests the OpenCL implementation of MonteCarloAnisotropicBarostat.
+ */
+
+#include "openmm/internal/AssertionUtilities.h"
+#include "openmm/CustomExternalForce.h"
+#include "openmm/MonteCarloBarostat.h"
+#include "openmm/MonteCarloAnisotropicBarostat.h"
+#include "openmm/Context.h"
+#include "OpenCLPlatform.h"
+#include "openmm/NonbondedForce.h"
+#include "openmm/System.h"
+#include "openmm/LangevinIntegrator.h"
+#include "openmm/VerletIntegrator.h"
+#include "sfmt/SFMT.h"
+#include "SimTKOpenMMRealType.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+OpenCLPlatform platform;
+
+void testChangingBoxSize() {
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(4, 0, 0), Vec3(0, 5, 0), Vec3(0, 0, 6));
+    system.addParticle(1.0);
+    NonbondedForce* nb = new NonbondedForce();
+    nb->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    nb->setCutoffDistance(2.0);
+    nb->addParticle(1, 0.5, 0.5);
+    system.addForce(nb);
+    LangevinIntegrator integrator(300.0, 1.0, 0.01);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions;
+    positions.push_back(Vec3());
+    context.setPositions(positions);
+    Vec3 x, y, z;
+    context.getState(State::Forces).getPeriodicBoxVectors(x, y, z);
+    ASSERT_EQUAL_VEC(Vec3(4, 0, 0), x, 0);
+    ASSERT_EQUAL_VEC(Vec3(0, 5, 0), y, 0);
+    ASSERT_EQUAL_VEC(Vec3(0, 0, 6), z, 0);
+    context.setPeriodicBoxVectors(Vec3(7, 0, 0), Vec3(0, 8, 0), Vec3(0, 0, 9));
+    context.getState(State::Forces).getPeriodicBoxVectors(x, y, z);
+    ASSERT_EQUAL_VEC(Vec3(7, 0, 0), x, 0);
+    ASSERT_EQUAL_VEC(Vec3(0, 8, 0), y, 0);
+    ASSERT_EQUAL_VEC(Vec3(0, 0, 9), z, 0);
+    
+    // Shrinking the box too small should produce an exception.
+    
+    context.setPeriodicBoxVectors(Vec3(7, 0, 0), Vec3(0, 3.9, 0), Vec3(0, 0, 9));
+    bool ok = true;
+    try {
+        context.getState(State::Forces).getPeriodicBoxVectors(x, y, z);
+        ok = false;
+    }
+    catch (exception& ex) {
+    }
+    ASSERT(ok);
+}
+
+void testIdealGas() {
+    const int numParticles = 64;
+    const int frequency = 10;
+    const int steps = 1000;
+    const double pressure = 1.5;
+    const double pressureInMD = pressure*(AVOGADRO*1e-25);
+    const double temp[] = {300.0, 600.0, 1000.0};
+    const double initialVolume = numParticles*BOLTZ*temp[1]/pressureInMD;
+    const double initialLength = std::pow(initialVolume, 1.0/3.0);
+    
+    // Create a gas of noninteracting particles.
+    
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(initialLength, 0, 0), Vec3(0, 0.5*initialLength, 0), Vec3(0, 0, 2*initialLength));
+    vector<Vec3> positions(numParticles);
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    for (int i = 0; i < numParticles; ++i) {
+        system.addParticle(1.0);
+        positions[i] = Vec3(initialLength*genrand_real2(sfmt), 0.5*initialLength*genrand_real2(sfmt), 2*initialLength*genrand_real2(sfmt));
+    }
+    MonteCarloAnisotropicBarostat* barostat = new MonteCarloAnisotropicBarostat(Vec3(pressure, pressure, pressure), temp[0], frequency);
+    system.addForce(barostat);
+    
+    // Test it for three different temperatures.
+    
+    for (int i = 0; i < 3; i++) {
+        barostat->setTemperature(temp[i]);
+        LangevinIntegrator integrator(temp[i], 0.1, 0.01);
+        Context context(system, integrator, platform);
+        context.setPositions(positions);
+        
+        // Let it equilibrate.
+        
+        integrator.step(10000);
+        
+        // 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;
+        double expected = (numParticles+1)*BOLTZ*temp[i]/pressureInMD;
+        ASSERT_USUALLY_EQUAL_TOL(expected, volume, 3/std::sqrt((double) steps));
+    }
+}
+
+void testIdealGasAxis(int axis) {
+    // Test scaling just one axis.
+    const int numParticles = 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, 600.0, 1000.0};
+    const double initialVolume = numParticles*BOLTZ*temp[1]/pressureInMD;
+    const double initialLength = std::pow(initialVolume, 1.0/3.0);
+    const bool scaleX = (axis == 0);
+    const bool scaleY = (axis == 1);
+    const bool scaleZ = (axis == 2);
+    double boxX;
+    double boxY;
+    double boxZ;
+    
+    // Create a gas of noninteracting particles.
+    
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(initialLength, 0, 0), Vec3(0, 0.5*initialLength, 0), Vec3(0, 0, 2*initialLength));
+    vector<Vec3> positions(numParticles);
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    for (int i = 0; i < numParticles; ++i) {
+        system.addParticle(1.0);
+        positions[i] = Vec3(initialLength*genrand_real2(sfmt), 0.5*initialLength*genrand_real2(sfmt), 2*initialLength*genrand_real2(sfmt));
+    }
+    MonteCarloAnisotropicBarostat* barostat = new MonteCarloAnisotropicBarostat(Vec3(pressure, pressure, pressure), temp[0], frequency, scaleX, scaleY, scaleZ);
+    system.addForce(barostat);
+    
+    // Test it for three different temperatures.
+    
+    for (int i = 0; i < 3; i++) {
+        barostat->setTemperature(temp[i]);
+        LangevinIntegrator integrator(temp[i], 0.1, 0.01);
+        Context context(system, integrator, platform);
+        context.setPositions(positions);
+        
+        // Let it equilibrate.
+        
+        integrator.step(10000);
+        
+        // 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]);
+            boxX = box[0][0];
+            boxY = box[1][1];
+            boxZ = box[2][2];
+            volume += box[0][0]*box[1][1]*box[2][2];
+            integrator.step(frequency);
+        }
+        volume /= steps;
+        double expected = (numParticles+1)*BOLTZ*temp[i]/pressureInMD;
+        ASSERT_USUALLY_EQUAL_TOL(expected, volume, 3/std::sqrt((double) steps));
+        if (!scaleX) {
+            ASSERT(boxX == initialLength);
+        }
+        if (!scaleY) {
+            ASSERT(boxY == 0.5*initialLength);
+        }
+        if (!scaleZ) {
+            ASSERT(boxZ == 2*initialLength);
+        }
+    }
+}
+
+void testRandomSeed() {
+    const int numParticles = 8;
+    const double temp = 100.0;
+    const double pressure = 1.5;
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(8, 0, 0), Vec3(0, 8, 0), Vec3(0, 0, 8));
+    VerletIntegrator integrator(0.01);
+    NonbondedForce* forceField = new NonbondedForce();
+    forceField->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    for (int i = 0; i < numParticles; ++i) {
+        system.addParticle(2.0);
+        forceField->addParticle((i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
+    }
+    system.addForce(forceField);
+    MonteCarloAnisotropicBarostat* barostat = new MonteCarloAnisotropicBarostat(Vec3(pressure, pressure, pressure), temp, 1);
+    system.addForce(barostat);
+    vector<Vec3> positions(numParticles);
+    vector<Vec3> velocities(numParticles);
+    for (int i = 0; i < numParticles; ++i) {
+        positions[i] = Vec3((i%2 == 0 ? 2 : -2), (i%4 < 2 ? 2 : -2), (i < 4 ? 2 : -2));
+        velocities[i] = Vec3(0, 0, 0);
+    }
+    
+    // Try twice with the same random seed.
+    
+    barostat->setRandomNumberSeed(5);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+    integrator.step(10);
+    State state1 = context.getState(State::Positions);
+    context.reinitialize();
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+    integrator.step(10);
+    State state2 = context.getState(State::Positions);
+    
+    // Try twice with a different random seed.
+    
+    barostat->setRandomNumberSeed(10);
+    context.reinitialize();
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+    integrator.step(10);
+    State state3 = context.getState(State::Positions);
+    context.reinitialize();
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+    integrator.step(10);
+    State state4 = context.getState(State::Positions);
+    
+    // Compare the results.
+    
+    for (int i = 0; i < numParticles; i++) {
+        for (int j = 0; j < 3; j++) {
+            ASSERT(state1.getPositions()[i][j] == state2.getPositions()[i][j]);
+            ASSERT(state3.getPositions()[i][j] == state4.getPositions()[i][j]);
+            ASSERT(state1.getPositions()[i][j] != state3.getPositions()[i][j]);
+        }
+    }
+}
+
+/**
+ * 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);
+    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(int argc, char* argv[]) {
+    try {
+        if (argc > 1)
+            platform.setPropertyDefaultValue("OpenCLPrecision", string(argv[1]));
+        testChangingBoxSize();
+        testIdealGas();
+        testIdealGasAxis(0);
+        testIdealGasAxis(1);
+        testIdealGasAxis(2);
+        testRandomSeed();
+        testEinsteinCrystal();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}

platforms/reference/include/ReferenceKernels.h

@@ -1182,17 +1182,20 @@ public:
      * @param system     the System this kernel will be applied to
      * @param barostat   the MonteCarloBarostat this kernel will be used for
      */
-    void initialize(const System& system, const MonteCarloBarostat& barostat);
+    void initialize(const System& system, const Force& barostat);
     /**
      * Attempt a Monte Carlo step, scaling particle positions (or cluster centers) by a specified value.
+     * This version scales the x, y, and z positions independently.
      * This is called BEFORE the periodic box size is modified.  It should begin by translating each particle
      * or cluster into the first periodic box, so that coordinates will still be correct after the box size
      * is changed.
      *
      * @param context    the context in which to execute this kernel
-     * @param scale      the scale factor by which to multiply particle positions
+     * @param scaleX     the scale factor by which to multiply particle x-coordinate
+     * @param scaleY     the scale factor by which to multiply particle y-coordinate
+     * @param scaleZ     the scale factor by which to multiply particle z-coordinate
      */
-    void scaleCoordinates(ContextImpl& context, double scale);
+    void scaleCoordinates(ContextImpl& context, double scaleX, double scaleY, double scaleZ);
     /**
      * Reject the most recent Monte Carlo step, restoring the particle positions to where they were before
      * scaleCoordinates() was last called.

platforms/reference/include/ReferenceMonteCarloBarostat.h

@@ -58,15 +58,17 @@ class ReferenceMonteCarloBarostat {
 
       /**---------------------------------------------------------------------------------------
 
-         Apply the barostat at the start of a time step.
+         Apply the barostat at the start of a time step, scaling x, y, and z coordinates independently.
 
          @param atomPositions      atom positions
          @param boxSize            the periodic box dimensions
-         @param scale              the factor by which to scale atom positions
+         @param scaleX             the factor by which to scale atomic x coordinates
+         @param scaleY             the factor by which to scale atomic y coordinates
+         @param scaleZ             the factor by which to scale atomic z coordinates
 
          --------------------------------------------------------------------------------------- */
 
-      void applyBarostat(std::vector<OpenMM::RealVec>& atomPositions, const OpenMM::RealVec& boxSize, RealOpenMM scale);
+      void applyBarostat(std::vector<OpenMM::RealVec>& atomPositions, const OpenMM::RealVec& boxSize, RealOpenMM scaleX, RealOpenMM scaleY, RealOpenMM scaleZ);
 
       /**---------------------------------------------------------------------------------------