[gmx-users] frozen ligand for free energy calculations

Ahmet Yildirim ahmedo047 at gmail.com
Thu Jul 13 15:59:31 CEST 2017


Dear users,

I come across with an issue when I try to do free energy calculations. The
issue is about the roto-translational motions of the ligand in the
decoupled state. I mean the ligand doesn't stay stable in the binding
pocket as in the coupled state.
It seems that the restraints that are applied on the atoms of the protein
and ligand under [ intermolecular_interactions ] part (an example of it is
below) in the compex top file aren't sufficient to keep the ligand from
repositioning/rotation with respect to the protein in the decoupled state.
Even one, two and three sets of restraints couldn't solve the issue.

[ intermolecular_interactions ]
[ bonds ]
; ai      aj    type   bA           kA       bB        kB
  629     3     6      0.597        0.0      0.597     4184.0

[ angles ]
; ai      aj    ak    type   thA          fcA      thB        fcB
  281   629     3     1       37.5        0.0       37.5      41.84
  629     3    21     1      121.5        0.0      121.5      41.84

[ dihedrals ]
; ai     aj     ak    al    type    thA         fcA       thB       fcB
  249   281   629     3     2     -147.4        0.0     -147.4      41.84
  281   629     3    21     2      -60.5        0.0      -60.5      41.84
  629     3    21    16     2     -153.9        0.0     -153.9      41.84

The mdp file used for production run with GROMACS2016.2 is here.

; RUN CONTROL
;----------------------------------------------------
integrator   = sd            ; leap-frog integrator
nsteps       = 20000000
dt           = 0.001         ; 2 fs
comm-mode    = Linear        ; remove center of mass translation
nstcomm      = 100           ; frequency for center of mass motion removal
comm-grps    = System

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

; BONDS
;----------------------------------------------------
constraint_algorithm   = lincs    ; holonomic constraints
constraints            = h-bonds  ; hydrogens only are constrained
lincs_iter             = 2        ; accuracy of LINCS (1 is default)
lincs_order            = 12        ; 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  ; 'group' is default in 4.6 - from GMX5 on
'Verlet' is default
ns-type          = grid   ; search neighboring grid cells
nstlist          = 10     ; 20 fs (default is 10)
pbc              = xyz    ; 3D PBC

; ELECTROSTATICS & EWALD
;----------------------------------------------------
coulombtype      = PME      ; Particle Mesh Ewald for long-range
electrostatics
rcoulomb         = 1.2      ; 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.12     ; grid spacing for FFT
ewald-rtol       = 1e-6     ; relative strength of the Ewald-shifted direct
potential at rcoulomb

; VAN DER WAALS
;----------------------------------------------------
vdwtype         = Cut-off
vdw_modifier    = Potential-switch
rvdw            = 1.1          ; short-range van der Waals cutoff (in nm)
rvdw-switch     = 1.0          ; where to start switching the LJ force
DispCorr        = EnerPres     ; apply long range dispersion corrections
for Energy and Pressure
; TEMPERATURE COUPLING (SD ==> Langevin dynamics)
;----------------------------------------------------
tc_grps    =  Protein_ligand Water_and_ions
tau_t      =  1.0 1.0
ref_t      =  300 300

; PRESSURE COUPLING
;----------------------------------------------------
pcoupl           = no
pcoupltype       = isotropic          ; uniform scaling of box vectors
tau_p            = 1.0                ; 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 on (if gen_vel is 'yes',
continuation should be 'no')
gen_seed     = -1       ; Use random seed
gen_temp     = 300

;----------------------------------------------------
; FREE ENERGY
free-energy              = yes
init-lambda              = 1
delta-lambda             = 0
sc-alpha                 = 0.3
sc-power                 = 1
sc-sigma                 = 0.25
sc-coul                  = yes
couple-moltype           = ligand
couple-intramol          = no
couple-lambda0           = vdw-q
couple-lambda1           = none
nstdhdl                  = 100

I would try to freeze the ligand in the decoupled state in the canonical
ensemle with the above restrains under [ intermolecular_interactions ] but
I am not sure whether that makes sense or not? Justin says (
https://mailman-1.sys.kth.se/pipermail/gromacs.org_gmx-users/2013-August/083647.html):

"...Anything that is frozen, by definition, never has its position
updated.  Under the influence of
pressure coupling, other particles around the frozen group can have their
positions scaled and thus collide with the frozen group, which has remained
in
its original location". I think I can use the frozen ligand in both coupled
and decoupled state? And I should take into consideration the effect of the
frozen ligand on the free energy calculation, right?

Any comment is appreciated.

-- 
Ahmet Yildirim


More information about the gromacs.org_gmx-users mailing list