Added anisotropic MC barostat option, allowing x, y, and z lattice vectors to change independently.

olla/include/openmm/kernels.h

@@ -1118,6 +1118,19 @@ public:
      * @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.
+     * 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 scaleCoordinatesXYZ(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

@@ -64,8 +64,9 @@ 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);
+    MonteCarloBarostat(double defaultPressure, double temperature, int frequency = 25, int anisotropic = 0);
     /**
      * Get the default pressure acting on the system (in bar).
      *
@@ -74,6 +75,18 @@ public:
     double getDefaultPressure() const {
         return defaultPressure;
     }
+    /**
+     * Get the switch which specifies whether anisotropic barostat is turned on. 1 = On, 0 = Off
+     */
+    int getAnisotropic() const {
+        return anisotropic;
+    }
+    /**
+     * Set the switch which specifies whether anisotropic barostat is turned on. 1 = On, 0 = Off
+     */
+    void setAnisotropic(int aniso) {
+        anisotropic = aniso;
+    }
     /**
      * 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.
@@ -122,7 +135,7 @@ protected:
     ForceImpl* createImpl() const;
 private:
     double defaultPressure, temperature;
-    int frequency, randomNumberSeed;
+    int frequency, anisotropic, randomNumberSeed;
 
         double GetTemperature() const {
             return temperature;

openmmapi/src/MonteCarloBarostat.cpp

@@ -35,9 +35,9 @@
 
 using namespace OpenMM;
 
-MonteCarloBarostat::MonteCarloBarostat(double defaultPressure, double temperature, int frequency) :
-        defaultPressure(defaultPressure), temperature(temperature), frequency(frequency) {
-    setRandomNumberSeed((int) time(NULL));
+MonteCarloBarostat::MonteCarloBarostat(double defaultPressure, double temperature, int frequency, int anisotropic) :
+  defaultPressure(defaultPressure), temperature(temperature), frequency(frequency), anisotropic(anisotropic) {
+  setRandomNumberSeed((int) time(NULL));
 }
 
 ForceImpl* MonteCarloBarostat::createImpl() const {

openmmapi/src/MonteCarloBarostatImpl.cpp

@@ -74,11 +74,36 @@ 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);
-    kernel.getAs<ApplyMonteCarloBarostatKernel>().scaleCoordinates(context, lengthScale);
-    context.getOwner().setPeriodicBoxVectors(box[0]*lengthScale, box[1]*lengthScale, box[2]*lengthScale);
+    if ((! owner.getAnisotropic()) || (genrand_real2(random) < 0.5)) {
+      kernel.getAs<ApplyMonteCarloBarostatKernel>().scaleCoordinates(context, 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.
     

platforms/cuda/src/CudaKernels.cpp

@@ -5314,7 +5314,7 @@ 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, "scalePositions");
+    kernel = cu.getKernel(module, "scalePositionsXYZ");
 }
 
 void CudaApplyMonteCarloBarostatKernel::scaleCoordinates(ContextImpl& context, double scale) {
@@ -5351,7 +5351,7 @@ void CudaApplyMonteCarloBarostatKernel::scaleCoordinates(ContextImpl& context, d
         throw OpenMMException(m.str());
     }        
     float scalef = (float) scale;
-    void* args[] = {&scalef, &numMolecules, cu.getPeriodicBoxSizePointer(), cu.getInvPeriodicBoxSizePointer(),
+    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++)
@@ -5359,6 +5359,50 @@ void CudaApplyMonteCarloBarostatKernel::scaleCoordinates(ContextImpl& context, d
     lastAtomOrder = cu.getAtomIndex();
 }
 
+void CudaApplyMonteCarloBarostatKernel::scaleCoordinatesXYZ(ContextImpl& context, double scaleX, double scaleY, double scaleZ) {
+    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 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++)
+        cu.getPosCellOffsets()[i] = make_int4(0, 0, 0, 0);
+    lastAtomOrder = cu.getAtomIndex();
+}
+
 void CudaApplyMonteCarloBarostatKernel::restoreCoordinates(ContextImpl& context) {
     cu.setAsCurrent();
     int bytesToCopy = cu.getPosq().getSize()*(cu.getUseDoublePrecision() ? sizeof(double4) : sizeof(float4));

platforms/cuda/src/CudaKernels.h

@@ -1265,6 +1265,19 @@ public:
      * @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.
+     * 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
+     */
+    void scaleCoordinatesXYZ(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

@@ -47,3 +47,53 @@ extern "C" __global__ void scalePositions(float scale, int numMolecules, real4 p
         }
     }
 }
+
+/**
+ * Scale the particle positions.
+ */
+
+extern "C" __global__ void scalePositionsXYZ(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];
+        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*(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;
+            pos.y += delta.y;
+            pos.z += delta.z;
+            posq[moleculeAtoms[atom]] = pos;
+        }
+    }
+}

platforms/cuda/tests/TestCudaMonteCarloBarostat.cpp

@@ -89,7 +89,7 @@ void testChangingBoxSize() {
     ASSERT(ok);
 }
 
-void testIdealGas() {
+void testIdealGas(int aniso) {
     const int numParticles = 64;
     const int frequency = 10;
     const int steps = 1000;
@@ -110,7 +110,7 @@ void testIdealGas() {
         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);
+    MonteCarloBarostat* barostat = new MonteCarloBarostat(pressure, temp[0], frequency, aniso);
     system.addForce(barostat);
 
     // Test it for three different temperatures.
@@ -132,8 +132,10 @@ void testIdealGas() {
             Vec3 box[3];
             context.getState(0).getPeriodicBoxVectors(box[0], box[1], box[2]);
             volume += box[0][0]*box[1][1]*box[2][2];
-            ASSERT_EQUAL_TOL(0.5*box[0][0], box[1][1], 1e-5);
-            ASSERT_EQUAL_TOL(2*box[0][0], box[2][2], 1e-5);
+	    if (!aniso) {
+	      ASSERT_EQUAL_TOL(0.5*box[0][0], box[1][1], 1e-5);
+	      ASSERT_EQUAL_TOL(2*box[0][0], box[2][2], 1e-5);
+	    }
             integrator.step(frequency);
         }
         volume /= steps;
@@ -204,7 +206,7 @@ void testRandomSeed() {
     }
 }
 
-void testWater() {
+void testWater(int aniso) {
     const int gridSize = 8;
     const int numMolecules = gridSize*gridSize*gridSize;
     const int frequency = 10;
@@ -250,7 +252,7 @@ void testWater() {
         }
     }
     system.addForce(nonbonded);
-    MonteCarloBarostat* barostat = new MonteCarloBarostat(pressure, temp, frequency);
+    MonteCarloBarostat* barostat = new MonteCarloBarostat(pressure, temp, frequency, aniso);
     system.addForce(barostat);
 
     // Simulate it and see if the density matches the expected value (1 g/mL).
@@ -276,9 +278,11 @@ int main(int argc, char* argv[]) {
         if (argc > 1)
             platform.setPropertyDefaultValue("CudaPrecision", string(argv[1]));
         testChangingBoxSize();
-        testIdealGas();
+        testIdealGas(0);
+        testIdealGas(1);
         testRandomSeed();
-        testWater();
+        testWater(0);
+        testWater(1);
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;

platforms/opencl/src/OpenCLKernels.cpp

@@ -5539,7 +5539,7 @@ 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, "scalePositions");
+    kernel = cl::Kernel(program, "scalePositionsXYZ");
 }
 
 void OpenCLApplyMonteCarloBarostatKernel::scaleCoordinates(ContextImpl& context, double scale) {
@@ -5566,15 +5566,58 @@ 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>(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) {
+    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) 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/OpenCLKernels.h

@@ -1279,6 +1279,19 @@ public:
      * @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.
+     * 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
+     */
+    void scaleCoordinatesXYZ(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

@@ -34,3 +34,41 @@ __kernel void scalePositions(float scale, int numMolecules, real4 periodicBoxSiz
         }
     }
 }
+
+/**
+ * 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,
+        __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);
+	real4 scaleXYZ = (real4) (scaleX, scaleY, scaleZ, 1);
+        center -= delta;
+
+        // Now scale the position of the molecule center.
+
+        delta = center*(scaleXYZ-1)-delta;
+        for (int atom = first; atom < last; atom++) {
+            real4 pos = posq[moleculeAtoms[atom]];
+            pos.xyz += delta.xyz;
+            posq[moleculeAtoms[atom]] = pos;
+        }
+    }
+}

platforms/opencl/tests/TestOpenCLMonteCarloBarostat.cpp

@@ -89,7 +89,7 @@ void testChangingBoxSize() {
     ASSERT(ok);
 }
 
-void testIdealGas() {
+void testIdealGas(int aniso) {
     const int numParticles = 64;
     const int frequency = 10;
     const int steps = 1000;
@@ -110,7 +110,7 @@ void testIdealGas() {
         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);
+    MonteCarloBarostat* barostat = new MonteCarloBarostat(pressure, temp[0], frequency, aniso);
     system.addForce(barostat);
 
     // Test it for three different temperatures.
@@ -132,8 +132,13 @@ void testIdealGas() {
             Vec3 box[3];
             context.getState(0).getPeriodicBoxVectors(box[0], box[1], box[2]);
             volume += box[0][0]*box[1][1]*box[2][2];
-            ASSERT_EQUAL_TOL(0.5*box[0][0], box[1][1], 1e-5);
-            ASSERT_EQUAL_TOL(2*box[0][0], box[2][2], 1e-5);
+	    
+	    // Ratios will only be correct if box deformations are isotropic.
+	    
+	    if (!aniso) {
+	      ASSERT_EQUAL_TOL(0.5*box[0][0], box[1][1], 1e-5);
+	      ASSERT_EQUAL_TOL(2*box[0][0], box[2][2], 1e-5);
+	    }
             integrator.step(frequency);
         }
         volume /= steps;
@@ -204,7 +209,7 @@ void testRandomSeed() {
     }
 }
 
-void testWater() {
+void testWater(int aniso) {
     const int gridSize = 8;
     const int numMolecules = gridSize*gridSize*gridSize;
     const int frequency = 10;
@@ -250,7 +255,7 @@ void testWater() {
         }
     }
     system.addForce(nonbonded);
-    MonteCarloBarostat* barostat = new MonteCarloBarostat(pressure, temp, frequency);
+    MonteCarloBarostat* barostat = new MonteCarloBarostat(pressure, temp, frequency, aniso);
     system.addForce(barostat);
 
     // Simulate it and see if the density matches the expected value (1 g/mL).
@@ -276,9 +281,11 @@ int main(int argc, char* argv[]) {
         if (argc > 1)
             platform.setPropertyDefaultValue("OpenCLPrecision", string(argv[1]));
         testChangingBoxSize();
-        testIdealGas();
+        testIdealGas(0);
+        testIdealGas(1);
         testRandomSeed();
-        testWater();
+        testWater(0);
+        testWater(1);
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;

platforms/reference/src/ReferenceKernels.cpp

@@ -2148,7 +2148,15 @@ void ReferenceApplyMonteCarloBarostatKernel::scaleCoordinates(ContextImpl& conte
         barostat = new ReferenceMonteCarloBarostat(context.getSystem().getNumParticles(), context.getMolecules());
     vector<RealVec>& posData = extractPositions(context);
     RealVec& boxSize = extractBoxSize(context);
-    barostat->applyBarostat(posData, boxSize, scale);
+    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);
 }
 
 void ReferenceApplyMonteCarloBarostatKernel::restoreCoordinates(ContextImpl& context) {

platforms/reference/src/ReferenceKernels.h

@@ -1193,6 +1193,19 @@ public:
      * @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.
+     * 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
+     */
+    void scaleCoordinatesXYZ(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

@@ -110,6 +110,54 @@ void ReferenceMonteCarloBarostat::applyBarostat(vector<RealVec>& atomPositions,
     }
 }
 
+void ReferenceMonteCarloBarostat::applyBarostatXYZ(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++)
+            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]*(scaleX-1)-dx;
+        dy = pos[1]*(scaleY-1)-dy;
+        dz = pos[2]*(scaleZ-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;
+        }
+    }
+}
+
 /**---------------------------------------------------------------------------------------
 
   Restore atom positions to what they were before applyBarostat() was called.

platforms/reference/src/SimTKReference/ReferenceMonteCarloBarostat.h

@@ -68,6 +68,20 @@ class ReferenceMonteCarloBarostat {
 
       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.
+
+         @param atomPositions      atom positions
+         @param boxSize            the periodic box dimensions
+         @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 applyBarostatXYZ(std::vector<OpenMM::RealVec>& atomPositions, const OpenMM::RealVec& boxSize, RealOpenMM scaleX, RealOpenMM scaleY, RealOpenMM scaleZ);
+
       /**---------------------------------------------------------------------------------------
 
          Restore atom positions to what they were before applyBarostat() was called.

platforms/reference/tests/TestReferenceMonteCarloBarostat.cpp

@@ -88,7 +88,7 @@ void testChangingBoxSize() {
     ASSERT(ok);
 }
 
-void testIdealGas() {
+void testIdealGas(int aniso) {
     const int numParticles = 64;
     const int frequency = 10;
     const int steps = 1000;
@@ -110,7 +110,7 @@ void testIdealGas() {
         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);
+    MonteCarloBarostat* barostat = new MonteCarloBarostat(pressure, temp[0], frequency, aniso);
     system.addForce(barostat);
 
     // Test it for three different temperatures.
@@ -132,8 +132,10 @@ void testIdealGas() {
             Vec3 box[3];
             context.getState(0).getPeriodicBoxVectors(box[0], box[1], box[2]);
             volume += box[0][0]*box[1][1]*box[2][2];
-            ASSERT_EQUAL_TOL(0.5*box[0][0], box[1][1], 1e-5);
-            ASSERT_EQUAL_TOL(2*box[0][0], box[2][2], 1e-5);
+	    if (!aniso) {
+	      ASSERT_EQUAL_TOL(0.5*box[0][0], box[1][1], 1e-5);
+	      ASSERT_EQUAL_TOL(2*box[0][0], box[2][2], 1e-5);
+	    }
             integrator.step(frequency);
         }
         volume /= steps;
@@ -208,7 +210,8 @@ void testRandomSeed() {
 int main() {
     try {
         testChangingBoxSize();
-        testIdealGas();
+        testIdealGas(0);
+        testIdealGas(1);
         testRandomSeed();
     }
     catch(const exception& e) {