[gmx-users] Huge acceleration needed to reproduce results!

Jennifer Williams Jennifer.Williams at ed.ac.uk
Mon Jul 27 20:00:50 CEST 2009

Hello users,

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
nm ps-2.

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
acceleration.

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?