[gmx-users] Huge acceleration needed to reproduce results!
Jennifer.Williams at ed.ac.uk
Mon Jul 27 20:00:50 CEST 2009
I am trying to reproduce a calculation that I carried out in DL_POLY.
It is to calculate the transport diffusion coefficient for CH4 in a
frozen mesoporous silica.
In DL_POLY I used an external force of 0.1 kJ mol-1 A-1. (0.1 KJ per
mole per angstrom). This equates to 10 dl_poly internal units which I
add in this way at each timestep;
Fsum = Fsum + Fex
In Gromacs, I want to apply the same force as I used in DL_POLY so I
calculated the required acceleration using F=ma. Where I took the mass
to be the mass of one molecule of methane (16 a.m.u).
The final value for acceleration that I came up with (which
corresponds to a force of 0.1kj mol-1 A-1 on each molecule) was 0.0625
The first hiccup was when I used this value, the MSD was negative
(though linear in the negative region of the graph). I assumed that
this had something to do with the orientation of the unit cell and
tried applying 0.0 0.0 -0.0625. The plot then looked much better.
The problem is when I calculate the Mean displacement of the CH4
molecules. (I do this using a slightly altered version of the g_msd
code). The Mean displacement from gromacs is very different to that
which I calculate using DL_POLY,
Gromacs gives 95.0, dl_poly 21347.0.
The MSD however (where I don?t add an acceleration are similar) so the
problem lies with the force I am adding.
To test that it wasn?t some bug in my code to calculate the Mean
displacement, I also looked at how the acceleration/force altered the
MSD in DL_POLY and gromacs.
In DL_POLY, adding an external force of 0.1kj mol-1 A-1 would change
the MSD of methane by 3 orders of magnitude compared to a run with no
My equivalent acceleration of -0.0625 in gromacs, in comparison,
barely changes the MSD from that of a run with no acceleration added.
In fact it takes an acceleration of -200 in the z direction to cause
such a difference in the Ds coefficients between runs with and without
Does anyone have any idea what is going on here? An acceleration of
200 ps nm-2 surely is not reasonable is it?. It seems very large
compared to the example in the manual of 0.1. This would then imply
that my back of the envelope calculation for relating force and
acceleration is wrong. Am I missing something? I?m quite sure that the
force I am adding in DL_POLY is equivalent to 0.1 kJ mol-1 A-1 so why
are my methane molecules moving so much less in gromacs in response to
the equivalent acceleration?
Also I noticed that although in the .mdp file I specify:
; Non-equilibrium MD stuff
acc-grps = CH4
accelerate = 0.0 0.0 -200.00
In the md.log file I get the following output
acc: 0 0 -154.549 0 0 45.4514
Can someone clarify what this means?
Any advice/comments appreciated
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.
More information about the gromacs.org_gmx-users