[gmx-users] Re: Energy Drift in Gromacs 4.0

[Mark Abraham]

Ilya Chorny wrote:
> Forgot column labels on energies. 
> 
>                                                          Average       
> RMSD     Fluct.      Drift  Tot-Drift
> 2 fs time step Total Energy                46497.2    36445.6     2630.7 
>   -148.833    -148833
> 4 fs time step Total Energy                41580.8    37693.4    2856.08 
>   -130.198    -130198
> 
> On Tue, Mar 10, 2009 at 7:17 PM, Ilya Chorny <ichorny at gmail.com 
> <mailto:ichorny at gmail.com>> wrote:
> 
>     Hello All,
> 
>     I am trying to run some calibration calculations with 2/4 fs time
>     steps. I am trying to reproduce the results in the Gromacs 4.0 paper
>     on a protein/water (not the same protein) system with ~100K atoms. I
>     ran 1 ns simulation in the NVE ensemble. My mdp params are shown
>     below. Topolgy was created using -vsite h. My system has 218441
>     degrees of freedom according to GMXDUMP. My P.E. energy is ~ -1E5
>     and has not full relaxed after 1ns but is close. 
> 
>     2 fs time step Total Energy                46497.2    36445.6    
>     2630.7   -148.833    -148833
>     4 fs time step Total Energy                41580.8    37693.4  
>      2856.08   -130.198    -130198
> 
>     Drift (kt/ns) 2fs/4fs = .34/.29
> 
>     The drift seems to be much larger for both the 2fs and 4fs time
>     steps then in the paper. Do I have a clear bug in my params? Should
>     I wait till the system fully relaxes before doing the drift calculation?

You should definitely equilibrate the system before trying to measure 
energy drift. I suggest doing so in NPT to fix your density and 
temperature, then switching to NVE, wait a bit, and only then start to 
collect your data.

>     coulombtype              = PME
>     rcoulomb                 = 1.0
>     fourierspacing           = 0.12
>     ; FFT grid size, when a value is 0 fourierspacing will be used
>     fourier_nx               = 0
>     fourier_ny               = 0
>     fourier_nz               = 0
>     ; EWALD/PME/PPPM parameters
>     pme_order                = 4
>     ewald_rtol               = 1.e-04

I would think the ewald_rtol value used in the paper was 1e-5, but it is 
unfortunately not stated. The different numbers will lead to different 
costs and accuracies for the direct and FFT parts of the PME algorithm.

rcoulomb, fourier_whatever, pme_order and ewald_rtol are all relevant if 
you're trying to reproduce results.

Mark