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
--
gmx-users mailing list gmx-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
