Compute heat transfer and work with stochastic rescaling 93/4693/3 release-2-13-beta-2
authorradakb <brian.radak@gmail.com>
Tue, 16 Oct 2018 15:17:08 +0000 (11:17 -0400)
committerJim Phillips <jim@ks.uiuc.edu>
Tue, 16 Oct 2018 16:32:45 +0000 (11:32 -0500)
The stochRescale thermostat is now extended to enable computation
of heat transfer and work. The added overhead is incredibly small,
but for now reporting (to stdout and Tcl callback) only occurs with
the new stochRescaleHeat keyword. Other thermostats may also be able
to accomodate this computation, so the generic HEAT and WORK fields
have been added.

Both computations are _cumulative_ wrt firstTimestep, so resetting
from Tcl causes the accumulation to start over. Note that this also
requires tracking the TOTAL energy change between the current step
and firstTimestep.

For true equilibrium simulations, the cumulative work should be
rigorously zero. In practice one instead gets a small number that is
polynomial in the timestep (roughly quartic?).

For neMD/MC simulations this should be an appropriate route to
computing work for a Metropolis criterion. While this approach
mixes protocol and shadow components, test comparisons with direct
computation of the protocol component (CUMALCHWORK) show no obvious
differences within statistical precision. There are minor debates
about asymmetric (used by NAMD) versus symmetric integrator
splittings, so this assertion rests on its empirical validity.

Change-Id: Iddbf54ccd6edb78594705d287681ba4159b4a2e1

src/Controller.C
src/Controller.h
src/SimParameters.C
src/SimParameters.h

index 89ca427..7908a19 100644 (file)
@@ -183,6 +183,7 @@ Controller::Controller(NamdState *s) :
     random = new Random(simParams->randomSeed);
     random->split(0,PatchMap::Object()->numPatches()+1);
 
+    heat = totalEnergy0 = 0.0;
     rescaleVelocities_sumTemps = 0;  
     rescaleVelocities_numTemps = 0;
     stochRescale_count = 0;
@@ -1370,6 +1371,7 @@ void Controller::stochRescaleVelocities(int step)
                   + 2*R1*sqrt(tempfactor*(1 - stochRescaleTimefactor)*stochRescaleTimefactor)); 
       }
       broadcast->stochRescaleCoefficient.publish(step,coefficient);
+      heat += 0.5*numDegFreedom*BOLTZMANN*temperature*(coefficient*coefficient - 1.0);
       stochRescale_count = 0;
     }
   }
@@ -2989,6 +2991,7 @@ void Controller::printEnergies(int step, int minimize)
     BigReal potentialEnergy;
     BigReal flatEnergy;
     BigReal smoothEnergy;
+    BigReal work;
 
     Vector momentum;
     Vector angularMomentum;
@@ -3104,6 +3107,12 @@ void Controller::printEnergies(int step, int minimize)
     smoothEnergy = (flatEnergy + smooth2_avg
         - (4.0/3.0)*(kineticEnergyHalfstep - kineticEnergyCentered));
 
+    // Reset values for accumulated heat and work.
+    if (step <= simParameters->firstTimestep &&
+        (simParameters->stochRescaleOn && simParameters->stochRescaleHeat)) {
+      heat = 0.0;
+      totalEnergy0 = totalEnergy;
+    }
     if ( simParameters->outputMomenta && ! minimize &&
          ! ( step % simParameters->outputMomenta ) )
     {
@@ -3209,6 +3218,11 @@ void Controller::printEnergies(int step, int minimize)
       GET_VECTOR(pairElectForce,reduction,REDUCTION_PAIR_ELECT_FORCE);
     }
 
+    // Compute cumulative nonequilibrium work (including shadow work).
+    if (simParameters->stochRescaleOn && simParameters->stochRescaleHeat) {
+      work = totalEnergy - totalEnergy0 - heat;
+    }
+
     // callback to Tcl with whatever we can
 #ifdef NAMD_TCL
 #define CALLBACKDATA(LABEL,VALUE) \
@@ -3253,6 +3267,10 @@ void Controller::printEnergies(int step, int minimize)
         CALLBACKLIST("VDW_FORCE",pairVDWForce);
         CALLBACKLIST("ELECT_FORCE",pairElectForce);
       }
+      if (simParameters->stochRescaleOn && simParameters->stochRescaleHeat) {
+        CALLBACKLIST("HEAT",heat);
+        CALLBACKLIST("WORK",work);
+      }
       if (simParameters->alchOn) {
         if (simParameters->alchThermIntOn) {
           CALLBACKLIST("BOND1", bondedEnergy_ti_1);
@@ -3356,6 +3374,11 @@ void Controller::printEnergies(int step, int minimize)
          iout << FORMAT("PRESSAVG");
          iout << FORMAT("GPRESSAVG");
        }
+        if (simParameters->stochRescaleOn && simParameters->stochRescaleHeat) {
+          iout << "     ";
+          iout << FORMAT("HEAT");
+          iout << FORMAT("WORK");
+        }
         if (simParameters->drudeOn) {
           iout << "     ";
          iout << FORMAT("DRUDEBOND");
@@ -3451,6 +3474,11 @@ void Controller::printEnergies(int step, int minimize)
        iout << FORMAT(pressure_avg*PRESSUREFACTOR/avg_count);
        iout << FORMAT(groupPressure_avg*PRESSUREFACTOR/avg_count);
     }
+    if (simParameters->stochRescaleOn && simParameters->stochRescaleHeat) {
+          iout << "     ";
+          iout << FORMAT(heat);
+          iout << FORMAT(work);
+    }
     if (simParameters->drudeOn) {
         iout << "     ";
        iout << FORMAT(drudeBondTemp);
index a63a88d..55a2c38 100644 (file)
@@ -158,6 +158,9 @@ protected:
       BigReal kineticEnergyHalfstep;
       BigReal kineticEnergyCentered;
       BigReal temperature;
+      BigReal heat;
+      /**< heat exchanged with the thermostat since firstTimestep */
+      BigReal totalEnergy0; /**< totalEnergy at firstTimestep */
       // BigReal smooth2_avg;
       BigReal smooth2_avg2;  // avoid internal compiler error
       Tensor pressure;
index 21aeeb2..e281623 100644 (file)
@@ -1315,6 +1315,8 @@ void SimParameters::config_parser_methods(ParseOptions &opts) {
        "Number of steps between stochastic rescalings",
         &stochRescaleFreq);
    opts.range("stochRescaleFreq", POSITIVE);
+   opts.optionalB("stochRescale", "stochRescaleHeat",
+       "Should heat transfer and work be computed?", &stochRescaleHeat, FALSE);
 
    opts.optional("main", "rescaleFreq", "Number of steps between "
     "velocity rescaling", &rescaleFreq);
index efdc029..c160da1 100644 (file)
@@ -595,6 +595,11 @@ public:
   int stochRescaleFreq;
   /**< How frequently (time steps) stochastic velocity rescaling occurs. */
 
+  Bool stochRescaleHeat;
+  /**< Flag TRUE enables calculation and reporting of heat transfer and work.
+    *  The computation is _cumulative_ from step = firstTimestep.
+    */
+
        int rescaleFreq;                //  Velocity rescale frequency
        BigReal rescaleTemp;            //  Temperature to rescale to