Protein-ligand atom groups coulombic interaction energies discrepancies between PME and cut-off

User discussions
1

GROMACS version:2020
GROMACS modification: No

Hello,

I am trying to compute the interaction energies between multiple protein and ligand atom groups in the most sensible possible to obtain information about which part of the ligand interacts the most with the protein. I have about 70 protein-ligand atom groups all stored in index files.

I ran my initial trajectory (see bellow: mdp-prod bellow) then I perform a re-run on that trajectory to extract the interaction energies for the different atom groups.

I am comparing two rerun protocols at a cutoff of 2.5 nm for the calculation of coulombic interactions. The first uses PME (see bellow: mdp-rerun-PME), and the second uses cut-off (see bellow: mdp-rerun-cutoff) as coulomb type. My understanding is that if the atoms groups are within the cutoff values (and they are) then the computed average energy values should be almost identical between PME and cut-off. However, this is not the case. In the PME protocol, the coulombic average potentials values are closer in magnitude to the ones of the LJ. The cut-off protocol leads to more extreme coulombic average potentials values when compared to their PME equivalent. The calculated average LJ potentials are almost identical between the two protocols. See table below for illustrative example values (tab1).

Thus, I am wondering what causes these discrepancies between the two protocols? And also, what approach makes the most sense in term of getting sensible values?

Please see the files below my name.

Many thanks,

Harold

mdp-prod
;====================================================
; Production simulation
;====================================================

;----------------------------------------------------
; RUN CONTROL
;----------------------------------------------------
integrator = sd
nsteps = 15000000
dt = 0.002
comm-mode = Linear
nstcomm = 100

;----------------------------------------------------
; OUTPUT CONTROL
;----------------------------------------------------
nstxout = 0
nstvout = 0
nstfout = 0
nstxout-compressed = 1000
compressed-x-precision = 1000
nstlog = 1000
nstenergy = 1000
nstcalcenergy = 100

;----------------------------------------------------
; BONDS
;----------------------------------------------------
constraint_algorithm = lincs
constraints = h-bonds
lincs_iter = 1
lincs_order = 4
lincs-warnangle = 30
continuation = yes

; NEIGHBOR SEARCHING
;----------------------------------------------------
cutoff-scheme = Verlet
ns-type = grid
nstlist = 10
rlist = 1.0
pbc = xyz

;----------------------------------------------------
; ELECTROSTATICS
;----------------------------------------------------
coulombtype = PME
rcoulomb = 1.0
ewald_geometry = 3d
pme-order = 4
fourierspacing = 0.10
ewald-rtol = 1e-6

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

;----------------------------------------------------
; TEMPERATURE & PRESSURE COUPL
;----------------------------------------------------
tcoupl = V-rescale
tc_grps = System
tau_t = 1.0
ref_t = 300
pcoupl = Parrinello-Rahman
pcoupltype = isotropic
tau_p = 2
ref_p = 1.0
compressibility = 4.5e-05

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

mdp-rerun-PME
;====================================================
; RERUN FOR Energy GROUPS
;====================================================

;----------------------------------------------------
; RUN CONTROL
;----------------------------------------------------
integrator = sd ; stochastic leap-frog integrator
nsteps = 15000000 ; 2 * 15,000,000 fs = 30,000 ps = 30 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 = 5000 ; xtc compressed trajectory output every 5000 steps (10 ps)
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
energygrps = SG0LIG SG0PROT
;----------------------------------------------------
; BONDS
;----------------------------------------------------
constraint_algorithm = lincs ; holonomic constraints
constraints = h-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 = 2.5 ; short-range neighborlist cutoff (in nm)
pbc = xyz ; 3D PBC

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

mdp-rerun-cutoff
;====================================================
; RERUN FOR Energy GROUPS
;====================================================

;----------------------------------------------------
; RUN CONTROL
;----------------------------------------------------
integrator = sd ; stochastic leap-frog integrator
nsteps = 15000000 ; 2 * 15,000,000 fs = 30,000 ps = 30 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 = 5000 ; xtc compressed trajectory output every 5000 steps (10 ps)
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
energygrps = SG0LIG SG0PROT
;----------------------------------------------------
; BONDS
;----------------------------------------------------
constraint_algorithm = lincs ; holonomic constraints
constraints = h-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 = 2.5 ; short-range neighborlist cutoff (in nm)
pbc = xyz ; 3D PBC

;----------------------------------------------------
; ELECTROSTATICS
;----------------------------------------------------
coulombtype = Cut-off ; Cut-off
rcoulomb = 2.5 ; short-range electrostatic cutoff (in nm)
coulomb-modifier = Potential-shift
ewald_geometry = 3d ; Ewald sum is performed in all three dimensions
ewald-rtol = 1e-6 ; relative strength of the Ewald-shifted direct potential at rcoulomb

;----------------------------------------------------
; VDW
;----------------------------------------------------
vdwtype = Cut-off
rvdw = 2.5
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

tab1
Coul_average_PME_2.5 Coul_average_cutoff_2.5
ID
SG0 3.726400 8.767360
SG1 -9.243990 44.363600
SG2 -13.026000 -33.384900
SG3 -13.784100 -36.331200
SG4 -0.438439 3.125040
SG5 -2.585790 4.864530
SG6 6.151790 12.554200
SG7 -2.251410 18.794800
SG8 -0.897619 -1.094030
SG9 17.354100 69.411700

2

PME calculates electrostatics assuming an infinite grid of periodic images, which is fine for dynamics but not quite realistic for inter-group interaction analysis, unless your system is a crystal. If you just want the interaction between two closest images, use the cut-off approach.

Also remember to always set nstlist to 1 for reruns (with large strides you cannot trust the neighbor list to stay similar between steps), although I believe Gromacs does that automatically now.

Best,
Miłosz

3

Hello,

many thanks for your response. I have been pushing this approach recently and I am still hitting problems when looking at the Coulombics. Calculation of non bonded interactions between certain atom groups (protein vs ligand) still yields abnormally large values.

Initially, I thought that this was due to a charge imbalance between the 2 groups. After looking more into this, the total charge does not seem to explain those extreme values. The idea was that if the two groups hold charges of the same sign then the sum of those charges (Qtot) would be different than 0 which could lead in large repulsions (or vice versa). But this does not seem to be the case (see figures here: Coulombic_interactions_gmx - Google Drive).

Would you have a clue about what could cause those large values? And if so, could you suggest an approach that normalizes the extreme values obtained or, better, makes the values similar to the LJ profile which is all in the same order of magnitude.

Please let me know if you require any further information.

Many thanks in advance for your help,

Harold

fig. 1 - Calculated non-bonded coulombic interactions between protein and ligand sub-atom groups for cutoff and PME schemes. The data points are colored by the total charge of the two groups (Qlig + Qprot = Qtot)

fig. 2 - Relationship between the total charge of the two groups (Qlig + Qprot = Qtot) and calculated non-bonded coulombic interactions for Cutoff scheme

fig. 3 - Relationship between the total charge of the two groups (Qlig + Qprot = Qtot) and calculated non-bonded coulombic interactions for PME scheme

Here is the relevant bit of my .mdp file

;----------------------------------------------------
; NEIGHBOR SEARCHING
;----------------------------------------------------
cutoff-scheme = Verlet
ns-type = grid ; search neighboring grid cells
nstlist = 1 ; 1 fs (set to 1 for re-runs)
rlist = 3.5 ; short-range neighborlist cutoff (in nm)
pbc = xyz ; 3D PBC
;----------------------------------------------------
; ELECTROSTATICS
;----------------------------------------------------
rcoulomb = 3.5 ; short-range electrostatic cutoff (in nm)
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

;Cutoff parameters
coulombtype = Cut-off
coulomb-modifier = Potential-shift

;PME parameters
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme-order = 4 ; interpolation order for PME (default is 4)
fourierspacing = 0.10 ; grid spacing for FFT