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!