Created MonteCarloBarostat, including reference implementation

olla/include/openmm/kernels.h

@@ -48,6 +48,7 @@
 #include "openmm/HarmonicBondForce.h"
 #include "openmm/KernelImpl.h"
 #include "openmm/LangevinIntegrator.h"
+#include "openmm/MonteCarloBarostat.h"
 #include "openmm/PeriodicTorsionForce.h"
 #include "openmm/RBTorsionForce.h"
 #include "openmm/NonbondedForce.h"
@@ -802,6 +803,42 @@ public:
     virtual void execute(ContextImpl& context) = 0;
 };
 
+/**
+ * This kernel is invoked by MonteCarloBarostat to adjust the periodic box volume
+ */
+class ApplyMonteCarloBarostatKernel : public KernelImpl {
+public:
+    static std::string Name() {
+        return "ApplyMonteCarloBarostat";
+    }
+    ApplyMonteCarloBarostatKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
+    }
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param barostat   the 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;
+    /**
+     * Reject the most recent Monte Carlo step, restoring the particle positions to where they were before
+     * scaleCoordinates() was last called.
+     *
+     * @param context    the context in which to execute this kernel
+     */
+    virtual void restoreCoordinates(ContextImpl& context) = 0;
+};
+
 /**
  * This kernel is invoked to calculate the kinetic energy of the system.
  */

openmmapi/include/openmm/MonteCarloBarostat.h

+#ifndef OPENMM_MONTECARLOBAROSTAT_H_
+#define OPENMM_MONTECARLOBAROSTAT_H_
+
+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2010 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "Force.h"
+#include <string>
+#include "internal/windowsExport.h"
+
+namespace OpenMM {
+
+/**
+ * This class uses a Monte Carlo algorithm to adjust the size of the periodic box, simulating the
+ * effect of constant pressure.
+ *
+ * This class assumes the simulation is also being run at constant temperature, and requires you
+ * to specify the system temperature (since it affects the acceptance probability for Monte Carlo
+ * moves).  It does not actually perform temperature regulation, however.  You must use another
+ * mechanism along with it to maintain the temperature, such as LangevinIntegrator or AndersenThermostat.
+ */
+
+class OPENMM_EXPORT MonteCarloBarostat : 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;
+    }
+    /**
+     * Create a MonteCarloBarostat.
+     *
+     * @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)
+     */
+    MonteCarloBarostat(double defaultPressure, double temperature, int frequency = 25);
+    /**
+     * Get the default pressure acting on the system (in bar).
+     *
+     * @return the default pressure acting on the system, measured in bar.
+     */
+    double getDefaultPressure() const {
+        return defaultPressure;
+    }
+    /**
+     * Get the frequency (in time steps) at which Monte Carlo pressure changes should be attempted.
+     */
+    int getFrequency() const {
+        return frequency;
+    }
+    /**
+     * Set the frequency (in time steps) at which Monte Carlo pressure changes should be attempted.
+     */
+    void setFrequency(int freq) {
+        frequency = freq;
+    }
+    /**
+     * Get the temperature at which the system is being maintained, measured in Kelvin.
+     */
+    double getTemperature() {
+        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();
+private:
+    double defaultPressure, temperature;
+    int frequency, randomNumberSeed;
+
+        double GetTemperature() const {
+            return temperature;
+        }
+
+        void SetTemperature(double temperature) {
+            this->temperature = temperature;
+        }
+    };
+
+} // namespace OpenMM
+
+#endif /*OPENMM_MONTECARLOBAROSTAT_H_*/

openmmapi/include/openmm/Vec3.h

@@ -92,6 +92,12 @@ public:
         return Vec3(lhs[0] - rhs[0], lhs[1] - rhs[1], lhs[2] - rhs[2]);
     }
 
+    // scalar product
+    Vec3 operator*(double rhs) const {
+        const Vec3& lhs = *this;
+        return Vec3(lhs[0]*rhs, lhs[1]*rhs, lhs[2]*rhs);
+    }
+    
     // dot product
     double dot(const Vec3& rhs) const {
         const Vec3& lhs = *this;

openmmapi/include/openmm/internal/MonteCarloBarostatImpl.h

+#ifndef OPENMM_MONTECARLOBAROSTATIMPL_H_
+#define OPENMM_MONTECARLOBAROSTATIMPL_H_
+
+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2010 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "ForceImpl.h"
+#include "openmm/MonteCarloBarostat.h"
+#include "openmm/Kernel.h"
+#include "sfmt/SFMT.h"
+#include <string>
+
+namespace OpenMM {
+
+/**
+ * This is the internal implementation of MonteCarloBarostat.
+ */
+
+class MonteCarloBarostatImpl : public ForceImpl {
+public:
+    MonteCarloBarostatImpl(MonteCarloBarostat& owner);
+    void initialize(ContextImpl& context);
+    MonteCarloBarostat& getOwner() {
+        return owner;
+    }
+    void updateContextState(ContextImpl& context);
+    void calcForces(ContextImpl& context) {
+        // This force doesn't apply forces to particles.
+    }
+    double calcEnergy(ContextImpl& context) {
+        return 0.0; // This force doesn't contribute to the potential energy.
+    }
+    std::map<std::string, double> getDefaultParameters();
+    std::vector<std::string> getKernelNames();
+private:
+    MonteCarloBarostat& owner;
+    int step;
+    double volumeScale;
+    OpenMM_SFMT::SFMT random;
+    Kernel kernel;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_MONTECARLOBAROSTATIMPL_H_*/

openmmapi/src/Context.cpp

@@ -37,18 +37,18 @@ using namespace OpenMM;
 using namespace std;
 
 Context::Context(System& system, Integrator& integrator) : properties(map<string, string>()) {
-    impl = new ContextImpl(*this, system, integrator, 0, properties);
     system.getDefaultPeriodicBoxVectors(periodicBoxVectors[0], periodicBoxVectors[1], periodicBoxVectors[2]);
+    impl = new ContextImpl(*this, system, integrator, 0, properties);
 }
 
 Context::Context(System& system, Integrator& integrator, Platform& platform) : properties(map<string, string>()) {
-    impl = new ContextImpl(*this, system, integrator, &platform, properties);
     system.getDefaultPeriodicBoxVectors(periodicBoxVectors[0], periodicBoxVectors[1], periodicBoxVectors[2]);
+    impl = new ContextImpl(*this, system, integrator, &platform, properties);
 }
 
 Context::Context(System& system, Integrator& integrator, Platform& platform, const map<string, string>& properties) : properties(properties) {
-    impl = new ContextImpl(*this, system, integrator, &platform, properties);
     system.getDefaultPeriodicBoxVectors(periodicBoxVectors[0], periodicBoxVectors[1], periodicBoxVectors[2]);
+    impl = new ContextImpl(*this, system, integrator, &platform, properties);
 }
 
 Context::~Context() {

openmmapi/src/MonteCarloBarostat.cpp

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

openmmapi/src/MonteCarloBarostatImpl.cpp

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

platforms/reference/src/ReferenceKernelFactory.cpp

@@ -84,6 +84,8 @@ KernelImpl* ReferenceKernelFactory::createKernelImpl(std::string name, const Pla
         return new ReferenceIntegrateVariableVerletStepKernel(name, platform, data);
     if (name == ApplyAndersenThermostatKernel::Name())
         return new ReferenceApplyAndersenThermostatKernel(name, platform);
+    if (name == ApplyMonteCarloBarostatKernel::Name())
+        return new ReferenceApplyMonteCarloBarostatKernel(name, platform);
     if (name == CalcKineticEnergyKernel::Name())
         return new ReferenceCalcKineticEnergyKernel(name, platform);
     if (name == RemoveCMMotionKernel::Name())

platforms/reference/src/ReferenceKernels.cpp

@@ -47,6 +47,7 @@
 #include "SimTKReference/ReferenceHarmonicBondIxn.h"
 #include "SimTKReference/ReferenceLJCoulomb14.h"
 #include "SimTKReference/ReferenceLJCoulombIxn.h"
+#include "SimTKReference/ReferenceMonteCarloBarostat.h"
 #include "SimTKReference/ReferenceProperDihedralBond.h"
 #include "SimTKReference/ReferenceRbDihedralBond.h"
 #include "SimTKReference/ReferenceStochasticDynamics.h"
@@ -1757,6 +1758,36 @@ void ReferenceApplyAndersenThermostatKernel::execute(ContextImpl& context) {
 			static_cast<RealOpenMM>(context.getIntegrator().getStepSize()) );
 }
 
+ReferenceApplyMonteCarloBarostatKernel::~ReferenceApplyMonteCarloBarostatKernel() {
+    if (barostat)
+        delete barostat;
+}
+
+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) {
+    RealOpenMM** posData = extractPositions(context);
+    Vec3 box[3];
+    context.getOwner().getPeriodicBoxVectors(box[0], box[1], box[2]);
+    RealOpenMM boxSize[] = {box[0][0], box[1][1], box[2][2]};
+    barostat->applyBarostat(posData, boxSize, scale);
+}
+
+void ReferenceApplyMonteCarloBarostatKernel::restoreCoordinates(ContextImpl& context) {
+    RealOpenMM** posData = extractPositions(context);
+    barostat->restorePositions(posData);
+}
+
 void ReferenceCalcKineticEnergyKernel::initialize(const System& system) {
     int numParticles = system.getNumParticles();
     masses.resize(numParticles);

platforms/reference/src/ReferenceKernels.h

@@ -45,6 +45,7 @@ class ReferenceCustomHbondIxn;
 class ReferenceBrownianDynamics;
 class ReferenceStochasticDynamics;
 class ReferenceConstraintAlgorithm;
+class ReferenceMonteCarloBarostat;
 class ReferenceVariableStochasticDynamics;
 class ReferenceVariableVerletDynamics;
 class ReferenceVerletDynamics;
@@ -871,6 +872,42 @@ private:
     RealOpenMM* masses;
 };
 
+/**
+ * This kernel is invoked by MonteCarloBarostat to adjust the periodic box volume
+ */
+class ReferenceApplyMonteCarloBarostatKernel : public ApplyMonteCarloBarostatKernel {
+public:
+    ReferenceApplyMonteCarloBarostatKernel(std::string name, const Platform& platform) : ApplyMonteCarloBarostatKernel(name, platform) {
+    }
+    ~ReferenceApplyMonteCarloBarostatKernel();
+    /**
+     * 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:
+    ReferenceMonteCarloBarostat* barostat;
+};
+
 /**
  * This kernel is invoked to calculate the kinetic energy of the system.
  */

platforms/reference/src/ReferencePlatform.cpp

@@ -62,6 +62,7 @@ ReferencePlatform::ReferencePlatform() {
     registerKernelFactory(IntegrateVariableLangevinStepKernel::Name(), factory);
     registerKernelFactory(IntegrateVariableVerletStepKernel::Name(), factory);
     registerKernelFactory(ApplyAndersenThermostatKernel::Name(), factory);
+    registerKernelFactory(ApplyMonteCarloBarostatKernel::Name(), factory);
     registerKernelFactory(CalcKineticEnergyKernel::Name(), factory);
     registerKernelFactory(RemoveCMMotionKernel::Name(), factory);
 }

platforms/reference/src/SimTKReference/ReferenceMonteCarloBarostat.cpp

+
+/* Portions copyright (c) 2010 Stanford University and Simbios.
+ * Contributors: Peter Eastman
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject
+ * to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+ * IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
+ * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+ * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+
+#include <string.h>
+#include <sstream>
+
+#include "../SimTKUtilities/SimTKOpenMMUtilities.h"
+#include "ReferenceMonteCarloBarostat.h"
+
+using namespace std;
+
+/**---------------------------------------------------------------------------------------
+
+  Constructor
+
+  --------------------------------------------------------------------------------------- */
+
+ReferenceMonteCarloBarostat::ReferenceMonteCarloBarostat(int numAtoms, const vector<pair<int, int> >& constraints) {
+    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);
+}
+
+/**---------------------------------------------------------------------------------------
+
+  Destructor
+
+  --------------------------------------------------------------------------------------- */
+
+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 ReferenceMonteCarloBarostat::applyBarostat(RealOpenMM** atomPositions, RealOpenMM* 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 clusters.
+
+    for (int i = 0; i < (int) clusters.size(); i++) {
+        // Find the cluster center.
+
+        RealOpenMM pos[3] = {0, 0, 0};
+        for (int j = 0; j < (int) clusters[i].size(); j++) {
+            RealOpenMM* atomPos = atomPositions[clusters[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();
+
+        // 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]);
+        float dx = xcell*boxSize[0];
+        float dy = ycell*boxSize[1];
+        float dz = zcell*boxSize[2];
+        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.
+
+        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]];
+            atomPos[0] += dx;
+            atomPos[1] += dy;
+            atomPos[2] += dz;
+        }
+    }
+}
+
+/**---------------------------------------------------------------------------------------
+
+  Restore atom positions to what they were before applyBarostat() was called.
+
+  @param atomPositions      atom positions
+
+  --------------------------------------------------------------------------------------- */
+
+void ReferenceMonteCarloBarostat::restorePositions(RealOpenMM** atomPositions) {
+    int numAtoms = savedAtomPositions[0].size();
+    for (int i = 0; i < numAtoms; i++)
+        for (int j = 0; j < 3; j++)
+            atomPositions[i][j] = savedAtomPositions[j][i];
+}

platforms/reference/src/SimTKReference/ReferenceMonteCarloBarostat.h

+
+/* Portions copyright (c) 2010 Stanford University and Simbios.
+ * Contributors: Peter Eastman
+ *
+ * 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.
+ */
+
+#ifndef __ReferenceMonteCarloBarostat_H__
+#define __ReferenceMonteCarloBarostat_H__
+
+#include "../SimTKUtilities/SimTKOpenMMCommon.h"
+#include <utility>
+#include <vector>
+
+// ---------------------------------------------------------------------------------------
+
+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);
+
+   public:
+
+      /**---------------------------------------------------------------------------------------
+
+         Constructor
+
+         --------------------------------------------------------------------------------------- */
+
+       ReferenceMonteCarloBarostat(int numAtoms, const std::vector<std::pair<int, int> >& constraints);
+
+      /**---------------------------------------------------------------------------------------
+
+         Destructor
+
+         --------------------------------------------------------------------------------------- */
+
+       ~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(RealOpenMM** atomPositions, RealOpenMM* boxSize, RealOpenMM scale);
+
+      /**---------------------------------------------------------------------------------------
+
+         Restore atom positions to what they were before applyBarostat() was called.
+
+         @param atomPositions      atom positions
+
+         --------------------------------------------------------------------------------------- */
+
+      void restorePositions(RealOpenMM** atomPositions);
+
+};
+
+// ---------------------------------------------------------------------------------------
+
+#endif // __ReferenceMonteCarloBarostat_H__

platforms/reference/tests/TestReferenceMonteCarloBarostat.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 reference implementation of MonteCarloBarostat.
+ */
+
+#include "../../../tests/AssertionUtilities.h"
+#include "openmm/MonteCarloBarostat.h"
+#include "openmm/Context.h"
+#include "ReferencePlatform.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() {
+    ReferencePlatform platform;
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(4, 0, 0), Vec3(0, 5, 0), Vec3(0, 0, 6));
+    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.
+
+    ReferencePlatform 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;
+    ReferencePlatform 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]);
+        }
+    }
+}
+
+int main() {
+    try {
+        testChangingBoxSize();
+        testIdealGas();
+        testRandomSeed();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}
+