[gmx-users] Re: Possible bug in the temperature calculation from rerun

[Mark Abraham]

On 21/09/2012 11:08 PM, Bastien Loubet wrote:
> Justin Lemkul wrote
>> On 9/21/12 8:29 AM, Bastien Loubet wrote:
>>> Dear gmx users,
>>>
>>> We recently got a problem with the rerun feature of mdrun, and we request
>>> your help in order to help to solve it.
>>>
>>> We have run a simulation of a large POPC membrane using the coarse
>>> grained
>>> Martini force fields. From these simulations we obtained a trr trajectory
>>> file which contains both the position and the velocity of each martini
>>> particle every 20 ps. From this trajectory we rerun the simulation using
>>> the
>>> -rerun option of the mdrun program but with a different topology for
>>> which
>>> we have set all the bonded and non bonded interaction to zero.
>>> Specifically
>>> we have set all the force constant to zero in the martini_v2.P.itp and
>>> the
>>> martini_v2.0_lipids.itp file, the aim is to calculate the virial due to
>>> the
>>> electrostatic forces alone.
>>> The problem is the following: from the edr file we get from the rerun we
>>> extract the mean value of the energies using g_energy. For the edr file
>>> obtained during the simulation (and not the rerun) we obtained for the
>>> temperature and the kinetic energy, using also a 20ps time interval:
>>>
>>> Energy                      Average   Err.Est.       RMSD  Tot-Drift
>>> -------------------------------------------------------------------------------
>>> Temperature                 324.998     0.0049   0.413669 -0.00743801
>>> (K)
>>> Kinetic En.              1.9525e+06         29    2485.21   -44.6815
>>> (kJ/mol)
>>>
>>> This is to be expected because we have a Nose-Hover temperature coupling
>>> of
>>> 325 K. The kinetic energy is also equals to the number of degree of
>>> freedom
>>> times half the Boltzmann constant times the temperature.
>>> However for the rerun, with no bonded and non-bonded interactions, we
>>> get:
>>>
>>> Energy                      Average   Err.Est.       RMSD  Tot-Drift
>>> -------------------------------------------------------------------------------
>>> Temperature                  327.35     0.0053   0.415978 -0.00514402
>>> (K)
>>> Kinetic En.              1.96663e+06         32    2499.08   -30.8973
>>> (kJ/mol)
>>>
>>> The kinetic energy is again equals to the number of degree of freedom
>>> times
>>> half the Boltzmann constant times the temperature, however the
>>> temperature
>>> is 2 K off ! This is surprising as only the velocities and the masses of
>>> the
>>> particle should enter the kinetic energy. The velocities are taken from
>>> the
>>> trr trajectory file and we have checked using the Linux diff utility that
>>> our topology have the same masses than the original one.
>>> Another cause can be that we are doing the rerun on a local machine while
>>> the original simulation was run on a cluster, however we believe that
>>> this
>>> kind of numerical error cannot be responsible for a 2 degree Kelvin
>>> mistake.
>>> Furthermore we have rerun the simulation using the original topology and
>>> we
>>> obtained:
>>>
>>> Energy                      Average   Err.Est.       RMSD  Tot-Drift
>>> -------------------------------------------------------------------------------
>>> Temperature                 324.998      0.005   0.411827 -0.00761929
>>> (K)
>>> Kinetic En.              1.9525e+06         30    2474.15   -45.7788
>>> (kJ/mol)
>>>
>>> The slightly different value are certainly due to the previously
>>> mentioned
>>> fact that the reruns were not run on the same machine as the simulation
>>> and
>>> numerical rounding errors. This point out that it is really the topology
>>> that changes the value of the temperature. However we really have no
>>> masses
>>> differences between the two topologies and the velocities are draw from
>>> the
>>> same trr trajectory file.
>>>
>>> Any thought on this problem would be welcome.
>>>
>> Another thing to consider is the fact that an .edr file is accurate over
>> all
>> simulation steps, while reruns only do calculations present in the frames
>> supplied.  If you are using frames every 20 ps, that may only represent
>> less
>> than 1% of the total data collected in the simulation, depending on the
>> value of dt.
>>
>> -Justin
>>
>> -- 
>> ========================================
>>
>> Justin A. Lemkul, Ph.D.
>> Research Scientist
>> Department of Biochemistry
>> Virginia Tech
>> Blacksburg, VA
>> jalemkul[at]vt.edu | (540) 231-9080
>> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
>>
>> ========================================
>> -- 
>> gmx-users mailing list
>> gmx-users@
>> http://lists.gromacs.org/mailman/listinfo/gmx-users
>> * Please search the archive at
>> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
>> * Please don't post (un)subscribe requests to the list. Use the
>> www interface or send it to
>> gmx-users-request@
>> .
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
> Alright, but we did make a rerun with the same trajectory and the original
> topology and it gives the third results in my original post: the 2 Kelvin
> difference is not there. So it is really the different topology file that
> give the temperature difference I believe.

The .log file (and grompp) report the number of degrees of freedom per 
T-coupling group. Search for nrdf in the .log file. Does that differ? 
What differences are reported for a single trajectory frame?

Mark