Merge branch 'main' into cookbook-flowchart

.github/workflows/ci.yaml

@@ -74,7 +74,7 @@ jobs:
           # Set the OFE_SLOW_TESTS to True if running a Cron job
           OFE_SLOW_TESTS: ${{ fromJSON('{"false":"false","true":"true"}')[github.event_name != 'pull_request'] }}
         run: |
-          pytest -n auto -v --cov=openfe --cov=openfecli --cov-report=xml
+          pytest -n auto -v --cov=openfe --cov=openfecli --cov-report=xml --durations=10
 
       - name: codecov
         if: ${{ github.repository == 'OpenFreeEnergy/openfe'

environment.yml

@@ -13,6 +13,7 @@ dependencies:
   - pytest-rerunfailures
   - pydantic<2.0
   - coverage
+  - cinnabar ==0.3.0
   - openff-toolkit>=0.13.0
   - openff-units==0.2.0
   - pint<0.22

openfe/protocols/openmm_rfe/_rfe_utils/relative.py

@@ -86,9 +86,6 @@ def __init__(self,
                  softcore_alpha=0.5,
                  softcore_LJ_v2=True,
                  softcore_LJ_v2_alpha=0.85,
-                 softcore_electrostatics=True,
-                 softcore_electrostatics_alpha=0.3,
-                 softcore_sigma_Q=1.0,
                  interpolate_old_and_new_14s=False,
                  flatten_torsions=False,
                  **kwargs):
@@ -125,12 +122,6 @@ def __init__(self,
             Implement the softcore LJ as defined by Gapsys et al. JCTC 2012.
         softcore_LJ_v2_alpha : float, default 0.85
             Softcore alpha parameter for LJ v2
-        softcore_electrostatics : bool, default True
-            Use softcore electrostatics as defined by Gapsys et al. JCTC 2021.
-        softcore_electrostatics_alpha : float, default 0.3
-            Softcore alpha parameter for softcore electrostatics.
-        softcore_sigma_Q : float, default 1.0
-            Softcore sigma parameter for softcore electrostatics.
         interpolate_old_and_new_14s : bool, default False
             Whether to turn off interactions for new exceptions (not just
             1,4s) at lambda = 0 and old exceptions at lambda = 1; if False,
@@ -172,14 +163,6 @@ def __init__(self,
             self._check_bounds(softcore_LJ_v2_alpha, "softcore_LJ_v2_alpha")
             self._softcore_LJ_v2_alpha = softcore_LJ_v2_alpha
 
-        self._softcore_electrostatics = softcore_electrostatics
-        if self._softcore_electrostatics:
-            self._softcore_electrostatics_alpha = softcore_electrostatics_alpha
-            self._check_bounds(softcore_electrostatics_alpha,
-                               "softcore_electrostatics_alpha")
-            self._softcore_sigma_Q = softcore_sigma_Q
-            self._check_bounds(softcore_sigma_Q, "softcore_sigma_Q")
-
         # TODO: end __init__ here and move everything else to
         # create_hybrid_system() or equivalent
 

openfe/protocols/openmm_rfe/equil_rfe_methods.py

@@ -489,9 +489,6 @@ def run(self, *, dry=False, verbose=True,
             softcore_alpha=alchem_settings.softcore_alpha,
             softcore_LJ_v2=alchem_settings.softcore_LJ_v2,
             softcore_LJ_v2_alpha=alchem_settings.softcore_alpha,
-            softcore_electrostatics=alchem_settings.softcore_electrostatics,
-            softcore_electrostatics_alpha=alchem_settings.softcore_electrostatics_alpha,
-            softcore_sigma_Q=alchem_settings.softcore_sigma_Q,
             interpolate_old_and_new_14s=alchem_settings.interpolate_old_and_new_14s,
             flatten_torsions=alchem_settings.flatten_torsions,
         )

openfe/protocols/openmm_rfe/equil_rfe_settings.py

@@ -26,6 +26,9 @@
 
 
 class AlchemicalSettings(SettingsBaseModel):
+    class Config:
+        extra = 'ignore'
+
     """Settings for the alchemical protocol
 
     This describes the lambda schedule and the creation of the
@@ -56,16 +59,8 @@ class AlchemicalSettings(SettingsBaseModel):
     Whether to use the LJ softcore function as defined by
     Gapsys et al. JCTC 2012 Default True.
     """
-    softcore_electrostatics = True
-    """Whether to use softcore electrostatics. Default True."""
     softcore_alpha = 0.85
     """Softcore alpha parameter. Default 0.85"""
-    softcore_electrostatics_alpha = 0.3
-    """Softcore alpha parameter for electrostatics. Default 0.3"""
-    softcore_sigma_Q = 1.0
-    """
-    Softcore sigma parameter for softcore electrostatics. Default 1.0.
-    """
     interpolate_old_and_new_14s = False
     """
     Whether to turn off interactions for new exceptions (not just 1,4s)

openfe/setup/atom_mapping/lomap_scorers.py

@@ -291,10 +291,11 @@ def creates_heterocyle(mol):
 
 def transmuting_methyl_into_ring_score(mapping: LigandAtomMapping,
                                        beta=0.1, penalty=6.0):
-    """Penalises ring forming
+    """Penalises having a non-mapped ring atoms become a non-ring
 
-    Check if any atoms transition to/from rings in the mapping, if so
-    returns a score of::
+    This score would for example penalise R-CH3 to R-Ph where R is the same
+    mapped atom and both CH3 and Ph are unmapped. Does not penalise R-H to R-Ph.
+    If any atoms trigger the rule returns a score of::
 
       1 - exp(-1 * beta * penalty)
 
@@ -310,14 +311,24 @@ def transmuting_methyl_into_ring_score(mapping: LigandAtomMapping,
     """
     molA = mapping.componentA.to_rdkit()
     molB = mapping.componentB.to_rdkit()
-    molA_to_molB = mapping.componentA_to_componentB
+    # filter out hydrogens from mapping
+    # this matches the current implementation in Lomap
+    molA_to_molB = {}
+    for i, j in mapping.componentA_to_componentB.items():
+        if molA.GetAtomWithIdx(i).GetAtomicNum() == 1:
+            continue
+        if molB.GetAtomWithIdx(j).GetAtomicNum() == 1:
+            continue
+        molA_to_molB[i] = j
 
     ringbreak = False
     for i, j in molA_to_molB.items():
         atomA = molA.GetAtomWithIdx(i)
 
         for bA in atomA.GetBonds():
             otherA = bA.GetOtherAtom(atomA)
+            if otherA.GetAtomicNum() == 1:
+                continue
             if otherA.GetIdx() in molA_to_molB:
                 # if other end of bond in core, ignore
                 continue
@@ -326,6 +337,8 @@ def transmuting_methyl_into_ring_score(mapping: LigandAtomMapping,
             atomB = molB.GetAtomWithIdx(j)
             for bB in atomB.GetBonds():
                 otherB = bB.GetOtherAtom(atomB)
+                if otherB.GetAtomicNum() == 1:
+                    continue
                 if otherB.GetIdx() in molA_to_molB.values():
                     continue
 

openfe/tests/protocols/test_openmm_equil_rfe_protocols.py

@@ -200,7 +200,7 @@ def test_dry_run_default_vacuum(benzene_vacuum_system, toluene_vacuum_system,
         htf = sampler._hybrid_factory
         # 16 atoms:
         # 11 common atoms, 1 extra hydrogen in benzene, 4 extra in toluene
-        # 12 bonds in benzene + 4 extra tolunee bonds
+        # 12 bonds in benzene + 4 extra toluene bonds
         assert len(list(htf.hybrid_topology.atoms)) == 16
         assert len(list(htf.omm_hybrid_topology.atoms())) == 16
         assert len(list(htf.hybrid_topology.bonds)) == 16
@@ -211,6 +211,10 @@ def test_dry_run_default_vacuum(benzene_vacuum_system, toluene_vacuum_system,
         assert len(list(ret_top.atoms())) == 16
         assert len(list(ret_top.bonds())) == 16
 
+        # check that our PDB has the right number of atoms
+        pdb = mdt.load_pdb('hybrid_system.pdb')
+        assert pdb.n_atoms == 16
+
 
 def test_dry_run_gaff_vacuum(benzene_vacuum_system, toluene_vacuum_system,
                              benzene_to_toluene_mapping, tmpdir):
@@ -373,6 +377,7 @@ def test_dry_run_ligand(benzene_system, toluene_system,
     settings = openmm_rfe.RelativeHybridTopologyProtocol.default_settings()
     settings.alchemical_sampler_settings.sampler_method = method
     settings.alchemical_sampler_settings.n_repeats = 1
+    settings.simulation_settings.output_indices = 'resname UNK'
 
     protocol = openmm_rfe.RelativeHybridTopologyProtocol(
             settings=settings,
@@ -392,6 +397,10 @@ def test_dry_run_ligand(benzene_system, toluene_system,
                           MonteCarloBarostat)
         assert sampler._thermodynamic_states[1].pressure == 1 * omm_unit.bar
 
+        # Check we have the right number of atoms in the PDB
+        pdb = mdt.load_pdb('hybrid_system.pdb')
+        assert pdb.n_atoms == 16
+
 
 def test_dry_run_ligand_tip4p(benzene_system, toluene_system,
                               benzene_to_toluene_mapping, tmpdir):
@@ -592,6 +601,7 @@ def test_dry_run_complex(benzene_complex_system, toluene_complex_system,
     settings = openmm_rfe.RelativeHybridTopologyProtocol.default_settings()
     settings.alchemical_sampler_settings.sampler_method = method
     settings.alchemical_sampler_settings.n_repeats = 1
+    settings.simulation_settings.output_indices = 'protein or resname  UNK'
 
     protocol = openmm_rfe.RelativeHybridTopologyProtocol(
             settings=settings,
@@ -611,6 +621,10 @@ def test_dry_run_complex(benzene_complex_system, toluene_complex_system,
                           MonteCarloBarostat)
         assert sampler._thermodynamic_states[1].pressure == 1 * omm_unit.bar
 
+        # Check we have the right number of atoms in the PDB
+        pdb = mdt.load_pdb('hybrid_system.pdb')
+        assert pdb.n_atoms == 2629
+
 
 def test_lambda_schedule_default():
     lambdas = openmm_rfe._rfe_utils.lambdaprotocol.LambdaProtocol(functions='default')

openfe/tests/protocols/test_rfe_tokenization.py

@@ -39,7 +39,7 @@ def instance(self):
 
 class TestRelativeHybridTopologyProtocol(GufeTokenizableTestsMixin):
     cls = openmm_rfe.RelativeHybridTopologyProtocol
-    key = "RelativeHybridTopologyProtocol-ba8bae8ffc0f48e0839312daa7e74e98"
+    key = "RelativeHybridTopologyProtocol-8fe3eb4c318673db7d57c1bffca602df"
     repr = f"<{key}>"
 
     @pytest.fixture()

openfe/tests/setup/atom_mapping/test_lomap_scorers.py

@@ -265,3 +265,34 @@ def test_lomap_regression(lomap_basic_test_files_dir,  # in a dir for lomap
         scores[i, i] = 0
 
     assert_allclose(matrix, scores)
+
+
+def test_transmuting_methyl_into_ring_score():
+    """
+    Sets up two mappings:
+      RC_to_RPh = [CCC]C -> [CCC]Ph
+      RH_to_RPh = [CCC]H -> [CCC]Ph
+    Where square brackets show mapped (core) region
+
+    The first mapping should trigger this rule, the second shouldn't
+    """
+    def makemol(smi):
+        m = Chem.MolFromSmiles(smi)
+        m = Chem.AddHs(m)
+        m.Compute2DCoords()
+
+        return openfe.SmallMoleculeComponent(m)
+
+    core = 'CCC{}'
+    RC = makemol(core.format('C'))
+    RPh = makemol(core.format('c1ccccc1'))
+    RH = makemol(core.format('[H]'))
+
+    RC_to_RPh = openfe.LigandAtomMapping(RC, RPh, {i: i for i in range(3)})
+    RH_to_RPh = openfe.LigandAtomMapping(RH, RPh, {i: i for i in range(3)})
+
+    score1 = lomap_scorers.transmuting_methyl_into_ring_score(RC_to_RPh)
+    score2 = lomap_scorers.transmuting_methyl_into_ring_score(RH_to_RPh)
+
+    assert score2 == 0.0
+    assert score1 == pytest.approx(1 - math.exp(-0.1 * 6.0))

openfe/tests/utils/test_system_probe.py

@@ -1,9 +1,11 @@
 import contextlib
 from collections import namedtuple
+import pathlib
 import sys
 from unittest.mock import Mock, patch
 
 import psutil
+from psutil._common import sdiskusage
 import pytest
 
 from openfe.utils.system_probe import (
@@ -38,6 +40,17 @@
 )
 
 
+def fake_disk_usage(path):
+    if str(path) == "/foo":
+        return sdiskusage(
+            total=1958854045696, used=1232985415680, free=626288726016, percent=66.3
+        )
+    if str(path) == "/bar":
+        return sdiskusage(
+            total=4000770252800, used=1678226952192, free=2322615496704, percent=41.9
+        )
+
+
 @contextlib.contextmanager
 def patch_system():
     # single patch to fully patch the system
@@ -91,6 +104,9 @@ def patch_system():
             )
         ),
     )
+    patch_psutil_disk_usage = patch(
+        "psutil.disk_usage", Mock(side_effect=fake_disk_usage)
+    )
 
     # assumes that each shell command is called in only one way
     cmd_to_output = {
@@ -111,10 +127,11 @@ def patch_system():
     )
     with contextlib.ExitStack() as stack:
         for ctx in [
+            patch_check_output,
             patch_hostname,
             patch_psutil_Process_as_dict,
-            patch_check_output,
             patch_psutil_Process_rlimit,
+            patch_psutil_disk_usage,
             patch_psutil_virtual_memory,
         ]:
             stack.enter_context(ctx)
@@ -229,6 +246,30 @@ def test_get_disk_usage():
         assert disk_info == expected_disk_info
 
 
+@pytest.mark.skipif(
+    sys.platform == "darwin", reason="test requires psutil.Process.rlimit"
+)
+def test_get_disk_usage_with_path():
+    with patch_system():
+        disk_info = _get_disk_usage(paths=[pathlib.Path("/foo"), pathlib.Path("/bar")])
+        expected_disk_info = {
+            "/bar": {
+                "available": "2.1T",
+                "percent_used": "42%",
+                "size": "3.6T",
+                "used": "1.5T",
+            },
+            "/foo": {
+                "available": "583.3G",
+                "percent_used": "66%",
+                "size": "1.8T",
+                "used": "1.1T",
+            },
+        }
+
+        assert disk_info == expected_disk_info
+
+
 @pytest.mark.skipif(
     sys.platform == "darwin", reason="test requires psutil.Process.rlimit"
 )
@@ -317,3 +358,4 @@ def test_probe_system():
 
 def test_probe_system_smoke_test():
     _probe_system()
+    _probe_system(paths=[pathlib.Path("/")])