Manual Calculation of Coulomb (SR) is inconsistent with Gromacs

User discussions
1

GROMACS version: 2024.1
GROMACS modification: No

I have been trying to manually calculate the interaction energy of a test system, a simple 5x5x5 nm box of TIP3 water in the CHARMM36 forcefield.

I am using PME with a 1.2nm cutoff and for a single frame, gromacs calculates a coulomb SR potential of around -177,000.

I have taken this frame using MDAnalysis to get the atomic coordinates and the charges and attempted to calculate the interaction energy using equation 168 in the 2024.1 documentation. As I understand this, the entire coulomb SR interaction should be a sum over atoms of charge products (qij) times coulomb constant times the complementary error function of the inter-atomic separation times the ewald splitting constant divided by the interatomic distance.

When I apply this formula to the system with the same cutoff and ewald splitting parameter, my energy is inconsistent with gromacs, I get -153,000. Am I missing any extra contribution towards the SR coulomb?

I want to use this to calculate interaction energies for use in a widom particle insertion analysis and the standard shifted coulomb potential up to now does not seem to produce the correct shape excess chemical potential profile for membranes, hence the need to try and reproduce PME.

2

I can think of two causes of the difference:

  1. GROMACS shifts the pair interaction by a constant such that the potential is zero at the cut-off
  2. Pairs up to 1-4 interactions are excluded.

Are you using the GROMACS TPI integrator for your test particle insertion?

3

The gromacs TPI integrator as far as I can see does not allow us to calculate how the excess chemical potential changes across the box so I am writing my own as a post-processing tool.

4

You can insert around a given location using TPIC.

It wouldn’t be much work to add a scan along z.

5

We have a cavity locating script already written and outputting cavity coordinates for each trajectory frame, do you have any guidance on adding our pre-defined cavity coordinates to the trajectory as required by TPIC?

6

I did this by spitting out pdb files using trjconv and pasting the extra line for the cavity location. That’s a bit cumbersome though. If you want to do a scan along z, the simplest way is to add a few line to tpi.cpp to read an environment variable that sets the z-location.

7

Thank you, I’ll take a look at modifying the .cpp file.

In the meantime, just to clarify, it is possible to do this by dumping the pdb or .gro files for each frame, adding the extra coordinate for TPIC, generate the tpr file using the tpic integrator and run this using gmx mdrun -v -deffnm tpic -rerun Frame_n.gro

This should allow us to perform test particle insertion on a particular pre-defined cavity, in a particular frame so we can calculate across the box?

8

Correct.

9

Here is the modification for a scan along z for GROMACS 2025. You can modify zRange to your liking. I guess that the location should be at least zRange from the box edge. I haven’t tried the code though.

diff --git a/src/gromacs/mdrun/tpi.cpp b/src/gromacs/mdrun/tpi.cpp
index 70af0a4fed..cb85efab7d 100644
--- a/src/gromacs/mdrun/tpi.cpp
+++ b/src/gromacs/mdrun/tpi.cpp
@@ -728,6 +728,15 @@ double TestParticleInsertion::insertIntoFrame(const double           t,
                                               gmx_wallcycle*         wallCycleCounters,
                                               t_nrnb*                nrnb)
 {
+    bool limitZ    = false;
+    real zRange    = 0.25;
+    real zLocation = 0;
+    if (char* result = std::getenv("GMX_TPI_Z"))
+    {
+        limitZ = true;
+        zLocation =  std::strtod(result, nullptr);
+    }
+
     /* Copy the coordinates from the input trajectory */
     auto x = makeArrayRef(stateGlobal->x);
     for (Index i = 0; i < rerunX.ssize(); i++)
@@ -769,7 +778,14 @@ double TestParticleInsertion::insertIntoFrame(const double           t,
                 /* Generate a random position in the box */
                 for (int d = 0; d < DIM; d++)
                 {
-                    xInit[d] = dist_(rng_) * box[d][d];
+                    if (d < ZZ || !limitZ)
+                    {
+                        xInit[d] = dist_(rng_) * box[d][d];
+                    }
+                    else
+                    {
+                        xInit[d] = zLocation + (dist_(rng_) * 2 - 1) * zRange;
+                    }
                 }
             }
         }
10

I also added the shift of the Ewald direct potential to the manual.

11

That is brilliant, thank you for your time!

If I could trouble you for one more thing, when gromacs is calculating the short range coulomb interaction using PME, is there any other contribution to the short range potential other than the direct space sum of equation 168 in the gromacs manual. I could get my plain and shifted coulomb potential python code to return consistent results with gromacs but the SR coulomb using PME varied greatly from that returned by gromacs.

12

There shouldn’t be. The long-range correction terms are added to the LR energy.

13

It’s quite strange, I have made sure I am using the same gaussian width (beta) parameter as gromacs and cutoff however the potential is consistently different.

14

I think I know the source of the difference. GROMACS includes the mesh energy of excluded pairs in the mesh energy. This energy is subtracted in the short range term.

15

I’ll take a look and see if this brings the energies in line, thank you!