[gmx-users] energy drift due to LINCS

[Roman Martonak (CSCS)]

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.
	 Thanks in advance,
				Roman Martonak.
-- 
Dr. Roman Martonak                         Tel:    +41(91) 610 8226    
Swiss Center for Scientific Computing      Fax:    +41(91) 610 8282    
Via Cantonale - Galleria 2                 E-mail: martonak at cscs.ch    
CH-6928 Manno, Switzerland                 WWW:  http://www.cscs.ch/