[gmx-users] Free energy calculation and Lincs warning

[Sadaf Rani]

Dear Gromacs users
I am doing free energy calculation of the protein-ligand system in gromacs
2020  in which I am annihilating ligand by removing charges and vdw
interactions. During Vdw removal I am facing lincs warning and my system
crashes by generating different pdb structures. I have tried to use
constraints on all bonds and also reduced time step to 1fs but still
getting same problem as below:-

Overriding thread affinity set outside gmx mdrun

WARNING: There are no atom pairs for dispersion correction
starting mdrun 'GROtesk MACabre and Sinister in water'
1000000 steps,   1000.0 ps.

Step 853159, time 853.159 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000569, max 0.041969 (between atoms 3430 and 3439)
bonds that rotated more than 30 degrees:
 atom 1 atom 2  angle  previous, current, constraint length
   3435   3437   38.0    0.1400   0.1415      0.1400
   3437   3438   52.2    0.1080   0.1068      0.1080
   3437   3439   39.0    0.1400   0.1402      0.1400

Step 853159, time 853.159 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000550, max 0.040414 (between atoms 3430 and 3439)
bonds that rotated more than 30 degrees:
 atom 1 atom 2  angle  previous, current, constraint length
   3435   3437   37.8    0.1400   0.1415      0.1400
   3437   3438   51.9    0.1080   0.1070      0.1080
   3437   3439   38.6    0.1400   0.1403      0.1400

Step 853160, time 853.16 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000833, max 0.054713 (between atoms 3437 and 3438)
bonds that rotated more than 30 degrees:
 atom 1 atom 2  angle  previous, current, constraint length
   3430   3439   30.5    0.1457   0.1428      0.1400

Please find the mdp file as below:-

; Run control
integrator               = sd       ; stochastic leap-frog integrator
dt                       = 0.001
nsteps                   = 1000000   ; 1 ns
comm-mode     = Linear        ; remove center of mass translation
nstcomm     = 100           ; frequency for center of mass motion removal


; Output control
nstxout                  = 5000
nstvout                  = 5000
nstfout                  = 0
nstlog                   = 5000
nstenergy                = 5000
nstxout-compressed       = 5000


; Bond parameters
continuation            = yes       ; formerly known as
'unconstrained-start' - useful for exact continuations and reruns
constraint_algorithm    = lincs     ; holonomic constraints
constraints             = all-bonds   ; bonds to H are constrained
lincs_iter              = 1         ; accuracy of LINCS
lincs-order             = 6    ; Highest order in the expansion of the
constraint coupling matrix


; Neighborsearching and short-range nonbonded interactions
cutoff-scheme            = verlet
nstlist                  = 10
ns_type                  = grid
rlist                    = 1.2
pbc                     = xyz       ; 3-D PBC ; Periodic boundary conditions

; Electrostatics
coulombtype              = PME
rcoulomb                 = 1.2
pme-order       = 6        ; interpolation order for PME (default is 4)
fourierspacing   = 0.10     ; grid spacing for FFT
ewald-rtol       = 1e-6     ; relative strength of the Ewald-shifted direct
potential at rcoulomb
ewald_geometry   = 3d       ; Ewald sum is performed in all three dimensions

; van der Waals
vdwtype                  = cutoff
vdw-modifier             = Potential-shift-Verlet
verlet-buffer-tolerance  = 0.005
rvdw                     = 1.2
rvdw-switch              = 1.0
DispCorr                 = EnerPres     ; apply long range dispersion
corrections for Energy and Pressure

; Temperature coupling
; tcoupl is implicitly handled by the sd integrator
tc-grps                 = system ; two coupling groups - more accurate
tau_t                   = 0.1                       ; time constant, in ps
ref_t                   = 298.15

; Pressure coupling is on for NPT
Pcoupl                  = Parrinello-Rahman
pcoupltype       = isotropic            ; uniform scaling of box vectors
tau_p           = 2.0                  ; time constant (ps)
ref_p             = 1.0                  ; reference pressure (bar)
compressibility   = 4.5e-05              ; isothermal compressibility of
water (bar^-1)

; velocities
gen_vel      = no      ; Velocity generation is on (if gen_vel is 'yes',
continuation should be 'no')
gen_seed     = -1       ; Use random seed
gen_temp     = 298.15

; Free energy control stuff
free_energy              = yes
init_lambda_state        = 17
delta_lambda             = 0
calc_lambda_neighbors    =-1        ; only immediate neighboring windows
; Vectors of lambda specified here
; Each combination is an index that is retrieved from init_lambda_state for
each simulation
 vdw_lambdas             = 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40
0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00
coul_lambdas             = 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
restraint_lambdas        = 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
; Options for the decoupling
sc-alpha                 = 0.5
sc-coul                  = no       ; linear interpolation of Coulomb (none
in this case)
sc-power                 = 1
sc-sigma                 = 0.3
nstdhdl                  = 100
dhdl-print-energy       = yes

Could you please suggest me how should I solve this problem?

Thanks in advance.

Sadaf