search

choderalab / espaloma

Extensible Surrogate Potential of Ab initio Learned and Optimized by Message-passing Algorithm 🍹https://arxiv.org/abs/2010.01196
https://docs.espaloma.org/en/latest/
MIT License
202 stars 23 forks source link

Espaloma cannot use constrained (default) OpenFF force fields #189

Open mattwthompson opened 9 months ago

mattwthompson commented 9 months ago 
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.