Inclusion of electrostatic 1-4 interactions when exclusions nR are set at 3

User discussions
1

GROMACS version: 2022
GROMACS modification: Yes
Hi, I have a question about the inclusion of electrostatic charge-to-charge interactions in a topology. We are building our own force field and are using a custom modification to gromacs that allows for LJ 6-10 rather than 6-12. Our original intention was to use the force that fully includes 1-4 interactions, with nR set to 2 so that 1-4 interactions are fully included, but I’ve found that 6-10 was a little too strong in the LJ 1-4 interactions. I decided to use a fudge factor on the LJ 1-4, so I switched to nR = 3 and included a fudge factor that is directly in the C6 and C10 constants I’ve placed with the [ pairs ] directive. My question is this: if I have gen-pairs set to “no” in the force field defaults, do the electrostatic 1-4 pairs get included? From some initial results, it looks like something is missing… it would make sense if the electrostatic 1-4 interactions have been missed.

Related question is this: if I set gen-pairs as “yes” in the force field defaults, do the custom constants included in the [ pairs ] topology get overridden? Or, is there a direct way to set electrostatic 1-4 pairs in the topology if setting nR as 3 omits the electrostatic pairs when gen-pairs is set as “no”?

Sorry, confusing question! thanks ahead for your help!

Greg Smith

2

I hope someone can help here: I think I have insight into why I’m having this problem. It looks like Gromacs must compute 1-4 interactions separately from the nbnxm kernels. The fudge factors do not seem to appear in the nbnxm code or the atomdata.cpp code, but maybe I’m missing something. My modification switching from 6-12 to 6-10 anticipated no use of fudge factors or the [ pairs ] directive; we used the 1-4 interactions at full strength.

Can someone please confirm. Does Gromacs calculate LJ 1-4 interactions as a bonded interaction? Where do I look in order to find the code which regulates this calculation? It might correct my discrepancy if I switch the 12 → 10 in these calculations. Please help!

Thanks!
Greg Smith

3

Quick update. Hope I’m heading in the right direction.

I’ve found some further evidence that 1-4 interactions are treated in listed_forces. Wishing I was more capable at programming than I am. Tracking through this, I’ve found a function called do_pairs under the pairs.h header. It looks like the functions in pairs.cpp are doing the necessary calculation. I need to work through this code on paper to make sure I understand the calculation, but it seems that my instinct is correct. The modified version of gromacs that we’re using has no included accommodations for 1-4 pairs calculated as LJ 6-10… the energies are naturally off because I’m giving it 6-10 constants and it’s calculating 6-12 values!

Greg

4

This is correct - you might get better feedback in the developer section of this forum.

Some of my notes on listed-forces that might be useful:

there are multiple types of listed_forces (bond, angle, dihedral, …) that involves different number of partner (2,3,4) and have different parameters and number of parameters.

A sketch of the call tree:
MD loop
do_force
ListedForces::calculate()
calc_listed()
…
do_pairs()
do_pairs_simple()

It is helpful to have an IDE that can track call to method across the codebase (e.g. CLion Find usages) to explore this easily, including finding where member of struct are used.

During mdrun, the structure containing the listed force data is InteractionDefinitions as contained in the ListedForces object.

Which particles interact is specified in an array of InteractionList, often abbreviated ilist or il in the code. This is a std::array<InteractionList*, F_NRE> with a pointer for each interaction type (e.g. F_IDIHS), with some empty as they are not computed using listed_force. The other important member of InteractionDefinitions is std::vector<t_iparams>& iparams.

The InteractionList is a linear list, with InteractionList.size() element. Say we have an interaction with 3 atoms i, j and k, then the list will be [index_param_1, i_1, j_1, k_1, index_param_2, i_2, j_2, k_2, …], where i_1 is the index of the atom in the local atom arrays (its position is found in x[i_1], same for velocity and force arrays), identically for j and k. The index_param_1 is the index of an entry in the iparams array. t_iparams, the element of the iparams array, is a polymorphic struct that contain the parameters for the interaction (an union), for instance it can be an harmonic if the interaction is an angle. There are parameters for the A and B state for free energy purpose.

Note that “normal” 1-4 are covered under the lj14 type (F_LJ14), to the best of my knowledge, not ljc14. The charges used are loaded from the general charge array using the atom index.

Depedending on how robust/clean you want your modification to be, you can either modify the way the 1-4 interaction is computed, which will result in a gmx that is only usable for your specific use, or try to introduce another interaction type (F_LJ14_612 ?), such that both LJ 1-4 type are supported (and similarly for the NB NxM code). While this is ultimately cleaner, this require extensive changes that might not make it the first things to do.

1 Like
5

Thank you very much! This helps enormously.

I think in the short term, I will modify for my group’s specific purpose, though I have some desire in the long term to do the “cleaner”, more thorough, more generally useful approach. At the moment at least, time is not on my side;-) I’ve been shying from calling myself a developer because my programming skill is very pidgin, but I’ll ask this sort of question in the developer section next time.

Again, thanks!
Greg Smith
University of Colorado