Code implemented and passes tests

olla/include/openmm/kernels.h

@@ -1130,6 +1130,45 @@ public:
     virtual void restoreCoordinates(ContextImpl& context) = 0;
 };
 
+/**
+ * This kernel is invoked by MonteCarloAnisotropicBarostat to adjust the periodic box volume
+ */
+class ApplyMonteCarloAnisotropicBarostatKernel : public KernelImpl {
+public:
+    static std::string Name() {
+        return "ApplyMonteCarloAnisotropicBarostat";
+    }
+    ApplyMonteCarloAnisotropicBarostatKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
+    }
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param barostat   the MonteCarloAnisotropicBarostat this kernel will be used for
+     */
+    virtual void initialize(const System& system, const MonteCarloAnisotropicBarostat& 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 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 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.
+     *
+     * @param context    the context in which to execute this kernel
+     */
+    virtual void restoreCoordinates(ContextImpl& context) = 0;
+};
+
 /**
  * This kernel is invoked to remove center of mass motion from the system.
  */

openmmapi/include/openmm/MonteCarloBarostat.h

@@ -151,16 +151,32 @@ class OPENMM_EXPORT MonteCarloAnisotropicBarostat : public Force {
 public:
     /**
      * This is the name of the parameter which stores the current pressure acting on
-     * the system (in bar).
+     * the X-axis (in bar).
      */
-    static const std::string& Pressure() {
-        static const std::string key = "MonteCarloPressure";
+    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 of the system (in bar)
+     * @param defaultPressure   3-element vector specifying 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
@@ -169,12 +185,52 @@ public:
      */
     MonteCarloAnisotropicBarostat(Vec3 defaultPressure, double temperature, int frequency = 25, bool scaleX = 1, bool scaleY = 1, bool scaleZ = 1);
     /**
-     * Get the default pressure acting on the system (in bar).
+     * Get the default pressure acting on the X-axis (in bar).
      *
      * @return the default pressure acting on the system, measured in bar.
      */
-    Vec3 getDefaultPressure() const {
-        return defaultPressure;
+    double getDefaultPressureX() const {
+        return defaultPressure[0];
+    }
+    /**
+     * Get the default pressure acting on the Y-axis (in bar).
+     *
+     * @return the default pressure acting on the system, measured in bar.
+     */
+    double getDefaultPressureY() const {
+        return defaultPressure[1];
+    }
+    /**
+     * Get the default pressure acting on the Z-axis (in bar).
+     *
+     * @return the default pressure acting on the system, measured in bar.
+     */
+    double getDefaultPressureZ() const {
+        return defaultPressure[2];
+    }
+    /**
+     * Get the true/false flag for scaling the X-axis.
+     *
+     * @return the true/false flag for scaling the X-axis.
+     */
+    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.
+     */
+    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.
+     */
+    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
@@ -225,6 +281,7 @@ protected:
 private:
     Vec3 defaultPressure;
     double temperature;
+    bool scaleX, scaleY, scaleZ;
     int frequency, randomNumberSeed;
 
         double GetTemperature() const {

openmmapi/include/openmm/internal/MonteCarloBarostatImpl.h

@@ -79,7 +79,7 @@ public:
         // This force doesn't apply forces to particles.
         return 0.0;
     }
-    std::map<std::string, Vec3> getDefaultParameters();
+    std::map<std::string, double> getDefaultParameters();
     std::vector<std::string> getKernelNames();
 private:
     const MonteCarloAnisotropicBarostat& owner;

openmmapi/src/MonteCarloBarostatImpl.cpp

@@ -122,12 +122,12 @@ std::vector<std::string> MonteCarloBarostatImpl::getKernelNames() {
     return names;
 }
 
-MonteCarloAnisotropicBarostatImpl::MonteCarloAnisotropicBarostatImpl(const MonteCarloBarostat& owner) : owner(owner), step(0) {
+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);
+    kernel = context.getPlatform().createKernel(ApplyMonteCarloAnisotropicBarostatKernel::Name(), context);
+    kernel.getAs<ApplyMonteCarloAnisotropicBarostatKernel>().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];
@@ -142,26 +142,30 @@ void MonteCarloAnisotropicBarostatImpl::initialize(ContextImpl& context) {
 void MonteCarloAnisotropicBarostatImpl::updateContextState(ContextImpl& context) {
     if (++step < owner.getFrequency() || owner.getFrequency() == 0)
         return;
-    if (scaleX == 0 && scaleY == 0 && scaleZ == 0)
+    if (owner.getScaleX() == 0 && owner.getScaleY() == 0 && owner.getScaleZ() == 0)
         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.
     double rnd = genrand_real2(random)*3.0;
     int axis;
     while (1) {
-        if (rnd < 1.0 && scaleX) {
+        if (rnd < 1.0 && owner.getScaleX()) {
             axis = 0;
+	    pressure = context.getParameter(MonteCarloAnisotropicBarostat::PressureX())*(AVOGADRO*1e-25);
             break;
-        } else if (rnd < 2.0 && scaleY) {
+        } else if (rnd < 2.0 && owner.getScaleY()) {
             axis = 1;
+	    pressure = context.getParameter(MonteCarloAnisotropicBarostat::PressureY())*(AVOGADRO*1e-25);
             break;
-        } else if (scaleZ) {
+        } else if (owner.getScaleZ()) {
             axis = 2;
+	    pressure = context.getParameter(MonteCarloAnisotropicBarostat::PressureZ())*(AVOGADRO*1e-25);
             break;
         }
     }
@@ -177,19 +181,18 @@ void MonteCarloAnisotropicBarostatImpl::updateContextState(ContextImpl& context)
     for (int i=0; i<3; i++)
       lengthScale[i] = 1.0;
     lengthScale[axis] = newVolume/volume;
-    kernel.getAs<ApplyMonteCarloBarostatKernel>().scaleCoordinates(context, lengthScale[0], lengthScale[1], lengthScale[2]);
+    kernel.getAs<ApplyMonteCarloAnisotropicBarostatKernel>().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 pressure = context.getParameter(MonteCarloBarostat::Pressure())[axis]*(AVOGADRO*1e-25);
     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);
+        kernel.getAs<ApplyMonteCarloAnisotropicBarostatKernel>().restoreCoordinates(context);
         context.getOwner().setPeriodicBoxVectors(box[0], box[1], box[2]);
         volume = newVolume;
     }
@@ -210,15 +213,17 @@ void MonteCarloAnisotropicBarostatImpl::updateContextState(ContextImpl& context)
     }
 }
 
-std::map<std::string, double[3]> MonteCarloAnisotropicBarostatImpl::getDefaultParameters() {
-    std::map<std::string, double[3]> parameters;
-    parameters[MonteCarloBarostat::Pressure()] = getOwner().getDefaultPressure();
+std::map<std::string, double> MonteCarloAnisotropicBarostatImpl::getDefaultParameters() {
+    std::map<std::string, double> parameters;
+    parameters[MonteCarloAnisotropicBarostat::PressureX()] = getOwner().getDefaultPressureX();
+    parameters[MonteCarloAnisotropicBarostat::PressureY()] = getOwner().getDefaultPressureY();
+    parameters[MonteCarloAnisotropicBarostat::PressureZ()] = getOwner().getDefaultPressureZ();
     return parameters;
 }
 
 std::vector<std::string> MonteCarloAnisotropicBarostatImpl::getKernelNames() {
     std::vector<std::string> names;
-    names.push_back(ApplyMonteCarloBarostatKernel::Name());
+    names.push_back(ApplyMonteCarloAnisotropicBarostatKernel::Name());
     return names;
 }
 

platforms/cuda/src/CudaKernels.cpp

@@ -5372,6 +5372,13 @@ void CudaApplyMonteCarloBarostatKernel::restoreCoordinates(ContextImpl& context)
     }
 }
 
+void CudaApplyMonteCarloAnisotropicBarostatKernel::initialize(const System& system, const MonteCarloAnisotropicBarostat& 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");
+}
+
 CudaRemoveCMMotionKernel::~CudaRemoveCMMotionKernel() {
     cu.setAsCurrent();
     if (cmMomentum != NULL)

platforms/cuda/src/CudaKernels.h

@@ -1286,6 +1286,29 @@ private:
     std::vector<int> lastAtomOrder;
 };
 
+/**
+ * This kernel is invoked by MonteCarloAnisotropicBarostat to adjust the periodic box volume
+ */
+class CudaApplyMonteCarloAnisotropicBarostatKernel : public CudaApplyMonteCarloBarostatKernel {
+public:
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param barostat   the MonteCarloAnisotropicBarostat this kernel will be used for
+     */
+    void initialize(const System& system, const MonteCarloAnisotropicBarostat& barostat);
+private:
+    CudaContext& cu;
+    bool hasInitializedKernels;
+    int numMolecules;
+    CudaArray* savedPositions;
+    CudaArray* moleculeAtoms;
+    CudaArray* moleculeStartIndex;
+    CUfunction kernel;
+    std::vector<int> lastAtomOrder;
+};
+
 /**
  * This kernel is invoked to remove center of mass motion from the system.
  */

platforms/cuda/tests/TestCudaMonteCarloBarostat.cpp

@@ -110,7 +110,9 @@ void testIdealGas(int aniso) {
         system.addParticle(1.0);
         positions[i] = Vec3(initialLength*genrand_real2(sfmt), 0.5*initialLength*genrand_real2(sfmt), 2*initialLength*genrand_real2(sfmt));
     }
-    MonteCarloBarostat* barostat = new MonteCarloBarostat(pressure, temp[0], frequency, aniso);
+    MonteCarloBarostat* barostat = new MonteCarloBarostat(pressure, temp[0], frequency);
+    if (aniso)
+      MonteCarloAnisotropicBarostat* barostat = new MonteCarloAnisotropicBarostat(Vec3(pressure, pressure, pressure), temp[0], frequency);
     system.addForce(barostat);
 
     // Test it for three different temperatures.
@@ -252,7 +254,9 @@ void testWater(int aniso) {
         }
     }
     system.addForce(nonbonded);
-    MonteCarloBarostat* barostat = new MonteCarloBarostat(pressure, temp, frequency, aniso);
+    MonteCarloBarostat* barostat = new MonteCarloBarostat(pressure, temp, frequency);
+    if (aniso)
+      MonteCarloAnisotropicBarostat* barostat = new MonteCarloAnisotropicBarostat(Vec3(pressure, pressure, pressure), temp, frequency);
     system.addForce(barostat);
 
     // Simulate it and see if the density matches the expected value (1 g/mL).

platforms/opencl/src/OpenCLKernels.h

@@ -1300,6 +1300,29 @@ private:
     std::vector<int> lastAtomOrder;
 };
 
+/**
+ * This kernel is invoked by MonteCarloAnisotropicBarostat to adjust the periodic box volume
+ */
+class OpenCLApplyMonteCarloAnisotropicBarostatKernel : public OpenCLApplyMonteCarloBarostatKernel {
+public:
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param barostat   the MonteCarloAnisotropicBarostat this kernel will be used for
+     */
+    void initialize(const System& system, const MonteCarloAnisotropicBarostat& barostat);
+private:
+    OpenCLContext& cl;
+    bool hasInitializedKernels;
+    int numMolecules;
+    OpenCLArray* savedPositions;
+    OpenCLArray* moleculeAtoms;
+    OpenCLArray* moleculeStartIndex;
+    cl::Kernel kernel;
+    std::vector<int> lastAtomOrder;
+};
+
 /**
  * This kernel is invoked to remove center of mass motion from the system.
  */

platforms/opencl/tests/TestOpenCLMonteCarloBarostat.cpp

@@ -110,7 +110,9 @@ void testIdealGas(int aniso) {
         system.addParticle(1.0);
         positions[i] = Vec3(initialLength*genrand_real2(sfmt), 0.5*initialLength*genrand_real2(sfmt), 2*initialLength*genrand_real2(sfmt));
     }
-    MonteCarloBarostat* barostat = new MonteCarloBarostat(pressure, temp[0], frequency, aniso);
+    MonteCarloBarostat* barostat = new MonteCarloBarostat(pressure, temp[0], frequency);
+    if (aniso)
+      MonteCarloAnisotropicBarostat* barostat = new MonteCarloAnisotropicBarostat(Vec3(pressure, pressure, pressure), temp[0], frequency);
     system.addForce(barostat);
 
     // Test it for three different temperatures.
@@ -255,7 +257,9 @@ void testWater(int aniso) {
         }
     }
     system.addForce(nonbonded);
-    MonteCarloBarostat* barostat = new MonteCarloBarostat(pressure, temp, frequency, aniso);
+    MonteCarloBarostat* barostat = new MonteCarloBarostat(pressure, temp, frequency);
+    if (aniso)
+      MonteCarloAnisotropicBarostat* barostat = new MonteCarloAnisotropicBarostat(Vec3(pressure, pressure, pressure), temp, frequency);
     system.addForce(barostat);
 
     // Simulate it and see if the density matches the expected value (1 g/mL).

platforms/reference/src/ReferenceKernels.cpp

@@ -2156,6 +2156,9 @@ void ReferenceApplyMonteCarloBarostatKernel::restoreCoordinates(ContextImpl& con
     barostat->restorePositions(posData);
 }
 
+void ReferenceApplyMonteCarloAnisotropicBarostatKernel::initialize(const System& system, const MonteCarloAnisotropicBarostat& barostat) {
+}
+
 void ReferenceRemoveCMMotionKernel::initialize(const System& system, const CMMotionRemover& force) {
     frequency = force.getFrequency();
     masses.resize(system.getNumParticles());

platforms/reference/src/ReferenceKernels.h

@@ -1207,6 +1207,22 @@ private:
     ReferenceMonteCarloBarostat* barostat;
 };
 
+/**
+ * This kernel is invoked by MonteCarloAnisotropicBarostat to adjust the periodic box volume
+ */
+class ReferenceApplyMonteCarloAnisotropicBarostatKernel : public ReferenceApplyMonteCarloBarostatKernel{
+public:
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param barostat   the MonteCarloAnisotropicBarostat this kernel will be used for
+     */
+    void initialize(const System& system, const MonteCarloAnisotropicBarostat& barostat);
+private:
+    ReferenceMonteCarloBarostat* barostat;
+};
+
 /**
  * This kernel is invoked to remove center of mass motion from the system.
  */

platforms/reference/tests/TestReferenceMonteCarloBarostat.cpp

@@ -110,7 +110,9 @@ void testIdealGas(int aniso) {
         system.addParticle(1.0);
         positions[i] = Vec3(initialLength*genrand_real2(sfmt), 0.5*initialLength*genrand_real2(sfmt), 2*initialLength*genrand_real2(sfmt));
     }
-    MonteCarloBarostat* barostat = new MonteCarloBarostat(pressure, temp[0], frequency, aniso);
+    MonteCarloBarostat* barostat = new MonteCarloBarostat(pressure, temp[0], frequency);
+    if (aniso)
+      MonteCarloAnisotropicBarostat* barostat = new MonteCarloAnisotropicBarostat(Vec3(pressure, pressure, pressure), temp[0], frequency);
     system.addForce(barostat);
 
     // Test it for three different temperatures.

wrappers/python/src/swig_doxygen/OpenMM.i

@@ -81,6 +81,7 @@ using namespace OpenMM;
     %template(XmlSerializer_serialize_HarmonicAngleForce) XmlSerializer::serialize<HarmonicAngleForce>;
     %template(XmlSerializer_serialize_HarmonicBondForce) XmlSerializer::serialize<HarmonicBondForce>;
     %template(XmlSerializer_serialize_MonteCarloBarostat) XmlSerializer::serialize<MonteCarloBarostat>;
+    %template(XmlSerializer_serialize_MonteCarloAnisotropicBarostat) XmlSerializer::serialize<MonteCarloAnisotropicBarostat>;
     %template(XmlSerializer_serialize_NonbondedForce) XmlSerializer::serialize<NonbondedForce>;
     %template(XmlSerializer_serialize_RBTorsionForce) XmlSerializer::serialize<RBTorsionForce>;
     %template(XmlSerializer_serialize_System) XmlSerializer::serialize<System>;
@@ -101,6 +102,7 @@ using namespace OpenMM;
     %template(XmlSerializer_deserialize_HarmonicAngleForce) XmlSerializer::deserialize<HarmonicAngleForce>;
     %template(XmlSerializer_deserialize_HarmonicBondForce) XmlSerializer::deserialize<HarmonicBondForce>;
     %template(XmlSerializer_deserialize_MonteCarloBarostat) XmlSerializer::deserialize<MonteCarloBarostat>;
+    %template(XmlSerializer_deserialize_MonteCarloAnisotropicBarostat) XmlSerializer::deserialize<MonteCarloAnisotropicBarostat>;
     %template(XmlSerializer_deserialize_NonbondedForce) XmlSerializer::deserialize<NonbondedForce>;
     %template(XmlSerializer_deserialize_RBTorsionForce) XmlSerializer::deserialize<RBTorsionForce>;
     %template(XmlSerializer_deserialize_System) XmlSerializer::deserialize<System>;

wrappers/python/src/swig_doxygen/swigInputConfig.py

@@ -32,6 +32,7 @@ SKIP_METHODS = [('State',),
                 ('ApplyAndersenThermostatKernel',),
                 ('ApplyConstraintsKernel',),
                 ('ApplyMonteCarloBarostatKernel',),
+                ('ApplyMonteCarloAnisotropicBarostatKernel',),
                 ('BondInfo',),
                 ('BondParameterInfo',),
                 ('CalcAmoebaGeneralizedKirkwoodForceKernel',),
@@ -98,6 +99,7 @@ SKIP_METHODS = [('State',),
                 ('KernelFactory',),
                 ('KernelImpl',),
                 ('MonteCarloBarostatImpl',),
+                ('MonteCarloAnisotropicBarostatImpl',),
                 ('MultipoleInfo',),
                 ('NonbondedForceImpl',),
                 ('OutOfPlaneBendInfo',),
@@ -181,6 +183,9 @@ UNITS = {
 ("*", "getDefaultPeriodicBoxVectors")
  : (None, ('unit.nanometer', 'unit.nanometer', 'unit.nanometer')),
 ("*", "getDefaultPressure") : ("unit.bar", ()),
+("*", "getDefaultPressureX") : ("unit.bar", ()),
+("*", "getDefaultPressureY") : ("unit.bar", ()),
+("*", "getDefaultPressureZ") : ("unit.bar", ()),
 ("*", "getDefaultTemperature") : ("unit.kelvin", ()),
 ("*", "getErrorTolerance") : (None, ()),
 ("*", "getEwaldErrorTolerance") : (None, ()),
@@ -421,6 +426,7 @@ UNITS = {
  : (None, (None, None, 'unit.nanometer',
            'unit.kilojoule_per_mole/(unit.nanometer*unit.nanometer)')),
 ("MonteCarloBarostat", "getFrequency") : (None, ()),
+("MonteCarloAnisotropicBarostat", "getFrequency") : (None, ()),
 ("NonbondedForce", "getExceptionParameters")
  : (None, (None, None,
            'unit.elementary_charge*unit.elementary_charge',