[WIP] Use PDBInf

environment.yml

@@ -10,6 +10,7 @@ dependencies:
   - openff-units ==0.2.0  # https://github.com/openforcefield/openff-units/issues/69
   - pint <0.22  # https://github.com/openforcefield/openff-units/issues/69
   - openff-models >=0.0.5
+  - pdbinf
   - pip
   - pydantic >1
   - pytest

gufe/components/explicitmoleculecomponent.py

@@ -2,12 +2,12 @@
 import numpy as np
 import warnings
 from rdkit import Chem
-from typing import Optional
+from typing import Any, Optional
 
 from .component import Component
 
-# typing
 from ..custom_typing import RDKitMol
+from ..molhashing import deserialize_numpy, serialize_numpy
 
 
 def _ensure_ofe_name(mol: RDKitMol, name: str) -> str:
@@ -37,6 +37,143 @@ def _ensure_ofe_name(mol: RDKitMol, name: str) -> str:
     return name
 
 
+_INT_TO_ATOMCHIRAL = {
+    0: Chem.rdchem.ChiralType.CHI_UNSPECIFIED,
+    1: Chem.rdchem.ChiralType.CHI_TETRAHEDRAL_CW,
+    2: Chem.rdchem.ChiralType.CHI_TETRAHEDRAL_CCW,
+    3: Chem.rdchem.ChiralType.CHI_OTHER,
+}
+# support for non-tetrahedral stereo requires rdkit 2022.09.1+
+if hasattr(Chem.rdchem.ChiralType, 'CHI_TETRAHEDRAL'):
+    _INT_TO_ATOMCHIRAL.update({
+        4: Chem.rdchem.ChiralType.CHI_TETRAHEDRAL,
+        5: Chem.rdchem.ChiralType.CHI_ALLENE,
+        6: Chem.rdchem.ChiralType.CHI_SQUAREPLANAR,
+        7: Chem.rdchem.ChiralType.CHI_TRIGONALBIPYRAMIDAL,
+        8: Chem.rdchem.ChiralType.CHI_OCTAHEDRAL,
+    })
+_ATOMCHIRAL_TO_INT = {v: k for k, v in _INT_TO_ATOMCHIRAL.items()}
+
+
+_INT_TO_BONDTYPE = {
+    0: Chem.rdchem.BondType.UNSPECIFIED,
+    1: Chem.rdchem.BondType.SINGLE,
+    2: Chem.rdchem.BondType.DOUBLE,
+    3: Chem.rdchem.BondType.TRIPLE,
+    4: Chem.rdchem.BondType.QUADRUPLE,
+    5: Chem.rdchem.BondType.QUINTUPLE,
+    6: Chem.rdchem.BondType.HEXTUPLE,
+    7: Chem.rdchem.BondType.ONEANDAHALF,
+    8: Chem.rdchem.BondType.TWOANDAHALF,
+    9: Chem.rdchem.BondType.THREEANDAHALF,
+    10: Chem.rdchem.BondType.FOURANDAHALF,
+    11: Chem.rdchem.BondType.FIVEANDAHALF,
+    12: Chem.rdchem.BondType.AROMATIC,
+    13: Chem.rdchem.BondType.IONIC,
+    14: Chem.rdchem.BondType.HYDROGEN,
+    15: Chem.rdchem.BondType.THREECENTER,
+    16: Chem.rdchem.BondType.DATIVEONE,
+    17: Chem.rdchem.BondType.DATIVE,
+    18: Chem.rdchem.BondType.DATIVEL,
+    19: Chem.rdchem.BondType.DATIVER,
+    20: Chem.rdchem.BondType.OTHER,
+    21: Chem.rdchem.BondType.ZERO}
+_BONDTYPE_TO_INT = {v: k for k, v in _INT_TO_BONDTYPE.items()}
+_INT_TO_BONDSTEREO = {
+    0: Chem.rdchem.BondStereo.STEREONONE,
+    1: Chem.rdchem.BondStereo.STEREOANY,
+    2: Chem.rdchem.BondStereo.STEREOZ,
+    3: Chem.rdchem.BondStereo.STEREOE,
+    4: Chem.rdchem.BondStereo.STEREOCIS,
+    5: Chem.rdchem.BondStereo.STEREOTRANS}
+_BONDSTEREO_TO_INT = {v: k for k, v in _INT_TO_BONDSTEREO.items()}
+
+
+def _setprops(obj, d: dict) -> None:
+    # add props onto rdkit "obj" (atom/bond/mol/conformer)
+    # props are guaranteed one of Bool, Int, Float or String type
+    for k, v in d.items():
+        if isinstance(v, bool):
+            obj.SetBoolProp(k, v)
+        elif isinstance(v, int):
+            obj.SetIntProp(k, v)
+        elif isinstance(v, float):
+            obj.SetDoubleProp(k, v)
+        else:  # isinstance(v, str):
+            obj.SetProp(k, v)
+
+
+def _mol_from_dict(d: dict) -> Chem.Mol:
+    m = Chem.Mol()
+    em = Chem.EditableMol(m)
+
+    for atom in d['atoms']:
+        a = Chem.Atom(atom[0])
+        a.SetIsotope(atom[1])
+        a.SetFormalCharge(atom[2])
+        a.SetIsAromatic(atom[3])
+        a.SetChiralTag(_INT_TO_ATOMCHIRAL[atom[4]])
+        a.SetAtomMapNum(atom[5])
+        _setprops(a, atom[6])
+        em.AddAtom(a)
+
+    for bond in d['bonds']:
+        em.AddBond(bond[0], bond[1], _INT_TO_BONDTYPE[bond[2]])
+        # other fields are applied onto the ROMol
+
+    m = em.GetMol()
+
+    for bond, b in zip(d['bonds'], m.GetBonds()):
+        b.SetStereo(_INT_TO_BONDSTEREO[bond[3]])
+        _setprops(b, bond[4])
+
+    pos = deserialize_numpy(d['conformer'][0])
+    c = Chem.Conformer(m.GetNumAtoms())
+    for i, p in enumerate(pos):
+        c.SetAtomPosition(i, p)
+    _setprops(c, d['conformer'][1])
+    m.AddConformer(c)
+
+    _setprops(m, d['molprops'])
+
+    m.UpdatePropertyCache()
+
+    return m
+
+
+def _mol_to_dict(m: Chem.Mol) -> dict[str, Any]:
+    # in a perfect world we'd use ToBinary()
+    # but this format slowly evolves, so the future hash of a SMC could change if rdkit were updated
+    # this is based on that method, with some irrelevant fields cut out
+
+    output: dict[str, Any] = {}
+
+    atoms = []
+    for atom in m.GetAtoms():
+        atoms.append((
+            atom.GetAtomicNum(), atom.GetIsotope(), atom.GetFormalCharge(), atom.GetIsAromatic(),
+            _ATOMCHIRAL_TO_INT[atom.GetChiralTag()], atom.GetAtomMapNum(),
+            atom.GetPropsAsDict(includePrivate=False),
+        ))
+    output['atoms'] = atoms
+
+    bonds = []
+    for bond in m.GetBonds():
+        bonds.append((
+            bond.GetBeginAtomIdx(), bond.GetEndAtomIdx(), _BONDTYPE_TO_INT[bond.GetBondType()],
+            _BONDSTEREO_TO_INT[bond.GetStereo()],
+            bond.GetPropsAsDict(includePrivate=False)
+        ))
+    output['bonds'] = bonds
+
+    conf = m.GetConformer()
+    output['conformer'] = (serialize_numpy(conf.GetPositions()), conf.GetPropsAsDict(includePrivate=False))
+
+    output['molprops'] = m.GetPropsAsDict(includePrivate=False)
+
+    return output
+
+
 def _check_partial_charges(mol: RDKitMol) -> None:
     """
     Checks for the presence of partial charges.
@@ -131,15 +268,13 @@ def smiles(self) -> str:
     def total_charge(self):
         return Chem.GetFormalCharge(self._rdkit)
 
-    # TO
     def to_rdkit(self) -> RDKitMol:
         """Return an RDKit copied representation of this molecule"""
         return Chem.Mol(self._rdkit)
 
     def to_json(self):
         return json.dumps(self.to_dict())
 
-    # From
     @classmethod
     def from_json(cls, json_str):
         dct = json.loads(json_str)