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

Atom elements which `espaloma` has been trained on? #206

Closed JSLJ23 closed 4 months ago

JSLJ23 commented 5 months ago

Dear espaloma devs,

I was trying to use espaloma to parameterise a few slightly more complex systems from PDBbind, specifically for the ligands and one system with a Heme prosthetic group is giving issues.

import torch
import espaloma as esp
from openff.toolkit.topology import Molecule
from rdkit import Chem

ligand_full_path = "./PDBbind_v2020_sdf/3csl_ligand.sdf"
rdkit_mol = next(Chem.SDMolSupplier(ligand_full_path, removeHs=False))
ligand_off_mol = Molecule.from_rdkit(rdkit_mol, allow_undefined_stereo=True)

molecule_graph = esp.Graph(ligand_off_mol)
espaloma_model = esp.get_model("latest")
espaloma_model(molecule_graph.heterograph)
openmm_system = esp.graphs.deploy.openmm_system_from_graph(molecule_graph)

This raises UnassignedBondError: BondHandler was not able to find parameters for the following valence terms:

---------------------------------------------------------------------------
UnassignedBondError                       Traceback (most recent call last)
Cell In[6], line 1
----> 1 openmm_system = esp.graphs.deploy.openmm_system_from_graph(molecule_graph)

File ~/mambaforge/envs/nnff_py310/lib/python3.10/site-packages/espaloma/graphs/deploy.py:116, in openmm_system_from_graph(g, forcefield, suffix, charge_method, create_system_kwargs)
    110 elif charge_method == "nn":
    111     g.mol.partial_charges = unit.elementary_charge * g.nodes["n1"].data[
    112         "q"
    113     ].flatten().detach().cpu().numpy().astype(
    114         np.float64,
    115     )
--> 116     sys = ff.create_openmm_system(
    117         g.mol.to_topology(),
    118         charge_from_molecules=[g.mol],
    119         allow_nonintegral_charges=True,
    120     )
    122 else:
    123     # create openmm system
    124     raise RuntimeError(
    125         "Charge method %s is not supported. " % charge_method
    126     )

File ~/mambaforge/envs/nnff_py310/lib/python3.10/site-packages/openff/utilities/utilities.py:80, in requires_package.<locals>.inner_decorator.<locals>.wrapper(*args, **kwargs)
     77 except Exception as e:
     78     raise e
---> 80 return function(*args, **kwargs)

File ~/mambaforge/envs/nnff_py310/lib/python3.10/site-packages/openff/toolkit/typing/engines/smirnoff/forcefield.py:1164, in ForceField.create_openmm_system(self, topology, toolkit_registry, charge_from_molecules, partial_bond_orders_from_molecules, allow_nonintegral_charges)
   1126 @requires_package("openmm")
   1127 def create_openmm_system(
   1128     self,
   (...)
   1134     allow_nonintegral_charges: bool = False,
   1135 ) -> "openmm.System":
   1136     """Create an OpenMM System from this ForceField and a Topology.
   1137 
   1138     Note that most force fields specify their own partial charges, and any
   (...)
   1162         Allow charges that do not sum to an integer.
   1163     """
-> 1164     return self.create_interchange(
   1165         topology,
   1166         toolkit_registry,
   1167         charge_from_molecules=charge_from_molecules,
   1168         partial_bond_orders_from_molecules=partial_bond_orders_from_molecules,
   1169         allow_nonintegral_charges=allow_nonintegral_charges,
   1170     ).to_openmm(combine_nonbonded_forces=True)

File ~/mambaforge/envs/nnff_py310/lib/python3.10/site-packages/openff/utilities/utilities.py:80, in requires_package.<locals>.inner_decorator.<locals>.wrapper(*args, **kwargs)
     77 except Exception as e:
     78     raise e
---> 80 return function(*args, **kwargs)

File ~/mambaforge/envs/nnff_py310/lib/python3.10/site-packages/openff/toolkit/typing/engines/smirnoff/forcefield.py:1223, in ForceField.create_interchange(self, topology, toolkit_registry, charge_from_molecules, partial_bond_orders_from_molecules, allow_nonintegral_charges)
   1220     used_registry = GLOBAL_TOOLKIT_REGISTRY
   1222 with toolkit_registry_manager(used_registry):
-> 1223     return Interchange.from_smirnoff(
   1224         force_field=self,
   1225         topology=topology,
   1226         charge_from_molecules=charge_from_molecules,
   1227         partial_bond_orders_from_molecules=partial_bond_orders_from_molecules,
   1228         allow_nonintegral_charges=allow_nonintegral_charges,
   1229     )

File ~/mambaforge/envs/nnff_py310/lib/python3.10/site-packages/openff/interchange/components/interchange.py:270, in Interchange.from_smirnoff(cls, force_field, topology, box, positions, charge_from_molecules, partial_bond_orders_from_molecules, allow_nonintegral_charges)
    221 """
    222 Create a new object by parameterizing a topology with a SMIRNOFF force field.
    223 
   (...)
    266 
    267 """
    268 from openff.interchange.smirnoff._create import _create_interchange
--> 270 return _create_interchange(
    271     force_field=force_field,
    272     topology=topology,
    273     box=box,
    274     positions=positions,
    275     charge_from_molecules=charge_from_molecules,
    276     partial_bond_orders_from_molecules=partial_bond_orders_from_molecules,
    277     allow_nonintegral_charges=allow_nonintegral_charges,
    278 )

File ~/mambaforge/envs/nnff_py310/lib/python3.10/site-packages/openff/interchange/smirnoff/_create.py:105, in _create_interchange(force_field, topology, box, positions, charge_from_molecules, partial_bond_orders_from_molecules, allow_nonintegral_charges)
    101 interchange.box = _topology.box_vectors if box is None else box
    103 _plugins(interchange, force_field, _topology)
--> 105 _bonds(interchange, force_field, _topology, partial_bond_orders_from_molecules)
    106 _constraints(
    107     interchange,
    108     force_field,
    109     _topology,
    110     bonds=interchange.collections.get("Bonds", None),  # type: ignore[arg-type]
    111 )
    112 _angles(interchange, force_field, _topology)

File ~/mambaforge/envs/nnff_py310/lib/python3.10/site-packages/openff/interchange/smirnoff/_create.py:149, in _bonds(interchange, force_field, _topology, partial_bond_orders_from_molecules)
    143     from openff.interchange.smirnoff._valence import _upconvert_bondhandler
    145     _upconvert_bondhandler(force_field["Bonds"])
    147 interchange.collections.update(
    148     {
--> 149         "Bonds": SMIRNOFFBondCollection.create(
    150             parameter_handler=force_field["Bonds"],
    151             topology=_topology,
    152             partial_bond_orders_from_molecules=partial_bond_orders_from_molecules,
    153         ),
    154     },
    155 )

File ~/mambaforge/envs/nnff_py310/lib/python3.10/site-packages/openff/interchange/smirnoff/_valence.py:277, in SMIRNOFFBondCollection.create(cls, parameter_handler, topology, partial_bond_orders_from_molecules)
    272         molecule.generate_conformers(n_conformers=1)
    273         molecule.assign_fractional_bond_orders(
    274             bond_order_model=handler.fractional_bond_order_method.lower(),
    275         )
--> 277 handler.store_matches(parameter_handler=parameter_handler, topology=topology)
    278 handler.store_potentials(parameter_handler=parameter_handler)
    280 return handler

File ~/mambaforge/envs/nnff_py310/lib/python3.10/site-packages/openff/interchange/smirnoff/_valence.py:169, in SMIRNOFFBondCollection.store_matches(self, parameter_handler, topology)
    165     self.key_map[topology_key] = potential_key
    167 valence_terms = self.valence_terms(topology)
--> 169 _check_all_valence_terms_assigned(
    170     handler=parameter_handler,
    171     topology=topology,
    172     assigned_terms=matches,
    173     valence_terms=valence_terms,
    174 )

File ~/mambaforge/envs/nnff_py310/lib/python3.10/site-packages/openff/interchange/smirnoff/_base.py:180, in _check_all_valence_terms_assigned(handler, assigned_terms, topology, valence_terms)
    178 exception.unassigned_topology_atom_tuples = unassigned_atom_tuples
    179 exception.handler_class = handler.__class__
--> 180 raise exception

UnassignedBondError: BondHandler was not able to find parameters for the following valence terms:

- Topology indices (40, 42): names and elements ( N), ( Fe), 
- Topology indices (41, 42): names and elements ( N), ( Fe), 
- Topology indices (39, 42): names and elements ( N), ( Fe), 
- Topology indices (38, 42): names and elements ( N), ( Fe), 

I am wondering if espaloma is having issues with the Fe2+ atom, and if so, may I know if there's a way to go about fixing this temporarily (i.e. like manually assigning the bonded term?). Also, is there documentation to point out which atom elements espaloma supports?

mattwthompson commented 4 months ago

Hey @JSLJ23 - I'm not on the Espaloma team but I can provide a couple of breadcrumbs

I am wondering if espaloma is having issues with the Fe2+ atom

Precisely - that's what the error message is trying to communicate out. Espaloma's non-bonded parameters come from Sage, which do not yet have parameters for coordinated metal ions.

https://github.com/choderalab/espaloma/blob/13af90b8158d50730c2c6888dacce8351021f797/espaloma/graphs/deploy.py#L41-L57

Additionally, probably because OpenFF doesn't have Fe2+ to begin with, Espaloma must not have been trained with this N-Fe bond. OpenFF (SMIRNOFF) force fields can be extended (here's an example) but I don't know how to extend Espaloma parameters. In this case you'd need to come up with both non-bonded and bonded parameters. Presumably if the SMIRNOFF force field has bonded parameters and Espaloma doesn't, it just won't override them ... I haven't tested this myself.

I did indulge myself to some code golf to see if I could get this from the OpenFF Toolkit's API:


In [1]: from openff.toolkit import ForceField

In [2]: from openff.units.elements import SYMBOLS

In [3]: from rdkit import Chem

In [4]: {
   ...:     SYMBOLS[Chem.MolFromSmarts(parameter.smirks).GetAtomWithIdx(0).G
   ...: etAtomicNum()]
   ...:     for parameter in ForceField("openff-2.0.0.offxml")["vdW"].parame
   ...: ters
   ...: }
Out[4]:
{'Br',
 'C',
 'Cl',
 'Cs',
 'F',
 'H',
 'I',
 'K',
 'Li',
 'N',
 'Na',
 'O',
 'P',
 'Rb',
 'S'}

Note that a couple of these elements are just aqueous ions pulled externally so that TIP3P can work in some solutions (don't recall the details).

JSLJ23 commented 4 months ago

This is really helpful, thank you. 👍