Coulomb (SR) calculation in PME

Developers discussions
,
1

Hi, I’m trying to understand the calculation of Coulomb (SR) in an alanine dipeptide molecule.

In my effort to do so, I read documentation and came across this statement ( Long Range Electrostatics - GROMACS 2026.1 documentation ):
”The term Coulomb (SR) contains the V_{dir} with the constant shift applied minus the reciprocal contribution for excluded atom pairs.”

When I ran it for the tpr file attached (as nve_pme.dat, just renamed from tpr because gromacs wouldn’t let me upload tpr files), (gmx mdrun -deffnm nve_pme), and took the energy dump of coulomb sr (gmx energy -f nve_pme.edr, select Coulomb SR), I got the following output:

”@ legend length 2
@ s0 legend “Coulomb (SR)”
0.000000 -340.237762”

I then took the human readable verbose version of the tpr, generated by (gmx dump -s nve_pme.tpr > nve_pme_verbose.tpr, attached as nve_pme_verbose.dat)

nve_pme_verbose.dat (43.5 KB)

, and passed it through the python script (attached) that calculates pairwise real part of the electrostatic energy, shifted with cutoff such that potential at cutoff is 0, and only includes pairs that have distance < rcutoff. It also calculates the total reciprocal energy for the excluded pairs, and subtracts it from the previously calculated sum. This final value is what I interpret to be Coulomb (SR), but the output of the script shows:
”Raw real-space sum : -73.27686428337606 kJ/mol
Exclusion correction : -190.07040972488696 kJ/mol
V_Coulomb_SR (total) : 116.7935454415109 kJ/mol”

Could someone please help in identifying which part of the calculation I did wrong? Thanks!

nve_pme.dat (11.4 KB)

calculate_coulomb_sr.dat (3.0 KB)

2

I think I have figured it out.

The documentation states that the self energy term (V0) is added to the term denoted by Coul. recip.

My experiments say that the self energy term is added to Coulomb (SR). I am currently working on verifying this with another test system of NaCl ions, but if someone could verify this, it would save me some time. Thanks!

You can use the same tpr files (nve_pme.tpr, uploaded as nve_pme.dat). Just uploading the updated python script to reproduce the exact (Coulomb (SR)). Also uploading the xvg file for Coulomb (SR) that I got.

calculate_coulomb_sr.dat (3.1 KB)

energy.xvg (1.0 KB)

3

Ah, the documentation is outdated then. Where do you read that the self term is added to Coul. recip.? The “Ewald summation” section says that it is added to SR.

4

Could you please share the link to the updated documentation?
Link to documentation I am referring to : Long Range Electrostatics - GROMACS 2026.1 documentation

For me, under “Ewald summation” section, it shows the following:
”The term Coulomb (SR) contains the V_{dir} with the constant shift applied minus the reciprocal contribution for excluded atom pairs. The term Coul. recip. contains the reciprocal sum $V_{rec}
$, including contributions of excluded atoms pairs, plus the charge correction $V_0
$ plus an optional surface correction term that depends on the system dipole.”

5

Ah, sorry. I was mixing things up while reading that part. I’ll fix it.

6

Thank you, I appreciate it.

I had a follow up question regarding exclusion lists and how they are handled.

For alanine dipeptide, I found the following exclusion list dump:

excls:
numLists=22
numElements=218
excls[0][num=7]={0, 1, 2, 3, 4, 5, 6}
excls[1][num=9]={0, 1, 2, 3, 4, 5, 6, 7, 8}
excls[2][num=7]={0, 1, 2, 3, 4, 5, 6}
excls[3][num=7]={0, 1, 2, 3, 4, 5, 6}
excls[4][num=12]={0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 14}
excls[5][num=9]={0, 1, 2, 3, 4, 5, 6, 7, 8}
excls[6][num=17]={0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,
14, 15, 16}
excls[7][num=9]={1, 4, 5, 6, 7, 8, 9, 10, 14}
excls[8][num=16]={1, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
16, 17, 18}
excls[9][num=12]={4, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16}
excls[10][num=12]={4, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16}
excls[11][num=8]={6, 8, 9, 10, 11, 12, 13, 14}
excls[12][num=8]={6, 8, 9, 10, 11, 12, 13, 14}
excls[13][num=8]={6, 8, 9, 10, 11, 12, 13, 14}
excls[14][num=17]={4, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,
17, 18, 19, 20, 21}
excls[15][num=9]={6, 8, 9, 10, 14, 15, 16, 17, 18}
excls[16][num=12]={6, 8, 9, 10, 14, 15, 16, 17, 18, 19, 20, 21}
excls[17][num=9]={8, 14, 15, 16, 17, 18, 19, 20, 21}
excls[18][num=9]={8, 14, 15, 16, 17, 18, 19, 20, 21}
excls[19][num=7]={14, 16, 17, 18, 19, 20, 21}
excls[20][num=7]={14, 16, 17, 18, 19, 20, 21}
excls[21][num=7]={14, 16, 17, 18, 19, 20, 21}

What I understood is that for each atom, it denotes the set of atoms it does NOT have coulombic interactions with. I used this set to calculate all values, and they tallied up with what gromacs told me.

For an NaCl system, with 256 Na+ ions and 256 Cl- ions, I got the following exclusion list:

moltype (0):
name=“NA”
atoms:
atom (1):
atom[ 0]={type= 0, typeB= 0, ptype= Atom, m= 2.29900e+01, q= 1.00000e+00, mB= 2.29900e+01, qB= 1.00000e+00, resind= 0, atomnumber= -1}
atom (1):
atom[0]={name=“Na”}
type (1):
type[0]={name=“Na”,nameB=“Na”}
residue (1):
residue[0]={name=“Na”, nr=1, ic=’ '}
excls:
numLists=1
numElements=1
excls[0][num=1]={0}

….

…

moltype (1):
name=“CL”
atoms:
atom (1):
atom[ 0]={type= 1, typeB= 1, ptype= Atom, m= 3.54530e+01, q=-1.00000e+00, mB= 3.54530e+01, qB=-1.00000e+00, resind= 0, atomnumber= -1}
atom (1):
atom[0]={name=“Cl”}
type (1):
type[0]={name=“Cl”,nameB=“Cl”}
residue (1):
residue[0]={name=“Cl”, nr=1, ic=’ '}
excls:
numLists=1
numElements=1
excls[0][num=1]={0}

Could you please explain how the exclusion lists work in this case? Does this mean Na+ ions don’t interact with each other coulombically? Similar for Cl-?

7

No. Exclusions are always within the same molecule. So each Na atom does not interact with itself.

8

Got it.

Just to confirm, atom 1 of type Na interacts with atom 2 of type Na coulombically, right? From what I understood, atom 1 of type Na does not interact with atom 1, i.e. itself, which is implicitly understood but explicitly stated in the file.

9

Exactly.

10

Great. Thank you so much for your help.