Created VariableLangevinIntegrator

olla/include/openmm/kernels.h

@@ -32,7 +32,6 @@
  * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
  * -------------------------------------------------------------------------- */
 
-#include "openmm/VariableVerletIntegrator.h"
 #include "openmm/AndersenThermostat.h"
 #include "openmm/BrownianIntegrator.h"
 #include "openmm/CMMotionRemover.h"
@@ -47,6 +46,8 @@
 #include "openmm/NonbondedForce.h"
 #include "openmm/Stream.h"
 #include "openmm/System.h"
+#include "openmm/VariableLangevinIntegrator.h"
+#include "openmm/VariableVerletIntegrator.h"
 #include "openmm/VerletIntegrator.h"
 #include <set>
 #include <string>
@@ -419,6 +420,33 @@ public:
     virtual void execute(OpenMMContextImpl& context, const BrownianIntegrator& integrator) = 0;
 };
 
+/**
+ * This kernel is invoked by VariableLangevinIntegrator to take one time step.
+ */
+class IntegrateVariableLangevinStepKernel : public KernelImpl {
+public:
+    static std::string Name() {
+        return "IntegrateVariableLangevinStep";
+    }
+    IntegrateVariableLangevinStepKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
+    }
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param integrator the VariableLangevinIntegrator this kernel will be used for
+     */
+    virtual void initialize(const System& system, const VariableLangevinIntegrator& integrator) = 0;
+    /**
+     * Execute the kernel.
+     *
+     * @param context    the context in which to execute this kernel
+     * @param integrator the LangevinIntegrator this kernel is being used for
+     * @param maxTime    the maximum time beyond which the simulation should not be advanced
+     */
+    virtual void execute(OpenMMContextImpl& context, const VariableLangevinIntegrator& integrator, double maxTime) = 0;
+};
+
 /**
  * This kernel is invoked by VariableVerletIntegrator to take one time step.
  */

openmmapi/include/OpenMM.h

@@ -32,7 +32,6 @@
  * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
  * -------------------------------------------------------------------------- */
 
-#include "openmm/VariableVerletIntegrator.h"
 #include "openmm/AndersenThermostat.h"
 #include "openmm/BrownianIntegrator.h"
 #include "openmm/CMMotionRemover.h"
@@ -51,6 +50,8 @@
 #include "openmm/State.h"
 #include "openmm/System.h"
 #include "openmm/Units.h"
+#include "openmm/VariableLangevinIntegrator.h"
+#include "openmm/VariableVerletIntegrator.h"
 #include "openmm/Vec3.h"
 #include "openmm/VerletIntegrator.h"
 #include "openmm/Platform.h"

openmmapi/include/openmm/BrownianIntegrator.h

@@ -66,14 +66,14 @@ public:
     }
     /**
      * Get the friction coefficient which determines how strongly the system is coupled to
-     * the heat bath.
+     * the heat bath (in inverse ps).
      */
     double getFriction() const {
         return friction;
     }
     /**
      * Set the friction coefficient which determines how strongly the system is coupled to
-     * the heat bath.
+     * the heat bath (in inverse ps).
      */
     void setFriction(double coeff) {
         friction = coeff;

openmmapi/include/openmm/LangevinIntegrator.h

@@ -66,14 +66,14 @@ public:
     }
     /**
      * Get the friction coefficient which determines how strongly the system is coupled to
-     * the heat bath.
+     * the heat bath (in inverse ps).
      */
     double getFriction() const {
         return friction;
     }
     /**
      * Set the friction coefficient which determines how strongly the system is coupled to
-     * the heat bath.
+     * the heat bath (in inverse ps).
      */
     void setFriction(double coeff) {
         friction = coeff;

openmmapi/include/openmm/VariableLangevinIntegrator.h

+#ifndef OPENMM_VARIABLELANGEVININTEGRATOR_H_
+#define OPENMM_VARIABLELANGEVININTEGRATOR_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) 2008-2009 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 "Integrator.h"
+#include "openmm/Kernel.h"
+#include "internal/windowsExport.h"
+
+namespace OpenMM {
+
+/**
+ * This is an error contolled, variable time step Integrator that simulates a System using Langevin
+ * dynamics.  It compares the result of the Langevin integrator to that of an
+ * explicit Euler integrator, takes the difference between the two as a measure of the integration
+ * error in each time step, and continuously adjusts the step size to keep the error below a
+ * specified tolerance.  This both improves the stability of the integrator and allows it to take
+ * larger steps on average, while still maintaining comparable accuracy to a fixed step size integrator.
+ *
+ * It is best not to think of the error tolerance as having any absolute meaning.  It is just an
+ * adjustable parameter that affects the step size and integration accuracy.  You
+ * should try different values to find the largest one that produces a trajectory sufficiently
+ * accurate for your purposes.  0.001 is often a good starting point.
+ */
+
+class OPENMM_EXPORT VariableLangevinIntegrator : public Integrator {
+public:
+    /**
+     * Create a VariableLangevinIntegrator.
+     *
+     * @param temperature    the temperature of the heat bath (in Kelvin)
+     * @param frictionCoeff  the friction coefficient which couples the system to the heat bath (in inverse picoseconds)
+     * @param errorTol       the error tolerance
+     */
+    VariableLangevinIntegrator(double temperature, double frictionCoeff, double errorTol);
+    /**
+     * Get the temperature of the heat bath (in Kelvin).
+     */
+    double getTemperature() const {
+        return temperature;
+    }
+    /**
+     * Set the temperature of the heat bath (in Kelvin).
+     */
+    void setTemperature(double temp) {
+        temperature = temp;
+    }
+    /**
+     * Get the friction coefficient which determines how strongly the system is coupled to
+     * the heat bath (in inverse ps).
+     */
+    double getFriction() const {
+        return friction;
+    }
+    /**
+     * Set the friction coefficient which determines how strongly the system is coupled to
+     * the heat bath (in inverse ps).
+     */
+    void setFriction(double coeff) {
+        friction = coeff;
+    }
+    /**
+     * Get the error tolerance.
+     */
+    double getErrorTolerance() const {
+        return errorTol;
+    }
+    /**
+     * Set the error tolerance.
+     */
+    void setErrorTolerance(double tol) {
+        errorTol = tol;
+    }
+    /**
+     * Get the random number seed.  See setRandomNumberSeed() for details.
+     */
+    int getRandomNumberSeed() const {
+        return randomNumberSeed;
+    }
+    /**
+     * Set the random number seed.  The precise meaning of this parameter is undefined, and is left up
+     * to each Platform to interpret in an appropriate way.  It is guaranteed that if two simulations
+     * are run with different random number seeds, the sequence of random forces 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;
+    }
+    /**
+     * Advance a simulation through time by taking a series of time steps.
+     *
+     * @param steps   the number of time steps to take
+     */
+    void step(int steps);
+    /**
+     * Advance a simulation through time by taking a series of steps until a specified time is
+     * reached.  When this method returns, the simulation time will exactly equal the time which
+     * was specified.  If you call this method and specify a time that is earlier than the
+     * current time, it will return without doing anything.
+     *
+     * @param time   the time to which the simulation should be advanced
+     */
+    void stepTo(double time);
+protected:
+    /**
+     * This will be called by the OpenMMContext when it is created.  It informs the Integrator
+     * of what context it will be integrating, and gives it a chance to do any necessary initialization.
+     * It will also get called again if the application calls reinitialize() on the OpenMMContext.
+     */
+    void initialize(OpenMMContextImpl& context);
+    /**
+     * Get the names of all Kernels used by this Integrator.
+     */
+    std::vector<std::string> getKernelNames();
+private:
+    double temperature, friction, errorTol;
+    int randomNumberSeed;
+    OpenMMContextImpl* context;
+    Kernel kernel;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_VARIABLELANGEVININTEGRATOR_H_*/

openmmapi/include/openmm/VariableVerletIntegrator.h

@@ -63,7 +63,7 @@ public:
     /**
      * Create a VariableVerletIntegrator.
      *
-     * @param tol      the error tolerance
+     * @param errorTol    the error tolerance
      */
     VariableVerletIntegrator(double errorTol);
     /**

openmmapi/src/VariableLangevinIntegrator.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-2009 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/VariableLangevinIntegrator.h"
+#include "openmm/OpenMMContext.h"
+#include "openmm/internal/OpenMMContextImpl.h"
+#include "openmm/kernels.h"
+#include <limits>
+#include <string>
+
+using namespace OpenMM;
+using std::string;
+using std::vector;
+
+VariableLangevinIntegrator::VariableLangevinIntegrator(double temperature, double frictionCoeff, double errorTol) {
+    setTemperature(temperature);
+    setFriction(frictionCoeff);
+    setErrorTolerance(errorTol);
+    setConstraintTolerance(1e-4);
+    setRandomNumberSeed((int) time(NULL));
+}
+
+void VariableLangevinIntegrator::initialize(OpenMMContextImpl& contextRef) {
+    context = &contextRef;
+    kernel = context->getPlatform().createKernel(IntegrateVariableLangevinStepKernel::Name(), contextRef);
+    dynamic_cast<IntegrateVariableLangevinStepKernel&>(kernel.getImpl()).initialize(contextRef.getSystem(), *this);
+}
+
+vector<string> VariableLangevinIntegrator::getKernelNames() {
+    std::vector<std::string> names;
+    names.push_back(IntegrateVariableLangevinStepKernel::Name());
+    return names;
+}
+
+void VariableLangevinIntegrator::step(int steps) {
+    for (int i = 0; i < steps; ++i) {
+        context->updateContextState();
+        context->calcForces();
+        dynamic_cast<IntegrateVariableLangevinStepKernel&>(kernel.getImpl()).execute(*context, *this, std::numeric_limits<double>::infinity());
+    }
+}
+
+void VariableLangevinIntegrator::stepTo(double time) {
+    while (time > context->getTime()) {
+        context->updateContextState();
+        context->calcForces();
+        dynamic_cast<IntegrateVariableLangevinStepKernel&>(kernel.getImpl()).execute(*context, *this, time);
+    }
+}

platforms/reference/src/ReferenceKernelFactory.cpp

@@ -64,6 +64,8 @@ KernelImpl* ReferenceKernelFactory::createKernelImpl(std::string name, const Pla
         return new ReferenceIntegrateLangevinStepKernel(name, platform, data);
     if (name == IntegrateBrownianStepKernel::Name())
         return new ReferenceIntegrateBrownianStepKernel(name, platform, data);
+    if (name == IntegrateVariableLangevinStepKernel::Name())
+        return new ReferenceIntegrateVariableLangevinStepKernel(name, platform, data);
     if (name == IntegrateVariableVerletStepKernel::Name())
         return new ReferenceIntegrateVariableVerletStepKernel(name, platform, data);
     if (name == ApplyAndersenThermostatKernel::Name())

platforms/reference/src/ReferenceKernels.cpp

@@ -33,18 +33,19 @@
 #include "ReferenceFloatStreamImpl.h"
 #include "gbsa/CpuObc.h"
 #include "gbsa/CpuGBVI.h"
-#include "SimTKReference/ReferenceVariableVerletDynamics.h"
 #include "SimTKReference/ReferenceAndersenThermostat.h"
 #include "SimTKReference/ReferenceAngleBondIxn.h"
 #include "SimTKReference/ReferenceBondForce.h"
 #include "SimTKReference/ReferenceBrownianDynamics.h"
+#include "SimTKReference/ReferenceCCMAAlgorithm.h"
 #include "SimTKReference/ReferenceHarmonicBondIxn.h"
 #include "SimTKReference/ReferenceLJCoulomb14.h"
 #include "SimTKReference/ReferenceLJCoulombIxn.h"
 #include "SimTKReference/ReferenceProperDihedralBond.h"
 #include "SimTKReference/ReferenceRbDihedralBond.h"
 #include "SimTKReference/ReferenceStochasticDynamics.h"
-#include "SimTKReference/ReferenceCCMAAlgorithm.h"
+#include "SimTKReference/ReferenceVariableStochasticDynamics.h"
+#include "SimTKReference/ReferenceVariableVerletDynamics.h"
 #include "SimTKReference/ReferenceVerletDynamics.h"
 #include "openmm/CMMotionRemover.h"
 #include "openmm/System.h"
@@ -787,6 +788,70 @@ void ReferenceIntegrateBrownianStepKernel::execute(OpenMMContextImpl& context, c
     data.stepCount++;
 }
 
+ReferenceIntegrateVariableLangevinStepKernel::~ReferenceIntegrateVariableLangevinStepKernel() {
+    if (dynamics)
+        delete dynamics;
+    if (constraints)
+        delete constraints;
+    if (masses)
+        delete[] masses;
+    if (constraintIndices)
+        disposeIntArray(constraintIndices, numConstraints);
+    if (constraintDistances)
+        delete[] constraintDistances;
+}
+
+void ReferenceIntegrateVariableLangevinStepKernel::initialize(const System& system, const VariableLangevinIntegrator& integrator) {
+    int numParticles = system.getNumParticles();
+    masses = new RealOpenMM[numParticles];
+    for (int i = 0; i < numParticles; ++i)
+        masses[i] = static_cast<RealOpenMM>(system.getParticleMass(i));
+    numConstraints = system.getNumConstraints();
+    constraintIndices = allocateIntArray(numConstraints, 2);
+    constraintDistances = new RealOpenMM[numConstraints];
+    for (int i = 0; i < numConstraints; ++i) {
+        int particle1, particle2;
+        double distance;
+        system.getConstraintParameters(i, particle1, particle2, distance);
+        constraintIndices[i][0] = particle1;
+        constraintIndices[i][1] = particle2;
+        constraintDistances[i] = static_cast<RealOpenMM>(distance);
+    }
+    SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) integrator.getRandomNumberSeed());
+}
+
+void ReferenceIntegrateVariableLangevinStepKernel::execute(OpenMMContextImpl& context, const VariableLangevinIntegrator& integrator, double maxTime) {
+    double temperature = integrator.getTemperature();
+    double friction = integrator.getFriction();
+    double errorTol = integrator.getErrorTolerance();
+    RealOpenMM** posData = ((ReferenceFloatStreamImpl&) context.getPositions().getImpl()).getData();
+    RealOpenMM** velData = ((ReferenceFloatStreamImpl&) context.getVelocities().getImpl()).getData();
+    RealOpenMM** forceData = const_cast<RealOpenMM**>(((ReferenceFloatStreamImpl&) context.getForces().getImpl()).getData()); // Reference code needs to be made const correct
+    if (dynamics == 0 || temperature != prevTemp || friction != prevFriction || errorTol != prevErrorTol) {
+        // Recreate the computation objects with the new parameters.
+
+        if (dynamics) {
+            delete dynamics;
+            delete constraints;
+        }
+        RealOpenMM tau = static_cast<RealOpenMM>( friction == 0.0 ? 0.0 : 1.0/friction );
+        dynamics = new ReferenceVariableStochasticDynamics(context.getSystem().getNumParticles(), (RealOpenMM) tau, (RealOpenMM) temperature, (RealOpenMM) errorTol);
+        vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
+        findAnglesForCCMA(context.getSystem(), angles);
+        constraints = new ReferenceCCMAAlgorithm(context.getSystem().getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, (RealOpenMM)integrator.getConstraintTolerance());
+        dynamics->setReferenceConstraintAlgorithm(constraints);
+        prevTemp = temperature;
+        prevFriction = friction;
+        prevErrorTol = errorTol;
+    }
+    RealOpenMM maxStepSize = (RealOpenMM) (maxTime-data.time);
+    dynamics->update(context.getSystem().getNumParticles(), posData, velData, forceData, masses, maxStepSize);
+    data.time += dynamics->getDeltaT();
+    if (dynamics->getDeltaT() == maxStepSize)
+        data.time = maxTime; // Avoid round-off error
+    data.stepCount++;
+}
+
 ReferenceIntegrateVariableVerletStepKernel::~ReferenceIntegrateVariableVerletStepKernel() {
     if (dynamics)
         delete dynamics;

platforms/reference/src/ReferenceKernels.h

@@ -43,6 +43,7 @@ class ReferenceAndersenThermostat;
 class ReferenceBrownianDynamics;
 class ReferenceStochasticDynamics;
 class ReferenceConstraintAlgorithm;
+class ReferenceVariableStochasticDynamics;
 class ReferenceVariableVerletDynamics;
 class ReferenceVerletDynamics;
 
@@ -442,6 +443,41 @@ private:
     double prevTemp, prevFriction, prevStepSize;
 };
 
+/**
+ * This kernel is invoked by VariableLangevinIntegrator to take one time step.
+ */
+class ReferenceIntegrateVariableLangevinStepKernel : public IntegrateVariableLangevinStepKernel {
+public:
+    ReferenceIntegrateVariableLangevinStepKernel(std::string name, const Platform& platform, ReferencePlatform::PlatformData& data) : IntegrateVariableLangevinStepKernel(name, platform),
+        data(data), dynamics(0), constraints(0), masses(0), constraintDistances(0), constraintIndices(0) {
+    }
+    ~ReferenceIntegrateVariableLangevinStepKernel();
+    /**
+     * Initialize the kernel.
+     *
+     * @param system     the System this kernel will be applied to
+     * @param integrator the LangevinIntegrator this kernel will be used for
+     */
+    void initialize(const System& system, const VariableLangevinIntegrator& integrator);
+    /**
+     * Execute the kernel.
+     *
+     * @param context    the context in which to execute this kernel
+     * @param integrator the LangevinIntegrator this kernel is being used for
+     * @param maxTime    the maximum time beyond which the simulation should not be advanced
+     */
+    void execute(OpenMMContextImpl& context, const VariableLangevinIntegrator& integrator, double maxTime);
+private:
+    ReferencePlatform::PlatformData& data;
+    ReferenceVariableStochasticDynamics* dynamics;
+    ReferenceConstraintAlgorithm* constraints;
+    RealOpenMM* masses;
+    RealOpenMM* constraintDistances;
+    int** constraintIndices;
+    int numConstraints;
+    double prevTemp, prevFriction, prevErrorTol;
+};
+
 /**
  * This kernel is invoked by VariableVerletIntegrator to take one time step.
  */

platforms/reference/src/ReferencePlatform.cpp

@@ -51,6 +51,7 @@ ReferencePlatform::ReferencePlatform() {
     registerKernelFactory(IntegrateVerletStepKernel::Name(), factory);
     registerKernelFactory(IntegrateLangevinStepKernel::Name(), factory);
     registerKernelFactory(IntegrateBrownianStepKernel::Name(), factory);
+    registerKernelFactory(IntegrateVariableLangevinStepKernel::Name(), factory);
     registerKernelFactory(IntegrateVariableVerletStepKernel::Name(), factory);
     registerKernelFactory(ApplyAndersenThermostatKernel::Name(), factory);
     registerKernelFactory(CalcKineticEnergyKernel::Name(), factory);

platforms/reference/src/SimTKReference/ReferenceStochasticDynamics.cpp

@@ -290,11 +290,11 @@ int ReferenceStochasticDynamics::updatePart1( int numberOfAtoms, RealOpenMM** at
 
            vVector[ii][jj]     = sqrtInvMass*fixedParameters[V]*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
           
-           velocities[ii][jj]  = oldVelocities[ii][jj]*fixedParameters[EM] +
+           RealOpenMM vPrime   = oldVelocities[ii][jj]*fixedParameters[EM] +
                                  inverseMasses[ii]*forces[ii][jj]*tau*( one - fixedParameters[EM]) +
                                  vVector[ii][jj] - fixedParameters[EM]*Vmh;
             
-           xPrime[ii][jj]      = atomCoordinates[ii][jj] + velocities[ii][jj]*fix1;
+           xPrime[ii][jj]      = atomCoordinates[ii][jj] + vPrime*fix1;
       }
    }
 

platforms/reference/src/SimTKReference/ReferenceVariableStochasticDynamics.cpp

+
+/* Portions copyright (c) 2006 Stanford University and Simbios.
+ * Contributors: Pande Group
+ *
+ * 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 <cstring>
+#include <sstream>
+
+#include "../SimTKUtilities/SimTKOpenMMCommon.h"
+#include "../SimTKUtilities/SimTKOpenMMLog.h"
+#include "../SimTKUtilities/SimTKOpenMMUtilities.h"
+#include "ReferenceVariableStochasticDynamics.h"
+
+#include <cstdio>
+
+/**---------------------------------------------------------------------------------------
+
+   ReferenceVariableStochasticDynamics constructor
+
+   @param numberOfAtoms  number of atoms
+   @param deltaT         delta t for dynamics
+   @param tau            viscosity(?)
+   @param temperature    temperature
+   @param accuracy       required accuracy
+
+   --------------------------------------------------------------------------------------- */
+
+ReferenceVariableStochasticDynamics::ReferenceVariableStochasticDynamics( int numberOfAtoms,
+                                                          RealOpenMM tau, RealOpenMM temperature,
+                                                          RealOpenMM accuracy ) :
+           ReferenceDynamics(numberOfAtoms, 0.0f, temperature), _tau(tau), _accuracy(accuracy) {
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName      = "\nReferenceVariableStochasticDynamics::ReferenceVariableStochasticDynamics";
+
+   static const RealOpenMM zero       =  0.0;
+   static const RealOpenMM one        =  1.0;
+
+   // ---------------------------------------------------------------------------------------
+
+   // insure tau is not zero -- if it is print warning message
+
+   if( _tau == zero ){
+
+      std::stringstream message;
+      message << methodName;
+      message << " input tau value=" << tau << " is invalid -- setting to 1.";
+      SimTKOpenMMLog::printError( message );
+
+      _tau = one;
+
+   }
+
+   allocate2DArrays( numberOfAtoms, 3, Max2DArrays );
+   allocate1DArrays( numberOfAtoms, Max1DArrays );
+
+}
+
+/**---------------------------------------------------------------------------------------
+
+   ReferenceVariableStochasticDynamics destructor
+
+   --------------------------------------------------------------------------------------- */
+
+ReferenceVariableStochasticDynamics::~ReferenceVariableStochasticDynamics( ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nReferenceVariableStochasticDynamics::~ReferenceVariableStochasticDynamics";
+
+   // ---------------------------------------------------------------------------------------
+
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Get the required accuracy
+
+   @return accuracy
+
+ --------------------------------------------------------------------------------------- */
+
+RealOpenMM ReferenceVariableStochasticDynamics::getAccuracy( void ) const {
+    return _accuracy;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Set the required accuracy
+
+ --------------------------------------------------------------------------------------- */
+
+void ReferenceVariableStochasticDynamics::setAccuracy( RealOpenMM accuracy ) {
+    _accuracy = accuracy;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Set fixed parameters
+
+   @return ReferenceDynamics::DefaultReturn
+
+   --------------------------------------------------------------------------------------- */
+
+int ReferenceVariableStochasticDynamics::_setFixedParameters( RealOpenMM timeStep, RealOpenMM prevTimeStep ){
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName = "\nReferenceVariableStochasticDynamics::_setFixedParameters";
+
+   static const RealOpenMM zero       =  0.0;
+   static const RealOpenMM one        =  1.0;
+   static const RealOpenMM two        =  2.0;
+   static const RealOpenMM three      =  3.0;
+   static const RealOpenMM four       =  4.0;
+   static const RealOpenMM half       =  0.5;
+
+   // ---------------------------------------------------------------------------------------
+
+   _fixedParameters[GDT]      = timeStep/getTau();
+   _fixedParameters[EPH]      = EXP(  half*_fixedParameters[GDT] );
+   _fixedParameters[EMH]      = EXP( -half*_fixedParameters[GDT] );
+   _fixedParameters[EM]       = EXP(      -_fixedParameters[GDT] );
+   _fixedParameters[EM_V]     = EXP(      -half*(timeStep+prevTimeStep)/getTau() );
+   _fixedParameters[EP]       = EXP(       _fixedParameters[GDT] );
+
+   if( _fixedParameters[GDT] >= (RealOpenMM) 0.1 ){
+
+      RealOpenMM term1        = _fixedParameters[EPH] - one;
+                 term1       *= term1;
+      _fixedParameters[B]     = _fixedParameters[GDT]*(_fixedParameters[EP] - one) - four*term1;
+
+      _fixedParameters[C]     = _fixedParameters[GDT] - three + four*_fixedParameters[EMH] - _fixedParameters[EM];
+      _fixedParameters[D]     = two - _fixedParameters[EPH] - _fixedParameters[EMH];
+
+    } else {
+
+      // this has not been debugged
+
+      RealOpenMM term1        = half*_fixedParameters[GDT];
+      RealOpenMM term2        = term1*term1;
+      RealOpenMM term4        = term2*term2;
+
+      RealOpenMM third        = (RealOpenMM) ( 1.0/3.0 );
+      RealOpenMM o7_9         = (RealOpenMM) ( 7.0/9.0 );
+      RealOpenMM o1_12        = (RealOpenMM) ( 1.0/12.0 );
+      RealOpenMM o17_90       = (RealOpenMM) ( 17.0/90.0 );
+      RealOpenMM o7_30        = (RealOpenMM) ( 7.0/30.0 );
+      RealOpenMM o31_1260     = (RealOpenMM) ( 31.0/1260.0 );
+      RealOpenMM o_360        = (RealOpenMM) ( 1.0/360.0 );
+
+      _fixedParameters[B]     = term4*( third  + term1*( third + term1*( o17_90 + term1*o7_9 )));
+      _fixedParameters[C]     = term2*term1*( two*third + term1*( -half + term1*( o7_30 + term1*(-o1_12 + term1*o31_1260 ))));
+      _fixedParameters[D]     = term2*( -one + term2*(-o1_12 - term2*o_360));
+   }
+
+   RealOpenMM kT        = ((RealOpenMM) BOLTZ)*getTemperature();
+
+   _fixedParameters[V]  = SQRT( kT*( one - _fixedParameters[EM]) );
+   _fixedParameters[X]  = getTau()*SQRT( kT*_fixedParameters[C] );
+   _fixedParameters[Yv] = SQRT( kT*_fixedParameters[B]/_fixedParameters[C] );
+   _fixedParameters[Yx] = getTau()*SQRT( kT*_fixedParameters[B]/(one - _fixedParameters[EM]) );
+
+   return ReferenceDynamics::DefaultReturn;
+
+};
+
+/**---------------------------------------------------------------------------------------
+
+   Get tau
+
+   @return tau
+
+   --------------------------------------------------------------------------------------- */
+
+RealOpenMM ReferenceVariableStochasticDynamics::getTau( void ) const {
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName  = "\nReferenceVariableStochasticDynamics::getTau";
+
+   // ---------------------------------------------------------------------------------------
+
+   return _tau;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Get array of fixed parameters indexed by 'FixedParameters' enums
+
+   @return array
+
+   --------------------------------------------------------------------------------------- */
+
+const RealOpenMM* ReferenceVariableStochasticDynamics::getFixedParameters( void ) const {
+
+   // ---------------------------------------------------------------------------------------
+
+   // static const char* methodName  = "\nReferenceVariableStochasticDynamics::getFixedParameters";
+
+   // ---------------------------------------------------------------------------------------
+
+   return _fixedParameters;
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Print parameters
+
+   @param message             message
+
+   @return ReferenceDynamics::DefaultReturn
+
+   --------------------------------------------------------------------------------------- */
+
+int ReferenceVariableStochasticDynamics::printParameters( std::stringstream& message ) const {
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName  = "\nReferenceVariableStochasticDynamics::printParameters";
+   static const char* parameterNames[MaxFixedParameters] = { "gdt", "ep", "eph", "emh", "em", "B", "C", "D",
+                                                             "V", "X", "Yv", "Yx" };
+
+   // ---------------------------------------------------------------------------------------
+
+   // print parameters
+
+   ReferenceDynamics::printParameters( message );
+   message << " tau=" << getTau();
+   message << " T=" << getTemperature();
+   int cut = 3;
+   for( int ii = 0; ii < MaxFixedParameters; ii++ ){
+      message << " " << parameterNames[ii] << "=" << _fixedParameters[ii];
+      if( cut++ > 5 ){
+         cut = 0;
+         message << std::endl;
+      }
+   }
+
+   return ReferenceDynamics::DefaultReturn;
+
+}
+
+/**---------------------------------------------------------------------------------------
+
+   First SD update; based on code in update.c do_update_sd() Gromacs 3.1.4
+
+   @param numberOfAtoms       number of atoms
+   @param atomCoordinates     atom coordinates
+   @param velocities          velocities
+   @param forces              forces
+   @param masses              atom masses
+   @param inverseMasses       inverse atom masses
+   @param xPrime              xPrime
+   @param oldVelocities       previous velocities
+   @param xVector             xVector
+   @param vVector             vVector
+   @param maxStepSize         maximum time step
+
+   @return ReferenceDynamics::DefaultReturn
+
+   --------------------------------------------------------------------------------------- */
+
+int ReferenceVariableStochasticDynamics::updatePart1( int numberOfAtoms, RealOpenMM** atomCoordinates,
+                                              RealOpenMM** velocities,
+                                              RealOpenMM** forces, RealOpenMM* masses, RealOpenMM* inverseMasses,
+                                              RealOpenMM** xPrime, RealOpenMM** oldVelocities,
+                                              RealOpenMM** xVector, RealOpenMM** vVector,
+                                              RealOpenMM maxStepSize ){
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName  = "\nReferenceVariableStochasticDynamics::updatePart1";
+
+   static const RealOpenMM zero   =  0.0;
+   static const RealOpenMM one    =  1.0;
+   static int debug               = 0;
+
+   // ---------------------------------------------------------------------------------------
+
+
+   // first-time-through initialization
+
+   if( getTimeStep() == 0 ){
+
+      std::stringstream message;
+      message << methodName;
+      int errors = 0;
+
+      // invert masses
+
+      for( int ii = 0; ii < numberOfAtoms; ii++ ){
+         if( masses[ii] <= zero ){
+            message << "mass at atom index=" << ii << " (" << masses[ii] << ") is <= 0" << std::endl;
+            errors++;
+         } else {
+            inverseMasses[ii] = one/masses[ii];
+         }
+      }
+
+      // exit if errors
+
+      if( errors ){
+         SimTKOpenMMLog::printError( message );
+      }
+   }
+
+   // Select the step size to use
+    RealOpenMM error = zero;
+    for (int i = 0; i < numberOfAtoms; ++i) {
+        for (int j = 0; j < 3; ++j) {
+            RealOpenMM xerror = inverseMasses[i]*forces[i][j];
+            error += xerror*xerror;
+        }
+    }
+    error = SQRT(error/(numberOfAtoms*3));
+    RealOpenMM newStepSize = SQRT(getAccuracy()/error);
+    if (getDeltaT() > 0.0f)
+        newStepSize = std::min(newStepSize, getDeltaT()*2.0f); // For safety, limit how quickly dt can increase.
+    if (newStepSize > getDeltaT() && newStepSize < 1.2f*getDeltaT())
+        newStepSize = getDeltaT(); // Keeping dt constant between steps improves the behavior of the integrator.
+    if (newStepSize > maxStepSize)
+        newStepSize = maxStepSize;
+   _setFixedParameters(newStepSize, getDeltaT());
+    setDeltaT(newStepSize);
+ 
+   if( getTimeStep() == 0 ){
+
+      // Initialize xVector
+
+      const RealOpenMM*  fixedParameters = getFixedParameters();
+      for( int ii = 0; ii < numberOfAtoms; ii++ ){
+         RealOpenMM sqrtInverseMass = SQRT( inverseMasses[ii] )*fixedParameters[X];
+         for( int jj = 0; jj < 3; jj++ ){
+            xVector[ii][jj] = sqrtInverseMass*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
+         }
+      }
+   }
+
+    // perform first update
+
+   const RealOpenMM*  fixedParameters = getFixedParameters();
+   RealOpenMM   tau                   = getTau();
+   RealOpenMM   fix1                  = tau*(fixedParameters[EPH] - fixedParameters[EMH]);
+
+   for( int ii = 0; ii < numberOfAtoms; ii++ ){
+
+      RealOpenMM sqrtInvMass = SQRT( inverseMasses[ii] );
+
+      for( int jj = 0; jj < 3; jj++ ){
+
+         oldVelocities[ii][jj] = velocities[ii][jj];
+
+         RealOpenMM Vmh        = xVector[ii][jj]*fixedParameters[D]/(tau*fixedParameters[C]) +
+                                 sqrtInvMass*fixedParameters[Yv]*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
+
+           vVector[ii][jj]     = sqrtInvMass*fixedParameters[V]*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
+
+           RealOpenMM vPrime   = oldVelocities[ii][jj]*fixedParameters[EM_V] +
+                                 inverseMasses[ii]*forces[ii][jj]*tau*(one - fixedParameters[EM_V]) +
+                                 vVector[ii][jj] - fixedParameters[EM_V]*Vmh;
+
+           xPrime[ii][jj]      = atomCoordinates[ii][jj] + vPrime*fix1;
+      }
+   }
+
+   return ReferenceDynamics::DefaultReturn;
+
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Second update; based on code in update.c do_update_sd() w/ bFirstHalf = false in Gromacs 3.1.4
+
+   @param numberOfAtoms       number of atoms
+   @param atomCoordinates     atom coordinates
+   @param velocities          velocities
+   @param forces              forces
+   @param masses              atom masses
+
+   @return ReferenceDynamics::DefaultReturn
+
+   --------------------------------------------------------------------------------------- */
+
+int ReferenceVariableStochasticDynamics::updatePart2( int numberOfAtoms, RealOpenMM** atomCoordinates,
+                                              RealOpenMM** velocities,
+                                              RealOpenMM** forces, RealOpenMM* inverseMasses,
+                                              RealOpenMM** xPrime, RealOpenMM** oldVelocities,
+                                              RealOpenMM** xVector, RealOpenMM** vVector ){
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName  = "\nReferenceVariableStochasticDynamics::updatePart2";
+
+   static const RealOpenMM one    =  1.0;
+   static int debug               = 0;
+
+   // ---------------------------------------------------------------------------------------
+
+   // perform second update
+
+   const RealOpenMM*  fixedParameters = getFixedParameters();
+   RealOpenMM   tau                   = getTau();
+   RealOpenMM   fix1                  = tau*(fixedParameters[EPH] - fixedParameters[EMH]);
+                fix1                  = one/fix1;
+
+   for( int ii = 0; ii < numberOfAtoms; ii++ ){
+
+      RealOpenMM sqrtInvMass = SQRT( inverseMasses[ii] );
+
+      for( int jj = 0; jj < 3; jj++ ){
+
+           velocities[ii][jj]  = (xPrime[ii][jj] - atomCoordinates[ii][jj])*fix1;
+         RealOpenMM Xmh        = vVector[ii][jj]*tau*fixedParameters[D]/(fixedParameters[EM] - one) +
+                                 sqrtInvMass*fixedParameters[Yx]*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
+
+           xVector[ii][jj]     = sqrtInvMass*fixedParameters[X]*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
+
+           xPrime[ii][jj]     += xVector[ii][jj] - Xmh;
+      }
+   }
+
+   return ReferenceDynamics::DefaultReturn;
+
+}
+
+/**---------------------------------------------------------------------------------------
+
+   Update -- driver routine for performing stochastic dynamics update of coordinates
+   and velocities
+
+   @param numberOfAtoms       number of atoms
+   @param atomCoordinates     atom coordinates
+   @param velocities          velocities
+   @param forces              forces
+   @param masses              atom masses
+
+   @return ReferenceDynamics::DefaultReturn
+
+   --------------------------------------------------------------------------------------- */
+
+int ReferenceVariableStochasticDynamics::update( int numberOfAtoms, RealOpenMM** atomCoordinates,
+                                          RealOpenMM** velocities,
+                                          RealOpenMM** forces, RealOpenMM* masses, RealOpenMM maxStepSize ){
+
+   // ---------------------------------------------------------------------------------------
+
+   static const char* methodName      = "\nReferenceVariableStochasticDynamics::update";
+
+   static const RealOpenMM zero       =  0.0;
+   static const RealOpenMM one        =  1.0;
+
+   static int debug                   =  0;
+
+   // ---------------------------------------------------------------------------------------
+
+   // get work arrays
+
+   RealOpenMM** xPrime          = get2DArrayAtIndex( xPrime2D );
+   RealOpenMM** oldVelocities   = get2DArrayAtIndex( OldV );
+   RealOpenMM** xVector         = get2DArrayAtIndex( X2D );
+   RealOpenMM** vVector         = get2DArrayAtIndex( V2D );
+
+   RealOpenMM* inverseMasses    = get1DArrayAtIndex( InverseMasses );
+
+   // 1st update
+
+   updatePart1( numberOfAtoms, atomCoordinates, velocities, forces, masses, inverseMasses,
+                xPrime, oldVelocities, xVector, vVector, maxStepSize );
+
+   ReferenceConstraintAlgorithm* referenceConstraintAlgorithm = getReferenceConstraintAlgorithm();
+   if( referenceConstraintAlgorithm ){
+      referenceConstraintAlgorithm->apply( numberOfAtoms, atomCoordinates, xPrime,
+                                           inverseMasses );
+   }
+
+   // 2nd update
+
+   updatePart2( numberOfAtoms, atomCoordinates, velocities, forces, inverseMasses,
+                xPrime, oldVelocities, xVector, vVector );
+
+   if( referenceConstraintAlgorithm ){
+      referenceConstraintAlgorithm->apply( numberOfAtoms, atomCoordinates, xPrime,
+                                           inverseMasses );
+   }
+
+   // copy xPrime -> atomCoordinates
+
+   for( int ii = 0; ii < numberOfAtoms; ii++ ){
+      atomCoordinates[ii][0] = xPrime[ii][0];
+      atomCoordinates[ii][1] = xPrime[ii][1];
+      atomCoordinates[ii][2] = xPrime[ii][2];
+   }
+
+   incrementTimeStep();
+
+   return ReferenceDynamics::DefaultReturn;
+
+}

platforms/reference/src/SimTKReference/ReferenceVariableStochasticDynamics.h

+
+/* Portions copyright (c) 2006 Stanford University and Simbios.
+ * Contributors: Pande Group
+ *
+ * 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 __ReferenceVariableStochasticDynamics_H__
+#define __ReferenceVariableStochasticDynamics_H__
+
+#include "ReferenceDynamics.h"
+
+// ---------------------------------------------------------------------------------------
+
+class ReferenceVariableStochasticDynamics : public ReferenceDynamics {
+
+   private:
+
+      enum FixedParameters   { GDT, EPH, EMH, EP, EM, EM_V, B, C, D, V, X, Yv, Yx, MaxFixedParameters };
+      enum TwoDArrayIndicies { X2D, V2D, OldV, xPrime2D, vPrime2D, Max2DArrays };
+      enum OneDArrayIndicies { InverseMasses, Max1DArrays };
+
+      RealOpenMM _tau, _accuracy;
+      RealOpenMM _fixedParameters[MaxFixedParameters];
+
+      /**---------------------------------------------------------------------------------------
+
+         Set fixed values
+
+         @return ReferenceDynamics::DefaultReturn
+
+         --------------------------------------------------------------------------------------- */
+
+      int _setFixedParameters( RealOpenMM timeStep, RealOpenMM prevTimeStep );
+
+   public:
+
+      /**---------------------------------------------------------------------------------------
+
+         Constructor
+
+         @param numberOfAtoms  number of atoms
+         @param tau            viscosity
+         @param temperature    temperature
+         @param accuracy       required accuracy
+
+         --------------------------------------------------------------------------------------- */
+
+       ReferenceVariableStochasticDynamics( int numberOfAtoms, RealOpenMM tau, RealOpenMM temperature, RealOpenMM accuracy );
+
+      /**---------------------------------------------------------------------------------------
+
+         Destructor
+
+         --------------------------------------------------------------------------------------- */
+
+       ~ReferenceVariableStochasticDynamics( );
+
+      /**---------------------------------------------------------------------------------------
+
+         Get tau
+
+         @return tau
+
+         --------------------------------------------------------------------------------------- */
+
+      RealOpenMM getTau( void ) const;
+      
+      /**---------------------------------------------------------------------------------------
+
+         Get the required accuracy
+
+         @return accuracy
+
+         --------------------------------------------------------------------------------------- */
+
+      RealOpenMM getAccuracy( void ) const;
+
+      /**---------------------------------------------------------------------------------------
+
+         Set the required accuracy
+
+         --------------------------------------------------------------------------------------- */
+
+      void setAccuracy( RealOpenMM accuracy );
+
+      /**---------------------------------------------------------------------------------------
+
+         Get array of fixed parameters indexed by 'FixedParameters' enums
+
+         @return array
+
+         --------------------------------------------------------------------------------------- */
+
+      const RealOpenMM* getFixedParameters( void ) const;
+
+      /**---------------------------------------------------------------------------------------
+
+         Print parameters
+
+         @param message message
+
+         @return ReferenceDynamics::DefaultReturn
+
+         --------------------------------------------------------------------------------------- */
+
+      int printParameters( std::stringstream& message ) const;
+
+      /**---------------------------------------------------------------------------------------
+
+         Update
+
+         @param numberOfAtoms       number of atoms
+         @param atomCoordinates     atom coordinates
+         @param velocities          velocities
+         @param forces              forces
+         @param masses              atom masses
+         @param maxStepSize         maximum time step
+
+         @return ReferenceDynamics::DefaultReturn
+
+         --------------------------------------------------------------------------------------- */
+
+      int update( int numberOfAtoms, RealOpenMM** atomCoordinates,
+                  RealOpenMM** velocities, RealOpenMM** forces, RealOpenMM* masses, RealOpenMM maxStepSize );
+
+      /**---------------------------------------------------------------------------------------
+
+         First update; based on code in update.c do_update_sd() Gromacs 3.1.4
+
+         @param numberOfAtoms       number of atoms
+         @param atomCoordinates     atom coordinates
+         @param velocities          velocities
+         @param forces              forces
+         @param masses              atom masses
+         @param inverseMasses       inverse atom masses
+         @param xPrime              xPrime
+         @param oldVelocities       previous velocities
+         @param xVector             xVector
+         @param vVector             vVector
+         @param maxStepSize         maximum time step
+
+         @return ReferenceDynamics::DefaultReturn
+
+         --------------------------------------------------------------------------------------- */
+
+      int updatePart1( int numberOfAtoms, RealOpenMM** atomCoordinates, RealOpenMM** velocities,
+                       RealOpenMM** forces, RealOpenMM* masses, RealOpenMM* inverseMasses,
+                       RealOpenMM** xPrime, RealOpenMM** oldVelocities,
+                       RealOpenMM** xVector, RealOpenMM** vVector, RealOpenMM maxStepSize );
+
+      /**---------------------------------------------------------------------------------------
+
+         Second update
+
+         @param numberOfAtoms       number of atoms
+         @param atomCoordinates     atom coordinates
+         @param velocities          velocities
+         @param forces              forces
+         @param masses              atom masses
+
+         @return ReferenceDynamics::DefaultReturn
+
+         --------------------------------------------------------------------------------------- */
+
+      int updatePart2( int numberOfAtoms, RealOpenMM** atomCoordinates, RealOpenMM** velocities,
+                       RealOpenMM** forces, RealOpenMM* inverseMasses,
+                       RealOpenMM** xPrime, RealOpenMM** oldVelocities,
+                       RealOpenMM** xVector, RealOpenMM** vVector );
+      
+};
+
+// ---------------------------------------------------------------------------------------
+
+#endif // __ReferenceVariableStochasticDynamics_H__

platforms/reference/tests/TestReferenceVariableLangevinIntegrator.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-2009 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 VariableLangevinIntegrator.
+ */
+
+#include "../../../tests/AssertionUtilities.h"
+#include "openmm/OpenMMContext.h"
+#include "ReferencePlatform.h"
+#include "openmm/HarmonicBondForce.h"
+#include "openmm/NonbondedForce.h"
+#include "openmm/System.h"
+#include "openmm/VariableLangevinIntegrator.h"
+#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
+#include "../src/sfmt/SFMT.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+const double TOL = 1e-5;
+
+void testSingleBond() {
+    ReferencePlatform platform;
+    System system;
+    system.addParticle(2.0);
+    system.addParticle(2.0);
+    VariableLangevinIntegrator integrator(0, 0.1, 1e-6);
+    HarmonicBondForce* forceField = new HarmonicBondForce();
+    forceField->addBond(0, 1, 1.5, 1);
+    system.addForce(forceField);
+    OpenMMContext context(system, integrator, platform);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(-1, 0, 0);
+    positions[1] = Vec3(1, 0, 0);
+    context.setPositions(positions);
+
+    // This is simply a damped harmonic oscillator, so compare it to the analytical solution.
+
+    double freq = std::sqrt(1-0.05*0.05);
+    for (int i = 0; i < 1000; ++i) {
+        State state = context.getState(State::Positions | State::Velocities);
+        double time = state.getTime();
+        double expectedDist = 1.5+0.5*std::exp(-0.05*time)*std::cos(freq*time);
+        ASSERT_EQUAL_VEC(Vec3(-0.5*expectedDist, 0, 0), state.getPositions()[0], 0.02);
+        ASSERT_EQUAL_VEC(Vec3(0.5*expectedDist, 0, 0), state.getPositions()[1], 0.02);
+        double expectedSpeed = -0.5*std::exp(-0.05*time)*(0.05*std::cos(freq*time)+freq*std::sin(freq*time));
+        ASSERT_EQUAL_VEC(Vec3(-0.5*expectedSpeed, 0, 0), state.getVelocities()[0], 0.02);
+        ASSERT_EQUAL_VEC(Vec3(0.5*expectedSpeed, 0, 0), state.getVelocities()[1], 0.02);
+        integrator.step(1);
+    }
+
+    // Now set the friction to a tiny value and see if it conserves energy.
+
+    integrator.setFriction(5e-5);
+    context.setPositions(positions);
+    State state = context.getState(State::Energy);
+    double initialEnergy = state.getKineticEnergy()+state.getPotentialEnergy();
+    for (int i = 0; i < 1000; ++i) {
+        state = context.getState(State::Energy);
+        double energy = state.getKineticEnergy()+state.getPotentialEnergy();
+        ASSERT_EQUAL_TOL(initialEnergy, energy, 0.05);
+        integrator.step(1);
+    }
+}
+
+void testTemperature() {
+    const int numParticles = 8;
+    const double temp = 100.0;
+    ReferencePlatform platform;
+    System system;
+    VariableLangevinIntegrator integrator(temp, 2.0, 1e-4);
+    NonbondedForce* forceField = new NonbondedForce();
+    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);
+    OpenMMContext context(system, integrator, platform);
+    vector<Vec3> positions(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));
+    context.setPositions(positions);
+
+    // Let it equilibrate.
+
+    integrator.step(10000);
+
+    // Now run it for a while and see if the temperature is correct.
+
+    double ke = 0.0;
+    for (int i = 0; i < 1000; ++i) {
+        State state = context.getState(State::Energy);
+        ke += state.getKineticEnergy();
+        integrator.step(1);
+    }
+    ke /= 1000;
+    double expected = 0.5*numParticles*3*BOLTZ*temp;
+    ASSERT_EQUAL_TOL(expected, ke, 3*expected/std::sqrt(1000.0));
+}
+
+void testConstraints() {
+    const int numParticles = 8;
+    const double temp = 100.0;
+    ReferencePlatform platform;
+    System system;
+    VariableLangevinIntegrator integrator(temp, 2.0, 1e-5);
+    integrator.setConstraintTolerance(1e-5);
+    integrator.setRandomNumberSeed(0);
+    NonbondedForce* forceField = new NonbondedForce();
+    for (int i = 0; i < numParticles; ++i) {
+        system.addParticle(10.0);
+        forceField->addParticle((i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
+    }
+    for (int i = 0; i < numParticles-1; ++i)
+        system.addConstraint(i, i+1, 1.0);
+    system.addForce(forceField);
+    OpenMMContext context(system, integrator, platform);
+    vector<Vec3> positions(numParticles);
+    vector<Vec3> velocities(numParticles);
+    init_gen_rand(0);
+    for (int i = 0; i < numParticles; ++i) {
+        positions[i] = Vec3(i/2, (i+1)/2, 0);
+        velocities[i] = Vec3(genrand_real2()-0.5, genrand_real2()-0.5, genrand_real2()-0.5);
+    }
+    context.setPositions(positions);
+    context.setVelocities(velocities);
+
+    // Simulate it and see whether the constraints remain satisfied.
+
+    for (int i = 0; i < 1000; ++i) {
+        State state = context.getState(State::Positions);
+        for (int j = 0; j < numParticles-1; ++j) {
+            Vec3 p1 = state.getPositions()[j];
+            Vec3 p2 = state.getPositions()[j+1];
+            double dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2]));
+            ASSERT_EQUAL_TOL(1.0, dist, 2e-5);
+        }
+        integrator.step(1);
+    }
+}
+
+void testRandomSeed() {
+    const int numParticles = 8;
+    const double temp = 100.0;
+    ReferencePlatform platform;
+    System system;
+    VariableLangevinIntegrator integrator(temp, 2.0, 1e-5);
+    NonbondedForce* forceField = new NonbondedForce();
+    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);
+    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.
+
+    integrator.setRandomNumberSeed(5);
+    OpenMMContext 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.
+
+    integrator.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 {
+        testSingleBond();
+        testTemperature();
+        testConstraints();
+        testRandomSeed();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}
+