Merge pull request #61 from Exscientia/feat__sync

Sync fix from OBS

Author: xiki-tempula

README.rst

@@ -1,8 +1,8 @@
 `BioSimSpace <https://biosimspace.openbiosim.org>`__
 ====================================================
 
-.. image:: https://github.com/openbiosim/biosimspace/workflows/Build/badge.svg
-   :target: https://github.com/openbiosim/biosimspace/actions?query=workflow%3ABuild)
+.. image:: https://github.com/openbiosim/biosimspace/actions/workflows/devel.yaml/badge.svg
+   :target: https://github.com/openbiosim/biosimspace/actions?query=workflow%3ARelease-Devel
    :alt: Build status
 
 .. image:: https://anaconda.org/openbiosim/biosimspace/badges/downloads.svg

doc/source/tutorials/alchemical_transfer.rst

@@ -88,8 +88,9 @@ more choices to be made before the system can be prepared. The most important of
 these is the choice of displacement vector, which defines the direction and
 distance at which the free ligand will be placed relative to the bound ligand.
 It is generally recommended that this displacement be at least 3 layers of water
-molecules (> 10 Å) thick. If no displacement is provided a default choice of
-[20Å, 20Å, 20Å] will be used.
+molecules (> 10 Å) thick. If no displacement is provided BioSimSpace will find a
+best guess translation vector based on the relative positions of the protein and
+bound ligands and translate the free ligand 20 Å along this vector.
 
 This is also the point at which a custom set of atoms can be chosen to define the
 centre of mass of both the ligands and the protein. In the majority of cases it

doc/source/tutorials/images/tyk2_prepared.png

@@ -1436,6 +1436,13 @@ def _analyse_amber(work_dir=None, estimator="MBAR", method="alchemlyb", **kwargs
                 if "lambda" in part:
                     lambdas.append(float(part.split("_")[-1]))
 
+        if len(lambdas) == 0:
+            raise ValueError(
+                "No lambda windows were detected from the output directory names! "
+                "Ensure that the directory being analysed contains sub-directories"
+                "name e.g `lambda_0`, `lambda_1` containing the output files."
+            )
+
         # Find the temperature for each lambda window.
         temperatures = []
         for file, lambda_ in zip(files, lambdas):
@@ -1517,6 +1524,13 @@ def _analyse_gromacs(work_dir=None, estimator="MBAR", method="alchemlyb", **kwar
                     if "lambda" in part:
                         lambdas.append(float(part.split("_")[-1]))
 
+            if len(lambdas) == 0:
+                raise ValueError(
+                    "No lambda windows were detected from the output directory names! "
+                    "Ensure that the directory being analysed contains sub-directories"
+                    "name e.g `lambda_0`, `lambda_1` containing the output files."
+                )
+
             # Find the temperature at each lambda window
             temperatures = []
             for file in files:
@@ -1686,6 +1700,13 @@ def _analyse_somd(work_dir=None, estimator="MBAR", method="alchemlyb", **kwargs)
                     if "lambda" in part:
                         lambdas.append(float(part.split("_")[-1]))
 
+            if len(lambdas) == 0:
+                raise ValueError(
+                    "No lambda windows were detected from the output directory names! "
+                    "Ensure that the directory being analysed contains sub-directories"
+                    "name e.g `lambda_0`, `lambda_1` containing the output files."
+                )
+
             temperatures = []
             for file in files:
                 found_temperature = False
@@ -1883,6 +1904,13 @@ def _analyse_somd2(work_dir=None, estimator="MBAR", method="alchemlyb", **kwargs
         temperatures = [x for _, x in sorted(zip(lambdas, temperatures))]
         lambdas = sorted(lambdas)
 
+        if len(lambdas) == 0:
+            raise ValueError(
+                "No lambda windows were detected from the output directory names! "
+                "Ensure that the directory being analysed contains sub-directories"
+                "name e.g `lambda_0`, `lambda_1` containing the output files."
+            )
+
         # Check that the temperatures at the end states match.
         if temperatures[0] != temperatures[-1]:
             raise ValueError("The temperatures at the endstates don't match!")

python/BioSimSpace/FreeEnergy/_relative.py

@@ -841,6 +841,7 @@ def __init__(
         net_charge=None,
         ensure_compatible=True,
         property_map={},
+        **kwargs,
     ):
         """
         Constructor.
@@ -873,6 +874,10 @@ def __init__(
             A dictionary that maps system "properties" to their user defined
             values. This allows the user to refer to properties with their
             own naming scheme, e.g. { "charge" : "my-charge" }
+
+        **kwargs: dict
+            Additional keyword arguments. These can be used to pass custom
+            parameters to the Antechamber program.
         """
 
         if type(version) is not int:
@@ -912,6 +917,11 @@ def __init__(
                     "'net_charge' must be of type 'int', or `BioSimSpace.Types.Charge'"
                 )
 
+        # Check the kwargs to see whether acdoctor is enabled.
+        self._acdoctor = kwargs.get("acdoctor", True)
+        if not isinstance(self._acdoctor, bool):
+            raise TypeError("'acdoctor' must be of type 'bool'")
+
         # Set the version.
         self._version = version
 
@@ -1086,6 +1096,10 @@ def run(self, molecule, work_dir=None, queue=None):
             charge,
         )
 
+        # Disable acdoctor if requested.
+        if not self._acdoctor:
+            command += " -dr no"
+
         with open(_os.path.join(str(work_dir), "README.txt"), "w") as file:
             # Write the command to file.
             file.write("# Antechamber was run with the following command:\n")

python/BioSimSpace/Parameters/_Protocol/_amber.py

@@ -362,6 +362,7 @@ def gaff(
         charge_method=charge_method,
         ensure_compatible=ensure_compatible,
         property_map=property_map,
+        **kwargs,
     )
 
     # Run the parameterisation protocol in the background and return
@@ -460,6 +461,7 @@ def gaff2(
         charge_method=charge_method,
         ensure_compatible=ensure_compatible,
         property_map=property_map,
+        **kwargs,
     )
 
     # Run the parameterisation protocol in the background and return

python/BioSimSpace/Parameters/_parameters.py

@@ -590,6 +590,14 @@ def analyse(
                 f"traj {type(traj)} must be of type 'BioSimSpace.Trajectory._trajectory.Trajectory'"
             )
 
+        n_frames = traj.nFrames()
+        if not n_frames >= 50:
+            _warnings.warn(
+                f"The trajectory for restraint selection has less than 50 frames ({n_frames} frames). "
+                "This may result in poor restraint selection with excessively high force constants. "
+                "Consider running a longer simulation or saving more frames."
+            )
+
         if not isinstance(temperature, _Temperature):
             raise ValueError(
                 f"temperature {type(temperature)} must be of type 'BioSimSpace.Types.Temperature'"
@@ -667,6 +675,18 @@ def analyse(
         decoupled_mol = _system.getDecoupledMolecules()[0]
         decoupled_resname = decoupled_mol.getResidues()[0].name()
 
+        # Warn the user if the decoupled molecule is water or a macromoelcule
+        if decoupled_mol.isWater():
+            _warnings.warn(
+                "The decoupled molecule is water. Ensure that this is intended. Using constrained water hydrogens as anchor "
+                "points may produce instabilities."
+            )
+
+        if decoupled_mol.nResidues() > 1:
+            _warnings.warn(
+                "The decoupled molecule has more than one residue and is likely a macromolecule. Ensure that this is intended."
+            )
+
         ligand_selection_str = f"((resname {decoupled_resname}) and (not name H*))"
         if append_to_ligand_selection:
             ligand_selection_str += " and "
@@ -1012,6 +1032,20 @@ def _findOrderedPairs(u, ligand_selection_str, receptor_selection_str, cutoff):
         """
 
         lig_selection = u.select_atoms(ligand_selection_str)
+
+        # Ensure that there are no shared atoms between the ligand selection and all possible receptor selections.
+        all_possible_receptor_selection = u.select_atoms(receptor_selection_str)
+        shared_atoms = [
+            atom for atom in lig_selection if atom in all_possible_receptor_selection
+        ]
+        if shared_atoms:
+            raise _AnalysisError(
+                "Shared atoms between ligand and receptor selections detected. "
+                "Please ensure that you are decoupling the intended molecule and that "
+                "the ligand and receptor selections are mutually exclusive. "
+                f"Shared atoms: {shared_atoms}"
+            )
+
         pair_variance_dict = {}
 
         # Get all receptor atoms within specified distance of cutoff
@@ -1488,7 +1522,9 @@ def _findOrderedBoresch(
                             dtheta = abs(val - circmean)
                             corrected_values.append(min(dtheta, 2 * _np.pi - dtheta))
                         corrected_values = _np.array(corrected_values)
-                        boresch_dof_data[pair][dof]["var"] = corrected_values.var()
+                        boresch_dof_data[pair][dof]["var"] = _np.mean(
+                            corrected_values**2
+                        )
 
                     # Assume Gaussian distributions and calculate force constants for harmonic potentials
                     # so as to reproduce these distributions at 298 K

python/BioSimSpace/Sandpit/Exscientia/FreeEnergy/_restraint_search.py

@@ -841,6 +841,7 @@ def __init__(
         net_charge=None,
         ensure_compatible=True,
         property_map={},
+        **kwargs,
     ):
         """
         Constructor.
@@ -873,6 +874,10 @@ def __init__(
             A dictionary that maps system "properties" to their user defined
             values. This allows the user to refer to properties with their
             own naming scheme, e.g. { "charge" : "my-charge" }
+
+        **kwargs: dict
+            Additional keyword arguments. These can be used to pass custom
+            parameters to the Antechamber program.
         """
 
         if type(version) is not int:
@@ -912,6 +917,11 @@ def __init__(
                     "'net_charge' must be of type 'int', or `BioSimSpace.Types.Charge'"
                 )
 
+        # Check the kwargs to see whether acdoctor is enabled.
+        self._acdoctor = kwargs.get("acdoctor", True)
+        if not isinstance(self._acdoctor, bool):
+            raise TypeError("'acdoctor' must be of type 'bool'")
+
         # Set the version.
         self._version = version
 
@@ -1086,6 +1096,10 @@ def run(self, molecule, work_dir=None, queue=None):
             charge,
         )
 
+        # Disable acdoctor if requested.
+        if not self._acdoctor:
+            command += " -dr no"
+
         with open(_os.path.join(str(work_dir), "README.txt"), "w") as file:
             # Write the command to file.
             file.write("# Antechamber was run with the following command:\n")

python/BioSimSpace/Sandpit/Exscientia/Parameters/_Protocol/_amber.py

@@ -362,6 +362,7 @@ def gaff(
         charge_method=charge_method,
         ensure_compatible=ensure_compatible,
         property_map=property_map,
+        **kwargs,
     )
 
     # Run the parameterisation protocol in the background and return
@@ -460,6 +461,7 @@ def gaff2(
         charge_method=charge_method,
         ensure_compatible=ensure_compatible,
         property_map=property_map,
+        **kwargs,
     )
 
     # Run the parameterisation protocol in the background and return

python/BioSimSpace/Sandpit/Exscientia/Parameters/_parameters.py

@@ -167,3 +167,24 @@ def test_smiles_stereo():
 
     # Make sure the SMILES strings are the same.
     assert rdmol0_smiles == rdmol1_smiles
+
+
+@pytest.mark.skipif(
+    has_antechamber is False or has_tleap is False,
+    reason="Requires AmberTools/antechamber and tLEaP to be installed.",
+)
+def test_acdoctor():
+    """
+    Test that parameterising negatively charged molecules works when acdoctor
+    is disabled.
+    """
+
+    # Load the molecule.
+    mol = BSS.IO.readMolecules(f"{url}/negative_charge.sdf")[0]
+
+    # Make sure parameterisation fails when acdoctor is enabled.
+    with pytest.raises(BSS._Exceptions.ParameterisationError):
+        BSS.Parameters.gaff(mol).getMolecule()
+
+    # Make sure parameterisation works when acdoctor is disabled.
+    mol = BSS.Parameters.gaff(mol, acdoctor=False).getMolecule()