I just wanted to confirm how Coulomb 1-4 energies are calculated. So I’m posting the pseudo code as I understand, please let me know if this is correct or not. (No PME or Ewald, only raw electrostatic force).
For Coulomb 1-4 energies, it considers all pairs in the [pairs] directive in the topology (topol.top) (In the verbose tpr file (output of gmx dump), this was stored under LJ-14 pairs. Is this expected behavior?).
For each pair, it calculates f q1 q2 / r, where f = 138.935 kJ/mol * nm / e^2
Once the calculation is done, it scales the sum of this value by (1/1.2) → I saw this in a research paper, but just wanted to confirm.
The Coulomb SR energies considers all atom pairs that are closer than rCutoff distance from each other, and excludes atom pairs that are not part of the 1-4 pairs list as discussed above. For each pair, it calculates the same (f q1 q2 / r), and then sums it up. The scaling factor (1/1.2) is not multiplied.
Thanks hess. Didn’t know about the factor already being present in the tpr file.
For Coulomb SR, I wasn’t able to find the formula in the manual. It mentions:
When Ewald summation or particle-mesh Ewald is used to calculate the long-range interactions, the short-range Coulomb potential must also be modified. Here the potential is switched to (nearly) zero at the cut-off, instead of the force. In this case the short range potential is given by:
It does talk about Ewald and PME, but not standard cutoff. When I say standard, here I’m just trying to understand the simple O(N^2) calculation (with cutoff) of the coulomb energies.
Could you perhaps guide me to the file in gromacs where this calculation occurs? I know that for bonded forces it happens in bonded.cpp under listed_forces module, but I wasn’t able to pinpoint where the calculation of LJ and SR energies and forces is taking place.
The code is very complex. Plain Coulomb, with the default potential-shift modifier, is qiqj/eps*(1/r - 1/rc). The shift of -1/rc is also applied for excluded pairs (which is an anomaly).
I think it’s under pairs.cpp? I’m willing to give this some time to truly understand it, so would appreciate any help in this regard.
Also, wanted to ask, in the calculation of forces for LJ and Coulomb (with / without PME), are the 1-4 pairs excluded completely? I see the 1-4 interaction energies in the energy files, but I was wondering if that is just information output for the user for their reference, or is the force for the neighbour pairs also calculated for MD?
pairs.cpp computes only the 1-4 interactions, or, more accurately, the pairs listed in the [ pairs ] section in the topology. These interactions are usually excluded from the non-bonded interactions. But strictly speaking the pairs in the [ pairs ] section doesn’t need to be a subset of the excluded pairs (even though I don’t see a reason to not do this in practice).
The non-bonded interactions are computed in many different places, depending on the back-end used.