# Protein-ligand atom groups coulombic interaction energies discrepancies between PME and cut-off **URL:** \ **Category:** User discussions\ **Created:** [February 12, 2021, 7:20pm UTC](https://gromacs.bioexcel.eu/t/protein-ligand-atom-groups-coulombic-interaction-energies-discrepancies-between-pme-and-cut-off/1587 "2021-02-12T19:20:05Z")\ **Posts on this page:** 3\ **Page:** 1
**Author:** ![harold.grosjean](https://avatars.discourse-cdn.com/v4/letter/h/e9a140/32.png) [@harold.grosjean](https://gromacs.bioexcel.eu/u/harold.grosjean)\ **Post date:** [February 12, 2021, 7:20pm UTC](https://gromacs.bioexcel.eu/t/protein-ligand-atom-groups-coulombic-interaction-energies-discrepancies-between-pme-and-cut-off/1587/1 "2021-02-12T19:20:06Z")
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 ---
**Author:** ![milosz.wieczor](https://avatars.discourse-cdn.com/v4/letter/m/b19c9b/32.png) [@milosz.wieczor](https://gromacs.bioexcel.eu/u/milosz.wieczor)\ **Post date:** [February 14, 2021, 5:26pm UTC](https://gromacs.bioexcel.eu/t/protein-ligand-atom-groups-coulombic-interaction-energies-discrepancies-between-pme-and-cut-off/1587/2 "2021-02-14T17:26:14Z")
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 ---
**Author:** ![harold.grosjean](https://avatars.discourse-cdn.com/v4/letter/h/e9a140/32.png) [@harold.grosjean](https://gromacs.bioexcel.eu/u/harold.grosjean)\ **Post date:** [June 23, 2021, 12:50pm UTC](https://gromacs.bioexcel.eu/t/protein-ligand-atom-groups-coulombic-interaction-energies-discrepancies-between-pme-and-cut-off/1587/3 "2021-06-23T12:50:06Z")
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](https://drive.google.com/drive/folders/1DP5Ci3ABKnVCOYy1PJjMfON15bdkBNMz?usp=sharing)). 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