Unable to read GROMACS amber14sb topology due to missing dihedral terms

[lohedges]

We are unable to load GROMACS topologies generated using pdb2gmx with the amber14sb force field. I am not yet sure how general this is, but here are some example files. You'll need to move the amber14sb.ff directory (an adaptation of the user contribution found here) to the topology file associated with your GROMACS install. This can be found be parsing the output of running gmx, e.g.:

gmx
                 :-) GROMACS - gmx, 2022.3-plumed_2.9.0_dev (-:

Executable:   /usr/local/gromacs/bin/gmx
Data prefix:  /usr/local/gromacs
Working dir:  /home/lester/Downloads
Command line:
...

Take the "Data prefix" and append /share/gromacs/top.

To reproduce the issue. First generate the topology and coordinates:

gmx pdb2gmx -f molecule.pdb -ff amber14sb -water tip3p

Now try to load them, which gives lots of errors related to missing dihedral parameters:

import sire as sr
sr.load(["conf.gro", "topol.top"])
...
OSError: SireError::io_error: Cannot load the molecules: Could not construct a molecule system from the information stored in this Gromacs topology file. Errors include:
Error constructing the molecule associated with template 'Protein' : Cannot find the dihedral parameters for the dihedral between atoms GroAtom( name() = HG, number() = 1450 )-GroAtom( name() = CG, number() = 1449 )-GroAtom( name() = CD1, number() = 1451 )-GroAtom( name() = HD11, number() = 1452 ) (atom types HC-3C-CT-HC, function type 9).
...

Looking at the contents of, in my case, /usr/local/gromacs/share/gromacs/top/amber14sb.ff/ffbonded.itp, I do see terms for the dihedral, but with the ordering inverted, e.g.:

 HC  CT  3C  HC    9         0.0      0.62760     3

(This is one example, but I've checked that others are also present, although not all.)

Is it possible that our parser isn't taking to account reverse ordering of dihedral parameters? I remember that it was a nightmare to match terms correctly when writing the Sire.IO.CharmmPSF parser, particular as the ordering of dihedral terms is ill-defined and changes between topology formats.

Alternatively, perhaps this user contributed file uses non-standard terms, or has failed to apply the correct ordering when converting from AMBER to GROMACS.

For reference, if I re-parameterise using one of the built-in force fields, e.g. amber99sb, then everything works:

gmx pdb2gmx -f molecule.pdb -ff amber99sb -water tip3p

import sire as sr
sr.load(["conf.gro", "topol.top"])
System( name=Protein num_molecules=1 num_residues=164 num_atoms=2472 )

(This makes me think that it could be an issue with the amber14sb.ff files, rather than our parser.)

[chryswoods]

On a quick look, I believe we do correctly search for both the forwards and inverted dihedral. The function get_dihedral_id creates a unique dihedral ID from "A B C D" or "D C B A" by inverting if the first is less than the last...

But, as I write this, I realise that this doesn't account for the problem, as we have here, where we have "A B C A" and "A C B A". In that case, "A" == "A" so we don't invert them properly.

The fix is to check of "A" == "A" in get_dihedral_id and then invert if "B" < "C".

[lohedges]

Great, I'll give this a quick test.

[chryswoods]

I think this is the fixed code

static QString get_dihedral_id(const QString &atm0, const QString &atm1,
                               const QString &atm2, const QString &atm3,
                               int func_type)
{
    if ((atm0 < atm3) or (atm0 == atm3 and atm1 <= atm2))
    {
        return QString("%1;%2;%3;%4;%5").arg(atm0,atm1,atm2,atm3).arg(func_type);
    }
    else
    {
        return QString("%1;%2;%3;%4;%5").arg(atm3,atm2,atm1,atm0).arg(func_type);
    }
}

[lohedges]

Yes, that fixes it. Many thanks for the quick solution!