Grompp handling bonded potentials

GROMACS version: 2023.4
GROMACS modification: Yes

Dear GROMACS users,

is an exactly double defined dihedral (e.g., A B C D 9 0 0.1234 3) counted twice during energy minimization/integration, or not? I did some single point energy calculations with different configurations of the dihedral and then compared the energy values with gmx check. For example, the energies were equal between A B C D 9 0 1000 3 and A B C D 9 0 700 3, A B C D 9 0 300 3, but not between A B C D 9 0 1000 3 and A B C D 9 0 500 3, A B C D 9 0 500 3. So am I right in my assumption that for exact copies of a dihedral only one copy is taken into account for further calculations?

Thanks in advance for any answers!

Best,
Marius

Do you mean in the [ dihedrals ] section or [ dihedraltypes ] section? In the [ dihedrals ] section every interaction is computed. So if you list a dihedral twice you get double the potential.

I am talking about the [ dihedraltypes ] section defined in a .prm file that is added to my .top file.
E.g.,

[ dihedraltypes ]
A B C D 9 0 500 3
A B C D 9 0 500 3

Dihedral type 9 is the only interaction type where multiple parameter entries generate multiple dihedrals. So you will get two dihedrals, so double the potential.

I see now that the reference manual is not very clear on this point.

I agree with the multiple parameter entries for the proper dihedral type. However, why do the energies differ then for two dihedral potentials with force constants of 500 kJ/mol each, compared to a single dihedral potential with 1000 kJ/mol? But are equal between two dihedral potentials with 700 kJ/mol and 300 kJ/mol compared to 1000 kJ/mol?

The manual states the following:

If parameters for a certain interaction are defined multiple times for the same combination of atom types the last definition is used; starting with GROMACS version 3.1.3 grompp generates a warning for parameter redefinitions with different values

Based on that I was thinking that grompp maybe ignores exact copies of a dihedral, and taking only the last definition.

Two of 500 should equal one of 1000.

As I said, the manual is not clear on this. I just uploaded a correction/explanation change for the next release.

However, they do not appear to be equal in GROMACS. When comparing the single-point energies of two identical structures—one with a .tpr file prepared using a single dihedral potential of 1000 kJ/mol and the other with two dihedral potentials of 500 kJ/mol each—I obtain different energy values.

Additionally, I compared two .tpr files using gmx check: one generated with a single dihedral potential of 500 kJ/mol and the other with two dihedral potentials, each with a force constant of 500 kJ/mol. The output was as follows:

gmx check -s1 1x_500/min.tpr -s2 2x_500/min.tpr

Note: When comparing run input files, default tolerances are reduced.
Reading file 1x_500/min.tpr, VERSION 2023.4 (single precision)
Reading file 2x_500/min.tpr, VERSION 2023.4 (single precision)
comparing inputrec
inputrec->ld_seed (-142628066 - -1141310215)
comparing mtop topology
comparing force field parameters
comparing cmap 0
comparing cmap 1
comparing cmap 2
comparing cmap 3
comparing cmap 4
comparing cmap 5
comparing cmap 6
comparing cmap 7
comparing cmap 8
comparing cmap 9
comparing cmap 10
comparing cmap 11
comparing cmap 12
comparing cmap 13
comparing cmap 14
comparing cmap 15
comparing cmap 16
comparing cmap 17
comparing cmap 18
comparing cmap 19
comparing cmap 20
comparing cmap 21
comparing cmap 22
comparing cmap 23
comparing molecule types
comparing atoms
comparing t_resinfo
comparing InteractionLists
comparing blocka excls[0]
comparing molecule blocks
comparing InteractionLists
comparing groups
comparing intermolecular exclusions
comparing moleculeBlockIndices
comparing flags
comparing box
comparing box_rel
comparing boxv
comparing x

If I understand the output correctly, the .tpr files are practical identical, or am I wrong?

Yes, they are identical.

Could you attach both top files, and both itp files if you used those?

I created a toy model using the CGenFF web server since the original data has not yet been published. However, the behavior remains the same. The molecule is butane, and I modified the dihedral between the first and last carbon atoms of the chain. I have attached the files with two dihedrals (each with a force constant of 500 kJ/mol) as well as the forcefield.itp file used for the toy model. (Note: GROMACS threw an error because the dihedral in the *.prm file was already defined in the ffbonded.itp file.)

Thank you for looking into this issue!

but_ini.pdb.txt (942 Bytes)
but.itp.txt (5.5 KB)
but.prm.txt (1.5 KB)
but.top.txt (190 Bytes)
single_point.mdp (387 Bytes)
forcefield.itp.txt (1005 Bytes)

Those two lines are, obviously, exactly identical. Does grompp warn about this? What happens when you use 600 and 400 instead?

No, grompp does not issue a warning about this (see 2x_500_grompp_output.txt ). I have also attached the output of gmx dump for the cases of 600 kJ/mol + 400 kJ/mol and 500 kJ/mol + 500 kJ/mol. In the latter case, only one PDHIS entry appears, with a force constant of 500 kJ/mol. Additionally, using gmx check , I compared the energies between the 600 kJ/mol + 400 kJ/mol and the 500 kJ/mol + 500 kJ/mol (or 1000 kJ/mol) cases. As you can see, the tool only returns differences between the 600 kJ/mol + 400 kJ/mol and the 500 kJ/mol + 500 kJ/mol cases.

2x_500_grompp_output.txt (968 Bytes)
but.prm.txt (1.5 KB)
1x_600_1x_400_tpr_dump.txt (518.6 KB)
2x_500_tpr_dump.txt (518.4 KB)
1x_600_1x_400_compare_1x_1000.txt (597 Bytes)
1x_600_1x_400_compare_2x_500.txt (691 Bytes)

You can see in the 2x500 dump that there is only one dihedral defined. I assume this must mean that grompp simply ignores a second, identical dihedraltype 9 parameter definition. We should have a warning for that.

Okay, so the PDHIS entries refer to any dihedral potential and not only to unique ones? A warning in such cases would certainly be helpful! :)

Yes, but only for dihedral type 9.

I think grompp should warn and keep both parameter entries (and thus generate two dihedrals).

Generating both dihedrals is likely a valid approach, as it aligns with the force field files (i.e., itp or prm). However, if a copy-paste error occurs—as it did in my case—a warning would be extremely useful to catch such mistakes.

Thank you very much for your help! :)