Fix 1-4 interpolation for core-dummy interactions

docs/guide/protocols/relativehybridtopology.rst

@@ -39,9 +39,10 @@ ligands, meaning that a single set of coordinates is used to represent the
 common core of the two ligands while the atoms that differ between the two
 ligands are represented separately. An atom map defines which atoms belong
 to the core (mapped atoms) and which atoms are unmapped and represented
-separately (see :ref:`Creating atom mappings <Creating Atom Mappings>`). During the alchemical transformation, mapped atoms are switched
-from the type in one ligand to the type in the other ligands, while unmapped
-atoms are switched on or off, depending on which ligand they belong to.
+separately (see :ref:`Creating atom mappings <Creating Atom Mappings>`). During the alchemical transformation, mapped atoms are interpolated
+from their type in ligand at state A to the type in the other ligand at state B, while unmapped
+atoms (also known as dummy atoms) are switched inserted or uncoupled, depending on which ligand they belong to. By default all nonbonded interactions between the
+dummy region and the core region are removed to avoid coupling their motion.
 
 .. note:: In this hybrid topology approach, all bonded interactions between the dummy region and the core region are kept. 
           As pointed out by Fleck et al. [1]_, this can lead to systematic errors if the contribution of the dummy group does not cancel out

src/openfe/protocols/openmm_rfe/_rfe_utils/relative.py

@@ -1772,7 +1772,8 @@ def _check_indices(idx1, idx2):
                 # with zeroed out parameters; add old charge to
                 # standard_nonbonded and zero sterics
                 check_index = self._hybrid_system_forces['standard_nonbonded_force'].addParticle(
-                    charge_old, 0.5*(sigma_old+sigma_new), 0.0)
+                    # this term is off due to epsilon = 0, but just set sigma to the initial value to not confuse things
+                    charge_old, sigma_old, 0.0)
                 _check_indices(particle_index, check_index)
 
                 # Charge is charge_old at lambda_electrostatics = 0,

src/openfe/protocols/openmm_rfe/equil_rfe_settings.py

@@ -76,11 +76,11 @@ class AlchemicalSettings(SettingsBaseModel):
     """
     softcore_alpha: float = 0.85
     """Softcore alpha parameter. Default 0.85"""
-    turn_off_core_unique_exceptions: bool = False
+    turn_off_core_unique_exceptions: bool = True
     """
     Whether to turn off interactions for new exceptions (not just 1,4s)
     at lambda 0 and old exceptions at lambda 1 between unique atoms and core
-    atoms. If False they are present in the nonbonded force. Default False.
+    atoms. If False they are present in the nonbonded force. Default True.
     """
     explicit_charge_correction: bool = False
     """

src/openfe/tests/conftest.py

@@ -124,11 +124,6 @@ def mol_from_smiles(smiles: str) -> Chem.Mol:
     return m
 
 
-@pytest.fixture(scope="session")
-def ethane():
-    return SmallMoleculeComponent(mol_from_smiles("CC"))
-
-
 @pytest.fixture(scope="session")
 def simple_mapping():
     """Disappearing oxygen on end
@@ -398,6 +393,34 @@ def fluorobenzene():
         yield SmallMoleculeComponent.from_sdf_file(f / "t4_lysozyme_data" / "fluorobenzene.sdf")
 
 
+@pytest.fixture(scope="module")
+def ethane():
+    """Load ethane with partial charges from sdf file."""
+    with resources.as_file(resources.files("openfe.tests.data.htf")) as f:
+        yield SmallMoleculeComponent.from_sdf_file(f / "ethane.sdf")
+
+
+@pytest.fixture(scope="module")
+def chloroethane():
+    """Load chloroethane with partial charges from sdf file."""
+    with resources.as_file(resources.files("openfe.tests.data.htf")) as f:
+        yield SmallMoleculeComponent.from_sdf_file(f / "chloroethane.sdf")
+
+
+@pytest.fixture(scope="module")
+def fluoroethane():
+    """Load fluoroethane with partial charges from sdf file."""
+    with resources.as_file(resources.files("openfe.tests.data.htf")) as f:
+        yield SmallMoleculeComponent.from_sdf_file(f / "fluoroethane.sdf")
+
+
+@pytest.fixture(scope="module")
+def benzene():
+    """Load fluoroethane with partial charges from sdf file."""
+    with resources.as_file(resources.files("openfe.tests.data.htf")) as f:
+        yield SmallMoleculeComponent.from_sdf_file(f / "t4_lysozyme_data" / "benzene.sdf")
+
+
 @pytest.fixture(scope="module")
 def chlorobenzene_to_fluorobenzene_mapping(chlorobenzene, fluorobenzene):
     """Return a mapping from chlorobenzene to fluorobenzene."""
@@ -422,6 +445,68 @@ def chlorobenzene_to_fluorobenzene_mapping(chlorobenzene, fluorobenzene):
     )
 
 
+@pytest.fixture(scope="module")
+def chloroethane_to_ethane_mapping(chloroethane, ethane):
+    """Return a mapping from chloroethane to ethane."""
+    return LigandAtomMapping(
+        componentA=chloroethane,
+        componentB=ethane,
+        componentA_to_componentB={
+            # Cl-H not mapped, all others one-to-one
+            1: 1,
+            2: 2,
+            3: 3,
+            4: 4,
+            5: 5,
+            6: 6,
+            7: 7,
+        },
+    )
+
+
+@pytest.fixture(scope="module")
+def chloroethane_to_fluoroethane_mapping(chloroethane, fluoroethane):
+    """Return a mapping from chloroethane to fluoroethane."""
+    return LigandAtomMapping(
+        componentA=chloroethane,
+        componentB=fluoroethane,
+        componentA_to_componentB={
+            # perfect one-to-one mapping
+            0: 0,
+            1: 1,
+            2: 2,
+            3: 3,
+            4: 4,
+            5: 5,
+            6: 6,
+            7: 7,
+        },
+    )
+
+
+@pytest.fixture(scope="module")
+def chlorobenzene_to_benzene_mapping(chlorobenzene, benzene):
+    """Return a mapping from chlorobenzene to benzene."""
+    return LigandAtomMapping(
+        componentA=chlorobenzene,
+        componentB=benzene,
+        componentA_to_componentB={
+            # Cl-H not mapped, all others one-to-one
+            1: 1,
+            2: 2,
+            3: 3,
+            4: 4,
+            5: 5,
+            6: 6,
+            7: 7,
+            8: 8,
+            9: 9,
+            10: 10,
+            11: 11,
+        },
+    )
+
+
 @pytest.fixture(scope="module")
 def t4_lysozyme_solvated():
     """Load the T4 lysozyme L99A structure and solvent from the pdb file."""
@@ -485,8 +570,6 @@ def htf_cmap_chlorobenzene_to_fluorobenzene(
 ):
     """Generate the htf for chlorobenzene to fluorobenzene with a CMAP term."""
     settings = RelativeHybridTopologyProtocol.default_settings()
-    # make sure we interpolate the 1-4 exceptions involving dummy atoms if present
-    settings.alchemical_settings.turn_off_core_unique_exceptions = True
     small_ff = settings.forcefield_settings.small_molecule_forcefield
     if ".offxml" not in small_ff:
         small_ff += ".offxml"
@@ -519,3 +602,161 @@ def htf_cmap_chlorobenzene_to_fluorobenzene(
         "vdW_scale": ff.get_parameter_handler("vdW").scale14,
         "force_field": ff,
     }
+
+
+@pytest.fixture(scope="module")
+def htf_chloro_fluoroethane(chloroethane, fluoroethane, chloroethane_to_fluoroethane_mapping):
+    """Generate the htf for chloroethane to fluoroethane in vacuum."""
+    settings = RelativeHybridTopologyProtocol.default_settings()
+    small_ff = settings.forcefield_settings.small_molecule_forcefield
+    if ".offxml" not in small_ff:
+        small_ff += ".offxml"
+    ff = ForceField(small_ff)
+    chloro_openff = chloroethane.to_openff()
+    chloro_charges = chloro_openff.partial_charges.m_as(offunit.elementary_charge)
+    chloro_labels = ff.label_molecules(chloro_openff.to_topology())[0]
+    fluoro_openff = fluoroethane.to_openff()
+    fluoro_charges = fluoro_openff.partial_charges.m_as(offunit.elementary_charge)
+    fluoro_labels = ff.label_molecules(fluoro_openff.to_topology())[0]
+    htf = make_htf(mapping=chloroethane_to_fluoroethane_mapping, settings=settings)
+    hybrid_system = htf.hybrid_system
+    forces = {force.getName(): force for force in hybrid_system.getForces()}
+
+    return {
+        "htf": htf,
+        "hybrid_system": hybrid_system,
+        "forces": forces,
+        "chloro_labels": chloro_labels,
+        "fluoro_labels": fluoro_labels,
+        "mapping": chloroethane_to_fluoroethane_mapping,
+        "chloroethane": chloroethane,
+        "fluoroethane": fluoroethane,
+        "chloro_charges": chloro_charges,
+        "fluoro_charges": fluoro_charges,
+        "electrostatic_scale": ff.get_parameter_handler("Electrostatics").scale14,
+        "vdW_scale": ff.get_parameter_handler("vdW").scale14,
+        "force_field": ff,
+    }
+
+
+@pytest.fixture(scope="module")
+def htf_chloro_ethane(chloroethane, ethane, chloroethane_to_ethane_mapping):
+    """Generate the htf for chloroethane to ethane in vacuum."""
+    settings = RelativeHybridTopologyProtocol.default_settings()
+    small_ff = settings.forcefield_settings.small_molecule_forcefield
+    if ".offxml" not in small_ff:
+        small_ff += ".offxml"
+    ff = ForceField(small_ff)
+
+    chloro_openff = chloroethane.to_openff()
+    chloro_charges = chloro_openff.partial_charges.m_as(offunit.elementary_charge)
+    chloro_labels = ff.label_molecules(chloro_openff.to_topology())[0]
+    ethane_openff = ethane.to_openff()
+    ethane_charges = ethane_openff.partial_charges.m_as(offunit.elementary_charge)
+    ethane_labels = ff.label_molecules(ethane_openff.to_topology())[0]
+    htf = make_htf(mapping=chloroethane_to_ethane_mapping, settings=settings)
+    hybrid_system = htf.hybrid_system
+    forces = {force.getName(): force for force in hybrid_system.getForces()}
+
+    return {
+        "htf": htf,
+        "hybrid_system": hybrid_system,
+        "forces": forces,
+        "chloro_labels": chloro_labels,
+        "ethane_labels": ethane_labels,
+        "mapping": chloroethane_to_ethane_mapping,
+        "chloroethane": chloroethane,
+        "ethane": ethane,
+        "chloro_charges": chloro_charges,
+        "ethane_charges": ethane_charges,
+        "electrostatic_scale": ff.get_parameter_handler("Electrostatics").scale14,
+        "vdW_scale": ff.get_parameter_handler("vdW").scale14,
+        "force_field": ff,
+    }
+
+
+@pytest.fixture(scope="module")
+def htf_chlorobenzene_fluorobenzene(
+    chlorobenzene, fluorobenzene, chlorobenzene_to_fluorobenzene_mapping, t4_lysozyme_solvated
+):
+    """Generate the htf for chlorobenzene to fluorobenzene in complex with t4-lysozyme solvated."""
+    settings = RelativeHybridTopologyProtocol.default_settings()
+    small_ff = settings.forcefield_settings.small_molecule_forcefield
+    if ".offxml" not in small_ff:
+        small_ff += ".offxml"
+    ff = ForceField(small_ff)
+    chloro_openff = chlorobenzene.to_openff()
+    chloro_charges = chloro_openff.partial_charges.m_as(offunit.elementary_charge)
+    chloro_labels = ff.label_molecules(chloro_openff.to_topology())[0]
+    fluoro_openff = fluorobenzene.to_openff()
+    fluoro_charges = fluoro_openff.partial_charges.m_as(offunit.elementary_charge)
+    fluoro_labels = ff.label_molecules(fluoro_openff.to_topology())[0]
+    htf = make_htf(
+        mapping=chlorobenzene_to_fluorobenzene_mapping,
+        settings=settings,
+        protein=t4_lysozyme_solvated,
+    )
+    hybrid_system = htf.hybrid_system
+
+    apply_box_vectors_and_fix_nb_force(hybrid_topology_factory=htf, force_field=ff)
+
+    forces = {force.getName(): force for force in hybrid_system.getForces()}
+
+    return {
+        "htf": htf,
+        "hybrid_system": hybrid_system,
+        "forces": forces,
+        "chloro_labels": chloro_labels,
+        "fluoro_labels": fluoro_labels,
+        "mapping": chlorobenzene_to_fluorobenzene_mapping,
+        "chlorobenzene": chlorobenzene,
+        "fluorobenzene": fluorobenzene,
+        "chloro_charges": chloro_charges,
+        "fluoro_charges": fluoro_charges,
+        "electrostatic_scale": ff.get_parameter_handler("Electrostatics").scale14,
+        "vdW_scale": ff.get_parameter_handler("vdW").scale14,
+        "force_field": ff,
+    }
+
+
+@pytest.fixture(scope="module")
+def htf_chlorobenzene_benzene(
+    chlorobenzene, benzene, chlorobenzene_to_benzene_mapping, t4_lysozyme_solvated
+):
+    """Generate the htf for chlorobenzene to benzene with interpolate 1-4s on in complex with t4-lysozyme solvated."""
+    settings = RelativeHybridTopologyProtocol.default_settings()
+    small_ff = settings.forcefield_settings.small_molecule_forcefield
+    if ".offxml" not in small_ff:
+        small_ff += ".offxml"
+    ff = ForceField(small_ff)
+
+    chloro_openff = chlorobenzene.to_openff()
+    chloro_charges = chloro_openff.partial_charges.m_as(offunit.elementary_charge)
+    chloro_labels = ff.label_molecules(chloro_openff.to_topology())[0]
+    benzene_openff = benzene.to_openff()
+    benzene_charges = benzene_openff.partial_charges.m_as(offunit.elementary_charge)
+    benzene_labels = ff.label_molecules(benzene_openff.to_topology())[0]
+    htf = make_htf(
+        mapping=chlorobenzene_to_benzene_mapping, settings=settings, protein=t4_lysozyme_solvated
+    )
+    hybrid_system = htf.hybrid_system
+
+    apply_box_vectors_and_fix_nb_force(hybrid_topology_factory=htf, force_field=ff)
+
+    forces = {force.getName(): force for force in hybrid_system.getForces()}
+
+    return {
+        "htf": htf,
+        "hybrid_system": hybrid_system,
+        "forces": forces,
+        "chloro_labels": chloro_labels,
+        "benzene_labels": benzene_labels,
+        "mapping": chlorobenzene_to_benzene_mapping,
+        "chlorobenzene": chlorobenzene,
+        "benzene": benzene,
+        "chloro_charges": chloro_charges,
+        "benzene_charges": benzene_charges,
+        "electrostatic_scale": ff.get_parameter_handler("Electrostatics").scale14,
+        "vdW_scale": ff.get_parameter_handler("vdW").scale14,
+        "force_field": ff,
+    }

src/openfe/tests/data/htf/chloroethane.sdf

@@ -0,0 +1,24 @@
+chloroethane
+     RDKit          3D
+
+  8  7  0  0  0  0  0  0  0  0999 V2000
+   -1.3775    1.4853   -0.3054 Cl  0  0  0  0  0  0  0  0  0  0  0  0
+   -0.7289   -0.1486    0.0141 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.7801   -0.0222    0.0530 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.0446   -0.7625   -0.8423 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.1589   -0.4787    0.9822 H   0  0  0  0  0  0  0  0  0  0  0  0
+    1.1368    0.9399   -0.3538 H   0  0  0  0  0  0  0  0  0  0  0  0
+    1.2037   -0.8083   -0.6046 H   0  0  0  0  0  0  0  0  0  0  0  0
+    1.1892   -0.2049    1.0569 H   0  0  0  0  0  0  0  0  0  0  0  0
+  1  2  1  0
+  2  3  1  0
+  2  4  1  0
+  2  5  1  0
+  3  6  1  0
+  3  7  1  0
+  3  8  1  0
+M  END
+>  <atom.dprop.PartialCharge>  (1)
+-0.20252500000000001 0.038874999999999993 -0.107225 0.062575000000000006 0.062575000000000006 0.048574999999999993 0.048574999999999993 0.048574999999999993
+
+$$$$

src/openfe/tests/data/htf/ethane.sdf

@@ -0,0 +1,24 @@
+ethane
+     RDKit          3D
+
+  8  7  0  0  0  0  0  0  0  0999 V2000
+   -1.1246    0.8482   -0.1808 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.7289   -0.1486    0.0141 C   0  0  0  0  0  0  0  0  0  0  0  0
+    0.7801   -0.0222    0.0530 C   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.0446   -0.7625   -0.8423 H   0  0  0  0  0  0  0  0  0  0  0  0
+   -1.1589   -0.4787    0.9822 H   0  0  0  0  0  0  0  0  0  0  0  0
+    1.1368    0.9399   -0.3538 H   0  0  0  0  0  0  0  0  0  0  0  0
+    1.2037   -0.8083   -0.6046 H   0  0  0  0  0  0  0  0  0  0  0  0
+    1.1892   -0.2049    1.0569 H   0  0  0  0  0  0  0  0  0  0  0  0
+  1  2  1  0
+  2  3  1  0
+  2  4  1  0
+  2  5  1  0
+  3  6  1  0
+  3  7  1  0
+  3  8  1  0
+M  END
+>  <atom.dprop.PartialCharge>  (1)
+0.031449999999999999 -0.094350000000000003 -0.094350000000000003 0.031449999999999999 0.031449999999999999 0.031449999999999999 0.031449999999999999 0.031449999999999999
+
+$$$$