[gmx-users] Energy drift problem

[xzp02 at mails.tsinghua.edu.cn]

Dear Anton:

>> Dear All:
>>     I ran a simulation for 5ns, no cut-off for VDW interaction, but 
>> finally I found my system(1200 atoms) have a total energy drift of 
>> 60kJ/mol! I think it's quite large. How large is the tolenrance 
>> typically? and which parameter can cause it?

>60kJ drift on a total energy of 13000kJ over 10ns sounds pretty fine to me.
  Because my system contains two carbon nanotube, no solvent, and the dynamics
begin at 0K, so the 60kJ/mol take ~4K kevin into my system, which is noisy to me.
But the time scale is very long, I do not know how to control that.

>>     I have read some posts before, and found that nstlist and vdw 
>> switch/shift can cause this, but I have a 1ps test, with no VDW 
>> cut-off, nstlist = 1and 10 give exactly the same total energy, but VDW 
>> cut-off may take into some not very large drift.

>Some issues to consider:
>1) timestep in relation to dynamics of the system. for 'protein/water' like
>systems,  2 fs is fine, but 1 fs will give you better energy conservation.
>    (See: J. Comput. Chem. 20 (8), pp 786-798, 1999)
   I use 1fs, and I found your paper yesterday by google.:)
>2) cut-off. You have none for VDW (does that mean you don't use PBC as well?).
>    If you have charged atoms, you need to consider coulomb treatment, e.g. PME
>    will also give you energy drift. (P.S. without cut-offs the value of nstlist
>    becomes irrelevant as well - as you have already seen.)
   No PBC, and another PBC give almost the same energy drift. So it should not be
the reason.
>3) precision of the code. Single is fine for nearly all applications, but will
>    cause considerable energy drift. E.g, for my JCC paper (above) I had to use
>    double precision to produce the drift vs. timestep plots.
   I use single precision for speed.

> For most of these things it holds that if it is really important, you have to 
check, e.g. by making a 'drift vs dt' plot. Note, that energy drift in principle
is a diffusive process, so you should average the drift over several independent
simulations if you want to determine it accurately. (There is also a discussion on
this issue in my JCC paper - there is a link to the pdf reprint on my homepage
www.few.vu.nl/~feenstra )
   I will try.

>>     In the output of g_energy, which is the unit of Drift, kJ/mol? and 
>> is there some critical value for it?

>kJ/mol/ps. There isn't really a 'hard' critical value - that depends on the type
>of system you simulate and what you want to do with it.
   Thanks a lot for your helpful explanation, so maybe the drift can not be
suppressed during such a long time? And will the drift always be positive? because
I have seen some cases with negative total energy drift (decreasing). 


ZP