Added Context::applyVelocityConstraints() and Context::setVelocitiesToTemperature()

olla/include/openmm/kernels.h

@@ -222,6 +222,13 @@ public:
      * @param tol        the distance tolerance within which constraints must be satisfied.
      */
     virtual void apply(ContextImpl& context, double tol) = 0;
+    /**
+     * Update particle velocities to enforce constraints.
+     *
+     * @param context    the context in which to execute this kernel
+     * @param tol        the velocity tolerance within which constraints must be satisfied.
+     */
+    virtual void applyToVelocities(ContextImpl& context, double tol) = 0;
 };
 
 /**

openmmapi/include/openmm/Context.h

@@ -35,6 +35,7 @@
 #include "Integrator.h"
 #include "State.h"
 #include "System.h"
+#include <ctime>
 #include <iosfwd>
 #include <map>
 #include <string>
@@ -148,6 +149,14 @@ public:
      * contains the velocity of the i'th particle.
      */
     void setVelocities(const std::vector<Vec3>& velocities);
+    /**
+     * Set the velocities of all particles in the System to random values chosen from a Boltzmann
+     * distribution at a given temperature.
+     * 
+     * @param temperature    the temperature for which to select the velocities (measured in Kelvin)
+     * @param randomSeed     the random number seed to use when selecting velocities
+     */
+    void setVelocitiesToTemperature(double temperature, int randomSeed=time(NULL));
     /**
      * Get the value of an adjustable parameter defined by a Force object in the System.
      * 
@@ -180,6 +189,12 @@ public:
      * @param tol    the distance tolerance within which constraints must be satisfied.
      */
     void applyConstraints(double tol);
+    /**
+     * Update the velocities of particles so the net velocity of each constrained distance is zero.
+     *
+     * @param tol    the velocity tolerance within which constraints must be satisfied.
+     */
+    void applyVelocityConstraints(double tol);
     /**
      * Recompute the locations of all virtual sites.  There is rarely a reason to call
      * this, since virtual sites are also updated by applyConstraints().  This is only

openmmapi/include/openmm/internal/ContextImpl.h

@@ -169,6 +169,12 @@ public:
      * @param tol    the distance tolerance within which constraints must be satisfied.
      */
     void applyConstraints(double tol);
+    /**
+     * Update the velocities of particles so the net velocity of each constrained distance is zero.
+     *
+     * @param tol    the velocity tolerance within which constraints must be satisfied.
+     */
+    void applyVelocityConstraints(double tol);
     /**
      * Recompute the locations of all virtual sites.  There is rarely a reason to call
      * this, since virtual sites are also updated by applyConstraints().  This is only

openmmapi/src/Context.cpp

@@ -33,6 +33,8 @@
 #include "openmm/internal/ContextImpl.h"
 #include "openmm/OpenMMException.h"
 #include "openmm/internal/ForceImpl.h"
+#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
+#include "sfmt/SFMT.h"
 #include <cmath>
 
 using namespace OpenMM;
@@ -151,6 +153,41 @@ void Context::setVelocities(const vector<Vec3>& velocities) {
     impl->setVelocities(velocities);
 }
 
+void Context::setVelocitiesToTemperature(double temperature, int randomSeed) {
+    System& system = impl->getSystem();
+    
+    // Generate the list of Gaussian random numbers.
+    
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(randomSeed, sfmt);
+    vector<double> randoms;
+    while (randoms.size() < system.getNumParticles()*3) {
+        double x, y, r2;
+        do {
+            x = 2.0*genrand_real2(sfmt)-1.0;
+            y = 2.0*genrand_real2(sfmt)-1.0;
+            r2 = x*x + y*y;
+        } while (r2 >= 1.0 || r2 == 0.0);
+        double multiplier = sqrt((-2.0*log(r2))/r2);
+        randoms.push_back(x*multiplier);
+        randoms.push_back(y*multiplier);
+    }
+    
+    // Assign the velocities.
+    
+    vector<Vec3> velocities(system.getNumParticles(), Vec3());
+    int nextRandom = 0;
+    for (int i = 0; i < system.getNumParticles(); i++) {
+        double mass = system.getParticleMass(i);
+        if (mass != 0) {
+            double velocityScale = sqrt(BOLTZ*temperature/mass);
+            velocities[i] = Vec3(randoms[nextRandom++], randoms[nextRandom++], randoms[nextRandom++])*velocityScale;
+        }
+    }
+    setVelocities(velocities);
+    impl->applyVelocityConstraints(1e-5);
+}
+
 double Context::getParameter(const string& name) {
     return impl->getParameter(name);
 }
@@ -167,6 +204,10 @@ void Context::applyConstraints(double tol) {
     impl->applyConstraints(tol);
 }
 
+void Context::applyVelocityConstraints(double tol) {
+    impl->applyVelocityConstraints(tol);
+}
+
 void Context::computeVirtualSites() {
     impl->computeVirtualSites();
 }

openmmapi/src/ContextImpl.cpp

@@ -192,6 +192,10 @@ void ContextImpl::applyConstraints(double tol) {
     applyConstraintsKernel.getAs<ApplyConstraintsKernel>().apply(*this, tol);
 }
 
+void ContextImpl::applyVelocityConstraints(double tol) {
+    applyConstraintsKernel.getAs<ApplyConstraintsKernel>().applyToVelocities(*this, tol);
+}
+
 void ContextImpl::computeVirtualSites() {
     virtualSitesKernel.getAs<VirtualSitesKernel>().computePositions(*this);
 }

platforms/cuda/src/CudaKernels.cpp

@@ -399,6 +399,10 @@ void CudaApplyConstraintsKernel::apply(ContextImpl& context, double tol) {
     integration.computeVirtualSites();
 }
 
+void CudaApplyConstraintsKernel::applyToVelocities(ContextImpl& context, double tol) {
+    cu.getIntegrationUtilities().applyVelocityConstraints(tol);
+}
+
 void CudaVirtualSitesKernel::initialize(const System& system) {
 }
 

platforms/cuda/src/CudaKernels.h

@@ -189,6 +189,13 @@ public:
      * @param tol        the distance tolerance within which constraints must be satisfied.
      */
     void apply(ContextImpl& context, double tol);
+    /**
+     * Update particle velocities to enforce constraints.
+     *
+     * @param context    the context in which to execute this kernel
+     * @param tol        the velocity tolerance within which constraints must be satisfied.
+     */
+    void applyToVelocities(ContextImpl& context, double tol);
 private:
     CudaContext& cu;
     bool hasInitializedKernel;

platforms/cuda/tests/TestCudaRandom.cpp

@@ -38,6 +38,10 @@
 #include "../src/CudaContext.h"
 #include "../src/CudaIntegrationUtilities.h"
 #include "openmm/System.h"
+#include "openmm/Context.h"
+#include "CudaPlatform.h"
+#include "openmm/VerletIntegrator.h"
+#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
 #include <iostream>
 
 using namespace OpenMM;
@@ -84,9 +88,53 @@ void testGaussian() {
     ASSERT_EQUAL_TOL(0.0, c4, 3.0/pow(numValues, 1.0/4.0));
 }
 
+void testRandomVelocities() {
+    // Create a system.
+    
+    const int numParticles = 10000;
+    const double temperture = 100.0;
+    CudaPlatform platform;
+    System system;
+    VerletIntegrator integrator(0.01);
+    for (int i = 0; i < numParticles; ++i)
+        system.addParticle(10.0+sin(0.1*i));
+    for (int i = 0; i < numParticles-1; ++i)
+        system.addConstraint(i, i+1, 1.0);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(numParticles);
+    for (int i = 0; i < numParticles; ++i)
+        positions[i] = Vec3(i/2, (i+1)/2, 0);
+    context.setPositions(positions);
+    
+    // Ask the context to generate random velocities.
+    
+    context.setVelocitiesToTemperature(temperture);
+    State state = context.getState(State::Velocities);
+    
+    // See if they respect constraints.
+    
+    for (int i = 1; i < numParticles; i++) {
+        Vec3 v1 = state.getVelocities()[i-1];
+        Vec3 v2 = state.getVelocities()[i];
+        double vel = (v1-v2).dot(positions[i-1]-positions[i]);
+        ASSERT_EQUAL_TOL(0.0, vel, 2e-5);
+    }
+    
+    // See if the temperature is correct.
+
+    double ke = 0;
+    for (int i = 0; i < numParticles; i++) {
+        Vec3 v = state.getVelocities()[i];
+        ke += 0.5*system.getParticleMass(i)*v.dot(v);
+    }
+    double expected = 0.5*(numParticles*3-system.getNumConstraints())*BOLTZ*temperture;
+    ASSERT_USUALLY_EQUAL_TOL(expected, ke, 4/sqrt(numParticles));
+}
+
 int main() {
     try {
         testGaussian();
+        testRandomVelocities();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;

platforms/opencl/src/OpenCLKernels.cpp

@@ -424,6 +424,10 @@ void OpenCLApplyConstraintsKernel::apply(ContextImpl& context, double tol) {
     integration.computeVirtualSites();
 }
 
+void OpenCLApplyConstraintsKernel::applyToVelocities(ContextImpl& context, double tol) {
+    cl.getIntegrationUtilities().applyVelocityConstraints(tol);
+}
+
 void OpenCLVirtualSitesKernel::initialize(const System& system) {
 }
 

platforms/opencl/src/OpenCLKernels.h

@@ -189,6 +189,13 @@ public:
      * @param tol        the distance tolerance within which constraints must be satisfied.
      */
     void apply(ContextImpl& context, double tol);
+    /**
+     * Update particle velocities to enforce constraints.
+     *
+     * @param context    the context in which to execute this kernel
+     * @param tol        the velocity tolerance within which constraints must be satisfied.
+     */
+    void applyToVelocities(ContextImpl& context, double tol);
 private:
     OpenCLContext& cl;
     bool hasInitializedKernel;

platforms/opencl/tests/TestOpenCLRandom.cpp

@@ -38,6 +38,10 @@
 #include "../src/OpenCLContext.h"
 #include "../src/OpenCLIntegrationUtilities.h"
 #include "openmm/System.h"
+#include "openmm/Context.h"
+#include "OpenCLPlatform.h"
+#include "openmm/VerletIntegrator.h"
+#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
 #include <iostream>
 
 using namespace OpenMM;
@@ -82,9 +86,53 @@ void testGaussian() {
     ASSERT_EQUAL_TOL(0.0, c4, 3.0/pow(numValues, 1.0/4.0));
 }
 
+void testRandomVelocities() {
+    // Create a system.
+    
+    const int numParticles = 10000;
+    const double temperture = 100.0;
+    OpenCLPlatform platform;
+    System system;
+    VerletIntegrator integrator(0.01);
+    for (int i = 0; i < numParticles; ++i)
+        system.addParticle(10.0+sin(0.1*i));
+    for (int i = 0; i < numParticles-1; ++i)
+        system.addConstraint(i, i+1, 1.0);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(numParticles);
+    for (int i = 0; i < numParticles; ++i)
+        positions[i] = Vec3(i/2, (i+1)/2, 0);
+    context.setPositions(positions);
+    
+    // Ask the context to generate random velocities.
+    
+    context.setVelocitiesToTemperature(temperture);
+    State state = context.getState(State::Velocities);
+    
+    // See if they respect constraints.
+    
+    for (int i = 1; i < numParticles; i++) {
+        Vec3 v1 = state.getVelocities()[i-1];
+        Vec3 v2 = state.getVelocities()[i];
+        double vel = (v1-v2).dot(positions[i-1]-positions[i]);
+        ASSERT_EQUAL_TOL(0.0, vel, 2e-5);
+    }
+    
+    // See if the temperature is correct.
+
+    double ke = 0;
+    for (int i = 0; i < numParticles; i++) {
+        Vec3 v = state.getVelocities()[i];
+        ke += 0.5*system.getParticleMass(i)*v.dot(v);
+    }
+    double expected = 0.5*(numParticles*3-system.getNumConstraints())*BOLTZ*temperture;
+    ASSERT_USUALLY_EQUAL_TOL(expected, ke, 4/sqrt(numParticles));
+}
+
 int main() {
     try {
         testGaussian();
+        testRandomVelocities();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;

platforms/reference/src/ReferenceKernels.cpp

@@ -326,6 +326,18 @@ void ReferenceApplyConstraintsKernel::apply(ContextImpl& context, double tol) {
     ReferenceVirtualSites::computePositions(context.getSystem(), positions);
 }
 
+void ReferenceApplyConstraintsKernel::applyToVelocities(ContextImpl& context, double tol) {
+    if (constraints == NULL) {
+        vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
+        findAnglesForCCMA(context.getSystem(), angles);
+        constraints = new ReferenceCCMAAlgorithm(context.getSystem().getNumParticles(), numConstraints, constraintIndices, constraintDistances, masses, angles, tol);
+    }
+    vector<RealVec>& positions = extractPositions(context);
+    vector<RealVec>& velocities = extractVelocities(context);
+    constraints->setTolerance(tol);
+    constraints->applyToVelocities(data.numParticles, positions, velocities, inverseMasses);
+}
+
 void ReferenceVirtualSitesKernel::initialize(const System& system) {
 }
 

platforms/reference/src/ReferenceKernels.h

@@ -206,6 +206,13 @@ public:
      * @param tol        the distance tolerance within which constraints must be satisfied.
      */
     void apply(ContextImpl& context, double tol);
+    /**
+     * Update particle velocities to enforce constraints.
+     *
+     * @param context    the context in which to execute this kernel
+     * @param tol        the velocity tolerance within which constraints must be satisfied.
+     */
+    void applyToVelocities(ContextImpl& context, double tol);
 private:
     ReferencePlatform::PlatformData& data;
     ReferenceConstraintAlgorithm* constraints;

platforms/reference/tests/TestReferenceRandom.cpp

@@ -35,6 +35,9 @@
 
 #include "openmm/internal/AssertionUtilities.h"
 #include "../src/SimTKUtilities/SimTKOpenMMUtilities.h"
+#include "openmm/Context.h"
+#include "ReferencePlatform.h"
+#include "openmm/VerletIntegrator.h"
 #include <iostream>
 
 using namespace OpenMM;
@@ -46,7 +49,6 @@ void testGaussian() {
     double var = 0.0;
     double skew = 0.0;
     double kurtosis = 0.0;
-    unsigned long jran = 12399103;
     for (int i = 0; i < numValues; i++) {
         double value = SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
         mean += value;
@@ -67,9 +69,53 @@ void testGaussian() {
     ASSERT_EQUAL_TOL(0.0, c4, 0.01);
 }
 
+void testRandomVelocities() {
+    // Create a system.
+    
+    const int numParticles = 10000;
+    const double temperture = 100.0;
+    ReferencePlatform platform;
+    System system;
+    VerletIntegrator integrator(0.01);
+    for (int i = 0; i < numParticles; ++i)
+        system.addParticle(10.0+sin(0.1*i));
+    for (int i = 0; i < numParticles-1; ++i)
+        system.addConstraint(i, i+1, 1.0);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(numParticles);
+    for (int i = 0; i < numParticles; ++i)
+        positions[i] = Vec3(i/2, (i+1)/2, 0);
+    context.setPositions(positions);
+    
+    // Ask the context to generate random velocities.
+    
+    context.setVelocitiesToTemperature(temperture);
+    State state = context.getState(State::Velocities);
+    
+    // See if they respect constraints.
+    
+    for (int i = 1; i < numParticles; i++) {
+        Vec3 v1 = state.getVelocities()[i-1];
+        Vec3 v2 = state.getVelocities()[i];
+        double vel = (v1-v2).dot(positions[i-1]-positions[i]);
+        ASSERT_EQUAL_TOL(0.0, vel, 2e-5);
+    }
+    
+    // See if the temperature is correct.
+
+    double ke = 0;
+    for (int i = 0; i < numParticles; i++) {
+        Vec3 v = state.getVelocities()[i];
+        ke += 0.5*system.getParticleMass(i)*v.dot(v);
+    }
+    double expected = 0.5*(numParticles*3-system.getNumConstraints())*BOLTZ*temperture;
+    ASSERT_USUALLY_EQUAL_TOL(expected, ke, 4/sqrt(numParticles));
+}
+
 int main() {
     try {
         testGaussian();
+        testRandomVelocities();
     }
     catch(const exception& e) {
         cout << "exception: " << e.what() << endl;