Absolute binding free energy calculation

[jjmllen]

GROMACS version:
GROMACS modification: Yes/No
Here post your question

Dear all,

I hope this message finds you well.

I am currently attempting to calculate Absolute Binding Free Energy (ABFE) using GROMACS, following the tutorial provided at the link below:

http://www.alchemistry.org/wiki/Absolute_Binding_Free_Energy_-_Gromacs_2016

However, I have encountered an issue where some of the expected dhdl files are not being generated, regardless of whether I use my system or the system described in the tutorial. I would greatly appreciate any guidance from anyone who has successfully run ABFE calculations in GROMACS and could assist me in resolving this.

Thank you for your time and assistance.

Kind regards,
Justice

[MagnusL]

When you say “some of the expected dhdl files are not being generated”, do you mean that only the dhdl files are missing? Did those simulations start at all (or was there a problem with the script)? Did they finish properly?

[jjmllen]

these are the dhdl files which are missing and below here is the part in slurm.out file when the lamba.xx is generating the files during the nvt step. and i have attached as well the whole slurm.out file

[MagnusL]

It seems like some simulations have crashed, i.e., they did not finish properly. So, you are not only missing the dhdl files.

I would suggest running the crashed simulations one at a time trying to understand what is going wrong. Are there any warnings from `gmx grompp´ already?

[jjmllen]

thank you for your response. I have tried to run the nvt.07. step separately( which is one of the simulation which was not completed by the script).

and I am getting below error. I have tried to increase the rlist from 1.0 to 1.5
Fatal error:
There are 1008 perturbed, excluded non-bonded pair interactions beyond the
pair-list cut-off, 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.

For more information and tips for troubleshooting, please check the GROMACS
website at Common Errors — GROMACS webpage https://www.gromacs.org documentation

and this is the .mdp file
NVT equilibration
;====================================================

;----------------------------------------------------
; RUN CONTROL
;----------------------------------------------------
define = -DPOSRES
integrator = sd ; stochastic leap-frog integrator
nsteps = 50000 ; 2 * 5,000 fs = 10 ps
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 = 5000 ; xtc compressed trajectory output every 5000 steps
compressed-x-precision = 1000 ; precision with which to write to the compressed trajectory file
nstlog = 5000 ; update log file every 10 ps
nstenergy = 5000 ; save energies every 10 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 = no ; 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.5 ; 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 = no

;----------------------------------------------------
; VELOCITY GENERATION
;----------------------------------------------------
gen_vel = yes ; Velocity generation is on (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 = 7
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
~

[MagnusL]

Is there a reason why you are decoupling/mutating bonded interactions at all? What is the b-state of those bonded interactions? Why not just leave them?
Does it work if you run with:

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

?

[jjmllen]

it has run for other lamda windows but in some cases the system breaks down.