[gmx-users] Re: Possible bug in the temperature calculation from rerun
Dear gmxusers,
Thank you for your answers so far. I have run a few more test rerun and here
is what I got:
-The number of degree of freedom are the same whatever the topology I run
with, both in the log file and from the grompp output.
-Comparing the trajectories created from rerun with gmxdump and diff show
only minor differences, such as atoms positioned on one side or the other of
the simulation box. The velocities are strictly the same however. Here is an
example:
8836288,8836290c8836288,8836290
< x[394728]={ 1.32183e+02, 3.15341e+00, 5.17869e+00}
< x[394729]={ 1.32054e+02, 3.10839e+00, 5.14599e+00}
< x[394730]={ 1.32314e+02, 3.17612e+00, 5.13366e+00}
---
> x[394728]={-7.17163e-04, 3.15341e+00, 5.17869e+00}
> x[394729]={-1.29196e-01, 3.10839e+00, 5.14599e+00}
> x[394730]={ 1.29883e-01, 3.17612e+00, 5.13366e+00}
-If I use another topology where I only have the interaction with water
not set to zero (still with electrostatic), I get a temperature of ~326K
intermediate between the desired value of 325K and the value I got with all
the interaction shutdown of ~327K.
-Turning the temperature/pressure coupling on or off do not change the
results.
I think there is a problem in the calculation of the kinetic energy that
translate to the temperature. However setting all the interaction but the
electrostatic to zero should not modify the kinetic energy calculated from
'mdrun -rerun'.
I do not know how to proceed from here, maybe I should post my mdp top and
log files here or on the developers forum ?
Have a good day,
Bastien Loubet
--
View this message in context:
http://gromacs.5086.n6.nabble.com/Possible-bug-in-the-temperature-calculation-from-rerun-tp5001194p5001508.html
Sent from the GROMACS Users Forum mailing list archive at Nabble.com.
--
gmx-users mailing listgmx-users@gromacs.org
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-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] Re: Possible bug in the temperature calculation from rerun
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 292485.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 322499.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 302474.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 -- g
[gmx-users] Re: Possible bug in the temperature calculation from rerun
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 292485.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 322499.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 302474.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 dif
