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

[Bastien Loubet]

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.



--
View this message in context: http://gromacs.5086.n6.nabble.com/Possible-bug-in-the-temperature-calculation-from-rerun-tp5001194p5001196.html
Sent from the GROMACS Users Forum mailing list archive at Nabble.com.