Merge pull request #2518 from andysim/inittemp

Modify temperature initialization mechanism

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) {

wrappers/python/src/swig_doxygen/swigInputConfig.py

@@ -470,8 +470,11 @@ UNITS = {
 ("System", "getConstraintParameters") : (None, (None, None, 'unit.nanometer')),
 ("System", "getForce") : (None, ()),
 ("System", "getVirtualSite") : (None, ()),
+("DrudeIntegrator", "getDrudeTemperature") : ("unit.kelvin", ()),
+("DrudeIntegrator", "getMaxDrudeDistance") : ("unit.nanometer", ()),
 ("DrudeLangevinIntegrator", "getDrudeTemperature") : ("unit.kelvin", ()),
 ("DrudeLangevinIntegrator", "getMaxDrudeDistance") : ("unit.nanometer", ()),
+("DrudeNoseHooverIntegrator", "getVelocitiesForTemperature") : ("unit.nanometers / unit.picosecond", ()),
 ("MonteCarloMembraneBarostat", "MonteCarloMembraneBarostat") : (None, ("unit.bar", "unit.bar*unit.nanometer", "unit.kelvin", None, None, None)),
 ("MonteCarloMembraneBarostat", "getXYMode") : (None, ()),
 ("MonteCarloMembraneBarostat", "getZMode") : (None, ()),