Currently in the middle of implementing new version but not finished yet.

olla/include/openmm/kernels.h

@@ -1108,16 +1108,6 @@ public:
      * @param barostat   the MonteCarloBarostat this kernel will be used for
      */
     virtual void initialize(const System& system, const MonteCarloBarostat& barostat) = 0;
-    /**
-     * Attempt a Monte Carlo step, scaling particle positions (or cluster centers) by a specified value.
-     * 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
-     */
-    virtual void scaleCoordinates(ContextImpl& context, double scale) = 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.
@@ -1130,7 +1120,7 @@ public:
      * @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 scaleCoordinatesXYZ(ContextImpl& context, double scaleX, double scaleY, double scaleZ) = 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/MonteCarloBarostat.h

@@ -35,6 +35,7 @@
 #include "Force.h"
 #include <string>
 #include "internal/windowsExport.h"
+#include "Vec3.h"
 
 namespace OpenMM {
 
@@ -64,9 +65,8 @@ public:
      * @param defaultPressure   the default pressure acting on the system (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 anisotropic       switch to determine whether anisotropic barostat (independent xyz-coordinate scaling) is turned on.
      */
-    MonteCarloBarostat(double defaultPressure, double temperature, int frequency = 25, bool anisotropic = 0);
+    MonteCarloBarostat(double defaultPressure, double temperature, int frequency = 25);
     /**
      * Get the default pressure acting on the system (in bar).
      *
@@ -76,16 +76,105 @@ public:
         return defaultPressure;
     }
     /**
-     * Get the switch which specifies whether anisotropic barostat is turned on. 1 = On, 0 = Off
+     * 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:
+    double defaultPressure, temperature;
+    int frequency, randomNumberSeed;
+
+        double GetTemperature() const {
+            return temperature;
+        }
+
+        void SetTemperature(double temperature) {
+            this->temperature = temperature;
+        }
+    };
+
+/**
+ * This class uses a Monte Carlo algorithm to adjust the size of the periodic box, simulating the
+ * effect of constant pressure.
+ *
+ * Compared to MonteCarloBarostat, this class scales the three axes of the simulation cell independently.
+ * The user supplies a Vec3 instead of a double to specify the pressure along each axis if desired.
+ *
+ * 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.
  */
-    int getAnisotropic() const {
-        return anisotropic;
+
+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).
+     */
+    static const std::string& Pressure() {
+        static const std::string key = "MonteCarloPressure";
+        return key;
     }
     /**
-     * Set the switch which specifies whether anisotropic barostat is turned on. 1 = On, 0 = Off
+     * Create a MonteCarloAnisotropicBarostat.
+     *
+     * @param defaultPressure   the default pressure acting on each axis of the system (in bar)
+     * @param temperature       the temperature at which the system is being maintained (in Kelvin)
+     * @param frequency         the frequency at which Monte Carlo pressure changes should be attempted (in time steps)
+     * @param scaleX            on/off switch for whether to scale the X axis
+     * @param scaleY            on/off switch for whether to scale the Y axis
+     * @param scaleZ            on/off switch for whether to scale the Z axis
+     */
+    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).
+     *
+     * @return the default pressure acting on the system, measured in bar.
      */
-    void setAnisotropic(bool aniso) {
-        anisotropic = aniso;
+    Vec3 getDefaultPressure() const {
+        return defaultPressure;
     }
     /**
      * Get the frequency (in time steps) at which Monte Carlo pressure changes should be attempted.  If this is set to
@@ -134,9 +223,9 @@ public:
 protected:
     ForceImpl* createImpl() const;
 private:
-    double defaultPressure, temperature;
+    Vec3 defaultPressure;
+    double temperature;
     int frequency, randomNumberSeed;
-    bool anisotropic;
 
         double GetTemperature() const {
             return temperature;

openmmapi/include/openmm/internal/MonteCarloBarostatImpl.h

@@ -35,6 +35,7 @@
 #include "ForceImpl.h"
 #include "openmm/MonteCarloBarostat.h"
 #include "openmm/Kernel.h"
+#include "openmm/Vec3.h"
 #include "sfmt/SFMT.h"
 #include <string>
 
@@ -66,6 +67,28 @@ private:
     Kernel kernel;
 };
 
+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, Vec3> 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_MONTECARLOBAROSTATIMPL_H_*/

openmmapi/src/MonteCarloBarostat.cpp

@@ -35,11 +35,20 @@
 
 using namespace OpenMM;
 
-MonteCarloBarostat::MonteCarloBarostat(double defaultPressure, double temperature, int frequency, bool anisotropic) :
-  defaultPressure(defaultPressure), temperature(temperature), frequency(frequency), anisotropic(anisotropic) {
+MonteCarloBarostat::MonteCarloBarostat(double defaultPressure, double temperature, int frequency) :
+        defaultPressure(defaultPressure), temperature(temperature), frequency(frequency) {
     setRandomNumberSeed((int) time(NULL));
 }
 
 ForceImpl* MonteCarloBarostat::createImpl() const {
     return new MonteCarloBarostatImpl(*this);
 }
+
+MonteCarloAnisotropicBarostat::MonteCarloAnisotropicBarostat(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/MonteCarloBarostatImpl.cpp

@@ -74,36 +74,11 @@ void MonteCarloBarostatImpl::updateContextState(ContextImpl& context) {
     Vec3 box[3];
     context.getPeriodicBoxVectors(box[0], box[1], box[2]);
     double volume = box[0][0]*box[1][1]*box[2][2];
-    
-    // Decide whether to scale coordinates or distort box.  
-    // If anisotropic parameter is set, scale the box with with 
-    // 50% chance (otherwise perform volume-preserving distortion).
     double deltaVolume = volumeScale*2*(genrand_real2(random)-0.5);
     double newVolume = volume+deltaVolume;
     double lengthScale = std::pow(newVolume/volume, 1.0/3.0);
-    if ((! owner.getAnisotropic()) || (genrand_real2(random) < 0.5)) {
-      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);
-    } 
-    else {
-      deltaVolume = 0.0;
-      newVolume = volume;
-      double lengthScaleBack = std::pow(lengthScale, -1.0/2.0);
-      double lengthScaleX = lengthScaleBack;
-      double lengthScaleY = lengthScaleBack;
-      double lengthScaleZ = lengthScaleBack;
-      double randXYZ = 3.0 * genrand_real2(random);
-      // Randomly choose an axis to lengthen and shorten the other two in a volume-preserving way.
-      if (randXYZ < 1.0) {
-	lengthScaleX = lengthScale;
-      } else if (randXYZ < 2.0) {
-	lengthScaleY = lengthScale;
-      } else {
-	lengthScaleZ = lengthScale;
-      }
-      kernel.getAs<ApplyMonteCarloBarostatKernel>().scaleCoordinatesXYZ(context, lengthScaleX, lengthScaleY, lengthScaleZ);
-      context.getOwner().setPeriodicBoxVectors(box[0]*lengthScaleX, box[1]*lengthScaleY, box[2]*lengthScaleZ);
-    }
 
     // Compute the energy of the modified system.
     
@@ -147,3 +122,103 @@ std::vector<std::string> MonteCarloBarostatImpl::getKernelNames() {
     return names;
 }
 
+MonteCarloAnisotropicBarostatImpl::MonteCarloAnisotropicBarostatImpl(const MonteCarloBarostat& 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 (scaleX == 0 && scaleY == 0 && scaleZ == 0)
+        return;
+    step = 0;
+
+    // Compute the current potential energy.
+
+    double initialEnergy = context.getOwner().getState(State::Energy).getPotentialEnergy();
+
+    // Choose which axis to modify at random.
+    double rnd = genrand_real2(random)*3.0;
+    int axis;
+    while (1) {
+        if (rnd < 1.0 && scaleX) {
+            axis = 0;
+            break;
+        } else if (rnd < 2.0 && scaleY) {
+            axis = 1;
+            break;
+        } else if (scaleZ) {
+            axis = 2;
+            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;
+    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]);
+    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);
+        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[3]> MonteCarloAnisotropicBarostatImpl::getDefaultParameters() {
+    std::map<std::string, double[3]> parameters;
+    parameters[MonteCarloBarostat::Pressure()] = getOwner().getDefaultPressure();
+    return parameters;
+}
+
+std::vector<std::string> MonteCarloAnisotropicBarostatImpl::getKernelNames() {
+    std::vector<std::string> names;
+    names.push_back(ApplyMonteCarloBarostatKernel::Name());
+    return names;
+}
+

platforms/cuda/src/CudaKernels.cpp

@@ -5314,52 +5314,10 @@ void CudaApplyMonteCarloBarostatKernel::initialize(const System& system, const M
     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, "scalePositionsXYZ");
+    kernel = cu.getKernel(module, "scalePositions");
 }
 
-void CudaApplyMonteCarloBarostatKernel::scaleCoordinates(ContextImpl& context, double scale) {
-    cu.setAsCurrent();
-    if (!hasInitializedKernels) {
-        hasInitializedKernels = true;
-
-        // Create the arrays with the molecule definitions.
-
-        vector<vector<int> > molecules = context.getMolecules();
-        numMolecules = molecules.size();
-        moleculeAtoms = CudaArray::create<int>(cu, cu.getNumAtoms(), "moleculeAtoms");
-        moleculeStartIndex = CudaArray::create<int>(cu, numMolecules+1, "moleculeStartIndex");
-        vector<int> atoms(moleculeAtoms->getSize());
-        vector<int> startIndex(moleculeStartIndex->getSize());
-        int index = 0;
-        for (int i = 0; i < numMolecules; i++) {
-            startIndex[i] = index;
-            for (int j = 0; j < (int) molecules[i].size(); j++)
-                atoms[index++] = molecules[i][j];
-        }
-        startIndex[numMolecules] = index;
-        moleculeAtoms->upload(atoms);
-        moleculeStartIndex->upload(startIndex);
-
-        // Initialize the kernel arguments.
-        
-    }
-    int bytesToCopy = cu.getPosq().getSize()*(cu.getUseDoublePrecision() ? sizeof(double4) : sizeof(float4));
-    CUresult result = cuMemcpyDtoD(savedPositions->getDevicePointer(), cu.getPosq().getDevicePointer(), bytesToCopy);
-    if (result != CUDA_SUCCESS) {
-        std::stringstream m;
-        m<<"Error saving positions for MC barostat: "<<cu.getErrorString(result)<<" ("<<result<<")";
-        throw OpenMMException(m.str());
-    }        
-    float scalef = (float) scale;
-    void* args[] = {&scalef, &scalef, &scalef, &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++)
-        cu.getPosCellOffsets()[i] = make_int4(0, 0, 0, 0);
-    lastAtomOrder = cu.getAtomIndex();
-}
-
-void CudaApplyMonteCarloBarostatKernel::scaleCoordinatesXYZ(ContextImpl& context, double scaleX, double scaleY, double scaleZ) {
+void CudaApplyMonteCarloBarostatKernel::scaleCoordinates(ContextImpl& context, double scaleX, double scaleY, double scaleZ) {
     cu.setAsCurrent();
     if (!hasInitializedKernels) {
         hasInitializedKernels = true;

platforms/cuda/src/CudaKernels.h

@@ -1255,16 +1255,6 @@ public:
      * @param barostat   the MonteCarloBarostat this kernel will be used for
      */
     void initialize(const System& system, const MonteCarloBarostat& barostat);
-    /**
-     * Attempt a Monte Carlo step, scaling particle positions (or cluster centers) by a specified value.
-     * 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
-     */
-    void scaleCoordinates(ContextImpl& context, double scale);
     /**
      * 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.
@@ -1277,7 +1267,7 @@ public:
      * @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 scaleCoordinatesXYZ(ContextImpl& context, double scaleX, double scaleY, double scaleZ);
+    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/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,
-        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];
-        int last = moleculeStartIndex[index+1];
-        int numAtoms = last-first;
-
-        // Find the center of each molecule.
-
-        real3 center = make_real3(0, 0, 0);
-        for (int atom = first; atom < last; atom++) {
-            real4 pos = posq[moleculeAtoms[atom]];
-            center.x += pos.x;
-            center.y += pos.y;
-            center.z += pos.z;
-        }
-        real invNumAtoms = RECIP(numAtoms);
-        center.x *= invNumAtoms;
-        center.y *= invNumAtoms;
-        center.z *= invNumAtoms;
-
-        // Move it into the first periodic box.
-
-        int xcell = (int) floor(center.x*invPeriodicBoxSize.x);
-        int ycell = (int) floor(center.y*invPeriodicBoxSize.y);
-        int zcell = (int) floor(center.z*invPeriodicBoxSize.z);
-        real3 delta = make_real3(xcell*periodicBoxSize.x, ycell*periodicBoxSize.y, zcell*periodicBoxSize.z);
-        center.x -= delta.x;
-        center.y -= delta.y;
-        center.z -= delta.z;
-
-        // 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;
-        for (int atom = first; atom < last; atom++) {
-            real4 pos = posq[moleculeAtoms[atom]];
-            pos.x += delta.x;
-            pos.y += delta.y;
-            pos.z += delta.z;
-            posq[moleculeAtoms[atom]] = pos;
-        }
-    }
-}
-
-/**
- * Scale the particle positions.
- */
-
-extern "C" __global__ void scalePositionsXYZ(float scaleX, float scaleY, float scaleZ, 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];

platforms/opencl/src/OpenCLKernels.cpp

@@ -5539,51 +5539,10 @@ OpenCLApplyMonteCarloBarostatKernel::~OpenCLApplyMonteCarloBarostatKernel() {
 void OpenCLApplyMonteCarloBarostatKernel::initialize(const System& system, const MonteCarloBarostat& thermostat) {
     savedPositions = OpenCLArray::create<mm_float4>(cl, cl.getPaddedNumAtoms(), "savedPositions");
     cl::Program program = cl.createProgram(OpenCLKernelSources::monteCarloBarostat);
-    kernel = cl::Kernel(program, "scalePositionsXYZ");
+    kernel = cl::Kernel(program, "scalePositions");
 }
 
-void OpenCLApplyMonteCarloBarostatKernel::scaleCoordinates(ContextImpl& context, double scale) {
-    if (!hasInitializedKernels) {
-        hasInitializedKernels = true;
-
-        // Create the arrays with the molecule definitions.
-
-        vector<vector<int> > molecules = context.getMolecules();
-        numMolecules = molecules.size();
-        moleculeAtoms = OpenCLArray::create<int>(cl, cl.getNumAtoms(), "moleculeAtoms");
-        moleculeStartIndex = OpenCLArray::create<int>(cl, numMolecules+1, "moleculeStartIndex");
-        vector<int> atoms(moleculeAtoms->getSize());
-        vector<int> startIndex(moleculeStartIndex->getSize());
-        int index = 0;
-        for (int i = 0; i < numMolecules; i++) {
-            startIndex[i] = index;
-            for (int j = 0; j < (int) molecules[i].size(); j++)
-                atoms[index++] = molecules[i][j];
-        }
-        startIndex[numMolecules] = index;
-        moleculeAtoms->upload(atoms);
-        moleculeStartIndex->upload(startIndex);
-
-        // Initialize the kernel arguments.
-        
-        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);
-    kernel.setArg<cl_float>(1, (cl_float) scale);
-    kernel.setArg<cl_float>(2, (cl_float) scale);
-    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);
-    lastAtomOrder = cl.getAtomIndex();
-}
-
-void OpenCLApplyMonteCarloBarostatKernel::scaleCoordinatesXYZ(ContextImpl& context, double scaleX, double scaleY, double scaleZ) {
+void OpenCLApplyMonteCarloBarostatKernel::scaleCoordinates(ContextImpl& context, double scaleX, double scaleY, double scaleZ) {
     if (!hasInitializedKernels) {
         hasInitializedKernels = true;
 

platforms/opencl/src/OpenCLKernels.h

@@ -1269,16 +1269,6 @@ public:
      * @param barostat   the MonteCarloBarostat this kernel will be used for
      */
     void initialize(const System& system, const MonteCarloBarostat& barostat);
-    /**
-     * Attempt a Monte Carlo step, scaling particle positions (or cluster centers) by a specified value.
-     * 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
-     */
-    void scaleCoordinates(ContextImpl& context, double scale);
     /**
      * 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.
@@ -1291,7 +1281,7 @@ public:
      * @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 scaleCoordinatesXYZ(ContextImpl& context, double scaleX, double scaleY, double scaleZ);
+    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/kernels/monteCarloBarostat.cl

-/**
- * Scale the particle positions.
- */
-
-__kernel void scalePositions(float scale, 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];
-        int last = moleculeStartIndex[index+1];
-        int numAtoms = last-first;
-
-        // Find the center of each molecule.
-
-        real4 center = (real4) 0;
-        for (int atom = first; atom < last; atom++)
-            center += posq[moleculeAtoms[atom]];
-        center /= (real) numAtoms;
-
-        // Move it into the first periodic box.
-
-        int xcell = (int) floor(center.x*invPeriodicBoxSize.x);
-        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);
-        center -= delta;
-
-        // Now scale the position of the molecule center.
-
-        delta = center*(scale-1)-delta;
-        for (int atom = first; atom < last; atom++) {
-            real4 pos = posq[moleculeAtoms[atom]];
-            pos.xyz += delta.xyz;
-            posq[moleculeAtoms[atom]] = pos;
-        }
-    }
-}
-
 /**
  * Scale the particle positions with each axis independent.
  */
 
-__kernel void scalePositionsXYZ(float scaleX, float scaleY, float scaleZ, 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];

platforms/reference/src/ReferenceKernels.cpp

@@ -2143,20 +2143,12 @@ ReferenceApplyMonteCarloBarostatKernel::~ReferenceApplyMonteCarloBarostatKernel(
 void ReferenceApplyMonteCarloBarostatKernel::initialize(const System& system, const MonteCarloBarostat& barostat) {
 }
 
-void ReferenceApplyMonteCarloBarostatKernel::scaleCoordinates(ContextImpl& context, double scale) {
+void ReferenceApplyMonteCarloBarostatKernel::scaleCoordinates(ContextImpl& context, double scaleX, double scaleY, double scaleZ) {
     if (barostat == NULL)
         barostat = new ReferenceMonteCarloBarostat(context.getSystem().getNumParticles(), context.getMolecules());
     vector<RealVec>& posData = extractPositions(context);
     RealVec& boxSize = extractBoxSize(context);
-    barostat->applyBarostatXYZ(posData, boxSize, scale, scale, scale);
-}
-
-void ReferenceApplyMonteCarloBarostatKernel::scaleCoordinatesXYZ(ContextImpl& context, double scaleX, double scaleY, double scaleZ) {
-    if (barostat == NULL)
-        barostat = new ReferenceMonteCarloBarostat(context.getSystem().getNumParticles(), context.getMolecules());
-    vector<RealVec>& posData = extractPositions(context);
-    RealVec& boxSize = extractBoxSize(context);
-    barostat->applyBarostatXYZ(posData, boxSize, scaleX, scaleY, scaleZ);
+    barostat->applyBarostat(posData, boxSize, scaleX, scaleY, scaleZ);
 }
 
 void ReferenceApplyMonteCarloBarostatKernel::restoreCoordinates(ContextImpl& context) {

platforms/reference/src/ReferenceKernels.h

@@ -1183,16 +1183,6 @@ public:
      * @param barostat   the MonteCarloBarostat this kernel will be used for
      */
     void initialize(const System& system, const MonteCarloBarostat& barostat);
-    /**
-     * Attempt a Monte Carlo step, scaling particle positions (or cluster centers) by a specified value.
-     * 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
-     */
-    void scaleCoordinates(ContextImpl& context, double scale);
     /**
      * 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.
@@ -1205,7 +1195,7 @@ public:
      * @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 scaleCoordinatesXYZ(ContextImpl& context, double scaleX, double scaleY, double scaleZ);
+    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/src/SimTKReference/ReferenceMonteCarloBarostat.cpp

@@ -58,59 +58,13 @@ ReferenceMonteCarloBarostat::~ReferenceMonteCarloBarostat( ) {
 
   @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 atom x-coordinates
+  @param scaleY             the factor by which to scale atom y-coordinates
+  @param scaleZ             the factor by which to scale atom z-coordinates
 
   --------------------------------------------------------------------------------------- */
 
-void ReferenceMonteCarloBarostat::applyBarostat(vector<RealVec>& atomPositions, const RealVec& boxSize, RealOpenMM scale) {
-    int numAtoms = savedAtomPositions[0].size();
-    for (int i = 0; i < numAtoms; i++)
-        for (int j = 0; j < 3; j++)
-            savedAtomPositions[j][i] = atomPositions[i][j];
-
-    // Loop over molecules.
-
-    for (int i = 0; i < (int) molecules.size(); i++) {
-        // Find the molecule center.
-
-        RealOpenMM pos[3] = {0, 0, 0};
-        for (int j = 0; j < (int) molecules[i].size(); j++) {
-            RealVec& atomPos = atomPositions[molecules[i][j]];
-            pos[0] += atomPos[0];
-            pos[1] += atomPos[1];
-            pos[2] += atomPos[2];
-        }
-        pos[0] /= molecules[i].size();
-        pos[1] /= molecules[i].size();
-        pos[2] /= molecules[i].size();
-
-        // Move it into the first periodic box.
-
-        int xcell = (int) floor(pos[0]/boxSize[0]);
-        int ycell = (int) floor(pos[1]/boxSize[1]);
-        int zcell = (int) floor(pos[2]/boxSize[2]);
-        RealOpenMM dx = xcell*boxSize[0];
-        RealOpenMM dy = ycell*boxSize[1];
-        RealOpenMM dz = zcell*boxSize[2];
-        pos[0] -= dx;
-        pos[1] -= dy;
-        pos[2] -= dz;
-
-        // Now scale the position of the molecule center.
-
-        dx = pos[0]*(scale-1)-dx;
-        dy = pos[1]*(scale-1)-dy;
-        dz = pos[2]*(scale-1)-dz;
-        for (int j = 0; j < (int) molecules[i].size(); j++) {
-            RealVec& atomPos = atomPositions[molecules[i][j]];
-            atomPos[0] += dx;
-            atomPos[1] += dy;
-            atomPos[2] += dz;
-        }
-    }
-}
-
-void ReferenceMonteCarloBarostat::applyBarostatXYZ(vector<RealVec>& atomPositions, const RealVec& boxSize, RealOpenMM scaleX, RealOpenMM scaleY, RealOpenMM scaleZ) {
+void ReferenceMonteCarloBarostat::applyBarostat(vector<RealVec>& atomPositions, const RealVec& boxSize, RealOpenMM scaleX, RealOpenMM scaleY, RealOpenMM scaleZ) {
     int numAtoms = savedAtomPositions[0].size();
     for (int i = 0; i < numAtoms; i++)
         for (int j = 0; j < 3; j++)

platforms/reference/src/SimTKReference/ReferenceMonteCarloBarostat.h

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