Espaloma cannot use constrained (default) OpenFF force fields

[mattwthompson]

import espaloma
from openff.toolkit.topology import Molecule

molecule = Molecule.from_smiles("CCO")
graph = espaloma.Graph(molecule)

espaloma_model = espaloma.get_model("latest")
espaloma_model(graph.heterograph)

espaloma.graphs.deploy.openmm_system_from_graph(
    graph,
    forcefield="openff-2.1.0",
)

AssertionError                            Traceback (most recent call last)
Cell In[1], line 10
      7 espaloma_model = espaloma.get_model("latest")
      8 espaloma_model(graph.heterograph)
---> 10 espaloma.graphs.deploy.openmm_system_from_graph(
     11     graph,
     12     forcefield="openff-2.1.0",
     13 )

File ~/mambaforge/envs/espaloma/lib/python3.11/site-packages/espaloma/graphs/deploy.py:131, in openmm_system_from_graph(g, forcefield, suffix, charge_method, create_system_kwargs)
    129 name = force.__class__.__name__
    130 if "HarmonicBondForce" in name:
--> 131     assert force.getNumBonds() * 2 == g.heterograph.number_of_nodes(
    132         "n2"
    133     )
    135     for idx in range(force.getNumBonds()):
    136         idx0, idx1, eq, k = force.getBondParameters(idx)

AssertionError:

This bubbles up because the number of bonds in the openmm.System is not guaranteed to match the number of bonds in the Molecule. Interchange currently does not add bonds to constrained atom pairs; it may or may not have in the past, I don't remember. In general this behavior shouldn't be considered stable since, to my knowledge, it doesn't affect how OpenMM simulations are run.

In [2]: from openff.toolkit import ForceField

In [3]: system = ForceField("openff-2.1.0.offxml").create_openmm_system(molecule.to_topology())

In [4]: system.getForce(3).getNumBonds(), molecule.n_bonds
Out[4]: (2, 8)

In the current stack, the number of bonds in the two objects will only match for toy molecules with no hydrogens.