setVelocitiesToTemperature() applies offset for leapfrog integrators

openmmapi/include/openmm/CompoundIntegrator.h

@@ -187,6 +187,13 @@ protected:
      * The implementation calls computeKineticEnergy() on whichever Integrator has been set as current.
      */
     double computeKineticEnergy();
+    /**
+     * Get the time interval by which velocities are offset from positions.  This is used to
+     * adjust velocities when setVelocitiesToTemperature() is called on a Context.
+     */
+    double getVelocityTimeOffset() const {
+        return getIntegrator(0).getVelocityTimeOffset();
+    }
 private:
     int currentIntegrator;
     std::vector<Integrator*> integrators;

openmmapi/include/openmm/Integrator.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2019 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2020 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -144,6 +144,13 @@ protected:
      * @param randomSeed the random number seed to use when selecting velocities 
      */
     virtual std::vector<Vec3> getVelocitiesForTemperature(const System &system, double temperature, int randomSeed) const;
+    /**
+     * Get the time interval by which velocities are offset from positions.  This is used to
+     * adjust velocities when setVelocitiesToTemperature() is called on a Context.
+     */
+    virtual double getVelocityTimeOffset() const {
+        return 0.0;
+    }
 private:
     double stepSize, constraintTol;
 };

openmmapi/include/openmm/LangevinIntegrator.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2012 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2020 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -133,6 +133,13 @@ protected:
      * Compute the kinetic energy of the system at the current time.
      */
     double computeKineticEnergy();
+    /**
+     * Get the time interval by which velocities are offset from positions.  This is used to
+     * adjust velocities when setVelocitiesToTemperature() is called on a Context.
+     */
+    double getVelocityTimeOffset() const {
+        return getStepSize()/2;
+    }
 private:
     double temperature, friction;
     int randomNumberSeed;

openmmapi/include/openmm/LangevinMiddleIntegrator.h

@@ -140,6 +140,13 @@ protected:
      * Computing kinetic energy for this integrator does not require forces.
      */
     bool kineticEnergyRequiresForce() const;
+    /**
+     * Get the time interval by which velocities are offset from positions.  This is used to
+     * adjust velocities when setVelocitiesToTemperature() is called on a Context.
+     */
+    double getVelocityTimeOffset() const {
+        return getStepSize()/2;
+    }
 private:
     double temperature, friction;
     int randomNumberSeed;

openmmapi/include/openmm/VerletIntegrator.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2012 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2020 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -76,6 +76,13 @@ protected:
      * Compute the kinetic energy of the system at the current time.
      */
     double computeKineticEnergy();
+    /**
+     * Get the time interval by which velocities are offset from positions.  This is used to
+     * adjust velocities when setVelocitiesToTemperature() is called on a Context.
+     */
+    double getVelocityTimeOffset() const {
+        return getStepSize()/2;
+    }
 private:
     Kernel kernel;
 };

openmmapi/src/Context.cpp

@@ -183,7 +183,19 @@ void Context::setVelocities(const vector<Vec3>& velocities) {
 void Context::setVelocitiesToTemperature(double temperature, int randomSeed) {
     const Integrator& integrator = impl->getIntegrator();
     const System& system = impl->getSystem();
-    setVelocities(integrator.getVelocitiesForTemperature(system, temperature, randomSeed));
+    vector<Vec3> velocities = integrator.getVelocitiesForTemperature(system, temperature, randomSeed);
+    double offset = integrator.getVelocityTimeOffset();
+    if (offset != 0.0) {
+        impl->calcForcesAndEnergy(true, false, -1);
+        vector<Vec3> forces;
+        impl->getForces(forces);
+        for (int i = 0; i < system.getNumParticles(); i++) {
+            double mass = system.getParticleMass(i);
+            if (mass != 0.0)
+                velocities[i] -= (offset/mass)*forces[i];
+        }
+    }
+    setVelocities(velocities);
     impl->applyVelocityConstraints(1e-5);
 }
 

plugins/drude/platforms/common/src/CommonDrudeKernels.cpp

@@ -121,12 +121,12 @@ void CommonCalcDrudeForceKernel::initialize(const System& system, const DrudeFor
             double k1 = ONE_4PI_EPS0*charge*charge/(polarizability*a1) - k3;
             double k2 = ONE_4PI_EPS0*charge*charge/(polarizability*a2) - k3;
             if (atoms[i][2] == -1) {
-                atoms[i][2] = 0;
+                atoms[i][2] = atoms[i][0];
                 k1 = 0;
             }
             if (atoms[i][3] == -1 || atoms[i][4] == -1) {
-                atoms[i][3] = 0;
-                atoms[i][4] = 0;
+                atoms[i][3] = atoms[i][0];
+                atoms[i][4] = atoms[i][0];
                 k2 = 0;
             }
             paramVector[i] = mm_float4((float) k1, (float) k2, (float) k3, 0.0f);

plugins/drude/tests/TestDrudeLangevinIntegrator.h

@@ -231,7 +231,7 @@ void testForceEnergyConsistency() {
 
 void testInitialTemperature() {
     // Check temperature initialization for a collection of randomly placed particles
-    const int numRealParticles = 500000;
+    const int numRealParticles = 50000;
     const int numParticles = 2 * numRealParticles;
     const int nDoF = 3 * numRealParticles;
     const double targetTemperature = 300;

plugins/drude/tests/TestDrudeNoseHoover.h

@@ -227,7 +227,7 @@ double testWaterBoxWithHardWallConstraint(double hardWallConstraint){
 
 void testInitialTemperature() {
     // Check temperature initialization for a collection of randomly placed particles
-    const int numRealParticles = 500000;
+    const int numRealParticles = 50000;
     const int numParticles = 2 * numRealParticles;
     const int nDoF = 3 * numRealParticles;
     const double targetTemperature = 300;

plugins/drude/tests/TestDrudeSCFIntegrator.h

@@ -135,7 +135,7 @@ void testWater() {
 
 void testInitialTemperature() {
     // Check temperature initialization for a collection of randomly placed particles
-    const int numRealParticles = 500000;
+    const int numRealParticles = 50000;
     const int numParticles = 2 * numRealParticles;
     const int nDoF = 3 * numRealParticles;
     const double targetTemperature = 300;

tests/TestBrownianIntegrator.h

@@ -253,7 +253,7 @@ void testRandomSeed() {
 
 void testInitialTemperature() {
     // Check temperature initialization for a collection of randomly placed particles
-    const int numParticles = 500000;
+    const int numParticles = 50000;
     const int nDoF = 3 * numParticles;
     const double targetTemperature = 300;
     System system;

tests/TestCustomIntegrator.h

@@ -1130,7 +1130,7 @@ void testRecordEnergy() {
 
 void testInitialTemperature() {
     // Check temperature initialization for a collection of randomly placed particles
-    const int numParticles = 500000;
+    const int numParticles = 50000;
     const int nDoF = 3 * numParticles;
     const double targetTemperature = 300;
     System system;

tests/TestLangevinIntegrator.h

@@ -260,7 +260,7 @@ void testRandomSeed() {
 
 void testInitialTemperature() {
     // Check temperature initialization for a collection of randomly placed particles
-    const int numParticles = 500000;
+    const int numParticles = 50000;
     const int nDoF = 3 * numParticles;
     const double targetTemperature = 300;
     System system;

tests/TestLangevinMiddleIntegrator.h

@@ -263,7 +263,7 @@ void testRandomSeed() {
 
 void testInitialTemperature() {
     // Check temperature initialization for a collection of randomly placed particles
-    const int numParticles = 500000;
+    const int numParticles = 50000;
     const int nDoF = 3 * numParticles;
     const double targetTemperature = 300;
     System system;

tests/TestNoseHooverIntegrator.h

@@ -271,7 +271,7 @@ void testThreeParticleVirtualSite() {
 
 void testInitialTemperature() {
     // Check temperature initialization for a collection of randomly placed particles
-    const int numParticles = 500000;
+    const int numParticles = 50000;
     const int nDoF = 3 * numParticles;
     const double targetTemperature = 300;
     System system;

tests/TestVariableLangevinIntegrator.h

@@ -332,7 +332,7 @@ void testArgonBox() {
 
 void testInitialTemperature() {
     // Check temperature initialization for a collection of randomly placed particles
-    const int numParticles = 500000;
+    const int numParticles = 50000;
     const int nDoF = 3 * numParticles;
     const double targetTemperature = 300;
     System system;

tests/TestVariableVerletIntegrator.h

@@ -311,7 +311,7 @@ void testArgonBox() {
 
 void testInitialTemperature() {
     // Check temperature initialization for a collection of randomly placed particles
-    const int numParticles = 500000;
+    const int numParticles = 50000;
     const int nDoF = 3 * numParticles;
     const double targetTemperature = 300;
     System system;

tests/TestVerletIntegrator.h

@@ -228,7 +228,7 @@ void testConstrainedMasslessParticles() {
 
 void testInitialTemperature() {
     // Check temperature initialization for a collection of randomly placed particles
-    const int numParticles = 500000;
+    const int numParticles = 50000;
     const int nDoF = 3 * numParticles;
     const double targetTemperature = 300;
     System system;