Merge pull request #301 from OpenBioSim/backport_295_298_300

Backport fixes from #295, #298, and #300

python/BioSimSpace/Align/_merge.py

@@ -693,6 +693,34 @@ def merge(
                     .molecule()
                 )
 
+    # Tolerance for zero sigma values.
+    null_lj_sigma = 1e-9
+
+    # Atoms with zero LJ sigma values need to have their sigma values set to the
+    # value from the other end state.
+    for atom in edit_mol.atoms():
+        # Get the end state LJ sigma values.
+        lj0 = atom.property("LJ0")
+        lj1 = atom.property("LJ1")
+
+        # Lambda = 0 state has a zero sigma value.
+        if abs(lj0.sigma().value()) <= null_lj_sigma:
+            # Use the sigma value from the lambda = 1 state.
+            edit_mol = (
+                edit_mol.atom(atom.index())
+                .set_property("LJ0", _SireMM.LJParameter(lj1.sigma(), lj0.epsilon()))
+                .molecule()
+            )
+
+        # Lambda = 1 state has a zero sigma value.
+        if abs(lj1.sigma().value()) <= null_lj_sigma:
+            # Use the sigma value from the lambda = 0 state.
+            edit_mol = (
+                edit_mol.atom(atom.index())
+                .set_property("LJ1", _SireMM.LJParameter(lj0.sigma(), lj1.epsilon()))
+                .molecule()
+            )
+
     # We now need to merge "bond", "angle", "dihedral", and "improper" parameters.
     # To do so, we extract the properties from molecule1, then add the additional
     # properties from molecule0, making sure to update the atom indices, and bond

python/BioSimSpace/IO/_io.py

@@ -764,7 +764,7 @@ def saveMolecules(
             continue
 
         # Add the file format to the property map.
-        _property_map["fileformat"] = _SireBase.wrap(format)
+        _property_map["fileformat"] = format
 
         # Warn the user if any molecules are parameterised with a force field
         # that uses geometric combining rules. While we can write this to file

python/BioSimSpace/Parameters/_parameters.py

@@ -111,8 +111,10 @@ def parameterise(
     Returns
     -------
 
-    molecule : :class:`Molecule <BioSimSpace._SireWrappers.Molecule>`
-        The parameterised molecule.
+    process : :class:`Process <BioSimSpace.Parameters._process.Process>`
+        A process to parameterise the molecule in the background. Call the
+        .getMolecule() method on the returned process to block until the
+        parameterisation is complete and get the parameterised molecule.
     """
 
     if not isinstance(forcefield, str):
@@ -208,8 +210,10 @@ def _parameterise_amber_protein(
     Returns
     -------
 
-    molecule : :class:`Molecule <BioSimSpace._SireWrappers.Molecule>`
-        The parameterised molecule.
+    process : :class:`Process <BioSimSpace.Parameters._process.Process>`
+        A process to parameterise the molecule in the background. Call the
+        .getMolecule() method on the returned process to block until the
+        parameterisation is complete and get the parameterised molecule.
     """
 
     if not isinstance(forcefield, str):
@@ -314,8 +318,10 @@ def gaff(
     Returns
     -------
 
-    molecule : :class:`Molecule <BioSimSpace._SireWrappers.Molecule>`
-        The parameterised molecule.
+    process : :class:`Process <BioSimSpace.Parameters._process.Process>`
+        A process to parameterise the molecule in the background. Call the
+        .getMolecule() method on the returned process to block until the
+        parameterisation is complete and get the parameterised molecule.
     """
 
     if _amber_home is None:
@@ -407,8 +413,10 @@ def gaff2(
     Returns
     -------
 
-    molecule : :class:`Molecule <BioSimSpace._SireWrappers.Molecule>`
-        The parameterised molecule.
+    process : :class:`Process <BioSimSpace.Parameters._process.Process>`
+        A process to parameterise the molecule in the background. Call the
+        .getMolecule() method on the returned process to block until the
+        parameterisation is complete and get the parameterised molecule.
     """
 
     if _amber_home is None:
@@ -506,8 +514,10 @@ def _parameterise_openff(
     Returns
     -------
 
-    molecule : :class:`Molecule <BioSimSpace._SireWrappers.Molecule>`
-        The parameterised molecule.
+    process : :class:`Process <BioSimSpace.Parameters._process.Process>`
+        A process to parameterise the molecule in the background. Call the
+        .getMolecule() method on the returned process to block until the
+        parameterisation is complete and get the parameterised molecule.
     """
 
     from sire.legacy.Base import findExe as _findExe
@@ -1068,8 +1078,10 @@ def _function(
         Returns
         -------
 
-        molecule : :class:`Molecule <BioSimSpace._SireWrappers.Molecule>`
-            The parameterised molecule.
+        process : :class:`Process <BioSimSpace.Parameters._process.Process>`
+            A process to parameterise the molecule in the background. Call the
+            .getMolecule() method on the returned process to block until the
+            parameterisation is complete and get the parameterised molecule.
         """
         return _parameterise_amber_protein(
             name,
@@ -1133,8 +1145,10 @@ def _function(
         Returns
         -------
 
-        molecule : :class:`Molecule <BioSimSpace._SireWrappers.Molecule>`
-            The parameterised molecule.
+        process : :class:`Process <BioSimSpace.Parameters._process.Process>`
+            A process to parameterise the molecule in the background. Call the
+            .getMolecule() method on the returned process to block until the
+            parameterisation is complete and get the parameterised molecule.
         """
         return _parameterise_openff(
             name,

python/BioSimSpace/Process/_somd.py

@@ -1129,6 +1129,10 @@ def _to_pert_file(
 
         # 1) Atoms.
 
+        # Store a null LJParameter object for dummy atoms. SOMD requires that
+        # both sigma and epsilon are set to zero for dummy atoms.
+        dummy_lj = _SireMM.LJParameter()
+
         def atom_sorting_criteria(atom):
             LJ0 = atom.property("LJ0")
             LJ1 = atom.property("LJ1")
@@ -1156,6 +1160,13 @@ def atom_sorting_criteria(atom):
                     LJ0 = atom.property("LJ0")
                     LJ1 = atom.property("LJ1")
 
+                    # Dummy atom at lambda = 0.
+                    if _is_dummy(mol, atom.index(), is_lambda1=False):
+                        LJ0 = dummy_lj
+                    # Dummy atom at lambda = 1.
+                    if _is_dummy(mol, atom.index(), is_lambda1=True):
+                        LJ1 = dummy_lj
+
                     # Atom data.
                     file.write("        name           %s\n" % atom.name().value())
                     file.write(
@@ -1199,6 +1210,13 @@ def atom_sorting_criteria(atom):
                     LJ0 = atom.property("LJ0")
                     LJ1 = atom.property("LJ1")
 
+                    # Dummy atom at lambda = 0.
+                    if _is_dummy(mol, idx, is_lambda1=False):
+                        LJ0 = dummy_lj
+                    # Dummy atom at lambda = 1.
+                    if _is_dummy(mol, idx, is_lambda1=True):
+                        LJ1 = dummy_lj
+
                     # Atom data.
                     file.write("        name           %s\n" % atom.name().value())
                     file.write(
@@ -1252,11 +1270,11 @@ def atom_sorting_criteria(atom):
 
                 # Set LJ/charge based on requested perturbed term.
                 if perturbation_type == "discharge_soft":
-                    if atom.property("element0") == _SireMol.Element(
-                        "X"
-                    ) or atom.property("element1") == _SireMol.Element("X"):
+                    if _is_dummy(mol, idx, is_lambda1=False) or _is_dummy(
+                        mol, idx, is_lambda1=True
+                    ):
                         # If perturbing TO dummy:
-                        if atom.property("element1") == _SireMol.Element("X"):
+                        if is_dummy(mol, idx, is_lambda1=True):
                             atom_type1 = atom_type0
 
                             # In this step, only remove charges from soft-core perturbations.
@@ -1270,7 +1288,7 @@ def atom_sorting_criteria(atom):
 
                         # If perturbing FROM dummy:
                         else:
-                            # All terms have already been perturbed in "5_grow_soft".
+                            # All terms have already been perturbed in "grow_soft".
                             atom_type1 = atom_type0
                             LJ0_value = LJ1_value = (
                                 LJ0.sigma().value(),
@@ -1290,11 +1308,11 @@ def atom_sorting_criteria(atom):
                         charge0_value = charge1_value = atom.property("charge0").value()
 
                 elif perturbation_type == "vanish_soft":
-                    if atom.property("element0") == _SireMol.Element(
-                        "X"
-                    ) or atom.property("element1") == _SireMol.Element("X"):
+                    if _is_dummy(mol, idx, is_lambda1=False) or _is_dummy(
+                        mol, idx, is_lambda1=True
+                    ):
                         # If perturbing TO dummy:
-                        if atom.property("element1") == _SireMol.Element("X"):
+                        if _is_dummy(mol, idx, is_lambda1=True):
                             # allow atom types to change.
                             atom_type0 = atom_type0
                             atom_type1 = atom_type1
@@ -1308,7 +1326,7 @@ def atom_sorting_criteria(atom):
 
                         # If perturbing FROM dummy:
                         else:
-                            # All terms have already been perturbed in "5_grow_soft".
+                            # All terms have already been perturbed in "grow_soft".
                             atom_type1 = atom_type0
                             LJ0_value = LJ1_value = (
                                 LJ0.sigma().value(),
@@ -1328,11 +1346,11 @@ def atom_sorting_criteria(atom):
                         charge0_value = charge1_value = atom.property("charge0").value()
 
                 elif perturbation_type == "flip":
-                    if atom.property("element0") == _SireMol.Element(
-                        "X"
-                    ) or atom.property("element1") == _SireMol.Element("X"):
+                    if _is_dummy(mol, idx, is_lambda1=False) or _is_dummy(
+                        mol, idx, is_lambda1=True
+                    ):
                         # If perturbing TO dummy:
-                        if atom.property("element1") == _SireMol.Element("X"):
+                        if _is_dummy(mol, idx, is_lambda1=True):
                             # atom types have already been changed.
                             atom_type0 = atom_type1
 
@@ -1342,7 +1360,7 @@ def atom_sorting_criteria(atom):
 
                         # If perturbing FROM dummy:
                         else:
-                            # All terms have already been perturbed in "5_grow_soft".
+                            # All terms have already been perturbed in "grow_soft".
                             atom_type1 = atom_type0
                             LJ0_value = LJ1_value = (
                                 LJ0.sigma().value(),
@@ -1362,11 +1380,11 @@ def atom_sorting_criteria(atom):
                         charge1_value = atom.property("charge1").value()
 
                 elif perturbation_type == "grow_soft":
-                    if atom.property("element0") == _SireMol.Element(
-                        "X"
-                    ) or atom.property("element1") == _SireMol.Element("X"):
+                    if _is_dummy(mol, idx, is_lambda1=False) or _is_dummy(
+                        mol, idx, is_lambda1=True
+                    ):
                         # If perturbing TO dummy:
-                        if atom.property("element1") == _SireMol.Element("X"):
+                        if _is_dummy(mol, idx, is_lambda1=True):
                             # atom types have already been changed.
                             atom_type0 = atom_type1
 
@@ -1398,11 +1416,11 @@ def atom_sorting_criteria(atom):
                         charge0_value = charge1_value = atom.property("charge1").value()
 
                 elif perturbation_type == "charge_soft":
-                    if atom.property("element0") == _SireMol.Element(
-                        "X"
-                    ) or atom.property("element1") == _SireMol.Element("X"):
+                    if _is_dummy(mol, idx, is_lambda1=False) or _is_dummy(
+                        mol, idx, is_lambda1=True
+                    ):
                         # If perturbing TO dummy:
-                        if atom.property("element1") == _SireMol.Element("X"):
+                        if _is_dummy(mol, idx, is_lambda1=True):
                             # atom types have already been changed.
                             atom_type0 = atom_type1
 
@@ -2971,11 +2989,13 @@ def _has_dummy(mol, idxs, is_lambda1=False):
     element_dummy = _SireMol.Element(0)
     ambertype_dummy = "du"
 
+    # Check that the molecule has the ambertype property.
+    has_ambertype = mol.hasProperty(ambertype_prop)
+
     # Check whether an of the atoms is a dummy.
     for idx in idxs:
-        if (
-            mol.atom(idx).property(element_prop) == element_dummy
-            or mol.atom(idx).property(ambertype_prop) == ambertype_dummy
+        if mol.atom(idx).property(element_prop) == element_dummy or (
+            has_ambertype and mol.atom(idx).property(ambertype_prop) == ambertype_dummy
         ):
             return True
 
@@ -2985,6 +3005,7 @@ def _has_dummy(mol, idxs, is_lambda1=False):
 def _is_dummy(mol, idxs, is_lambda1=False):
     """
     Internal function to return whether each atom is a dummy.
+    If a single index is pased, then a single boolean is returned.
 
     Parameters
     ----------
@@ -3027,17 +3048,30 @@ def _is_dummy(mol, idxs, is_lambda1=False):
     element_dummy = _SireMol.Element(0)
     ambertype_dummy = "du"
 
+    # Check that the molecule has the ambertype property.
+    has_ambertype = mol.hasProperty(ambertype_prop)
+
     # Initialise a list to store the state of each atom.
     is_dummy = []
 
+    # Convert single index to list.
+    if not isinstance(idxs, list):
+        idxs = [idxs]
+
     # Check whether each of the atoms is a dummy.
     for idx in idxs:
         is_dummy.append(
             mol.atom(idx).property(element_prop) == element_dummy
-            or mol.atom(idx).property(ambertype_prop) == ambertype_dummy
+            or (
+                has_ambertype
+                and mol.atom(idx).property(ambertype_prop) == ambertype_dummy
+            )
         )
 
-    return is_dummy
+    if len(is_dummy) == 1:
+        return is_dummy[0]
+    else:
+        return is_dummy
 
 
 def _random_suffix(basename, size=4, chars=_string.ascii_uppercase + _string.digits):

python/BioSimSpace/Sandpit/Exscientia/Align/_merge.py

@@ -739,6 +739,34 @@ def merge(
                     edit_mol.atom(idx).setProperty(name, atom.property(prop)).molecule()
                 )
 
+    # Tolerance for zero sigma values.
+    null_lj_sigma = 1e-9
+
+    # Atoms with zero LJ sigma values need to have their sigma values set to the
+    # value from the other end state.
+    for atom in edit_mol.atoms():
+        # Get the end state LJ sigma values.
+        lj0 = atom.property("LJ0")
+        lj1 = atom.property("LJ1")
+
+        # Lambda = 0 state has a zero sigma value.
+        if abs(lj0.sigma().value()) <= null_lj_sigma:
+            # Use the sigma value from the lambda = 1 state.
+            edit_mol = (
+                edit_mol.atom(atom.index())
+                .set_property("LJ0", _SireMM.LJParameter(lj1.sigma(), lj0.epsilon()))
+                .molecule()
+            )
+
+        # Lambda = 1 state has a zero sigma value.
+        if abs(lj1.sigma().value()) <= null_lj_sigma:
+            # Use the sigma value from the lambda = 0 state.
+            edit_mol = (
+                edit_mol.atom(atom.index())
+                .set_property("LJ1", _SireMM.LJParameter(lj0.sigma(), lj1.epsilon()))
+                .molecule()
+            )
+
     # We now need to merge "bond", "angle", "dihedral", and "improper" parameters.
     # To do so, we extract the properties from molecule1, then add the additional
     # properties from molecule0, making sure to update the atom indices, and bond

python/BioSimSpace/Sandpit/Exscientia/IO/_io.py

@@ -764,7 +764,7 @@ def saveMolecules(
             continue
 
         # Add the file format to the property map.
-        _property_map["fileformat"] = _SireBase.wrap(format)
+        _property_map["fileformat"] = format
 
         # Warn the user if any molecules are parameterised with a force field
         # that uses geometric combining rules. While we can write this to file