diff --git a/.coveragerc b/.coveragerc index d93622c1..7a8d0b2e 100644 --- a/.coveragerc +++ b/.coveragerc @@ -9,3 +9,4 @@ exclude_lines = pragma: no-cover -no-cov raise NotImplementedError + ... diff --git a/gufe/__init__.py b/gufe/__init__.py index f860a85b..87a8fec3 100644 --- a/gufe/__init__.py +++ b/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 diff --git a/gufe/chemicalstate.py b/gufe/chemicalsystem.py similarity index 88% rename from gufe/chemicalstate.py rename to gufe/chemicalsystem.py index 106bf189..a80ffb66 100644 --- a/gufe/chemicalstate.py +++ b/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): diff --git a/gufe/component.py b/gufe/component.py index 62109e87..fd82e879 100644 --- a/gufe/component.py +++ b/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""" + ... diff --git a/gufe/proteincomponent.py b/gufe/proteincomponent.py index 6eb9dc6a..79eeee1c 100644 --- a/gufe/proteincomponent.py +++ b/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) diff --git a/gufe/smallmoleculecomponent.py b/gufe/smallmoleculecomponent.py index f6599240..37b6daf4 100644 --- a/gufe/smallmoleculecomponent.py +++ b/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) diff --git a/gufe/solventcomponent.py b/gufe/solventcomponent.py index cb161a51..9567e287 100644 --- a/gufe/solventcomponent.py +++ b/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 diff --git a/gufe/tests/test_chemicalstate.py b/gufe/tests/test_chemicalsystem.py similarity index 69% rename from gufe/tests/test_chemicalstate.py rename to gufe/tests/test_chemicalsystem.py index ced1d8ed..290d1264 100644 --- a/gufe/tests/test_chemicalstate.py +++ b/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,7 +40,7 @@ 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}, ) @@ -48,7 +48,7 @@ def solvated_ligand(solv_comp, toluene_ligand_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 diff --git a/gufe/tests/test_protein_molecule.py b/gufe/tests/test_protein_molecule.py index f4ae2e57..ab6d8310 100644 --- a/gufe/tests/test_protein_molecule.py +++ b/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 diff --git a/gufe/tests/test_smallmoleculecomponent.py b/gufe/tests/test_smallmoleculecomponent.py index bde71cd5..af19927d 100644 --- a/gufe/tests/test_smallmoleculecomponent.py +++ b/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 diff --git a/gufe/tests/test_solvents.py b/gufe/tests/test_solvents.py index 07d60069..8c8f09f0 100644 --- a/gufe/tests/test_solvents.py +++ b/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