Merge pull request #10 from OpenFreeEnergy/formal-charge

Added formal charge to components, chemicalstate

.coveragerc

@@ -9,3 +9,4 @@ exclude_lines =
     pragma: no-cover
     -no-cov
     raise NotImplementedError
+    ...

gufe/__init__.py

@@ -7,4 +7,4 @@
 from .smallmoleculecomponent import SmallMoleculeComponent
 from .proteincomponent import ProteinComponent
 from .solventcomponent import SolventComponent
-from .chemicalstate import ChemicalState
+from .chemicalsystem import ChemicalSystem

gufe/chemicalstate.py → gufe/chemicalsystem.py

@@ -9,7 +9,7 @@
 from .component import Component
 
 
-class ChemicalState(Serializable):
+class ChemicalSystem(Serializable):
     """A node of an alchemical network.
 
     Attributes
@@ -114,9 +114,18 @@ def identifier(self):
         return self._identifier
 
     @property
-    def charge(self):
-        """Total charge for the ChemicalState."""
-        return sum([component.charge for component in self._components.values()])
+    def total_charge(self):
+        """Formal charge for the ChemicalSystem."""
+        # This might evaluate the property twice?
+        #return sum(component.total_charge
+        #           for component in self._components.values()
+        #           if component.total_charge is not None)
+        total_charge = 0
+        for c in self._components.values():
+            fc = c.total_charge
+            if fc is not None:
+                total_charge += fc
+        return total_charge
 
     @classmethod
     def as_protein_smallmolecule_solvent(cls):

gufe/component.py

@@ -2,10 +2,11 @@
 # For details, see https://github.com/OpenFreeEnergy/gufe
 
 import abc
+from typing import Union
 
 
 class Component(abc.ABC):
-    """Base class for members of a ChemicalState"""
+    """Base class for members of a ChemicalSystem"""
     @abc.abstractmethod
     def __hash__(self):
         pass
@@ -22,3 +23,9 @@ def to_dict(self) -> dict:
     @abc.abstractmethod
     def from_dict(cls, d: dict):
         pass
+
+    @property
+    @abc.abstractmethod
+    def total_charge(self) -> Union[int, None]:
+        """Net formal charge for the Component if defined"""
+        ...

gufe/proteincomponent.py

@@ -7,7 +7,7 @@
 
 
 class ProteinComponent(Component):
-    """Wrapper around a Protein representation
+    """Wrapper around a Protein representation.
 
     This representation is immutable.  If you want to make any modifications,
     do this in an appropriate toolkit then remake this class.
@@ -65,3 +65,7 @@ def __hash__(self):
 
     def __eq__(self, other):
         return hash(self) == hash(other)
+
+    @property
+    def total_charge(self):
+        return Chem.GetFormalCharge(self._rdkit)

gufe/smallmoleculecomponent.py

@@ -51,7 +51,7 @@ class SmallMoleculeComponent(Component):
     .. note::
        This class is a read-only representation of a molecule, if you want to
        edit the molecule do this in an appropriate toolkit **before** creating
-       this class.
+       an instance from this class.
 
     This wrapper uses SMILES as the primary hash, so is best suited to smaller
     molecules.  It also supports reading/writing to .sdf format, which again
@@ -205,3 +205,7 @@ def _from_sdf_supplier(cls, supp):
             raise RuntimeError(f"SDF contains more than 1 molecule")
 
         return cls(rdkit=mol)  # name is obtained automatically
+
+    @property
+    def total_charge(self):
+        return Chem.GetFormalCharge(self._rdkit)

gufe/solventcomponent.py

@@ -1,6 +1,7 @@
 # This code is part of OpenFE and is licensed under the MIT license.
 # For details, see https://github.com/OpenFreeEnergy/gufe
 from typing import Tuple
+import math
 from openff.units import unit
 
 from gufe import Component
@@ -65,6 +66,11 @@ def concentration(self) -> unit.Quantity:
         """Concentration of ions in the solvent state"""
         return self._concentration
 
+    @property
+    def total_charge(self):
+        """Solvents don't have a formal charge defined so this returns None"""
+        return None
+
     def __eq__(self, other):
         try:
             return (self.smiles == other.smiles and

gufe/tests/test_chemicalstate.py → gufe/tests/test_chemicalsystem.py

@@ -31,7 +31,7 @@ def phenol_ligand_comp(benzene_modifications):
 
 @pytest.fixture
 def solvated_complex(prot_comp, solv_comp, toluene_ligand_comp):
-    return gufe.ChemicalState(
+    return gufe.ChemicalSystem(
         {'protein': prot_comp,
          'solvent': solv_comp,
          'ligand': toluene_ligand_comp},
@@ -40,15 +40,15 @@ def solvated_complex(prot_comp, solv_comp, toluene_ligand_comp):
 
 @pytest.fixture
 def solvated_ligand(solv_comp, toluene_ligand_comp):
-    return gufe.ChemicalState(
+    return gufe.ChemicalSystem(
         {'ligand': toluene_ligand_comp, 'solvent': solv_comp},
     )
 
 
 def test_ligand_construction(solv_comp, toluene_ligand_comp):
     # sanity checks on construction
 
-    state = gufe.ChemicalState(
+    state = gufe.ChemicalSystem(
         {'solvent': solv_comp,
          'ligand': toluene_ligand_comp},
     )
@@ -62,7 +62,7 @@ def test_ligand_construction(solv_comp, toluene_ligand_comp):
 def test_complex_construction(prot_comp, solv_comp, toluene_ligand_comp):
     # sanity checks on construction
 
-    state = gufe.ChemicalState(
+    state = gufe.ChemicalSystem(
         {'protein': prot_comp,
          'solvent': solv_comp,
          'ligand': toluene_ligand_comp},
@@ -76,28 +76,28 @@ def test_complex_construction(prot_comp, solv_comp, toluene_ligand_comp):
 
 
 def test_hash_and_eq(prot_comp, solv_comp, toluene_ligand_comp):
-    c1 = gufe.ChemicalState({'protein': prot_comp,
-                             'solvent': solv_comp,
-                             'ligand': toluene_ligand_comp})
+    c1 = gufe.ChemicalSystem({'protein': prot_comp,
+                              'solvent': solv_comp,
+                              'ligand': toluene_ligand_comp})
 
-    c2 = gufe.ChemicalState({'solvent': solv_comp,
-                             'ligand': toluene_ligand_comp,
-                             'protein': prot_comp})
+    c2 = gufe.ChemicalSystem({'solvent': solv_comp,
+                              'ligand': toluene_ligand_comp,
+                              'protein': prot_comp})
 
     assert c1 == c2
     assert hash(c1) == hash(c2)
 
 
-def test_chemical_state_neq_1(solvated_complex, prot_comp):
+def test_chemical_system_neq_1(solvated_complex, prot_comp):
     # wrong class
     assert solvated_complex != prot_comp
     assert hash(solvated_complex) != hash(prot_comp)
 
 
-def test_chemical_state_neq_2(solvated_complex, prot_comp, solv_comp,
-                              toluene_ligand_comp):
+def test_chemical_system_neq_2(solvated_complex, prot_comp, solv_comp,
+                               toluene_ligand_comp):
     # identifiers are different
-    complex2 = gufe.ChemicalState(
+    complex2 = gufe.ChemicalSystem(
         {'protein': prot_comp,
          'solvent': solv_comp,
          'ligand': toluene_ligand_comp},
@@ -108,10 +108,10 @@ def test_chemical_state_neq_2(solvated_complex, prot_comp, solv_comp,
     assert hash(solvated_complex) != hash(complex2)
 
 
-def test_chemical_state_neq_3(solvated_complex, prot_comp, solv_comp,
-                              toluene_ligand_comp):
+def test_chemical_system_neq_3(solvated_complex, prot_comp, solv_comp,
+                               toluene_ligand_comp):
     # different unit cell size
-    complex2 = gufe.ChemicalState(
+    complex2 = gufe.ChemicalSystem(
         {'protein': prot_comp,
          'solvent': solv_comp,
          'ligand': toluene_ligand_comp},
@@ -121,19 +121,28 @@ def test_chemical_state_neq_3(solvated_complex, prot_comp, solv_comp,
     assert hash(solvated_complex) != hash(complex2)
 
 
-def test_chemical_state_neq_4(solvated_complex, solvated_ligand):
+def test_chemical_system_neq_4(solvated_complex, solvated_ligand):
     # different component keys
     assert solvated_complex != solvated_ligand
     assert hash(solvated_complex) != hash(solvated_ligand)
 
 
-def test_chemical_state_neq_5(solvated_complex, prot_comp, solv_comp,
-                              phenol_ligand_comp):
+def test_chemical_system_neq_5(solvated_complex, prot_comp, solv_comp,
+                               phenol_ligand_comp):
     # same component keys, but different components
-    complex2 = gufe.ChemicalState(
+    complex2 = gufe.ChemicalSystem(
         {'protein': prot_comp,
          'solvent': solv_comp,
          'ligand': phenol_ligand_comp},
     )
     assert solvated_complex != complex2
     assert hash(solvated_complex) != hash(complex2)
+
+
+def test_complex_system_charge(solvated_complex):
+    # protein = 22, ligand = 0, solvent = 0
+    assert solvated_complex.total_charge == 22
+
+
+def test_ligand_system_charge(solvated_ligand):
+    assert solvated_ligand.total_charge == 0

gufe/tests/test_protein_molecule.py

@@ -88,3 +88,8 @@ def test_hash_neq_name(PDB_181L_path):
 
     assert hash(m1) != hash(m2)
 
+
+def test_protein_total_charge(PDB_181L_path):
+    m1 = ProteinComponent.from_pdbfile(PDB_181L_path)
+
+    assert m1.total_charge == 22

gufe/tests/test_smallmoleculecomponent.py

@@ -166,3 +166,12 @@ def test_to_oechem(self, ethane):
         oec_ethane = ethane.to_openeye()
 
         assert isinstance(oec_ethane, oechem.OEMol)
+
+
+@pytest.mark.parametrize('mol, charge', [
+    ('CC', 0), ('CC[O-]', -1),
+])
+def test_total_charge_neutral(mol, charge):
+    sm = SmallMoleculeComponent.from_rdkit(Chem.MolFromSmiles(mol))
+
+    assert sm.total_charge == charge

gufe/tests/test_solvents.py

@@ -51,3 +51,9 @@ def test_conc():
     s = SolventComponent(ions=('Na', 'Cl'), concentration=1.75 * unit.molar)
 
     assert s.concentration == unit.Quantity('1.75 M')
+
+
+def test_solvent_charge():
+    s = SolventComponent(ions=('Na', 'Cl'), concentration=1.75 * unit.molar)
+
+    assert s.total_charge is None