DatasetBuilder() fails to correctly read in Urey-Bradley bonds from CHARMM generated topology

[arnosschnizler]

Hi, I am trying to create a grappa training dataset using the DatasetBuilder() using orca calculations and custom topology files generated from charmm. The custom topology files parameterise bonds according to Urey-Bradley, i.e. the forcefield.itp file reads:

For reference: "normal" topology file:

Using the DatasetBuilder(), I try to add the nonbonded parameters calling:

add_nonbonded_from_gmx_top(mol_id=mol_id, top_file=mol_topol)

Resulting in the error:

This error arises in from_openmm_system due to the UB-angles parameterised in the .itp files being interpreted as additional harmonic bonds when generating the OpenMM System object leading to a topology mismatch. A similar problem has been documented before.

How to best deal with it? Is this something that could be addressed in the DatasetBuilder()? Thank you for your help!

[arnosschnizler]

Just to add some more context: this is the workflow I am trying to establish: I prepare the topology files (from the MD simulations I have run to generate conformations I have labeled with orca) so they are grappa conform: prep_topol_toppar then I run grappa_ds_build to generate the dataset object combining the data from both qm-labeling and md-setup/runs.

Unfortunately these scripts are very setup specific (myself relying on an sql query and specific paths etc.) but maybe they make the problem more clear.

grappa_ds_build.py

prep_topol_toppar.py

[ehhartmann]

Hi,
Glad to hear that you are trying out grappa!
I've never worked with Urey-Bradley terms but the point here is that there is no explicit bond definition in the topology, right?
If I understood you correctly, you just want to get the nonbonded parameters from the topology file. For this, you could try to add the bond terms explicitly and change the angle type to harmonic so that grappa can read the top file. Alternatively, you could extract the nonbonded parameters and feed them to grappa as list/dict.
Of course, we could also implement proper parsing of UB terms in grappa but I don't know if @emerbitz has time for it/it is useful for his project, so the first options might give you quicker results;)
Best,
Eric

[emerbitz]

Hey:), does the molecule from your example consist of 52 or 106 bonds or neither? If your molecule has 52 bonds, we could maybe just skip the validation of the number of bonds (e.g., by adding a conditional). But not sure whether this won't lead to other downstream problems.

[arnosschnizler]

Thank you for the suggestions: Being lazy (and – as far as I could see – the mismatching bonds only being constructed and used within the validate_bonds scope) – I just modified grappa method inplace:

_orig_from_openmm_system = Molecule.from_openmm_system.__func__
@classmethod
def _patched_from_openmm_system(cls, *args, validate_bonds=False, **kwargs):
    return _orig_from_openmm_system(cls, *args, validate_bonds=validate_bonds, **kwargs)
Molecule.from_openmm_system = _patched_from_openmm_system

downstream topology comparisons of the qm-calculation derived molecular graph constructed with ASE and the graph constructed from the topology again agree :) thanks for your help!

[ehhartmann]

As Eddie mentioned before, this approach should lead to issues during training because the topology objects are missing bonds and the corresponding bond parameters won't be used for training. Could you confirm that your system has 106 bonds?
Do you have impropers in your system? If not, we should have some methods to regenerate the bonded terms, which could quickly fix your issue.

[arnosschnizler]

At the moment something else has come up we need to deal with before we continue here. Also for this particular system we could confirm we have 52 bonds and we use proper dihedrals (funct. 9).

[ehhartmann]

Okay, your approach should be fine, then!
I'll leave this issue open in case you run into any subsequent problems.