[gmx-users] energy drift in nve with conservative parameters

Sun Mar 9 16:37:11 CET 2003

The problem with Ewald sums is that, although it has other great virtues, 
it is not strictly conservative, and with the approximations of PME that 
becomes a serious problem and pretty well precludes the use of NVE boundary 
conditions (at least that is our experience).  We considered some of this 
using just pure water as a test system in:

Chiu, S.-W., Clark, M., Subramaniam, S., and E. Jakobsson. 2000. Collective 
motion artefacts arising in long-duration molecular dynamics 
simulations.   J. Comp. Chem. 21:121-131

Eric

At 01:49 PM 3/4/2003 -0500, you wrote:
>Hi, I am trying to simulate a protein in flexible spc water and I believe 
>I have very conservative parameters. I however see a drift in the total 
>energy, and although it's not too large it will become large as I need the 
>simulation to be several ns long. If you plot the total energy you see 
>that the system is consistently heating up. I have equilibrated the system 
>for more than 100 ps doing npt and before that using position restraints,
>so everything should be fine ...
>
>time step is 1 fs ,  1nm for rvdw, fourierspacing is 0.12 for pme , i 
>update the neighbor list every step nstlist =1, etc ...
>the mdp file is attached.
>Do you think it's the PME? if so,  what should i use for fourier spacing 
>and tolerance to make it not drift? I am running single precision though.
>Thanks!!!
>Claudio.
>
>here is the output from g_energy statistics over 380 ps
>
>Energy                      Average       RMSD     Fluct.      Drift Tot-Drift
>-------------------------------------------------------------------------------
>Bond                        15717.5    254.829    252.943  -0.282489 -107.205
>Angle                         10366    194.383    193.899   0.125137  47.4897
>Proper Dih.                 122.033    12.9951    12.9539 0.00944353   3.58383
>Ryckaert-Bell.              1424.82    52.2058    52.0629 -0.0352289 -13.3694
>LJ-14                       1983.58    41.7575    41.6305 -0.0297046  -11.2729
>Coulomb-14                  11277.9    55.3076    53.6154   0.123927  47.0305
>LJ (SR)                     23216.1    415.996    415.932 -0.0663132  -25.1659
>Disper. corr.              -776.671          0          0          0    0
>Coulomb (SR)                -161371    627.813    626.229   0.406744   154.36
>Coulomb 
>(LR)                 -36995    51.9416    50.1115  -0.124746   -47.3413
>Potential                   -135035    251.587    251.203   0.126766 
>48.1079
>Kinetic 
>En.                 37008.7    245.956    245.244   0.170733    64.7933
>Total 
>Energy               -98026.1    35.0014    12.7625   0.297498    112.901
>
>
>--
>-----------------------------------------------
>Claudio Javier Margulis, Ph.D.
>Chemistry Department
>Columbia University
>Web http://www.chem.columbia.edu/~claudiom/
>Work phone (212) 854-8470
>-----------------------------------------------
>
>
>
>;;here no constraints and
>; nslist=1 nstcomm removed as parameter
>;notice that by default gromacs makes center of mass velocity removal
>; here we are trying flexible spc to see if we get better energy consevation
>;       User spoel (236)
>;       Wed Nov  3 17:12:44 1993
>;       Input file
>;
>title               =  Yo
>cpp                 =  /lib/cpp
>define              =  -DFLEX_SPC
>constraints         =  none
>integrator          =  md
>dt                  =  0.001    ; ps !
>nsteps              =  50000000 ; 50 ns
>;nstcomm             =  1 should not be there?
>nstxout             =  50000
>nstxtcout            =  1500
>nstvout             =  50000
>nstfout             =  0
>nstlog              =  5000
>nstenergy           =  1500
>nstlist             =  1
>ns_type             =  grid ;n-list parameter
>rlist               =  1.0 ;cutoff for the short range neighbor list
>
>
>; OPTIONS FOR ELECTROSTATICS AND VDW =
>; Method for doing electrostatics =
>coulombtype              = PME
>; Method for doing Van der Waals =
>vdw-type                 = Cut-off
>; cut-off lengths        =
>rvdw-switch              = 0
>rvdw                     = 1.0
>; Apply long range dispersion corrections for Energy and Pressure =
>DispCorr                 = EnerPres
>; Spacing for the PME/PPPM FFT grid =
>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               = 1e-05
>ewald_geometry           = 3d
>epsilon_surface          = 0
>optimize_fft             = yes
>; Nose-Hoover temperature coupling is off
>Tcoupl              = no
>tc-grps             =
>tau_t               =
>ref_t               =
>  Energy monitoring
>energygrps          =  Protein SOL Na
>; Pressure coupling is off
>Pcoupl              = no
>Pcoupltype          =
>tau_p               =
>compressibility     =
>ref_p               =
>; Generate velocites is on at 300 K.
>gen_vel             =  yes
>gen_temp            =  300.0
>gen_seed            =  173529

---------------------------------
Eric Jakobsson, Ph.D.
Professor, Department of Molecular and Integrative Physiology, and of 
Biochemistry
Senior Research Scientist, National Center for Supercomputing Applications
Professor, Beckman Institute for Advanced Science and Technology
4021 Beckman Institute, mc251
405 N. Mathews Avenue
University of Illinois, Urbana, IL 61801
ph. 217-244-2896       fax 217-244-2909

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20030309/b53ff12b/attachment.html>