computeDrudeTemperature (#3512)

* computeDrudeTemperature

* cleanup

* revert drude temperature computation in DrudeNoseHoover

plugins/drude/openmmapi/include/openmm/DrudeLangevinIntegrator.h

@@ -126,6 +126,12 @@ public:
      * getTemperature().
      */
     double computeSystemTemperature();
+    /**
+     * Compute the instantaneous temperature of the Drude system, measured in Kelvin.
+     * This is calculated based on the kinetic energy of the internal motion of Drude pairs
+     * and should remain close to the prescribed Drude temperature.
+     */
+    double computeDrudeTemperature();
 protected:
     /**
      * This will be called by the Context when it is created.  It informs the Integrator

plugins/drude/openmmapi/include/openmm/DrudeNoseHooverIntegrator.h

@@ -106,6 +106,12 @@ public:
      * getTemperature().
      */
     double computeSystemTemperature();
+    /**
+     * Compute the instantaneous temperature of the Drude system, measured in Kelvin.
+     * This is calculated based on the kinetic energy of the internal motion of Drude pairs
+     * and should remain close to the prescribed Drude temperature.
+     */ 
+     double computeDrudeTemperature();
 protected:
     /**
      * Return a list of velocities normally distributed around a target temperature, with the Drude

plugins/drude/openmmapi/src/DrudeHelpers.cpp

@@ -127,12 +127,16 @@ vector<Vec3> assignDrudeVelocities(const System &system, double temperature, dou
     return velocities;
 }
 
-double computeSystemTemperatureFromVelocities(const System& system, const vector<Vec3>& velocities) {
+
+/**
+ * Computes the instantaneous temperatures of the system and the internal Drude motion and returns a pair (T_system, T_drude)
+ */
+pair<double, double> computeTemperaturesFromVelocities(const System& system, const vector<Vec3>& velocities) {
     vector<int> normalParticles;
     vector<pair<int, int> > pairParticles;
     findParticlesAndPairs(system, normalParticles, pairParticles);
-    double energy = 0.0;
-    int dof = 0;
+    double energy = 0.0, drudeEnergy = 0;
+    int dof = 0, drudeDof = 0;
     
     // Kinetic energy and degrees of freedom from normal particles.
     
@@ -151,10 +155,14 @@ double computeSystemTemperatureFromVelocities(const System& system, const vector
         int p2 = pair.second;
         double mass1 = system.getParticleMass(p1);
         double mass2 = system.getParticleMass(p2);
+	double reducedMass = (mass1*mass2)/(mass1+mass2);
         if (mass1 != 0 || mass2 != 0) {
             Vec3 momentum = mass1*velocities[p1] + mass2*velocities[p2];
             energy += momentum.dot(momentum)/(mass1+mass2);
+	    Vec3 drudeVelocity = velocities[p1] - velocities[p2];
+            drudeEnergy += reducedMass*drudeVelocity.dot(drudeVelocity);
 	    dof += 3;
+	    drudeDof += 3;
         }
     }
     
@@ -176,7 +184,10 @@ double computeSystemTemperatureFromVelocities(const System& system, const vector
             break;
         }
     energy *= 0.5;
-    return 2*energy/(dof*BOLTZ);
+    drudeEnergy *= 0.5;
+    if (drudeDof == 0)
+        drudeDof = 1; // so that the drude temperature is reported as 0
+    return make_pair<double, double>(2*energy/(dof*BOLTZ), 2*drudeEnergy/(drudeDof*BOLTZ));
 }
 
 }

plugins/drude/openmmapi/src/DrudeLangevinIntegrator.cpp

@@ -42,9 +42,10 @@
 using namespace OpenMM;
 using std::string;
 using std::vector;
+using std::pair;
 
 namespace OpenMM {
-    double computeSystemTemperatureFromVelocities(const System& system, const vector<Vec3>& velocities);
+	pair<double, double> computeTemperaturesFromVelocities(const System& system, const vector<Vec3>& velocities);
 }
 
 DrudeLangevinIntegrator::DrudeLangevinIntegrator(double temperature, double frictionCoeff, double drudeTemperature, double drudeFrictionCoeff, double stepSize) : DrudeIntegrator(stepSize) {
@@ -126,5 +127,15 @@ double DrudeLangevinIntegrator::computeSystemTemperature() {
     context->calcForcesAndEnergy(true, false, getIntegrationForceGroups());
     vector<Vec3> velocities;
     context->computeShiftedVelocities(getVelocityTimeOffset(), velocities);
-    return computeSystemTemperatureFromVelocities(context->getSystem(), velocities);
+    return computeTemperaturesFromVelocities(context->getSystem(), velocities).first;
 }
+
+double DrudeLangevinIntegrator::computeDrudeTemperature() {
+    if (context == NULL)
+        throw OpenMMException("This Integrator is not bound to a context!");  
+    context->calcForcesAndEnergy(true, false, getIntegrationForceGroups());
+    vector<Vec3> velocities;
+    context->computeShiftedVelocities(getVelocityTimeOffset(), velocities);
+    return computeTemperaturesFromVelocities(context->getSystem(), velocities).second;
+} 
+

plugins/drude/openmmapi/src/DrudeNoseHooverIntegrator.cpp

@@ -30,6 +30,7 @@
  * -------------------------------------------------------------------------- */
 
 #include "openmm/DrudeNoseHooverIntegrator.h"
+#include "SimTKOpenMMRealType.h"
 #include "openmm/Context.h"
 #include "openmm/OpenMMException.h"
 #include "openmm/internal/ContextImpl.h"
@@ -45,12 +46,14 @@
 using namespace OpenMM;
 using std::string;
 using std::vector;
+using std::pair;
 
 namespace OpenMM {
     extern std::vector<Vec3> assignDrudeVelocities(const System &system, double temperature, double drudeTemperature, int randomSeed);
-    double computeSystemTemperatureFromVelocities(const System& system, const vector<Vec3>& velocities);
+    pair<double, double> computeTemperaturesFromVelocities(const System& system, const vector<Vec3>& velocities);
 }
 
+
 DrudeNoseHooverIntegrator::DrudeNoseHooverIntegrator(double temperature, double collisionFrequency, 
                                                      double drudeTemperature, double drudeCollisionFrequency,
                                                      double stepSize, int chainLength, int numMTS, int numYoshidaSuzuki) :
@@ -151,8 +154,20 @@ double DrudeNoseHooverIntegrator::computeSystemTemperature() {
     context->calcForcesAndEnergy(true, false, getIntegrationForceGroups());
     vector<Vec3> velocities;
     context->computeShiftedVelocities(getVelocityTimeOffset(), velocities);
-    return computeSystemTemperatureFromVelocities(context->getSystem(), velocities);
+    return computeTemperaturesFromVelocities(context->getSystem(), velocities).first;
+}
+
+double DrudeNoseHooverIntegrator::computeDrudeTemperature() {
+    if (context == NULL)
+        throw OpenMMException("This Integrator is not bound to a context!");  
+    double kE = computeDrudeKineticEnergy();
+    size_t num_dofs = 0;
+    for (const auto &nhc: noseHooverChains){
+       num_dofs += 3 * nhc.getThermostatedPairs().size();
+    }
+    return kE/(0.5 * num_dofs * BOLTZ);
 }
+
 std::vector<Vec3> DrudeNoseHooverIntegrator::getVelocitiesForTemperature(const System &system, double temperature,
                                                                          int randomSeedIn) const {
     return assignDrudeVelocities(system, temperature, drudeTemperature, randomSeedIn);

plugins/drude/tests/TestDrudeLangevinIntegrator.h

@@ -84,7 +84,9 @@ void testSinglePair() {
         Vec3 velCM = vel[0]*(mass1/totalMass) + vel[1]*(mass2/totalMass);
         keCM += 0.5*totalMass*velCM.dot(velCM);
         Vec3 velInternal = vel[0]-vel[1];
-        keInternal += 0.5*reducedMass*velInternal.dot(velInternal);
+        double _keInternal = 0.5*reducedMass*velInternal.dot(velInternal);
+        keInternal += _keInternal;
+	ASSERT_EQUAL_TOL(integ.computeDrudeTemperature(), _keInternal/(3*0.5*BOLTZ), 1e-6);
         Vec3 delta = state.getPositions()[0]-state.getPositions()[1];
         double distance = sqrt(delta.dot(delta));
         ASSERT(distance <= maxDistance*(1+1e-6));

plugins/drude/tests/TestDrudeNoseHoover.h

@@ -167,6 +167,7 @@ void testWaterBox() {
         double drudeKE = integ.computeDrudeKineticEnergy();
         double temp = KE/(0.5*numStandardDof*BOLTZ);
         double drudeTemp = drudeKE/(0.5*numDrudeDof*BOLTZ);
+	ASSERT_EQUAL_TOL(drudeTemp, integ.computeDrudeTemperature(), 1e-6);
         meanTemp = (i*meanTemp + temp)/(i+1);
         meanDrudeTemp = (i*meanDrudeTemp + drudeTemp)/(i+1);
         double heatBathEnergy = integ.computeHeatBathEnergy();

wrappers/python/src/swig_doxygen/swigInputConfig.py

@@ -226,6 +226,7 @@ UNITS = {
 ("*", "getGlobalParameterDefaultValue") : (None, ()),
 ("*", "getPermutationMode") : (None, ()),
 ("*", "computeSystemTemperature") : ("unit.kelvin", ()),
+("*", "computeDrudeTemperature") : ("unit.kelvin", ()),
 ("LocalCoordinatesSite", "getOriginWeights") : (None, ()),
 ("LocalCoordinatesSite", "getXWeights") : (None, ()),
 ("LocalCoordinatesSite", "getYWeights") : (None, ()),