Fixed problems in MonteCarloBarostat and created OpenCL implementation

openmmapi/include/openmm/internal/ContextImpl.h

@@ -163,13 +163,21 @@ public:
      * Set the platform-specific data stored in this context.
      */
     void setPlatformData(void* data);
+    /**
+     * Get a list of the particles in each molecules in the system.  Two particles are in the
+     * same molecule if they are connected by constraints or bonds.
+     */
+    const std::vector<std::vector<int> >& getMolecules() const;
 private:
     friend class Context;
+    static void tagParticlesInMolecule(int particle, int molecule, std::vector<int>& particleMolecule, std::vector<std::vector<int> >& particleBonds);
     Context& owner;
     System& system;
     Integrator& integrator;
     std::vector<ForceImpl*> forceImpls;
     std::map<std::string, double> parameters;
+    mutable std::vector<std::vector<int> > molecules;
+    bool hasInitializedForces;
     Platform* platform;
     Kernel initializeForcesKernel, kineticEnergyKernel, updateStateDataKernel;
     void* platformData;

openmmapi/include/openmm/internal/CustomBondForceImpl.h

@@ -60,6 +60,7 @@ public:
     double calcEnergy(ContextImpl& context);
     std::map<std::string, double> getDefaultParameters();
     std::vector<std::string> getKernelNames();
+    std::vector<std::pair<int, int> > getBondedParticles();
 private:
     CustomBondForce& owner;
     Kernel kernel;

openmmapi/include/openmm/internal/ForceImpl.h

@@ -35,6 +35,7 @@
 #include "openmm/internal/windowsExport.h"
 #include <map>
 #include <string>
+#include <utility>
 #include <vector>
 
 namespace OpenMM {
@@ -100,6 +101,13 @@ public:
      * Get the names of all Kernels used by this Force.
      */
     virtual std::vector<std::string> getKernelNames() = 0;
+    /**
+     * Get pairs of particles connected by bonds by this force.  This is used to determine which particles
+     * are part of the same molecule.
+     */
+    std::vector<std::pair<int, int> > getBondedParticles() {
+        return std::vector<std::pair<int, int> >(0);
+    }
 };
 
 } // namespace OpenMM

openmmapi/include/openmm/internal/HarmonicBondForceImpl.h

@@ -62,6 +62,7 @@ public:
         return std::map<std::string, double>(); // This force field doesn't define any parameters.
     }
     std::vector<std::string> getKernelNames();
+    std::vector<std::pair<int, int> > getBondedParticles();
 private:
     HarmonicBondForce& owner;
     Kernel kernel;

openmmapi/src/ContextImpl.cpp

@@ -37,15 +37,17 @@
 #include "openmm/internal/ForceImpl.h"
 #include "openmm/internal/ContextImpl.h"
 #include <map>
+#include <utility>
 #include <vector>
 
 using namespace OpenMM;
 using std::map;
+using std::pair;
 using std::vector;
 using std::string;
 
 ContextImpl::ContextImpl(Context& owner, System& system, Integrator& integrator, Platform* platform, const map<string, string>& properties) :
-         owner(owner), system(system), integrator(integrator), platform(platform), platformData(NULL) {
+         owner(owner), system(system), integrator(integrator), hasInitializedForces(false), platform(platform), platformData(NULL) {
     vector<string> kernelNames;
     kernelNames.push_back(CalcKineticEnergyKernel::Name());
     kernelNames.push_back(CalcForcesAndEnergyKernel::Name());
@@ -57,6 +59,7 @@ ContextImpl::ContextImpl(Context& owner, System& system, Integrator& integrator,
         vector<string> forceKernels = forceImpls[forceImpls.size()-1]->getKernelNames();
         kernelNames.insert(kernelNames.begin(), forceKernels.begin(), forceKernels.end());
     }
+    hasInitializedForces = true;
     vector<string> integratorKernels = integrator.getKernelNames();
     kernelNames.insert(kernelNames.begin(), integratorKernels.begin(), integratorKernels.end());
     if (platform == 0)
@@ -160,3 +163,54 @@ const void* ContextImpl::getPlatformData() const {
 void ContextImpl::setPlatformData(void* data) {
     platformData = data;
 }
+
+const vector<vector<int> >& ContextImpl::getMolecules() const {
+    if (!hasInitializedForces)
+        throw OpenMMException("ContextImpl: getMolecules() cannot be called until all ForceImpls have been initialized");
+    if (molecules.size() > 0 || system.getNumParticles() == 0)
+        return molecules;
+
+    // First make a list of bonds and constraints.
+
+    vector<pair<int, int> > bonds;
+    for (int i = 0; i < system.getNumConstraints(); i++) {
+        int particle1, particle2;
+        double distance;
+        system.getConstraintParameters(i, particle1, particle2, distance);
+        bonds.push_back(std::make_pair(particle1, particle2));
+    }
+    for (int i = 0; i < (int) forceImpls.size(); i++) {
+        vector<pair<int, int> > forceBonds = forceImpls[i]->getBondedParticles();
+        bonds.insert(bonds.end(), forceBonds.begin(), forceBonds.end());
+    }
+
+    // Make a list of every other particle to which each particle is connected
+
+    int numParticles = system.getNumParticles();
+    vector<vector<int> > particleBonds(numParticles);
+    for (int i = 0; i < (int) bonds.size(); i++) {
+        particleBonds[bonds[i].first].push_back(bonds[i].second);
+        particleBonds[bonds[i].second].push_back(bonds[i].first);
+    }
+
+    // Now tag particles by which molecule they belong to.
+
+    vector<int> particleMolecule(numParticles, -1);
+    int numMolecules = 0;
+    for (int i = 0; i < numParticles; i++)
+        if (particleMolecule[i] == -1)
+            tagParticlesInMolecule(i, numMolecules++, particleMolecule, particleBonds);
+    molecules.resize(numMolecules);
+    for (int i = 0; i < numParticles; i++)
+        molecules[particleMolecule[i]].push_back(i);
+    return molecules;
+}
+
+void ContextImpl::tagParticlesInMolecule(int particle, int molecule, vector<int>& particleMolecule, vector<vector<int> >& particleBonds) {
+    // Recursively tag particles as belonging to a particular molecule.
+
+    particleMolecule[particle] = molecule;
+    for (int i = 0; i < (int) particleBonds[particle].size(); i++)
+        if (particleMolecule[particleBonds[particle][i]] == -1)
+            tagParticlesInMolecule(particleBonds[particle][i], molecule, particleMolecule, particleBonds);
+}

openmmapi/src/CustomBondForceImpl.cpp

@@ -103,3 +103,12 @@ map<string, double> CustomBondForceImpl::getDefaultParameters() {
     return parameters;
 }
 
+vector<pair<int, int> > CustomBondForceImpl::getBondedParticles() {
+    int numBonds = owner.getNumBonds();
+    vector<pair<int, int> > bonds(numBonds);
+    for (int i = 0; i < numBonds; i++) {
+        vector<double> parameters;
+        owner.getBondParameters(i, bonds[i].first, bonds[i].second, parameters);
+    }
+    return bonds;
+}

openmmapi/src/HarmonicBondForceImpl.cpp

@@ -63,3 +63,12 @@ std::vector<std::string> HarmonicBondForceImpl::getKernelNames() {
     return names;
 }
 
+vector<pair<int, int> > HarmonicBondForceImpl::getBondedParticles() {
+    int numBonds = owner.getNumBonds();
+    vector<pair<int, int> > bonds(numBonds);
+    for (int i = 0; i < numBonds; i++) {
+        double length, k;
+        owner.getBondParameters(i, bonds[i].first, bonds[i].second, length, k);
+    }
+    return bonds;
+}

openmmapi/src/MonteCarloBarostatImpl.cpp

@@ -81,9 +81,9 @@ void MonteCarloBarostatImpl::updateContextState(ContextImpl& context) {
     // Compute the energy of the modified system.
     
     double finalEnergy = context.getOwner().getState(State::Energy).getPotentialEnergy();
-    double pressure = context.getParameter(MonteCarloBarostat::Pressure())*(AVOGADRO*1e-25);
+    double pressure = context.getParameter(MonteCarloBarostat::Pressure())/(AVOGADRO*1e-25);
     double kT = BOLTZ*owner.getTemperature();
-    double w = finalEnergy-initialEnergy + pressure*deltaVolume - context.getSystem().getNumParticles()*kT*std::log(newVolume/volume);
+    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.
 

platforms/opencl/src/OpenCLKernelFactory.cpp

@@ -76,6 +76,8 @@ KernelImpl* OpenCLKernelFactory::createKernelImpl(std::string name, const Platfo
         return new OpenCLIntegrateVariableLangevinStepKernel(name, platform, cl);
     if (name == ApplyAndersenThermostatKernel::Name())
         return new OpenCLApplyAndersenThermostatKernel(name, platform, cl);
+    if (name == ApplyMonteCarloBarostatKernel::Name())
+        return new OpenCLApplyMonteCarloBarostatKernel(name, platform, cl);
     if (name == CalcKineticEnergyKernel::Name())
         return new OpenCLCalcKineticEnergyKernel(name, platform, cl);
     if (name == RemoveCMMotionKernel::Name())

platforms/opencl/src/OpenCLKernels.cpp

@@ -3425,6 +3425,61 @@ void OpenCLApplyAndersenThermostatKernel::execute(ContextImpl& context) {
     cl.executeKernel(kernel, cl.getNumAtoms());
 }
 
+OpenCLApplyMonteCarloBarostatKernel::~OpenCLApplyMonteCarloBarostatKernel() {
+    if (savedPositions != NULL)
+        delete savedPositions;
+    if (moleculeAtoms != NULL)
+        delete moleculeAtoms;
+    if (moleculeStartIndex != NULL)
+        delete moleculeStartIndex;
+}
+
+void OpenCLApplyMonteCarloBarostatKernel::initialize(const System& system, const MonteCarloBarostat& thermostat) {
+    savedPositions = new OpenCLArray<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) {
+    if (!hasInitializedKernels) {
+        hasInitializedKernels = true;
+
+        // Create the arrays with the molecule definitions.
+
+        vector<vector<int> > molecules = context.getMolecules();
+        numMolecules = molecules.size();
+        moleculeAtoms = new OpenCLArray<int>(cl, cl.getNumAtoms(), "moleculeAtoms");
+        moleculeStartIndex = new OpenCLArray<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>(1, numMolecules);
+        kernel.setArg<cl::Buffer>(4, cl.getPosq().getDeviceBuffer());
+        kernel.setArg<cl::Buffer>(5, moleculeAtoms->getDeviceBuffer());
+        kernel.setArg<cl::Buffer>(6, 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<mm_float4>(2, cl.getPeriodicBoxSize());
+    kernel.setArg<mm_float4>(3, cl.getInvPeriodicBoxSize());
+    cl.executeKernel(kernel, cl.getNumAtoms());
+}
+
+void OpenCLApplyMonteCarloBarostatKernel::restoreCoordinates(ContextImpl& context) {
+    cl.getQueue().enqueueCopyBuffer(savedPositions->getDeviceBuffer(), cl.getPosq().getDeviceBuffer(), 0, 0, cl.getPosq().getSize()*sizeof(mm_float4));
+}
+
 void OpenCLCalcKineticEnergyKernel::initialize(const System& system) {
     int numParticles = system.getNumParticles();
     masses.resize(numParticles);

platforms/opencl/src/OpenCLKernels.h

@@ -902,7 +902,49 @@ private:
     bool hasInitializedKernels;
     int randomSeed;
     cl::Kernel kernel;
-    double prevTemp, prevFriction, prevStepSize;
+};
+
+/**
+ * This kernel is invoked by MonteCarloBarostat to adjust the periodic box volume
+ */
+class OpenCLApplyMonteCarloBarostatKernel : public ApplyMonteCarloBarostatKernel {
+public:
+    OpenCLApplyMonteCarloBarostatKernel(std::string name, const Platform& platform, OpenCLContext& cl) : ApplyMonteCarloBarostatKernel(name, platform), cl(cl),
+            savedPositions(NULL), moleculeAtoms(NULL), moleculeStartIndex(NULL) {
+    }
+    ~OpenCLApplyMonteCarloBarostatKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @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);
+    /**
+     * 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);
+    /**
+     * 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
+     */
+    void restoreCoordinates(ContextImpl& context);
+private:
+    OpenCLContext& cl;
+    bool hasInitializedKernels;
+    int numMolecules;
+    OpenCLArray<mm_float4>* savedPositions;
+    OpenCLArray<cl_int>* moleculeAtoms;
+    OpenCLArray<cl_int>* moleculeStartIndex;
+    cl::Kernel kernel;
 };
 
 /**

platforms/opencl/src/OpenCLPlatform.cpp

@@ -65,6 +65,7 @@ OpenCLPlatform::OpenCLPlatform() {
     registerKernelFactory(IntegrateVariableVerletStepKernel::Name(), factory);
     registerKernelFactory(IntegrateVariableLangevinStepKernel::Name(), factory);
     registerKernelFactory(ApplyAndersenThermostatKernel::Name(), factory);
+    registerKernelFactory(ApplyMonteCarloBarostatKernel::Name(), factory);
     registerKernelFactory(CalcKineticEnergyKernel::Name(), factory);
     registerKernelFactory(RemoveCMMotionKernel::Name(), factory);
     platformProperties.push_back(OpenCLDeviceIndex());

platforms/opencl/src/kernels/monteCarloBarostat.cl

+/**
+ * Scale the particle positions.
+ */
+
+__kernel void scalePositions(float scale, int numMolecules, float4 periodicBoxSize, float4 invPeriodicBoxSize, __global float4* posq,
+        __global int* moleculeAtoms, __global int* 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.
+
+        float4 center = (float4) 0.0f;
+        for (int atom = first; atom < last; atom++)
+            center += posq[moleculeAtoms[atom]];
+        center /= (float) 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);
+        float4 delta = (float) (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++) {
+            float4 pos = posq[moleculeAtoms[atom]];
+            pos.xyz += delta.xyz;
+            posq[moleculeAtoms[atom]] = pos;
+        }
+    }
+}

platforms/opencl/tests/TestOpenCLMonteCarloBarostat.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-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.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests the OpenCL implementation of MonteCarloBarostat.
+ */
+
+#include "../../../tests/AssertionUtilities.h"
+#include "openmm/MonteCarloBarostat.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 "../src/SimTKUtilities/SimTKOpenMMRealType.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+void testChangingBoxSize() {
+    OpenCLPlatform platform;
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(4, 0, 0), Vec3(0, 5, 0), Vec3(0, 0, 6));
+    system.addParticle(1.0);
+    LangevinIntegrator integrator(300.0, 1.0, 0.01);
+    Context context(system, integrator, platform);
+    Vec3 x, y, z;
+    context.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.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);
+}
+
+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.
+
+    OpenCLPlatform platform;
+    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));
+    }
+    MonteCarloBarostat* barostat = new MonteCarloBarostat(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.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);
+            integrator.step(frequency);
+        }
+        volume /= steps;
+        double expected = numParticles*BOLTZ*temp[i]/pressureInMD;//+numParticles*(4*M_PI/3)*sigma*sigma*sigma;
+        ASSERT_USUALLY_EQUAL_TOL(expected, volume, 3/std::sqrt(steps));
+    }
+}
+
+void testRandomSeed() {
+    const int numParticles = 8;
+    const double temp = 100.0;
+    const double pressure = 1.5;
+    OpenCLPlatform platform;
+    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);
+    MonteCarloBarostat* barostat = new MonteCarloBarostat(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]);
+        }
+    }
+}
+
+void testWater() {
+    const int gridSize = 8;
+    const int numMolecules = gridSize*gridSize*gridSize;
+    const int frequency = 10;
+    const int steps = 400;
+    const double temp = 273.15;
+    const double pressure = 3;
+    const double spacing = 0.32;
+    const double angle = 109.47*M_PI/180;
+    const double dOH = 0.1;
+    const double dHH = dOH*2*std::sin(0.5*angle);
+
+    // Create a box of SPC water molecules.
+
+    OpenCLPlatform platform;
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(gridSize*spacing, 0, 0), Vec3(0, gridSize*spacing, 0), Vec3(0, 0, gridSize*spacing));
+    NonbondedForce* nonbonded = new NonbondedForce();
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    vector<Vec3> positions;
+    Vec3 offset1(dOH, 0, 0);
+    Vec3 offset2(dOH*std::cos(angle), dOH*std::sin(angle), 0);
+    for (int i = 0; i < gridSize; ++i) {
+        for (int j = 0; j < gridSize; ++j) {
+            for (int k = 0; k < gridSize; ++k) {
+                int firstParticle = system.getNumParticles();
+                system.addParticle(16.0);
+                system.addParticle(1.0);
+                system.addParticle(1.0);
+                nonbonded->addParticle(-0.82, 0.316557, 0.650194);
+                nonbonded->addParticle(0.41, 1, 0);
+                nonbonded->addParticle(0.41, 1, 0);
+                Vec3 pos = Vec3(spacing*i, spacing*j, spacing*k);
+                positions.push_back(pos);
+                positions.push_back(pos+offset1);
+                positions.push_back(pos+offset2);
+                system.addConstraint(firstParticle, firstParticle+1, dOH);
+                system.addConstraint(firstParticle, firstParticle+2, dOH);
+                system.addConstraint(firstParticle+1, firstParticle+2, dHH);
+                nonbonded->addException(firstParticle, firstParticle+1, 0, 1, 0);
+                nonbonded->addException(firstParticle, firstParticle+2, 0, 1, 0);
+                nonbonded->addException(firstParticle+1, firstParticle+2, 0, 1, 0);
+            }
+        }
+    }
+    system.addForce(nonbonded);
+    MonteCarloBarostat* barostat = new MonteCarloBarostat(pressure, temp, frequency);
+    system.addForce(barostat);
+
+    // Simulate it and see if the density matches the expected value (1 g/mL).
+
+    LangevinIntegrator integrator(temp, 1.0, 0.002);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    integrator.step(2000);
+    double volume = 0.0;
+    for (int j = 0; j < steps; ++j) {
+        Vec3 box[3];
+        context.getPeriodicBoxVectors(box[0], box[1], box[2]);
+        volume += box[0][0]*box[1][1]*box[2][2];
+        integrator.step(frequency);
+    }
+    volume /= steps;
+    double density = numMolecules*18/(AVOGADRO*volume*1e-21);
+    ASSERT_USUALLY_EQUAL_TOL(1, density, 0.1);
+}
+
+int main() {
+    try {
+        testChangingBoxSize();
+        testIdealGas();
+        testRandomSeed();
+        testWater();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}
+
+

platforms/reference/src/ReferenceKernels.cpp

@@ -1764,18 +1764,11 @@ ReferenceApplyMonteCarloBarostatKernel::~ReferenceApplyMonteCarloBarostatKernel(
 }
 
 void ReferenceApplyMonteCarloBarostatKernel::initialize(const System& system, const MonteCarloBarostat& barostat) {
-    int numConstraints = system.getNumConstraints();
-    vector<pair<int, int> > constraints;
-    for (int i = 0; i < numConstraints; i++) {
-        int particle1, particle2;
-        double distance;
-        system.getConstraintParameters(i, particle1, particle2, distance);
-        constraints.push_back(make_pair(particle1, particle2));
-    }
-    this->barostat = new ReferenceMonteCarloBarostat(system.getNumParticles(), constraints);
 }
 
 void ReferenceApplyMonteCarloBarostatKernel::scaleCoordinates(ContextImpl& context, double scale) {
+    if (barostat == NULL)
+        barostat = new ReferenceMonteCarloBarostat(context.getSystem().getNumParticles(), context.getMolecules());
     RealOpenMM** posData = extractPositions(context);
     Vec3 box[3];
     context.getOwner().getPeriodicBoxVectors(box[0], box[1], box[2]);

platforms/reference/src/ReferenceKernels.h

@@ -877,7 +877,7 @@ private:
  */
 class ReferenceApplyMonteCarloBarostatKernel : public ApplyMonteCarloBarostatKernel {
 public:
-    ReferenceApplyMonteCarloBarostatKernel(std::string name, const Platform& platform) : ApplyMonteCarloBarostatKernel(name, platform) {
+    ReferenceApplyMonteCarloBarostatKernel(std::string name, const Platform& platform) : ApplyMonteCarloBarostatKernel(name, platform), barostat(NULL) {
     }
     ~ReferenceApplyMonteCarloBarostatKernel();
     /**

platforms/reference/src/SimTKReference/ReferenceMonteCarloBarostat.cpp

@@ -36,38 +36,10 @@ using namespace std;
 
   --------------------------------------------------------------------------------------- */
 
-ReferenceMonteCarloBarostat::ReferenceMonteCarloBarostat(int numAtoms, const vector<pair<int, int> >& constraints) {
+ReferenceMonteCarloBarostat::ReferenceMonteCarloBarostat(int numAtoms, const vector<vector<int> >& molecules) : molecules(molecules) {
     savedAtomPositions[0].resize(numAtoms);
     savedAtomPositions[1].resize(numAtoms);
     savedAtomPositions[2].resize(numAtoms);
-
-    // First make a list of every other atom to which each atom is connect by a constraint.
-
-    vector<vector<int> > atomBonds(numAtoms);
-    for (int i = 0; i < (int) constraints.size(); i++) {
-        atomBonds[constraints[i].first].push_back(constraints[i].second);
-        atomBonds[constraints[i].second].push_back(constraints[i].first);
-    }
-
-    // Now tag atoms by which cluster they belong to.
-
-    vector<int> atomCluster(numAtoms, -1);
-    int numClusters = 0;
-    for (int i = 0; i < numAtoms; i++)
-        if (atomCluster[i] == -1)
-            tagAtomsInCluster(i, numClusters++, atomCluster, atomBonds);
-    clusters.resize(numClusters);
-    for (int i = 0; i < numAtoms; i++)
-        clusters[atomCluster[i]].push_back(i);
-}
-
-void ReferenceMonteCarloBarostat::tagAtomsInCluster(int atom, int cluster, vector<int>& atomCluster, vector<vector<int> >& atomBonds) {
-    // Recursively tag atoms as belonging to a particular cluster.
-
-    atomCluster[atom] = cluster;
-    for (int i = 0; i < (int) atomBonds[atom].size(); i++)
-        if (atomCluster[atomBonds[atom][i]] == -1)
-            tagAtomsInCluster(atomBonds[atom][i], cluster, atomCluster, atomBonds);
 }
 
 /**---------------------------------------------------------------------------------------
@@ -95,21 +67,21 @@ void ReferenceMonteCarloBarostat::applyBarostat(RealOpenMM** atomPositions, Real
         for (int j = 0; j < 3; j++)
             savedAtomPositions[j][i] = atomPositions[i][j];
 
-    // Loop over clusters.
+    // Loop over molecules.
 
-    for (int i = 0; i < (int) clusters.size(); i++) {
-        // Find the cluster center.
+    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) clusters[i].size(); j++) {
-            RealOpenMM* atomPos = atomPositions[clusters[i][j]];
+        for (int j = 0; j < (int) molecules[i].size(); j++) {
+            RealOpenMM* atomPos = atomPositions[molecules[i][j]];
             pos[0] += atomPos[0];
             pos[1] += atomPos[1];
             pos[2] += atomPos[2];
         }
-        pos[0] /= clusters[i].size();
-        pos[1] /= clusters[i].size();
-        pos[2] /= clusters[i].size();
+        pos[0] /= molecules[i].size();
+        pos[1] /= molecules[i].size();
+        pos[2] /= molecules[i].size();
 
         // Move it into the first periodic box.
 
@@ -122,20 +94,14 @@ void ReferenceMonteCarloBarostat::applyBarostat(RealOpenMM** atomPositions, Real
         pos[0] -= dx;
         pos[1] -= dy;
         pos[2] -= dz;
-        for (int j = 0; j < (int) clusters[i].size(); j++) {
-            RealOpenMM* atomPos = atomPositions[clusters[i][j]];
-            atomPos[0] -= dx;
-            atomPos[1] -= dy;
-            atomPos[2] -= dz;
-        }
 
-        // Now scale the position of the cluster center.
+        // Now scale the position of the molecule center.
 
-        dx = pos[0]*(scale-1);
-        dy = pos[1]*(scale-1);
-        dz = pos[2]*(scale-1);
-        for (int j = 0; j < (int) clusters[i].size(); j++) {
-            RealOpenMM* atomPos = atomPositions[clusters[i][j]];
+        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++) {
+            RealOpenMM* atomPos = atomPositions[molecules[i][j]];
             atomPos[0] += dx;
             atomPos[1] += dy;
             atomPos[2] += dz;

platforms/reference/src/SimTKReference/ReferenceMonteCarloBarostat.h

@@ -36,9 +36,7 @@ class ReferenceMonteCarloBarostat {
    private:
 
        std::vector<RealOpenMM> savedAtomPositions[3];
-       std::vector<std::vector<int> > clusters;
-       
-       void tagAtomsInCluster(int atom, int cluster, std::vector<int>& atomCluster, std::vector<std::vector<int> >& atomBonds);
+       std::vector<std::vector<int> > molecules;
 
    public:
 
@@ -48,7 +46,7 @@ class ReferenceMonteCarloBarostat {
 
          --------------------------------------------------------------------------------------- */
 
-       ReferenceMonteCarloBarostat(int numAtoms, const std::vector<std::pair<int, int> >& constraints);
+       ReferenceMonteCarloBarostat(int numAtoms, const std::vector<std::vector<int> >& molecules);
 
       /**---------------------------------------------------------------------------------------
 

platforms/reference/tests/TestReferenceMonteCarloBarostat.cpp

@@ -72,7 +72,7 @@ void testIdealGas() {
     const int frequency = 10;
     const int steps = 1000;
     const double pressure = 1.5;
-    const double pressureInMD = pressure*AVOGADRO*1e-25;
+    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);