Error met in free energy calculation

User discussions
,
1

GROMACS version:2022.5
GROMACS modification: No
I am trying to do free energy calculation on a protein-ligand system. I met such error when I doing the ligand part:

There are perturbed non-bonded pair interactions beyond the pair-list cutoff
of 1.26317 nm, which is not supported. This can happen because the system is
unstable or because intra-molecular interactions at long distances are
excluded. If the latter is the case, you can try to increase nstlist or rlist
to avoid this.The error is likely triggered by the use of couple-intramol=no
and the maximal distance in the decoupled molecule exceeding rlist

Here is my mdp file for the particular lambda window:
;====================================================
; Production simulation
;====================================================

;----------------------------------------------------
; RUN CONTROL
;----------------------------------------------------
integrator = sd ; stochastic leap-frog integrator
nsteps = 5000000 ; 5000 ps
dt = 0.001 ; 1 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 = none ; 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 = 90 ; maximum angle that a bond can rotate before LINCS will complain (30 is default)
continuation = yes

;----------------------------------------------------
; 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 (if gen_vel is ‘yes’, continuation should be ‘no’)
gen_seed = -1 ; Use random seed
gen_temp = 300

;----------------------------------------------------
; FREE ENERGY CALCULATIONS
;----------------------------------------------------
free-energy = yes
couple-moltype = ligand
couple-lambda0 = vdw-q
couple-lambda1 = none
couple-intramol = yes
separate-dhdl-file = yes
sc-alpha = 0.5
sc-power = 1
sc-sigma = 0.3
init-lambda-state = 12
coul-lambdas = 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.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
dhdl-print-energy = potential
dhdl-derivatives = yes

I removed the LINCS since it always cause another error. One thing that is really confusing for me is that when I turn on the LINCS, the mdrun for lambda from 0 to 11 and 19 are all okay. The LINCS and such cut-off problem only occur in lambda 12 to 18. I don’t know how to solve such problem. Hope to get some suggestions. Thanks!

2

You need to have constraints on h-bonds. Integration with all bonds flexible is unstable with dt=0.001. And the force field is likely parameterized with constraints on h-bonds.

What error do you get with constraints?

3

Thanks for replying.

The problem I met when I turn on LINCS with ds=0.002 is “bond that rotates more than 30/60/90 degree”. I build this ligand by using LigParGen and this model performed well in a typical simulation. So I think the ligand model I build is still okay. That’s the reason I tried turn off the LINCS.

4

You should never turn of constraints (except for problematic energy minimization).

What do you mean with “performed well in a typical simulation”? Does the same topology run without LINCS warnings in a non free-energy simulation?

5

The same topology run well in non free-energy simulation, no matter it is with or without LINCS. And, of course, I adopt the simulation with LINCS in that case.

This cutoff problem can also appear in free-energy simulation with LINCS and ds=0.002, based on my latest simulation task.

6

Some ligand force fields require constraints on all bonds. But I don’t know if that is the case for LigParGen.

7

Some updates:

For an identical system, I successfully ran the entire decoupling process in the 2018.8 version of gmx, also using totally the same mdp settings. But I met the same problem with this: How to solve Lincs warning in free energy calculation in 2023.1 version. (particularly from init-lamda-state 11 to init-lambda 18, simulation can run well in the rest of lambda states.)

I successfully ran the example of gmx tutorial (the methane-water Free Energy Calculations) in gromacs version 2023.1. So I guess it’s the topology of the ligand. I tried to use LigParGen (OPLS/AA) and charmm-gui (CHARMM36) to build the topology, but they both encountered the above problem (How to solve Lincs warning in free energy calculation) in such calculation.