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

[Jennifer Williams]

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  
force added.

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?

Any advice/comments appreciated



-- 
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.