How to solve Lincs warning in free energy calculation

GROMACS version:2020
GROMACS modification: Yes/No
Here post your question Dear Gromacs users
I am trying to calculate the free energy of binding of a ligand to protein in a protein-ligand complex using dual topology approach. I am using Amer-99SB-ILDN force field for protein and GAFF for the ligand and constraints for h-bonds and a time step of 1 fs. In order to keep the ligand in place, I put a few restraints on some atoms of ligand and protein. During the annihilation of ligand, my system crashes while generating different PDB structures with the following error message:-
Step 341571, time 341.571 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000002, max 0.000007 (between atoms 4477 and 4478)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
3437 3438 52.4 0.1080 0.1080 0.1080

Step 341571, time 341.571 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000002, max 0.000007 (between atoms 3459 and 3461)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
3437 3438 52.1 0.1080 0.1080 0.1080

Step 341572, time 341.572 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000002, max 0.000007 (between atoms 5248 and 5250)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
3437 3438 46.5 0.1080 0.1080 0.1080

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 = h-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 = yesmd_17.mdp (3.6 KB)
Could you please help me to fix this problem?
Thanks.

I see that all your coulomb lambda points are 1.0, you are not using fep-lambda for other interactions, meaning that you will keep the coulomb and all other (e.g. bonded, dihedral) interactions while only annihilating the van-der-Walls interactions - could it be that this a reason for the instability in your system?

Thanks for your reply.
I am annihilating the ligand in two steps firstly I changed coumbic lambda from 0 to 1 keeping the Vdw switched on. Then in step two I am removing Vdw while keeping chargeg off.
I am not sure if this could be the reason.
Any further suggestion will help me to fix this issue.
Thanks in advance.
Sadaf

I am not very much sure about the fep-lambda. Can you please give a more detail.
Thanks

You have not specified coulomb-lambda0 or coulomb-lambda1 in your .mdp file, so these will both default to vdw-q, meaning charges are on in the B-state (lambda = 1).

Dear Justin
Thank you for your reply.

I am treating the whole of the system in a single molecule type so that I can define restraints between ligand and protein atoms. For this, I assume that I should not include these two parameters as I have defined a dual topology for my ligand in the topology file. In transforming ligand from lambda 0 to 1 means that the ligand is being transformed from state A (fully interaction) to a dummy entity. Please correct me if I am wrong. Please find the topology file attached.

Thanks in advance.

Sadaf

topol.top (2.18 MB)

Dear Gromacs users
Any comment on this please. It will really help me to solve the problem.
Thanks.
Sadaf