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
[ 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
[ 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
[ 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
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
[ 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
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.
[ 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.