[gmx-users] energy drift due to LINCS

[David]

On Wed, 2003-09-24 at 17:49, Roman Martonak (CSCS) wrote:
> I am simulating ice with TIP4P using GROMACS and optimized order of atoms
> employing the tip4p.itp file posted here by Dr. van der Spoel some time
> ago. I found that in an NVE simulation the energy is very strongly
> non-conserved and the system cools down rapidly. While there are obviously
> several more subtle sources of energy non-conservation related to cutoff
> treatment, PME etc., here it seems that the main origin of the drift is
> LINCS. I tried to make a test using a non-optimized order of atoms where
> both SETTLE and LINCS can be used and the results were dramatically
> different (see below).
> 
> delta t = 2 fs
> run started from well equilibrated sample of Ih ice at T = 80 K
> 
> settle
> ******
> Statistics over 100001 steps [ 0.0000 thru 200.0000 ps ], 3 data sets
> 
> Energy                      Average       RMSD     Fluct.      Drift  Tot-Drift
> -------------------------------------------------------------------------------
> Potential                  -19664.6    18.6506    18.6502 -0.00197686  -0.395376
> Total Energy               -18826.5   0.965364   0.938974 -0.00388277  -0.776561
> Temperature                 93.5755    2.05744     2.0574 -0.000212522 -0.0425048
> 
> 
> lincs
> *****
> 
> order4
> ******
> Statistics over 100001 steps [ 0.0000 thru 200.0000 ps ], 3 data sets
> 
> Energy                      Average       RMSD     Fluct.      Drift  Tot-Drift
> -------------------------------------------------------------------------------
> Potential                  -20345.4    180.325    46.3176   -3.01851   -603.708
> Total Energy               -19906.4    348.354    85.6164   -5.84854   -1169.72
> Temperature                 48.9448    18.8072    4.65787  -0.315599   -63.1204
> 
> order8
> ******
> Statistics over 100001 steps [ 0.0000 thru 200.0000 ps ], 3 data sets
> 
> Energy                      Average       RMSD     Fluct.      Drift  Tot-Drift
> -------------------------------------------------------------------------------
> Potential                  -19778.6    44.8999    17.3167  -0.717516   -143.505
> Total Energy               -19018.9    79.8277    3.57521   -1.38126   -276.254
> Temperature                 84.7156    4.66053    1.85948 -0.0740186   -14.8039
> 
> While it can be seen that the effect can be reduced by using lincs_order =
> 8 instead of the default value of 4, I still find it way too large if the
> temperature drops by 15 % during 200 ps of an NVE run. My question thus is
> whether it is possible to use SETTLE together with the optimized order of
> atoms.

Not right now, maybe in the future. It's just a technical issue of
course. Maybe you could try to use shake with different tolerance as
well?

> 	 Thanks in advance,
> 				Roman Martonak.
-- 
David.
________________________________________________________________________
David van der Spoel, PhD, Assist. Prof., Molecular Biophysics group,
Dept. of Cell and Molecular Biology, Uppsala University.
Husargatan 3, Box 596,  	75124 Uppsala, Sweden
phone:	46 18 471 4205		fax: 46 18 511 755
spoel at xray.bmc.uu.se	spoel at gromacs.org   http://xray.bmc.uu.se/~spoel
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++