Coulomb force calculation

GROMACS version: 2021
GROMACS modification: Yes


I want to calculate a Coulomb force contribution from my simulations for every atom in simulation to use it later in visualization. What is the easiest way to extract the force components?

I started to do it myself by applying Coulomb law in coordinate processing and then comparing the electrostatic forces with total forces which I calculated from an xvg dump file. As a result, I got some values that I decided to double-check. After all, simplifying the system step-by-step, I ended up with a simple synthetic example system with just one H and one Cl ions with turned off LJ interactions. I hoped that now I would get a 100% match so I can make sure that my calculations are correct. But it appears to be that they still mismatch quite significantly:
H ion: total force (by gromacs): 9.4910180592178843 kJ/mol/nm vs coulomb (my utility) 11.3263893059850318 kJ/mol/nm
Cl ion: total force (by gromacs): 9.4894741077024918 kJ/mol/nm vs coulomb (my utility) 11.3263893059850318 kJ/mol/nm

What is strange to me is that:

  1. Why are the force values for H ion and Cl ion are still different (9.491 vs 9.489 kJ/mol/nm)?
  2. Why is there such a discrepancy between my and gromacs values (11.326 vs 9.491 kJ/mol/nm for NPT and ~11 vs 13 kJ/mol/nm for NVT)? I understand that I apply a slightly different approach (PME vs Coulomb law), but why is there such a difference? For more complex examples the difference can be up to 50%. Are there different forces there except electrostatics? Does the thermostat/barostat/anything else modify the forces?

Here is the code on how I calculate the interactions by Coulomb law (C++):

double k_e_si = 8.9875517923e9; // Coulomb constant; N * m^2 / C^2
double N_a = 6.0221409e23;
double e_to_coulomb = 1.602176565e-19; // e to C
// force in SI -> gromacs: N / mol = J / (mol * m) = 10^-3 kJ / (mol * 10^9 nm) = 10^-12 kJ / (mol * nm)
const double N_per_m_to_kJ_per_nm = 1e-12;
for (int i = 0; i < in->N_atoms; ++i) {
    f_coul[i][0] = 0;
    f_coul[i][1] = 0;
    f_coul[i][2] = 0;

    for (int j = 0; j < in->N_atoms; ++j) {
        if (i != j) {
            // position delta vector components between particles i and j
            double delta_x = pbc(coords[i][0] - coords[j][0], in->box_x); // Angstrom
            double delta_y = pbc(coords[i][1] - coords[j][1], in->box_y); // Angstrom
            double delta_z = pbc(coords[i][2] - coords[j][2], in->box_z); // Angstrom
            // square distance between particles i and j in Angstrom and meters 
            double delta_r_sq_A2 = delta_x * delta_x + delta_y * delta_y + delta_z * delta_z; // Angstrom^2
            double delta_r_sq_m2 = 1e-20 * delta_r_sq_A2; // m^2

            // force absolute value between particles i and j
            double force_ij_abs = N_per_m_to_kJ_per_nm * N_a * e_to_coulomb * e_to_coulomb
                                  * k_e_si * charges[i] * charges[j] / delta_r_sq_m2;

            // force vector components from given delta r vector and the force modulus
            double delta_r = sqrt(delta_r_sq_A2);
            double fx = force_ij_abs * delta_x / delta_r;
            double fy = force_ij_abs * delta_y / delta_r;
            double fz = force_ij_abs * delta_z / delta_r;

            // add to the total vector
            f_coul[i][0] += fx;
            f_coul[i][1] += fy;
            f_coul[i][2] += fz;

The param file (a slightly modified file from a gromacs tutorial on lysozyme in water simulation; I also tried NVT with a similar result):

title = OPLS Lysozyme NPT equilibration
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 5000 ; 2 * 500000000 = 1000000 ps = 1000 ns = 1us
dt = 0.001 ; 2 fs
; Output control
nstxout = 1000
nstvout = 1000
nstfout = 1000
nstenergy = 1000
nstlog = 1000
nstxout-compressed = 1000 ; save compressed coordinates
compressed-x-grps = System ; save the whole system
; Bond parameters
continuation = yes ; Restarting after NPT
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds involving H are constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Neighborsearching
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet scheme
rcoulomb = 4.99 ; short-range electrostatic cutoff (in nm)
rvdw = 4.99 ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
;coulombtype = Cut-off
; Temperature coupling is on
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = Water ; two coupling groups - more accurate
tau_t = 0.1 ; time constant, in ps
ref_t = 300 ; reference temperature, one for each group, in K
; Pressure coupling is on
pcoupl = Parrinello-Rahman ; Pressure coupling on in NPT
pcoupltype = isotropic ; uniform scaling of box vectors
tau_p = 2.0 ; time constant, in ps
ref_p = 1.0 ; reference pressure, in bar
compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction
DispCorr = EnerPres ; account for cut-off vdW scheme
; Velocity generation
gen_vel = no ; Velocity generation is off

Here is how I get the xvg file ( file):

echo 0 | gmx trjconv -s "${OUTPREFIX}.gro" -f "${OUTPREFIX}.xtc" -o "${OUTPREFIX}_movie.gro" -dt 1
echo 0 | gmx traj -f "${OUTPREFIX}.trr" -s "${OUTPREFIX}.tpr" -of "${OUTPREFIX}_force_dump.xvg" -dt 1
(then I call my utility)

The rest of the files, including the xyz files (can be opened in text editor and Ovito), are also attached. (749 Bytes)
npt.mdp (2.5 KB)
The whole archive with xyz files

1 Like


the bigger difference is most likely due to PME versus your calculation using only Coulomb’s law. The PME calculation accounts for interactions with the infinite periodic images of all charges. That’s vastly different than just calculating within the simulation box.

I believe the smaller difference between the two atoms is due to floating point precision rounding when adding the forces.


1 Like

Hi Petter,

Yes, you are right, thank you! Once I switched my calculations to “coulombtype = Cut-off” and turned on cut-off in my tool as well, the values started to match with high precision. Before that, I was also trying to use that scheme but got zero force but only now I realized it was due to the fact the diagonal is sqrt(3) times longer than box size so the zero force value was indeed possible:)

Also, I found the difference between ensembles: in the NVT and NPT configs, the cut-off values were different and the cut-off value is used in the PME mesh. I’m sorry for the confusion.

All in all, it looks like my calculations were correct up to the precision of the PME method in comparison to Coulomb law. I will also try to continue my production simulations with a small duration with a modified param file where electrostatics calculated in the expensive way (Coulomb) so I can make a fair analysis of the contribution of electrostatics to the total force.

Thank you so much, Petter!