Merge branch 'main' into shared-object-v2

docs/api.rst

@@ -8,5 +8,6 @@ API docs
 
    api/components
    api/protocols
+   api/transformation
    api/settings
 

docs/api/transformation.rst

@@ -0,0 +1,12 @@
+GUFE Transformation API
+-----------------------
+
+Alchemical network edges
+========================
+
+.. autoclass:: gufe.transformations.Transformation
+           :members:
+
+.. autoclass:: gufe.transformations.NonTransformation
+           :members:
+           :show-inheritance:

environment.yml

@@ -11,7 +11,7 @@ dependencies:
   - pint <0.22  # https://github.com/openforcefield/openff-units/issues/69
   - openff-models >=0.0.5
   - pip
-  - pydantic
+  - pydantic <2.0.0
   - pytest
   - pytest-cov
   - pytest-xdist
@@ -20,5 +20,6 @@ dependencies:
   # docs
   - pydata-sphinx-theme
   - sphinx-jsonschema==1.15
+  - sphinx <7.1.2
   - pip:
-    - autodoc_pydantic
+    - autodoc_pydantic<2.0.0

gufe/components/explicitmoleculecomponent.py

@@ -1,4 +1,5 @@
 import json
+import numpy as np
 import warnings
 from rdkit import Chem
 from typing import Optional
@@ -36,6 +37,48 @@ def _ensure_ofe_name(mol: RDKitMol, name: str) -> str:
     return name
 
 
+def _check_partial_charges(mol: RDKitMol) -> None:
+    """
+    Checks for the presence of partial charges.
+
+    Raises
+    ------
+    ValueError
+      * If the partial charges are not of length atoms.
+      * If the sum of partial charges is not equal to the
+        formal charge.
+    UserWarning
+      * If partial charges are found.
+      * If the partial charges are near 0 for all atoms.
+    """
+    if 'atom.dprop.PartialCharge' not in mol.GetPropNames():
+        return
+
+    p_chgs = np.array(
+        mol.GetProp('atom.dprop.PartialCharge').split(), dtype=float
+    )
+
+    if len(p_chgs) != mol.GetNumAtoms():
+        errmsg = (f"Incorrect number of partial charges: {len(p_chgs)} "
+                  f" were provided for {mol.GetNumAtoms()} atoms")
+        raise ValueError(errmsg)
+
+    if (sum(p_chgs) - Chem.GetFormalCharge(mol)) > 0.01:
+        errmsg = (f"Sum of partial charges {sum(p_chgs)} differs from "
+                  f"RDKit formal charge {Chem.GetFormalCharge(mol)}")
+        raise ValueError(errmsg)
+
+    if np.all(np.isclose(p_chgs, 0.0)):
+        wmsg = (f"Partial charges provided all equal to "
+                "zero. These may be ignored by some Protocols.")
+        warnings.warn(wmsg)
+    else:
+        wmsg = ("Partial charges have been provided, these will "
+                "preferentially be used instead of generating new "
+                "partial charges")
+        warnings.warn(wmsg)
+
+
 class ExplicitMoleculeComponent(Component):
     """Base class for explicit molecules.
 
@@ -48,6 +91,7 @@ class ExplicitMoleculeComponent(Component):
 
     def __init__(self, rdkit: RDKitMol, name: str = ""):
         name = _ensure_ofe_name(rdkit, name)
+        _check_partial_charges(rdkit)
         conformers = list(rdkit.GetConformers())
         if not conformers:
             raise ValueError("Molecule was provided with no conformers.")

gufe/components/proteincomponent.py

@@ -27,6 +27,10 @@
     1: BondType.SINGLE,
     2: BondType.DOUBLE,
     3: BondType.TRIPLE,
+    app.Single: BondType.SINGLE,
+    app.Double: BondType.DOUBLE,
+    app.Triple: BondType.TRIPLE,
+    app.Aromatic: BondType.AROMATIC,
     None: BondType.UNSPECIFIED,
 }
 _BONDORDERS_RDKIT_TO_OPENMM = {
@@ -361,7 +365,8 @@ def chainkey(m):
             a1 = atom_lookup[bond.GetBeginAtomIdx()]
             a2 = atom_lookup[bond.GetEndAtomIdx()]
             top.addBond(a1, a2,
-                        order=_BONDORDERS_RDKIT_TO_OPENMM[bond.GetBondType()])
+                        order=_BONDORDERS_RDKIT_TO_OPENMM.get(
+                            bond.GetBondType(), None))
 
         return top
 

gufe/ligandnetwork.py

@@ -2,9 +2,11 @@
 # For details, see https://github.com/OpenFreeEnergy/gufe
 from __future__ import annotations
 
+from itertools import chain
 import json
 import networkx as nx
 from typing import FrozenSet, Iterable, Optional
+import gufe
 
 from gufe import SmallMoleculeComponent
 from .mapping import LigandAtomMapping
@@ -30,7 +32,8 @@ def __init__(
             nodes = []
 
         self._edges = frozenset(edges)
-        edge_nodes = set.union(*[{edge.componentA, edge.componentB} for edge in edges])
+        edge_nodes = set(chain.from_iterable((e.componentA, e.componentB)
+                                             for e in edges))
         self._nodes = frozenset(edge_nodes) | frozenset(nodes)
         self._graph = None
 
@@ -46,7 +49,7 @@ def _from_dict(cls, dct: dict):
         return cls.from_graphml(dct['graphml'])
 
     @property
-    def graph(self) -> nx.Graph:
+    def graph(self) -> nx.MultiDiGraph:
         """NetworkX graph for this network"""
         if self._graph is None:
             graph = nx.MultiDiGraph()
@@ -140,7 +143,6 @@ def to_graphml(self) -> str:
     @classmethod
     def from_graphml(cls, graphml_str: str):
         """Create ``Network`` from GraphML string.
-
         This is the primary deserialization mechanism for this class.
 
         Parameters
@@ -178,3 +180,137 @@ def enlarge_graph(self, *, edges=None, nodes=None) -> LigandNetwork:
             nodes = set([])
 
         return LigandNetwork(self.edges | set(edges), self.nodes | set(nodes))
+
+    def _to_rfe_alchemical_network(
+        self,
+        components: dict[str, gufe.Component],
+        leg_labels: dict[str, list[str]],
+        protocol: gufe.Protocol,
+        *,
+        alchemical_label: str = "ligand",
+        autoname=True,
+        autoname_prefix=""
+    ) -> gufe.AlchemicalNetwork:
+        """
+        Parameters
+        ----------
+        components: dict[str, :class:`.Component`]
+            non-alchemical components (components that will be on both sides
+            of a transformation)
+        leg_label: dict[str, list[str]]
+            mapping of the names for legs (the keys of this dict) to a list
+            of the component names. The componnent names must be the same as
+            as used in the ``componentns`` dict.
+        protocol: :class:`.Protocol`
+            the protocol to apply
+        alchemical_label: str
+            the label for the component undergoing an alchemical
+            transformation (default ``'ligand'``)
+        """
+        transformations = []
+        for edge in self.edges:
+            for leg_name, labels in leg_labels.items():
+
+                # define a helper func to avoid repeated code
+                def sys_from_dict(component):
+                    """
+                    Input component alchemically changing. Other info taken
+                    from the outer scope.
+                    """
+                    syscomps = {alchemical_label: component}
+                    other_labels = set(labels) - {alchemical_label}
+                    syscomps.update({label: components[label]
+                                     for label in other_labels})
+
+                    if autoname:
+                        name = f"{component.name}_{leg_name}"
+                    else:
+                        name = ""
+
+                    return gufe.ChemicalSystem(syscomps, name=name)
+
+                sysA = sys_from_dict(edge.componentA)
+                sysB = sys_from_dict(edge.componentB)
+                if autoname:
+                    prefix = f"{autoname_prefix}_" if autoname_prefix else ""
+                    name = f"{prefix}{sysA.name}_{sysB.name}"
+                else:
+                    name = ""
+
+                mapping: dict[str, gufe.ComponentMapping] = {
+                    alchemical_label: edge,
+                }
+
+                transformation = gufe.Transformation(sysA, sysB, protocol,
+                                                     mapping, name)
+
+                transformations.append(transformation)
+
+        return gufe.AlchemicalNetwork(transformations)
+
+    def to_rbfe_alchemical_network(
+        self,
+        solvent: gufe.SolventComponent,
+        protein: gufe.ProteinComponent,
+        protocol: gufe.Protocol,
+        *,
+        autoname: bool = True,
+        autoname_prefix: str = "easy_rbfe",
+        **other_components
+    ):
+        """
+        Parameters
+        ----------
+        protocol: :class:Protocol
+        autoname: bool
+            whether to automatically name objects by the ligand name and
+            state label
+        autoname_prefix: str
+            prefix for the autonaming; only used if autonaming is True
+        other_components:
+            additional non-alchemical components, keyword will be the string
+            label for the component
+        """
+        components = {
+            'protein': protein,
+            'solvent': solvent,
+            **other_components
+        }
+        leg_labels = {
+            "solvent": ["ligand", "solvent"],
+            "complex": (["ligand", "solvent", "protein"]
+                        + list(other_components)),
+        }
+        return self._to_rfe_alchemical_network(
+            components=components,
+            leg_labels=leg_labels,
+            protocol=protocol,
+            autoname=autoname,
+            autoname_prefix=autoname_prefix
+        )
+
+    # on hold until we figure out how to best hack in the PME/NoCutoff
+    # switch
+    # def to_rhfe_alchemical_network(self, *, solvent, protocol,
+    #                                autoname=True,
+    #                                autoname_prefix="easy_rhfe",
+    #                                **other_components):
+    #     leg_labels = {
+    #         "solvent": ["ligand", "solvent"] + list(other_components),
+    #         "vacuum": ["ligand"] + list(other_components),
+    #     }
+    #     return self._to_rfe_alchemical_network(
+    #         components={"solvent": solvent, **other_components},
+    #         leg_labels=leg_labels,
+    #         protocol=protocol,
+    #         autoname=autoname,
+    #         autoname_prefix=autoname_prefix
+    #     )
+
+    def is_connected(self) -> bool:
+        """Are all ligands in the network (indirectly) connected to each other
+
+        A "False" value indicates that either some ligands have no edges or that
+        there are separate networks that do not link to each other.
+        """
+        return nx.is_weakly_connected(self.graph)