# Role of \[ pairs \], \[nonbonded \], and combination rules **URL:** \ **Category:** User discussions\ **Tags:** forcefield\ **Created:** [August 16, 2024, 9:49am UTC](https://gromacs.bioexcel.eu/t/role-of-pairs-nonbonded-and-combination-rules/9899 "2024-08-16T09:49:00Z")\ **Posts on this page:** 3\ **Page:** 1
**Author:** ![obZehn](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/obzehn/32/2629_2.png) [@obZehn](https://gromacs.bioexcel.eu/u/obZehn)\ **Post date:** [August 16, 2024, 9:49am UTC](https://gromacs.bioexcel.eu/t/role-of-pairs-nonbonded-and-combination-rules/9899/1 "2024-08-16T09:49:01Z")
GROMACS version: 2023 & 2024 GROMACS modification: No Dear all, I am trying to port a small force field (~ ten atom types and 5 molecules) from the original implementation (DL\_POLY & LAMMPS) to GROMACS. Now, apart from a ton of technical limitations due to the difference in the engines, I am facing a few key problems in the definition of the force fields in GROMACS, specifically how non-bonded interactions are computed. I’ll try to split my question with numbers, as I (sadly) have a few of them. A couple of key points are a) the chemical species are small (carbonate ion, water, bicarbonate) and should not need any `[pair]` entry, as in the original implementation the neighbors are excluded and the intra-molecular forces are all due to `[bonds]`, `[dihedrals]`, and `[angles]`. As such, I have no `[pair]` entry and I do not want the force field to generate any pairs at all, and in the topology of the molecules I exclude all the nearby atoms by properly setting `nrexcl`. My understanding is that, formally, pair interactions are part of the “bonded” term of the force field (even if they are passed to non-bonded kernels), as these are calculated for atoms that are within the same molecule, not for all the atoms within a certain cut-off as for Coulomb and van der Waals (those would be the inter-molecular non-bonded terms). I can’t have `[pairs]` interactions between atoms that do not belong to the same molecule, those will be defined purely by non-bonded terms (coulomb + LJ). b) the Lennard-Jones interaction of the atoms in the box is defined for many of the possible pairs of atoms, but for a few of them they are omitted, which means (for the ff developer) that they are zero, i.e., some of them have only a pure Coulombic interaction as non-bonded term. Now, in the definition of the force field, I start with something like this ```auto [defaults] ; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ 1 2 no 0.000000 1.000000 ``` where I am setting the interaction as Lennard-Jones (`nbfunc = 1`) with `comb-rule = 2`, that is, the columns in the `[atomtypes]` and in the `[non-bonded]` section will report sigma and epsilon of a LJ 12-6 potential in nm and kJ/mol, respectively. Here I start getting a bit confused. I set `gen-pairs = no`, I do not want any pair within molecules and I want GROMACS to complain if it need to generate pairs, as in principle I do not have any. From the manual, I see that > Note that gen-pairs, fudgeLJ, fudgeQQ, and are optional. fudgeLJ is only used when generate pairs is set to ‘yes’, and fudgeQQ is always used. So, my first question is i) I set `gen-pairs = no` and I have an empty pair list. Whatever is the value of `fudgeLJ`, there won’t be any (1-4) pair calculated. However, `fudgeQQ` is always used. Does this refer to the pair generation again? That is, whenever I calculate a pair interaction, if it is calculated with a combination rule, then the LJ part will be calculated with the sigma and epsilon of the `[atomtypes]` and will be rescaled by `fudgeLJ`, while the coulomb part will be calculated with charges and rescaled by `fudgeQQ`. On the other hand, if it is tabulated in `[pairs]`, then the LJ part will be calculated with the parameters reported in the `[pairs]` section and NOT rescaled by `fudgeLJ` while the coulomb part again with the charges and also rescaled by `fudgeQQ`? And, as such, since I do not want and I do not have any pairs, are both `fudgeLJ` and `fudgeQQ` irrelevant in my case? Then, I have the `[atomtypes]` section. Here, I have something like this ```auto [atomtypes] ; name at.num mass charge ptype sigma epsilon H 1 1.0100 0.410000 A 0.000000 0.000000 O 8 16.0000 -0.820000 A 0.000000 0.000000 ``` where I set the the LJ sigma and epsilon to zero. Then I have the `[nonbond_params]` that looks something like this ```auto [nonbond_params] ; i j func sigma (nm) epsilon (kJ/mol) O O 1 5.98839454115e-01 3.15558573222e-04 ``` This takes me to the second question, which is ii) I want for some non-bonded interaction to be purely coulombic. In the example here, I have two atom types, so I have three possible combinations of non-bonded pairs, that is, O-O, H-H, and O-H. However, I define the interactions only for O-O. My understanding is that the non-bonded term for O-O will be coulomb (due to charges) and LJ (from the tabulated sigma and epsilon in `[nonbond_params]`). Here, `fudgeLJ` and `fudgeQQ` are used NOWHERE, they are irrelevant. This leaves me with two undefined interaction, namely H-H and O-H. The coulomb part is clearly defined (I have charges). How is the LJ contribution computed? My understanding is that the corresponding sigma and epsilon will be taken from `[atomtypes]`, a process similar to pair generation but without `fudgeLJ` and `fudgeQQ` rescaling. Since I set them all to zero, I expect that for the couple O-O there is a LJ term coming from `[nonbond_params]`, while for H-H and O-H the LJ term is calculated from `[atomtypes]` and, as such, will be zero due to the values being zero in `[atomtypes]`. This will implicitly leave me with pure coulomb interaction between O-H and H-H. In all these calculations, both `fudgeLJ` and `fudgeQQ` are irrelevant. Am I mistaken? Related to the question before, when I test such a simple system, `grompp` goes smooth and does not complain, which is okay. In the output I see something like ```auto Generated 2 of the 3 non-bonded parameter combinations ``` which I guess are H-H and O-H. Instead, if I add explicitly the non-bonded parameters as in the following ```auto [nonbond_params] ; i j func sigma (nm) epsilon (kJ/mol) O O 1 5.98839454115e-01 3.15558573222e-04 H H 1 0.0 0.0 H O 1 0.0 0.0 ``` I get ```auto Generated 0 of the 3 non-bonded parameter combinations ``` My understanding is that, since I have all the interactions listed in the `[nonbond_params]` section, there is no need for any generation of non-bonded combinations. iii) Supposing I am right, my question then is, are the two implementations identical, that is, explicitly setting to zero the interaction in `[nonbond_params]` is as setting to zero the parameters in `[atomtypes]`? Because, since I have no pairs, in my case would be easier to just set sigma=epsilon=0 for all atoms in `[atomtypes]` rather than generating all the couples in `[nonbond_params]`. Since I have no pairs and I do not generate them, I shouldn’t have any problem. Am I mistaken? Are the (sigma,espilon) values in `[atomtypes]` used for something else and I am messing up in some way with the interactions? And finally, I have a more general question. Here I report the topology for a TIP3 water molecule from CHARMM-GUI. ```auto [moleculetype] ; name nrexcl TIP3 2 [atoms] ; nr type resnr residu atom cgnr charge mass 1 OT 1 TIP3 OH2 1 -0.834000 15.9994 ; qtot -0.834 2 HT 1 TIP3 H1 2 0.417000 1.0080 ; qtot -0.417 3 HT 1 TIP3 H2 3 0.417000 1.0080 ; qtot 0.000 [settles] ; OW funct doh dhh 1 1 9.572000e-02 1.513900e-01 [exclusions] 1 2 3 2 1 3 3 1 2 ``` iv) Are `nrexcl = 2` and the `[exclusions]` redundant? Because excluding the two next bonded atoms is the same as listing all the combinations below `[exclusions]`. I understand that `[exclusions]` gives me more freedom, as I can selectively kill non-bonded terms while `nrexcl` applies everywhere in the molecule, but here within this simple example they seem to be doing the same thing. Am I wrong and the `[exclusions]` is also used somewhere else? Sorry for the long post! I hope that I was exhaustive enough. ---
**Author:** ![milosz.wieczor](https://avatars.discourse-cdn.com/v4/letter/m/b19c9b/32.png) [@milosz.wieczor](https://gromacs.bioexcel.eu/u/milosz.wieczor)\ **Post date:** [August 16, 2024, 10:34am UTC](https://gromacs.bioexcel.eu/t/role-of-pairs-nonbonded-and-combination-rules/9899/2 "2024-08-16T10:34:29Z")
Hi @obZehn, Your explanations match my intuitions about the features you discuss. I understand that (i) fudgeXX is only applied to interactions listed in `[pairs]`, (ii) setting epsilon to 0 leaves a coulomb-only interaction (with `[nonbond_params]` always taking precedence over `[atomtypes]`), (iii) sigma/epsilon do not have other uses than what you listed, and (iv) the two definitions of exclusions should be equivalent (although I don’t quite know what the technical reasons are for always listing them explicitly in the `.itp`). Still, one reasonable thing to do would be to test your assumptions using this topology comparison function (that uses Gromacs’ `gmx energy` under the hood) so that you can verify that what you assume is equivalent is indeed seen as such by Gromacs: > **[Checking whether topologies yield identical energies · Wiki · KBM /...](https://gitlab.com/KomBioMol/gromologist/-/wikis/Checking-whether-topologies-yield-identical-energies)** > > Package to handle GROMACS topology stuff Just remember to choose a compatible structure where any potential differences will show up. ---
**Author:** ![obZehn](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/obzehn/32/2629_2.png) [@obZehn](https://gromacs.bioexcel.eu/u/obZehn)\ **Post date:** [August 18, 2024, 12:26pm UTC](https://gromacs.bioexcel.eu/t/role-of-pairs-nonbonded-and-combination-rules/9899/3 "2024-08-18T12:26:12Z")
Dear @milosz.wieczor Thanks for the insight, and for the tool! :)