[gmx-users] Berendsen Temperature Coupling and Heat

[Daniel S. Han]

Hello,
I am running a steered MD simulation in an NPT ensemble using Berendsen 
thermostats.  I expect that my system could heat up due to the external 
forces, and this heat would eventually be absorbed by the heat bath.

I thought I could get an estimate of the heat released by looking at the 
Temperature (T) and lambda values (velocity scaling factors) available 
in the gromacs energy output files, but I wasn't sure I was calculating 
it correctly.

Here's the algorithm I thought might work:
at every time step, we have a Temperature (T) and a Lambda (L) taken 
from the edr output file.

Where we can calculate the kinetic energy at any step as:
Kinetic Energy (KE) = .5 * N_degrees_of_freedom * K_boltzman * T
(which I'll call KE_output)

If I understand correctly, KE_output is the result after scaling by 
Lambda^2.  So if there were no temperature scaling, then the kinetic 
energy would be equal to:
KE_no_scaling = KE_output / (lambda*lambda)

Thus the heat exchanged between the system and the bath is:
Heat = KE_output - KE_no_scaling
where if heat is negative, it means that heat was released by the system 
and absorbed by the bath.

Using this algorithm I found some odd results, where I found that both 
my Steered MD simulations and a simulation of the same system without 
external forces tended to absorb heat.
I can post more details of that later, but for now, if anyone has any 
corrections/comments to this line of thought, it would be much appreciated.

Thanks in advance,
Danny.