Force function for Dispersive PME

Developers discussions
1

Hi,

Does anyone know which scheme is used for determining the forces for dispersive PME (r^-6) non bonded interactions? As far as I know, for PME electrostatics, you can use ik scheme and spline differentiation schemes.

The ik scheme works well for dispersive but is time intensive as it requires 3+ 3d fast fourier transforms. The spline differentiation is faster, but I didn’t see any dispersive PME papers mention the use of spline differentiation for force calculation. In my numerical experiments as well, the spline differentiation for dispersive was WAY off from the ewald and ik force function values, which seemed to reconcile.

In summary, I have the following 2 questions:

  1. Which force calculation scheme does GROMACS use for non bonded electrostatic force calculations?

2. Which force calculation scheme does GROMACS use for non bonded dispersive PME calculations?

Thanks in advance!

2

GROMACS uses spline differentiation for all. The main reason is that we require a conservative potential / forces. But it is also much faster, in particular in parallel. There is nothing special about r^-6. Just that the functions have a slightly different form.

3

Got it. I had a few follow up questions / clarifications :

I’m trying to understand how the term in LJ (SR) and LJ (Recip) is calculated.

  1. Is LJ (SR) = V_dir + V_0, and LJ (recip) = V_recip from this link : Long Range Van der Waals interactions - GROMACS 2026.1 documentation?
  2. What is the significance of LJ_14 parameters and interactions? If PME is used, LJ 14 is not considered, and if vdw type = cutoff is used, only then does that energy get included. Is this correct?
  3. From what I understand, the c6ii parameters are determined from the atom “types”. the atom type is determined, and the corresponding index (0, 1, 2, …. ) is used in the lookup table for LJ_SR. Those c6 and c12 values are c6i and c12i, for which geometric rules can be applied (c12ij = sqrt(c12ii * c12jj)). Is this correct? And what is the unit of c12i and c6i parameters? Couldn’t find them specified in any documentation, or else I must not have found the right page for it.
4
  1. No, LJ (SR) = V_dr + V_0 - V_recip_for_excluded_pairs
  2. LJ_14 are listed interactions and are applied independently from the non-bonded interactions. Usually 1-4 atom pairs are excluded from the non-bonded pair interactions
  3. Yes. Simple dimension analysis tells you that units have to be energy/distance^(6 or 12), this is also listed in the topology parameter table. Note that the non-bonded machinery internally stores 6*c6 *and 12 **c12.
5

Got it. I was able to reproduce the exact LJ recip. values for Argon, which doesn’t have any exclusion lists, but not able to reproduce the exact LJ values for alanine dipeptide, which does have exclusion lists. Also I’m having a hard time reconciling the V_dir values from the formula you gave me.

  1. I was able to deduce from a 2 particle Argon system with distance between the particles > rC, that the negative of V0 is added to LJ (SR), i.e. as per my understanding, LJ (SR) = V_dir - V_0 - V_recip_for_excluded_pairs. Could you please check this, or let me know if I am wrong?
  2. Do exclusion lists play any part in calculating the reciprocal part from PME, or are they not considered?
  3. What happens to the C12 / r^12 terms in LJ-PME? Do we neglect them in V_dir, or are they included? I couldn’t find any papers that mention why they neglect that term, or if the dispersive function (r^-6) is a subset of LJ (r^-12 with r^-6), and that we have to add the r^-12 terms as well.

Any clarification would be greatly appreciated.

6
  1. Why the minus sign in front of V_0? Otherwise correct.
  2. No, reciprocal part and energy contains all interactions, including excluded and self.
  3. 1/r^12 decays so fast that there is never a need for long-range treatment.
7

I was able to calculate the LJ (SR) value for argon system of any number of particles as follows:

V_{rep-sr} = \frac{1}{2}\sum_{n_x}\sum_{n_y}\sum_{n_z*}\sum_{i,j}^{N} \frac{C^{(12)}_{ij}}{(r_{ijn})^{12}}

Where V_{rep-sr} is the repulsive part of the short term interactions.

For Argon, I got :

LJ (SR) = V_{rep-sr} + V_{dir} - V_0 - V_{recip-excluded-pairs}

I am attaching a tpr file here (argon_42_tpr). I did the following:

argon_42_tpr.dat (4.7 KB)

Run the simulation, extract energy:

#move file from dat to tpr
mv argon_42_tpr.dat argon_42.tpr
# run simulation
gmx_mpi mdrun -deffnm argon_42
gmx_mpi energy -f argon_42.edr -o lj_sr_argon_42.xvg
head -n 30 lj_sr_argon_42.xvg
# note the lj_sr value at timestep 0, that is the expected value for LJ (SR)

Use my python script to calculate LJ (SR) as per my understanding:

calculate_lj_sr_py.dat (9.5 KB)

# convert dat to py
mv calculate_lj_sr_py.dat calculate_lj_sr.py
# Port tpr file to human readable format for python
gmx_mpi dump -s argon_42.tpr > argon_42_verbose.tpr
# Run the python script with the verbose file
python calculate_lj_sr.py argon_42_verbose.tpr

This script only gave the correct LJ (SR) for configurations of argon system. For NaCl, and alanine dipeptide, it did not reconcile with energy values from gromacs.

So there’s something I am doing wrong.

8

I misunderstood how c6 values are interpreted from the verbose tpr file. I have fixed the script, and LJ (SR) calculated with the script now reconciles with argon as well as nacl tpr files. I am attaching the nacl tpr file here:

nvt_nacl_tpr.dat (16.3 KB)

Updated script:

calculate_lj_sr_py.dat (9.5 KB)

However, I am still unable to reconcile it with alanine dipeptide :

aladi_nve_tpr.dat (11.4 KB)

The difference between them is that V_recip_excluded_pairs != 0 for aladi.

  1. Is this formula correct for V_recip_excluded_pairs?

    V_{recip-excl-pairs} = -\frac{1}{2} \sum_{i,j, i != j}^N \frac{C_{ij}^{6}(1 - \phi(\beta . r_{ij}))}{r_{ij}^6} - \frac{1}{2} . \sum_{i}^N C_{ii}^{(6)}.\frac{\beta^6}{6} \\ \phi(x) = (1 + x^2 + \frac{1}{2} x^4) e^{-x^2}

  2. I’m assuming that V_dir also excludes terms of (i,j) which are part of the excluded pair. Is this correct?

9
  1. Looks correct.
  2. I don’t completely understand your question. You mean that term of 1. is included in V_dir? The answer to that is yes.

Note that things are more complicated when the parameters do not follow geometric combinations rules.

10

Let me clarify a little more on the question 2 above.

We calculate V_{dir} as follows:

V_{dir} = -\frac{1}{2}\sum_{i,j}^N \sum_{n_x}\sum_{n_y}\sum_{n_z*} \frac{C_{ij}^{(6)} \phi(\beta r_{ij, n})}{r_{ij,n}^6}

In the above summation, we are summing over all possible (i, j) pairs. When exclusion lists are not empty, do we exclude the pairs from the above sum that belong to the exclusion list?

Is this the correct formula then for V_dir? :

V_{dir} = -\frac{1}{2}\sum_{i,j, (i,j) \not\in EXCL}^N \sum_{n_x}\sum_{n_y}\sum_{n_z*} \frac{C_{ij}^{(6)} \phi(\beta r_{ij, n})}{r_{ij,n}^6}

Yes, for now I am only discussing about geometric combination rules.

I did read up about how things can get a little more complex with Lorentz Berthelot combination rules.

11

I got all of my answers with numerical experiments.

Thank you Dr. Hess.

I still think this is the formula for LJ (SR):

𝐿𝐽(𝑆𝑅)=𝑉𝑟𝑒𝑝−𝑠𝑟 + 𝑉𝑑𝑖𝑟 − 𝑉0 − 𝑉𝑟𝑒𝑐𝑖𝑝−𝑒𝑥𝑐𝑙𝑢𝑑𝑒𝑑−𝑝𝑎𝑖𝑟𝑠

as I mentioned above. All of my numerical experiments say that, but please let me know if you think I am wrong.