How to solve Lincs warning in free energy calculation

User discussions
1

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.

2

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?

1 Like
3

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

4

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

5

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).

6

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)

7

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

8

Hello all,
I have run into similar issues where for lambda state 11,12 and 13 I encountered LINCS warning and the NVT equilibration did not go through. These states will correspond to the coul-lambdas switching themselves on. I am trying to compute absolute binding free energy of a protein ligand pair. The ligand is not charged (at least that is what antechamber/acpype said). I am pasting an example MDP file. I have been following the alchemical tutorial to reach this point as described here.

##############################
;====================================================
; Production simulation
;====================================================

;----------------------------------------------------
; RUN CONTROL
;----------------------------------------------------
integrator = sd ; stochastic leap-frog integrator
nsteps = 500000 ; 2 * 500,000 fs = 1000 ps = 1 ns
dt = 0.002 ; 2 fs
comm-mode = Linear ; remove center of mass translation
nstcomm = 100 ; frequency for center of mass motion removal

;----------------------------------------------------
; OUTPUT CONTROL
;----------------------------------------------------
nstxout = 0 ; donā€™t save coordinates to .trr
nstvout = 0 ; donā€™t save velocities to .trr
nstfout = 0 ; donā€™t save forces to .trr
nstxout-compressed = 1000 ; xtc compressed trajectory output every 1000 steps (2 ps)
compressed-x-precision = 1000 ; precision with which to write to the compressed trajectory file
nstlog = 1000 ; update log file every 2 ps
nstenergy = 1000 ; save energies every 2 ps
nstcalcenergy = 100 ; calculate energies every 100 steps

;----------------------------------------------------
; BONDS
;----------------------------------------------------
constraint_algorithm = lincs ; holonomic constraints
constraints = all-bonds ; hydrogens only are constrained
lincs_iter = 1 ; accuracy of LINCS (1 is default)
lincs_order = 4 ; also related to accuracy (4 is default)
lincs-warnangle = 30 ; maximum angle that a bond can rotate before LINCS will complain (30 is default)
continuation = yes ; formerly known as ā€˜unconstrained-startā€™ - useful for exact continuations and reruns

;----------------------------------------------------
; NEIGHBOR SEARCHING
;----------------------------------------------------
cutoff-scheme = Verlet
ns-type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs (default is 10)
rlist = 1.0 ; short-range neighborlist cutoff (in nm)
pbc = xyz ; 3D PBC

;----------------------------------------------------
; ELECTROSTATICS
;----------------------------------------------------
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm)
ewald_geometry = 3d ; Ewald sum is performed in all three dimensions
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

;----------------------------------------------------
; VDW
;----------------------------------------------------
vdw-type = PME
rvdw = 1.0
vdw-modifier = Potential-Shift
ewald-rtol-lj = 1e-3
lj-pme-comb-rule = Geometric
DispCorr = EnerPres

;----------------------------------------------------
; TEMPERATURE & PRESSURE COUPL
;----------------------------------------------------
tc_grps = System
tau_t = 1.0
ref_t = 300
pcoupl = Parrinello-Rahman
pcoupltype = isotropic ; uniform scaling of box vectors
tau_p = 2 ; time constant (ps)
ref_p = 1.0 ; reference pressure (bar)
compressibility = 4.5e-05 ; isothermal compressibility of water (bar^-1)

;----------------------------------------------------
; VELOCITY GENERATION
;----------------------------------------------------
gen_vel = no ; Velocity generation is off

;----------------------------------------------------
; FREE ENERGY CALCULATIONS
;----------------------------------------------------
free-energy = yes
couple-moltype = ligand
couple-lambda0 = vdw-q
couple-lambda1 = none
couple-intramol = no
separate-dhdl-file = yes
sc-alpha = 0.5
sc-power = 1
sc-sigma = 0.3
init-lambda-state = 11
bonded-lambdas = 0.0 0.01 0.025 0.05 0.075 0.1 0.2 0.35 0.5 0.75 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.0 1.0 1.0 1.0 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0
coul-lambdas = 0.0 0.00 0.000 0.00 0.000 0.0 0.0 0.00 0.0 0.00 0.0 0.25 0.5 0.75 1.0 1.00 1.0 1.0 1.0 1.0 1.0 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0
vdw-lambdas = 0.0 0.00 0.000 0.00 0.000 0.0 0.0 0.00 0.0 0.00 0.0 0.00 0.0 0.00 0.0 0.05 0.1 0.2 0.3 0.4 0.5 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1.0
nstdhdl = 100
calc-lambda-neighbors = -1
########################################

Any help will be very much appreciated as to how to solve this critical problem.

Best,
Amit

9

Ok ā€¦ I think may be I had the wrong idea. Please correct me if I am wrong.
The direction from lamba=1 to lambda=30 is of anihilation/decoupling and not the other way around (lambda=30 to lambda=1). If that is the case as soon as the electrostatics are switched off (and not switched on) from lambda=11 and higher, the major force that is keeping the complex intact (and is electrostatic in nature) is no more and the ligand falls out and starts to rotate and therefore the LINCS warning. If that is the case how do we solve this problem of calculating the free energy? ā€¦or there is something else that I am not aware of.

Update:
Ok I looked in the log (something I should have done first) and checked the bonds that had rotated by more than the margin. There is an sp3 carbon and all four atoms had rotated as per warnings. Also to this central atom are two rings attached and all the bonds involving hydrogen and carbon/nitrogen had rotated by huge values. As I was suspecting this region is also the region that is involved in hydrogen bonds. However I am not sure if it has to do anything with the increasing coul-lambdas values.

Even though more informed about the cause of the errors/warnings, I am still not sure how to go about solving this problem.

10

@amit7urmc , @sadaf , Have you solved this? Iā€™ve encountered a smiliar problem.

11

@kexul I could not solve it. Every time I get the hint that my topology may be wrong. I tried to rerun few lambda windows worked but for few I get my pdb structure crashed.
If anyone can have a suggestion please guide us.
Regards
Sadaf

12

@kexul I am sorry to hear that you are facing the same problemā€¦ to some extent I was able to solve it by reducing the dt by half. However it came back later for other lambdas. May be the ligand was in strain, since I was dealing with a docked structure and may be I did not generate correct topology of the ligand ā€¦ but this is only a guess ā€¦ An expert can tell you better about this reason.

13

@kexul @amit7urmc
I also reduced the dt to half but got the same as @amit7urmc
Didnā€™t get any further suggestion in the mailing list.
Regards
Sadaf

14

Digging back into this thread, I do not even understand what is going on with @sadafā€™s system. There is no [moleculetype] that has been specified for transformation and the ligands from GAFF are merged into the protein [molecueltype]. This isnā€™t a viable approach for transforming a ligand, which must be defined separately in the topology and coupled to Ī».

In general, for randomly crashing systems, the ligand topology is to blame. The troubleshooting protocol is to run the ligand in water by itself (without free energy options) to see if itā€™s stable, then attempt the transformation. If the first simulation fails, the topology is unsound and needs to be refined. If the free energy simulations fail, then investigate the minimization and equilibration protocol, and potentially refine the Ī» interval.

15

@sadaf @amit7urmc Iā€™ve raised an issue in the third party software Iā€™ve used to do free energy calculation, and got some help from experts, hope this give some clues.

16

@jalemkul,

Do we have to actually do all the energy minimization, NVT, and NPT equilibration for each init_lambda_state before run production dynamics? i.e. if I have for instance 20 states and I want to compare similar equilibrations that I am doing in Tinker such as EM, NVT @ 50K, NVT @ 100K, NVT @ 200K, NVT @ 300K, and NPT @ 300K, then I assumed I could do this at my initial state where all my lambdas are 0 and then run production on that NPT equilibrated structure for all my lambda states.

@sadaf posted his/her NVT.mdp file and it was for state 158 so this makes me second guess my assumption.

I seem to only be getting these errors for states where vdw-lambdas > 0 and coul-lambdas == 0. For instance I donā€™t get LINCs or ā€œunsettled water moleculeā€ warnings for vdw-lambdas=1 and coul-lambdas=0.1 but I do get these warnings for the following states:

vdw-lambdas             = 0.45 1.00
coul-lambdas            = 0.00 0.00

I also just tested from switching off both vdw and electrostatics vdw-q from couple-lambda0 to couple-lambda1 = none to only switching off vdw I donā€™t not get these warnings and can run dynamics. I would like to know why changing this matters in the sense of how GROMACs runs. Also should my couple-lambda0 and couple-lambda1 be changed to accommodate each lambda state? And if so how are they related and what should they be for:

vdw-lambdas             = 0.00 0.45 0.52 0.56 0.58 0.60 0.62 0.64 0.67 0.70 0.75 0.80 0.85 0.90 0.95 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
coul-lambdas            = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00

If my coul-lambdas == 0.00 for N different states should I only be turning off vdw interactions for these N states. And likewise if vdw-lambdas == 1 for M different states should I only be turning off coulomb interactions for these M states?