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