[gmx-users] Increase in kinetic energy in LINCS with smaller time step

[Nakamura, Hideya]

Dear Gromacs users,

I have a question about an effect of time-step on results of
simulation with LINCS.

Now, I am running on a MD simulation of a rigid molecule in vacuum
under NVE ensemble.
No temperature coupling was used. No inter-atomic potential was used.
The rigid molecule is consisting of just three atoms and the geometry
of the molecule is triangle.
We are using LINCS to constrain all bond length.

We have run simulations with different time steps.
The initial velocities of each atom was defined in .gro and the
initial velocities were identical at different times steps.

In calculation results, the bond lengths were successfully constrained
at each time step.
However, when the smaller time step, the kinetic energy of the rigid
molecule increased with an increase in simulation time.
(At other time steps, the kinetic energies were almost constant.)

I cannot understand why the kinetic energy was increased at "smaller" time step.
Is this derived from LINCS algorithm, or not?

Please give me a suggestion.

Thanks in advance for your help.

Hideya Nakamura


------------------------------------------------------------------------------------------------
The .mdp file and .top file were shown as follows:
------------------------------------------------------------------------------------------------

; RUN CONTROL PARAMETERS =
integrator               = md
; start time and timestep in ps =
tinit                    = 0
dt                       = 0.00002 ;[ps]
nsteps                   = 250000000 ;[STEP]  ;
;nsteps                   = 1000 ;[STEP]
; number of steps for center of mass motion removal =
;nstcomm = 1  ;(1)
; mode for center of mass motion removal
comm-mode                = None

; NEIGHBORSEARCHING PARAMETERS =
; nblist update frequency =
nstlist                  = 10 ;[step]
; ns algorithm (simple or grid) =
ns_type                  = grid
; nblist cut-off         = [nm]
rlist                    = 1.5
; Periodic boundary conditions: xyz or none =
pbc                      = xyz

;OPTIONS FOR TEMPERATURE COUPLING
tcoupl = no
tc_grps                  = TRI
tau_t                    = 0 ;[ps]
ref_t                    = 350 ;[K]
;OPTIONS FOR PRESSURE COUPLING
Pcoupl                   = no
;Pcoupl                   = Parrinello-Rahman
Pcoupltype               = isotropic   ;useful for membrane
tau_p                    = 1.5        ;[ps]
compressibility          = 4.5e-05 ;[bar^-1]
ref_p                    = 1.01325        ;[bar]

; GENERATE VELOCITIES FOR STARTUP RUN =
gen_vel                  = no  ;yes or no
gen_temp                 = 323 ;[K]
gen_seed                 = 173529 ;(173529)

; OPTIONS FOR BONDS     =
constraints              = none  ;all-bonds/none
; Type of constraint algorithm =
constraint-algorithm     = LINCS
; Do not constrain the start configuration =
unconstrained-start      = no
; Highest order in the expansion of the constraint coupling matrix =
lincs-order              = 4   ;(4)
; Lincs will write a warning to the stderr if in one step a bond =
; rotates over more degrees than =
lincs-warnangle          = 30  ;(30)
------------------------------------------------------------------------------------------------------------------------
[ system ]
; name
1 triangle in Vacuum

[ moleculetype ]
; name  nrexcl
TRI  1

[ defaults ]
1 1 yes 0.5 0.5

[ atomtypes ]
;name  at.num      mass        charge   ptype       c6           c12
   AA   79     196.96655       0.000       A   0.00000E+00   0.00000E+00

[ atoms ]
;   nr   type   resnr  residu  atom    cgnr    charge
    1    AA        1    TRI     1AA        1        0.000
    2    AA        1    TRI     2AA        2        0.000
    3    AA        1    TRI     3AA        3        0.000

[ constraints ]
;  ai    aj funct        b0
    1    2     1     0.408
    1    3     1     0.408
    2    3     1     0.408

[ exclusions ]
1  2  3
2  1  3
3  1  2

[ molecules ]
; name number
TRI     1

------------------------------------------------------------------------------------------------------------------------