update

python/BioSimSpace/Sandpit/Exscientia/Protocol/_config.py

@@ -565,34 +565,35 @@ def generateGromacsConfig(
                 [
                     mol,
                 ] = self.system.getDecoupledMolecules()
-                decouple_dict = mol._sire_object.property("decouple")
-                protocol_dict["couple-moltype"] = mol._sire_object.name().value()
-
-                def tranform(charge, LJ):
-                    if charge and LJ:
-                        return "vdw-q"
-                    elif charge and not LJ:
-                        return "q"
-                    elif not charge and LJ:
-                        return "vdw"
-                    else:
-                        return "none"
+                if not mol.isPerturbable():
+                    decouple_dict = mol._sire_object.property("decouple")
+                    protocol_dict["couple-moltype"] = mol._sire_object.name().value()
+
+                    def tranform(charge, LJ):
+                        if charge and LJ:
+                            return "vdw-q"
+                        elif charge and not LJ:
+                            return "q"
+                        elif not charge and LJ:
+                            return "vdw"
+                        else:
+                            return "none"
 
-                protocol_dict["couple-lambda0"] = tranform(
-                    decouple_dict["charge"][0], decouple_dict["LJ"][0]
-                )
-                protocol_dict["couple-lambda1"] = tranform(
-                    decouple_dict["charge"][1], decouple_dict["LJ"][1]
-                )
+                    protocol_dict["couple-lambda0"] = tranform(
+                        decouple_dict["charge"][0], decouple_dict["LJ"][0]
+                    )
+                    protocol_dict["couple-lambda1"] = tranform(
+                        decouple_dict["charge"][1], decouple_dict["LJ"][1]
+                    )
+                    if decouple_dict["intramol"].value():
+                        # The intramol is being coupled to the lambda change and thus being annihilated.
+                        protocol_dict["couple-intramol"] = "yes"
+                    else:
+                        protocol_dict["couple-intramol"] = "no"
                 # Add the soft-core parameters for the ABFE
                 protocol_dict["sc-alpha"] = 0.5
                 protocol_dict["sc-power"] = 1
                 protocol_dict["sc-sigma"] = 0.3
-                if decouple_dict["intramol"].value():
-                    # The intramol is being coupled to the lambda change and thus being annihilated.
-                    protocol_dict["couple-intramol"] = "yes"
-                else:
-                    protocol_dict["couple-intramol"] = "no"
             elif nDecoupledMolecules > 1:
                 raise ValueError(
                     "Gromacs cannot handle more than one decoupled molecule."

tests/Sandpit/Exscientia/Protocol/test_config.py

@@ -18,6 +18,7 @@
 from BioSimSpace.Sandpit.Exscientia.Units.Energy import kcal_per_mol
 from BioSimSpace.Sandpit.Exscientia.Units.Temperature import kelvin
 from BioSimSpace.Sandpit.Exscientia.FreeEnergy import Restraint
+from BioSimSpace.Sandpit.Exscientia._SireWrappers import Molecule
 from BioSimSpace.Sandpit.Exscientia._Utils import _try_import, _have_imported
 
 
@@ -301,6 +302,49 @@ def test_decouple_vdw_q(self, system):
             assert "couple-lambda1 = none" in mdp_text
             assert "couple-intramol = yes" in mdp_text
 
+
+    def test_decouple_perturbable(self, system):
+        m, protocol = system
+        mol = decouple(m)
+        sire_mol = mol._sire_object
+        c = sire_mol.cursor()
+        for key in [
+            "charge",
+            "LJ",
+            "bond",
+            "angle",
+            "dihedral",
+            "improper",
+            "forcefield",
+            "intrascale",
+            "mass",
+            "element",
+            "atomtype",
+            "coordinates",
+            "velocity",
+            "ambertype"
+        ]:
+            if f"{key}1" not in c and key in c:
+                c[f"{key}0"] = c[key]
+                c[f"{key}1"] = c[key]
+
+        c["is_perturbable"] = True
+        sire_mol = c.commit()
+        mol = Molecule(sire_mol)
+
+        freenrg = BSS.FreeEnergy.AlchemicalFreeEnergy(
+            mol.toSystem(),
+            protocol,
+            engine="GROMACS",
+        )
+        with open(f"{freenrg._work_dir}/lambda_6/gromacs.mdp", "r") as f:
+            mdp_text = f.read()
+            assert "couple-moltype = LIG" not in mdp_text
+            assert "couple-lambda0 = vdw-q" not in mdp_text
+            assert "couple-lambda1 = none" not in mdp_text
+            assert "couple-intramol = yes" not in mdp_text
+
+
     @pytest.mark.skipif(
         has_gromacs is False, reason="Requires GROMACS to be installed."
     )