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

[Erik Lindahl]

Hi Claudio,

If you use PME the problem is probably the constraints combined with  
single precision. It should conserve energy in double, and with a  
little luck the next version should conserve energy in single even with  
constraints (I'm working on it, but hasn't been trivial to find what  
the accuracy-sensitive operation is - probably a combination of the  
integration and constraint).

Cheers,

Erik


On Tuesday, Mar 4, 2003, at 10:49 US/Pacific, Claudio J. Margulis 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