Delegate temperature initialization to integrators

openmmapi/include/openmm/Integrator.h

@@ -42,6 +42,7 @@ namespace OpenMM {
 
 class Context;
 class ContextImpl;
+class System;
 
 /**
  * An Integrator defines a method for simulating a System by integrating the equations of motion.
@@ -133,6 +134,16 @@ protected:
     virtual bool kineticEnergyRequiresForce() const {
         return true;
     }
+    /**
+     * Return a list of velocities normally distributed around a target temperature.  This may be
+     * overridden by Drude integrators to ensure that Drude pairs have their center of mass velocity
+     * assigned as a single entity, rather than treating both particles as being independent.
+     *
+     * @param system the system whose velocities are to be initialized.
+     * @param temperature the target temperature in Kelvin.
+     * @param randomSeed the random number seed to use when selecting velocities 
+     */
+    virtual std::vector<Vec3> getVelocitiesForTemperature(const System &system, double temperature, int randomSeed) const;
 private:
     double stepSize, constraintTol;
 };

openmmapi/src/Context.cpp

@@ -33,8 +33,6 @@
 #include "openmm/internal/ContextImpl.h"
 #include "openmm/OpenMMException.h"
 #include "openmm/internal/ForceImpl.h"
-#include "SimTKOpenMMRealType.h"
-#include "sfmt/SFMT.h"
 #include <cmath>
 #include <iostream>
 #include <sstream>
@@ -183,37 +181,9 @@ void Context::setVelocities(const vector<Vec3>& velocities) {
 }
 
 void Context::setVelocitiesToTemperature(double temperature, int randomSeed) {
+    const Integrator& integrator = impl->getIntegrator();
     const 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);
+    setVelocities(integrator.getVelocitiesForTemperature(system, temperature, randomSeed));
     impl->applyVelocityConstraints(1e-5);
 }
 

openmmapi/src/Integrator.cpp

@@ -29,9 +29,14 @@
  * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
  * -------------------------------------------------------------------------- */
 
+#include "sfmt/SFMT.h"
+#include "SimTKOpenMMRealType.h"
 #include "openmm/Integrator.h"
+#include "openmm/System.h"
 #include "openmm/internal/ContextImpl.h"
 
+#include <cmath>
+
 using namespace OpenMM;
 
 Integrator::Integrator() : owner(NULL), context(NULL) {
@@ -63,3 +68,34 @@ double Integrator::getConstraintTolerance() const {
 void Integrator::setConstraintTolerance(double tol) {
     constraintTol = tol;
 }
+
+std::vector<Vec3> Integrator::getVelocitiesForTemperature(const System &system, double temperature, int randomSeed) const {
+    // Generate the list of Gaussian random numbers.
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(randomSeed, sfmt);
+    std::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*std::log(r2))/r2);
+        randoms.push_back(x*multiplier);
+        randoms.push_back(y*multiplier);
+    }
+
+    // Assign the velocities.
+    std::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;
+        }
+    }
+    return velocities;
+}
+

tests/TestBrownianIntegrator.h

@@ -251,6 +251,34 @@ void testRandomSeed() {
     }
 }
 
+void testInitialTemperature() {
+    // Check temperature initialization for a collection of randomly placed particles
+    const int numParticles = 500000;
+    const int nDoF = 3 * numParticles;
+    const double targetTemperature = 300;
+    System system;
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    std::vector<Vec3> positions(numParticles);
+
+    for (int i = 0; i < numParticles; i++) {
+        system.addParticle(1.0);
+        positions[i][0] = genrand_real2(sfmt);
+        positions[i][1] = genrand_real2(sfmt);
+        positions[i][2] = genrand_real2(sfmt);
+    }
+
+    BrownianIntegrator integrator(300, 2, 0.001);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    context.setVelocitiesToTemperature(targetTemperature);
+    auto velocities = context.getState(State::Velocities).getVelocities();
+    double kineticEnergy = 0;
+    for(const auto &v : velocities) kineticEnergy += 0.5 * v.dot(v);
+    double temperature = (2*kineticEnergy / (nDoF*BOLTZ));
+    ASSERT_USUALLY_EQUAL_TOL(targetTemperature, temperature, 0.01);
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
@@ -261,6 +289,7 @@ int main(int argc, char* argv[]) {
         testConstraints();
         testConstrainedMasslessParticles();
         testRandomSeed();
+        testInitialTemperature();
         runPlatformTests();
     }
     catch(const exception& e) {

tests/TestCustomIntegrator.h

@@ -1128,6 +1128,34 @@ void testRecordEnergy() {
     }
 }
 
+void testInitialTemperature() {
+    // Check temperature initialization for a collection of randomly placed particles
+    const int numParticles = 500000;
+    const int nDoF = 3 * numParticles;
+    const double targetTemperature = 300;
+    System system;
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    std::vector<Vec3> positions(numParticles);
+
+    for (int i = 0; i < numParticles; i++) {
+        system.addParticle(1.0);
+        positions[i][0] = genrand_real2(sfmt);
+        positions[i][1] = genrand_real2(sfmt);
+        positions[i][2] = genrand_real2(sfmt);
+    }
+
+    CustomIntegrator integrator(0.001);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    context.setVelocitiesToTemperature(targetTemperature);
+    auto velocities = context.getState(State::Velocities).getVelocities();
+    double kineticEnergy = 0;
+    for(const auto &v : velocities) kineticEnergy += 0.5 * v.dot(v);
+    double temperature = (2*kineticEnergy / (nDoF*BOLTZ));
+    ASSERT_USUALLY_EQUAL_TOL(targetTemperature, temperature, 0.01);
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
@@ -1155,6 +1183,7 @@ int main(int argc, char* argv[]) {
         testUpdateContextState();
         testVectorFunctions();
         testRecordEnergy();
+        testInitialTemperature();
         runPlatformTests();
     }
     catch(const exception& e) {

tests/TestLangevinIntegrator.h

@@ -258,6 +258,34 @@ void testRandomSeed() {
     }
 }
 
+void testInitialTemperature() {
+    // Check temperature initialization for a collection of randomly placed particles
+    const int numParticles = 500000;
+    const int nDoF = 3 * numParticles;
+    const double targetTemperature = 300;
+    System system;
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    std::vector<Vec3> positions(numParticles);
+
+    for (int i = 0; i < numParticles; i++) {
+        system.addParticle(1.0);
+        positions[i][0] = genrand_real2(sfmt);
+        positions[i][1] = genrand_real2(sfmt);
+        positions[i][2] = genrand_real2(sfmt);
+    }
+
+    LangevinIntegrator integrator(300, 25, 0.001);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    context.setVelocitiesToTemperature(targetTemperature);
+    auto velocities = context.getState(State::Velocities).getVelocities();
+    double kineticEnergy = 0;
+    for(const auto &v : velocities) kineticEnergy += 0.5 * v.dot(v);
+    double temperature = (2*kineticEnergy / (nDoF*BOLTZ));
+    ASSERT_USUALLY_EQUAL_TOL(targetTemperature, temperature, 0.01);
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
@@ -268,6 +296,7 @@ int main(int argc, char* argv[]) {
         testConstraints();
         testConstrainedMasslessParticles();
         testRandomSeed();
+        testInitialTemperature();
         runPlatformTests();
     }
     catch(const exception& e) {

tests/TestLangevinMiddleIntegrator.h

@@ -261,6 +261,34 @@ void testRandomSeed() {
     }
 }
 
+void testInitialTemperature() {
+    // Check temperature initialization for a collection of randomly placed particles
+    const int numParticles = 500000;
+    const int nDoF = 3 * numParticles;
+    const double targetTemperature = 300;
+    System system;
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    std::vector<Vec3> positions(numParticles);
+
+    for (int i = 0; i < numParticles; i++) {
+        system.addParticle(1.0);
+        positions[i][0] = genrand_real2(sfmt);
+        positions[i][1] = genrand_real2(sfmt);
+        positions[i][2] = genrand_real2(sfmt);
+    }
+
+    BAOABLangevinIntegrator integrator(300, 25, 0.001);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    context.setVelocitiesToTemperature(targetTemperature);
+    auto velocities = context.getState(State::Velocities).getVelocities();
+    double kineticEnergy = 0;
+    for(const auto &v : velocities) kineticEnergy += 0.5 * v.dot(v);
+    double temperature = (2*kineticEnergy / (nDoF*BOLTZ));
+    ASSERT_USUALLY_EQUAL_TOL(targetTemperature, temperature, 0.01);
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
@@ -271,6 +299,7 @@ int main(int argc, char* argv[]) {
         testConstraints();
         testConstrainedMasslessParticles();
         testRandomSeed();
+        testInitialTemperature();
         runPlatformTests();
     }
     catch(const exception& e) {

tests/TestNoseHooverIntegrator.h

@@ -269,6 +269,34 @@ void testThreeParticleVirtualSite() {
     }
 }
 
+void testInitialTemperature() {
+    // Check temperature initialization for a collection of randomly placed particles
+    const int numParticles = 500000;
+    const int nDoF = 3 * numParticles;
+    const double targetTemperature = 300;
+    System system;
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    std::vector<Vec3> positions(numParticles);
+
+    for (int i = 0; i < numParticles; i++) {
+        system.addParticle(1.0);
+        positions[i][0] = genrand_real2(sfmt);
+        positions[i][1] = genrand_real2(sfmt);
+        positions[i][2] = genrand_real2(sfmt);
+    }
+
+    NoseHooverIntegrator integrator(300, 25, 0.001);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    context.setVelocitiesToTemperature(targetTemperature);
+    auto velocities = context.getState(State::Velocities).getVelocities();
+    double kineticEnergy = 0;
+    for(const auto &v : velocities) kineticEnergy += 0.5 * v.dot(v);
+    double temperature = (2*kineticEnergy / (nDoF*BOLTZ));
+    ASSERT_USUALLY_EQUAL_TOL(targetTemperature, temperature, 0.01);
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
@@ -279,6 +307,7 @@ int main(int argc, char* argv[]) {
         testVVConstrainedClusters();
         testVVConstrainedMasslessParticles();
         testThreeParticleVirtualSite();
+        testInitialTemperature();
         runPlatformTests();
     }
     catch(const exception& e) {

tests/TestVariableLangevinIntegrator.h

@@ -330,6 +330,34 @@ void testArgonBox() {
     }
 }
 
+void testInitialTemperature() {
+    // Check temperature initialization for a collection of randomly placed particles
+    const int numParticles = 500000;
+    const int nDoF = 3 * numParticles;
+    const double targetTemperature = 300;
+    System system;
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    std::vector<Vec3> positions(numParticles);
+
+    for (int i = 0; i < numParticles; i++) {
+        system.addParticle(1.0);
+        positions[i][0] = genrand_real2(sfmt);
+        positions[i][1] = genrand_real2(sfmt);
+        positions[i][2] = genrand_real2(sfmt);
+    }
+
+    VariableLangevinIntegrator integrator(300, 25, 1e-5);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    context.setVelocitiesToTemperature(targetTemperature);
+    auto velocities = context.getState(State::Velocities).getVelocities();
+    double kineticEnergy = 0;
+    for(const auto &v : velocities) kineticEnergy += 0.5 * v.dot(v);
+    double temperature = (2*kineticEnergy / (nDoF*BOLTZ));
+    ASSERT_USUALLY_EQUAL_TOL(targetTemperature, temperature, 0.01);
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
@@ -341,6 +369,7 @@ int main(int argc, char* argv[]) {
         testConstrainedMasslessParticles();
         testRandomSeed();
         testArgonBox();
+        testInitialTemperature();
         runPlatformTests();
     }
     catch(const exception& e) {

tests/TestVariableVerletIntegrator.h

@@ -309,6 +309,34 @@ void testArgonBox() {
     }
 }
 
+void testInitialTemperature() {
+    // Check temperature initialization for a collection of randomly placed particles
+    const int numParticles = 500000;
+    const int nDoF = 3 * numParticles;
+    const double targetTemperature = 300;
+    System system;
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    std::vector<Vec3> positions(numParticles);
+
+    for (int i = 0; i < numParticles; i++) {
+        system.addParticle(1.0);
+        positions[i][0] = genrand_real2(sfmt);
+        positions[i][1] = genrand_real2(sfmt);
+        positions[i][2] = genrand_real2(sfmt);
+    }
+
+    VariableVerletIntegrator integrator(0.001);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    context.setVelocitiesToTemperature(targetTemperature);
+    auto velocities = context.getState(State::Velocities).getVelocities();
+    double kineticEnergy = 0;
+    for(const auto &v : velocities) kineticEnergy += 0.5 * v.dot(v);
+    double temperature = (2*kineticEnergy / (nDoF*BOLTZ));
+    ASSERT_USUALLY_EQUAL_TOL(targetTemperature, temperature, 0.01);
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
@@ -319,6 +347,7 @@ int main(int argc, char* argv[]) {
         testConstrainedClusters();
         testConstrainedMasslessParticles();
         testArgonBox();
+        testInitialTemperature();
         runPlatformTests();
     }
     catch(const exception& e) {

tests/TestVerletIntegrator.h

@@ -226,6 +226,34 @@ void testConstrainedMasslessParticles() {
     ASSERT_EQUAL(0.0, state.getVelocities()[0][0]);
 }
 
+void testInitialTemperature() {
+    // Check temperature initialization for a collection of randomly placed particles
+    const int numParticles = 500000;
+    const int nDoF = 3 * numParticles;
+    const double targetTemperature = 300;
+    System system;
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    std::vector<Vec3> positions(numParticles);
+
+    for (int i = 0; i < numParticles; i++) {
+        system.addParticle(1.0);
+        positions[i][0] = genrand_real2(sfmt);
+        positions[i][1] = genrand_real2(sfmt);
+        positions[i][2] = genrand_real2(sfmt);
+    }
+
+    VerletIntegrator integrator(0.001);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    context.setVelocitiesToTemperature(targetTemperature);
+    auto velocities = context.getState(State::Velocities).getVelocities();
+    double kineticEnergy = 0;
+    for(const auto &v : velocities) kineticEnergy += 0.5 * v.dot(v);
+    double temperature = (2*kineticEnergy / (nDoF*BOLTZ));
+    ASSERT_USUALLY_EQUAL_TOL(targetTemperature, temperature, 0.01);
+}
+
 void runPlatformTests();
 
 int main(int argc, char* argv[]) {
@@ -235,6 +263,7 @@ int main(int argc, char* argv[]) {
         testConstraints();
         testConstrainedClusters();
         testConstrainedMasslessParticles();
+        testInitialTemperature();
         runPlatformTests();
     }
     catch(const exception& e) {