Incomplete GROMACS top file when using openff with parmed or interchange for GAFF

[pablo-osmo]

I am trying to use OpenFF and either parmed or Interchange to produce gromacs input files with GAFF 2 parameters for small molecules (example below is cis-2-butene), however something is wrong with the bond parameterization using gaff-2.11 (I have also tried other version of gaff for the same problem). How would I get all the necessary bond parameters to run my simulations?

The final GROMACS .top file does not have bond coefficients for all the bond types, which causes the simulations to not work. I have attached an example of this incomplete TOP file from using parmed. I see the correct number of bonds in the original topology object.

[ bonds ]
; ai aj funct c0 c1 c2 c3
1 2 1 0.15100 213886.080000
2 3 1 0.13340 403170.240000
3 4 1 0.15100 213886.080000
1 5 1
1 6 1
[...]

I originally tried using Interchange, but the "Interchange.from_openmm" says it is experimental only and fails at "bond conversion" because it cannot find the correct bond information.

Failed to find parameters for bond with topology indices (0, 4)

Please let me know if any additional information is needed to help. I have attached a Jupyter notebook with the Parmed and interchange examples that are causing the errors above.

openff_incorrect_gromacs_topfile.zip

[mattwthompson]

This system was created with hydrogen-containing bonds constrained:

forcefield_kwargs={
    "constraints": HBonds,
    ...
},

so the bond force does not include information about the force constants of those bonds. There are only three bonds in that force, one for each heavy atom-heavy atom bond in the molecule:

import openmm

bond_force = [force for force in system.getForces() if isinstance(force, openmm.HarmonicBondForce)][0]

assert bond_force.getNumBonds() < molecule.n_bonds  # True; 3 < 11

OpenMM doesn't need force constants and lengths for constrained bonds and Amber has some way to work around it, I think. GROMACS needs those missing parameters, so converting through either ParmEd or Interchange is going to fail for the same reason. Interchange attempts to make this error descriptive; the bond between atoms 0 and 4 is a C-H bond for which it is not provided parameters.

Dropping the "constraints": HBonds, from the kwargs dict produces files that at least pass the eye check; I didn't see if they compile or not. If you only care about GROMACS simulations, I figure this should be fine as you'll just set the constraints in your MDP file.

[mattwthompson]

@pablo-osmo does dropping "constraints": HBonds fix this for you, or are there other issues you're running into?