# [gmx-users] PME constant shifts.

Sergio Perez sperezconesa at gmail.com
Fri Mar 29 12:34:10 CET 2019

```Dear gmx comunity,

I am trying to calculate the electrostatic interaction of my system for
force field development. My system consists of 7 positively charged
particles interacting (system A) with another positively charged particle
(system B). I calculate the electrostatic interaction by doing reruns of
the trajectory in which some of them have been removed.

E_int(electrostatic) =
E_AB(coul-SR)+E_AB(coul-LR)-E_A(coul-SR)-E_A(coul-LR)-E_B(coul-SR)-E_B(coul-LR)

And to my surprise I get negative energies. I know the PME method has its
energy shifted:

"With the Verlet cut-off scheme, the PME direct space potential is shifted
by a constant such that
the potential is zero at the cut-off. This shift is small and since the net
system charge is close to
zero, the total shift is very small, unlike in the case of the
Lennard-Jones potential where all shifts
add up. We apply the shift anyhow, such that the potential is the exact
integral of the force."

Since for the system with just one charge I get a possitive E_1q(coul-LR),
I imagine the LR-coul term is shifted (eventhough I have to use the group
cut-off scheme and not the verlet). The problem is why don't all these
shifts cancel? Is there a way to not add the shifts?

I leave bellow the .mdp parameters. Note that for the VDW part I need to
use user-tables and therefore I need to use the group scheme.

Thank you very much!
Sergio Pérez-Conesa

; Integrator stuff
integrator = md-vv       ; steep or md
dt = 0.001
nsteps = 300000000 ; -1 no maximum
init-step = 30305000 ;
ld-seed = -2018962629
continuation = yes
emtol       = 10.0
emstep      = 0.01

; neighbour search
cutoff-scheme = group
nst-list = 10
verlet-buffer-tolerance = 0.005
ns-type = grid
; Pbc
pbc = xyz ; xyz or xy

; Vdw
rvdw = 1.4
rlist = 1.4
vdwtype = user  ; PME or User to look for table /user
DispCorr = No   ; corrections to energy and pressure or No

;Electrostatic
coulombtype = PME   ; User if look for table
rcoulomb = 1.4
pme-order =  4
fourierspacing = 0.2
ewald-geometry = 3d ; 3d or 3dc

;IMD-group                =

; T and P
tcoupl =no  ; none, nose-hoover, v-rescale
nh-chain-length = 10
nsttcouple = 10 ; 10 if not used
pcoupl = no ;
pcoupltype = semiisotropic ;
compressibility     = 4.5e-5 4.5e-5      ; isothermal compressibility of
water, bar^-1
refcoord-scaling = com
tau-p = 5.0  ; ps
ref-p = 1.0 1.0 ; bar

; Constraints
constraints = none      ; constraints  distances
constraint-algorithm = LINCS ; or SHAKE
lincs_iter                  = 1
lincs_order                 = 4

;Generate velocities
gen_vel         = no            ; assign velocities from Maxwell
distribution
gen_temp        = 300           ; temperature for Maxwell distribution
gen_seed        = 180611                ; generate a random seed

; OUTPUT
; Output control
nstenergy       = 1             ; save energies
nstlog          = 10 ; update log file every
;nstxout-compressed  = 500      ; save compressed coordinates every 10.0 ps
;nstxout = 20
;nstvout = 20
;nstfout = 20

; COM removal . There are more options like removing from groups or the
frequency
;comm-grps = System
comm-mode = none    ; linear, angular (removes linear and angular motion)
or none
nstcomm = 0  ; frequency of removal of com motion
freezegrps =
energygrp-excl =
;freezedim = Y Y Y
nstcalcenergy =
```