diff --git a/.github/workflows/Sandpit_exs.yml b/.github/workflows/Sandpit_exs.yml index f5b0fc0d3..58243a58a 100644 --- a/.github/workflows/Sandpit_exs.yml +++ b/.github/workflows/Sandpit_exs.yml @@ -18,7 +18,7 @@ jobs: strategy: fail-fast: false matrix: - os: ["ubuntu-latest", "macOS-latest",] + os: ["ubuntu-latest", ] python-version: ["3.10",] steps: @@ -35,7 +35,7 @@ jobs: - name: Install dependency run: | - mamba install -c conda-forge -c openbiosim/label/main biosimspace python=3.10 ambertools gromacs "sire=2023.5" "alchemlyb>=2.1" pytest openff-interchange pint=0.21 rdkit "jaxlib>0.3.7" tqdm + mamba install -c conda-forge -c openbiosim/label/main biosimspace python=3.10 ambertools gromacs "sire=2024.1.0" "alchemlyb>=2.1" pytest openff-interchange pint=0.21 rdkit "jaxlib>0.3.7" tqdm python -m pip install git+https://github.com/Exscientia/MDRestraintsGenerator.git # For the testing of BSS.FreeEnergy.AlchemicalFreeEnergy.analysis python -m pip install https://github.com/alchemistry/alchemtest/archive/master.zip diff --git a/python/BioSimSpace/Sandpit/Exscientia/FreeEnergy/_restraint_search.py b/python/BioSimSpace/Sandpit/Exscientia/FreeEnergy/_restraint_search.py index 31a53e956..2fa2c1b90 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/FreeEnergy/_restraint_search.py +++ b/python/BioSimSpace/Sandpit/Exscientia/FreeEnergy/_restraint_search.py @@ -92,6 +92,9 @@ install_command="pip install MDRestraintsGenerator", ) if _have_imported(_MDRestraintsGenerator): + from Bio import BiopythonDeprecationWarning + + _warnings.simplefilter("ignore", BiopythonDeprecationWarning) from MDRestraintsGenerator import search as _search from MDRestraintsGenerator.restraints import ( FindBoreschRestraint as _FindBoreschRestraint, diff --git a/python/BioSimSpace/Sandpit/Exscientia/Process/_amber.py b/python/BioSimSpace/Sandpit/Exscientia/Process/_amber.py index 64c7a173f..8b557d51f 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/Process/_amber.py +++ b/python/BioSimSpace/Sandpit/Exscientia/Process/_amber.py @@ -145,6 +145,14 @@ def __init__( property_map, ) + if ( + not isinstance(protocol, _Protocol._FreeEnergyMixin) + and system.nPerturbableMolecules() + ): + raise _IncompatibleError( + "Perturbable system is not compatible with none free energy protocol." + ) + # Set the package name. self._package_name = "AMBER" @@ -2936,7 +2944,7 @@ def _init_stdout_dict(self): for file in _Path(self.workDir()).glob("*.out.offset"): os.remove(file) - def saveMetric( + def _saveMetric( self, filename="metric.parquet", u_nk="u_nk.parquet", dHdl="dHdl.parquet" ): """ @@ -2945,41 +2953,54 @@ def saveMetric( is Free Energy protocol, the dHdl and the u_nk data will be saved in the same parquet format as well. """ - - _assert_imported(_alchemlyb) - - self._init_stdout_dict() - datadict = dict() - if isinstance(self._protocol, _Protocol.Minimisation): - datadict_keys = [ - ("Time (ps)", None, "getStep"), - ( - "PotentialEnergy (kJ/mol)", - _Units.Energy.kj_per_mol, - "getTotalEnergy", - ), - ] - else: - datadict_keys = [ - ("Time (ps)", _Units.Time.picosecond, "getTime"), - ( - "PotentialEnergy (kJ/mol)", - _Units.Energy.kj_per_mol, - "getPotentialEnergy", - ), - ("Volume (nm^3)", _Units.Volume.nanometer3, "getVolume"), - ("Pressure (bar)", _Units.Pressure.bar, "getPressure"), - ("Temperature (kelvin)", _Units.Temperature.kelvin, "getTemperature"), - ] - # # Disable this now - df = self._convert_datadict_keys(datadict_keys) - df.to_parquet(path=f"{self.workDir()}/{filename}", index=True) - if isinstance(self._protocol, _Protocol.FreeEnergy): - energy = _extract( - f"{self.workDir()}/{self._name}.out", - T=self._protocol.getTemperature() / _Units.Temperature.kelvin, - ) - if "u_nk" in energy and energy["u_nk"] is not None: - energy["u_nk"].to_parquet(path=f"{self.workDir()}/{u_nk}", index=True) - if "dHdl" in energy and energy["dHdl"] is not None: - energy["dHdl"].to_parquet(path=f"{self.workDir()}/{dHdl}", index=True) + if filename is not None: + self._init_stdout_dict() + if isinstance(self._protocol, _Protocol.Minimisation): + datadict_keys = [ + ("Time (ps)", None, "getStep"), + ( + "PotentialEnergy (kJ/mol)", + _Units.Energy.kj_per_mol, + "getTotalEnergy", + ), + ] + else: + datadict_keys = [ + ("Time (ps)", _Units.Time.picosecond, "getTime"), + ( + "PotentialEnergy (kJ/mol)", + _Units.Energy.kj_per_mol, + "getPotentialEnergy", + ), + ("Volume (nm^3)", _Units.Volume.nanometer3, "getVolume"), + ("Pressure (bar)", _Units.Pressure.bar, "getPressure"), + ( + "Temperature (kelvin)", + _Units.Temperature.kelvin, + "getTemperature", + ), + ] + # # Disable this now + df = self._convert_datadict_keys(datadict_keys) + df.to_parquet(path=f"{self.workDir()}/{filename}", index=True) + if u_nk is not None or dHdl is not None: + _assert_imported(_alchemlyb) + + if isinstance(self._protocol, _Protocol.FreeEnergy): + energy = _extract( + f"{self.workDir()}/{self._name}.out", + T=self._protocol.getTemperature() / _Units.Temperature.kelvin, + ) + with _warnings.catch_warnings(): + _warnings.filterwarnings( + "ignore", + message="The DataFrame has column names of mixed type.", + ) + if "u_nk" in energy and energy["u_nk"] is not None: + energy["u_nk"].to_parquet( + path=f"{self.workDir()}/{u_nk}", index=True + ) + if "dHdl" in energy and energy["dHdl"] is not None: + energy["dHdl"].to_parquet( + path=f"{self.workDir()}/{dHdl}", index=True + ) diff --git a/python/BioSimSpace/Sandpit/Exscientia/Process/_gromacs.py b/python/BioSimSpace/Sandpit/Exscientia/Process/_gromacs.py index 788450c84..b6e02992f 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/Process/_gromacs.py +++ b/python/BioSimSpace/Sandpit/Exscientia/Process/_gromacs.py @@ -54,7 +54,8 @@ _alchemlyb = _try_import("alchemlyb") if _have_imported(_alchemlyb): - from alchemlyb.parsing.gmx import extract as _extract + from alchemlyb.parsing.gmx import extract_u_nk as _extract_u_nk + from alchemlyb.parsing.gmx import extract_dHdl as _extract_dHdl from .. import _gmx_exe from .. import _isVerbose @@ -400,18 +401,19 @@ def _write_system(self, system, coord_file=None, topol_file=None, ref_file=None) property_map=self._property_map, ) - # Write the restraint to the topology file - if self._restraint: - with open(topol_file, "a") as f: - f.write("\n") - f.write( - self._restraint.toString( - engine="GROMACS", - perturbation_type=self._protocol.getPerturbationType(), - restraint_lambda="restraint" - in self._protocol.getLambda(type="series"), - ) + def _apply_ABFE_restraint(self): + # Write the restraint to the topology file + if self._restraint: + with open(self._top_file, "a") as f: + f.write("\n") + f.write( + self._restraint.toString( + engine="GROMACS", + perturbation_type=self._protocol.getPerturbationType(), + restraint_lambda="restraint" + in self._protocol.getLambda(type="series"), ) + ) def _generate_config(self): """Generate GROMACS configuration file strings.""" @@ -488,6 +490,7 @@ def _generate_config(self): if isinstance(self._protocol, _Protocol._FreeEnergyMixin) else None ) + self._apply_ABFE_restraint() self.addToConfig( config.generateGromacsConfig( extra_options={**config_options, **self._extra_options}, @@ -841,7 +844,6 @@ def getSystem(self, block="AUTO"): system : :class:`System ` The latest molecular system. """ - # Wait for the process to finish. if block is True: self.wait() @@ -856,6 +858,10 @@ def getSystem(self, block="AUTO"): if block is True and not self.isError(): return self._getFinalFrame() else: + raise ValueError( + "getSystem should not use `trjconv` to get the frame in production settings." + "Trigger an exception here to show the traceback when that happens." + ) # Minimisation trajectories have a single frame, i.e. the final state. if isinstance(self._protocol, _Protocol.Minimisation): time = 0 * _Units.Time.nanosecond @@ -2396,22 +2402,27 @@ def _parse_energy_terms(text): U-B, Proper-Dih.. Note that this order is absolute and will not be changed by the input to `gmx energy`. """ - sections = text.split("---") + # Get rid of all the nasty message from argo + content = text.split("End your selection with an empty line or a zero.")[1] + sections = content.split("---") # Remove the empty sections - sections = [section for section in sections if section] + sections = [section for section in sections if section.strip()] # Concatenate the lines - section = sections[1].replace("\n", "") + section = sections[0].replace("\n", "") terms = section.split() # Remove the possible '-' from the separation line terms = [term for term in terms if term != "-"] # Check if the index order is correct - indexes = [int(term) for term in terms[::2]] + try: + indexes = [int(term) for term in terms[::2]] + except ValueError: + raise ValueError("Cannot parse the terms: {}".format("\n".join(terms))) energy_names = terms[1::2] length_nomatch = len(indexes) != len(energy_names) # -1 as the index is 1-based. index_nomatch = (_np.arange(len(indexes)) != _np.array(indexes) - 1).any() if length_nomatch or index_nomatch: - raise ValueError(f"Cannot parse the energy terms in the {edr_file} file.") + raise ValueError(f"Cannot parse the energy terms in the {text} file.") else: return energy_names @@ -2858,7 +2869,7 @@ def _find_trajectory_file(self): else: return self._traj_file - def saveMetric( + def _saveMetric( self, filename="metric.parquet", u_nk="u_nk.parquet", dHdl="dHdl.parquet" ): """ @@ -2867,41 +2878,56 @@ def saveMetric( is Free Energy protocol, the dHdl and the u_nk data will be saved in the same parquet format as well. """ - - _assert_imported(_alchemlyb) - - self._update_energy_dict(initialise=True) - datadict_keys = [ - ("Time (ps)", _Units.Time.picosecond, "getTime"), - ( - "PotentialEnergy (kJ/mol)", - _Units.Energy.kj_per_mol, - "getPotentialEnergy", - ), - ] - if not isinstance(self._protocol, _Protocol.Minimisation): - datadict_keys.extend( - [ - ("Volume (nm^3)", _Units.Volume.nanometer3, "getVolume"), - ("Pressure (bar)", _Units.Pressure.bar, "getPressure"), - ( - "Temperature (kelvin)", - _Units.Temperature.kelvin, - "getTemperature", - ), - ] - ) - df = self._convert_datadict_keys(datadict_keys) - df.to_parquet(path=f"{self.workDir()}/{filename}", index=True) - if isinstance(self._protocol, _Protocol.FreeEnergy): - energy = _extract( - f"{self.workDir()}/{self._name}.xvg", - T=self._protocol.getTemperature() / _Units.Temperature.kelvin, - ) - if "u_nk" in energy: - energy["u_nk"].to_parquet(path=f"{self.workDir()}/{u_nk}", index=True) - if "dHdl" in energy: - energy["dHdl"].to_parquet(path=f"{self.workDir()}/{dHdl}", index=True) + if filename is not None: + self._update_energy_dict(initialise=True) + datadict_keys = [ + ("Time (ps)", _Units.Time.picosecond, "getTime"), + ( + "PotentialEnergy (kJ/mol)", + _Units.Energy.kj_per_mol, + "getPotentialEnergy", + ), + ] + if not isinstance(self._protocol, _Protocol.Minimisation): + datadict_keys.extend( + [ + ("Volume (nm^3)", _Units.Volume.nanometer3, "getVolume"), + ("Pressure (bar)", _Units.Pressure.bar, "getPressure"), + ( + "Temperature (kelvin)", + _Units.Temperature.kelvin, + "getTemperature", + ), + ] + ) + df = self._convert_datadict_keys(datadict_keys) + df.to_parquet(path=f"{self.workDir()}/{filename}", index=True) + if u_nk is not None: + _assert_imported(_alchemlyb) + if isinstance(self._protocol, _Protocol.FreeEnergy): + energy = _extract_u_nk( + f"{self.workDir()}/{self._name}.xvg", + T=self._protocol.getTemperature() / _Units.Temperature.kelvin, + ) + with _warnings.catch_warnings(): + _warnings.filterwarnings( + "ignore", + message="The DataFrame has column names of mixed type.", + ) + energy.to_parquet(path=f"{self.workDir()}/{u_nk}", index=True) + if dHdl is not None: + _assert_imported(_alchemlyb) + if isinstance(self._protocol, _Protocol.FreeEnergy): + energy = _extract_dHdl( + f"{self.workDir()}/{self._name}.xvg", + T=self._protocol.getTemperature() / _Units.Temperature.kelvin, + ) + with _warnings.catch_warnings(): + _warnings.filterwarnings( + "ignore", + message="The DataFrame has column names of mixed type.", + ) + energy.to_parquet(path=f"{self.workDir()}/{dHdl}", index=True) def _is_minimisation(config): diff --git a/python/BioSimSpace/Sandpit/Exscientia/Process/_process.py b/python/BioSimSpace/Sandpit/Exscientia/Process/_process.py index 1228a11ed..45579e0c2 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/Process/_process.py +++ b/python/BioSimSpace/Sandpit/Exscientia/Process/_process.py @@ -911,16 +911,6 @@ def wait(self, max_time=None): # The process isn't running. if not self.isRunning(): - try: - self.saveMetric() - except Exception: - exception_info = traceback.format_exc() - with open(f"{self.workDir()}/{self._name}.err", "a+") as f: - f.write("Exception Information during saveMetric():\n") - f.write("======================\n") - f.write(exception_info) - f.write("\n\n") - return if max_time is not None: @@ -1650,12 +1640,25 @@ def _generate_args(self): """Generate the dictionary of command-line arguments.""" self.clearArgs() + def _saveMetric(self, filename, u_nk, dHdl): + """The abstract function to save the metric and free energy data. Need to be + defined for each MD engine.""" + pass + def saveMetric( self, filename="metric.parquet", u_nk="u_nk.parquet", dHdl="dHdl.parquet" ): """The abstract function to save the metric and free energy data. Need to be defined for each MD engine.""" - pass + try: + self._saveMetric(filename, u_nk, dHdl) + except Exception: + exception_info = traceback.format_exc() + with open(f"{self.workDir()}/{self._name}.err", "a+") as f: + f.write("Exception Information during saveMetric():\n") + f.write("======================\n") + f.write(exception_info) + f.write("\n\n") def _convert_datadict_keys(self, datadict_keys): """This function is a helper function for saveMetric that converts a diff --git a/python/BioSimSpace/Sandpit/Exscientia/Protocol/_config.py b/python/BioSimSpace/Sandpit/Exscientia/Protocol/_config.py index 99446e14b..3daebf40d 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/Protocol/_config.py +++ b/python/BioSimSpace/Sandpit/Exscientia/Protocol/_config.py @@ -247,10 +247,24 @@ def generateAmberConfig(self, extra_options=None, extra_lines=None): # Restrain the backbone. restraint = self.protocol.getRestraint() + # Convert to a squashed representation, if needed + if isinstance(self.protocol, _Protocol._FreeEnergyMixin): + atom_mapping0 = _squashed_atom_mapping(self.system, is_lambda1=False) + atom_mapping1 = _squashed_atom_mapping(self.system, is_lambda1=True) + if self.system.getAlchemicalIon(): alchem_ion_idx = self.system.getAlchemicalIonIdx() protein_com_idx = _get_protein_com_idx(self.system) - alchemical_ion_mask = f"@{alchem_ion_idx} | @{protein_com_idx}" + atom_idxs = [alchem_ion_idx, protein_com_idx] + if isinstance(self.protocol, _Protocol._FreeEnergyMixin): + atom_idxs = sorted( + {atom_mapping0[x] for x in atom_idxs if x in atom_mapping0} + | {atom_mapping1[x] for x in atom_idxs if x in atom_mapping1} + ) + # Convert to 1-based index + alchemical_ion_mask = "@" + ",".join( + [str(atom_idx + 1) for atom_idx in atom_idxs] + ) else: alchemical_ion_mask = None @@ -263,10 +277,6 @@ def generateAmberConfig(self, extra_options=None, extra_lines=None): # Convert to a squashed representation, if needed if isinstance(self.protocol, _Protocol._FreeEnergyMixin): - atom_mapping0 = _squashed_atom_mapping( - self.system, is_lambda1=False - ) - atom_mapping1 = _squashed_atom_mapping(self.system, is_lambda1=True) atom_idxs = sorted( {atom_mapping0[x] for x in atom_idxs if x in atom_mapping0} | {atom_mapping1[x] for x in atom_idxs if x in atom_mapping1} diff --git a/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py b/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py index 64426775e..f5c578ece 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py +++ b/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py @@ -1359,7 +1359,11 @@ def getAlchemicalIon(self): The Alchemical Ion or None if there isn't any. """ try: - return self.search("mols with property AlchemicalIon").molecules()[0] + return ( + self.search("mols with property AlchemicalIon") + .molecules()[0] + .getAtoms()[0] + ) except: return None diff --git a/tests/Sandpit/Exscientia/Align/test_alchemical_ion.py b/tests/Sandpit/Exscientia/Align/test_alchemical_ion.py index a3320e5da..f06471e6a 100644 --- a/tests/Sandpit/Exscientia/Align/test_alchemical_ion.py +++ b/tests/Sandpit/Exscientia/Align/test_alchemical_ion.py @@ -5,7 +5,7 @@ _get_protein_com_idx, _mark_alchemical_ion, ) -from BioSimSpace.Sandpit.Exscientia._SireWrappers import Molecule +from BioSimSpace.Sandpit.Exscientia._SireWrappers import Atom from tests.conftest import has_gromacs, root_fp @@ -45,13 +45,13 @@ def test_getAlchemicalIon(input_system, isalchem, request): if isalchem is None: assert ion is None else: - assert isinstance(ion, Molecule) + assert isinstance(ion, Atom) @pytest.mark.skipif(has_gromacs is False, reason="Requires GROMACS to be installed.") def test_getAlchemicalIonIdx(alchemical_ion_system): index = alchemical_ion_system.getAlchemicalIonIdx() - assert index == 680 + assert index == 2053 @pytest.mark.skipif(has_gromacs is False, reason="Requires GROMACS to be installed.") diff --git a/tests/Sandpit/Exscientia/FreeEnergy/test_restraint.py b/tests/Sandpit/Exscientia/FreeEnergy/test_restraint.py index 0149a64e3..a82c89be7 100644 --- a/tests/Sandpit/Exscientia/FreeEnergy/test_restraint.py +++ b/tests/Sandpit/Exscientia/FreeEnergy/test_restraint.py @@ -1,20 +1,22 @@ -import pytest - import numpy as np - -from sire.legacy.Units import angstrom3 as _Sire_angstrom3 -from sire.legacy.Units import k_boltz as _k_boltz -from sire.legacy.Units import meter3 as _Sire_meter3 -from sire.legacy.Units import mole as _Sire_mole +import pytest +from sire.legacy.Units import ( + angstrom3 as _Sire_angstrom3, + k_boltz as _k_boltz, + meter3 as _Sire_meter3, + mole as _Sire_mole, +) import BioSimSpace.Sandpit.Exscientia as BSS from BioSimSpace.Sandpit.Exscientia.Align import decouple from BioSimSpace.Sandpit.Exscientia.FreeEnergy import Restraint -from BioSimSpace.Sandpit.Exscientia.Units.Length import angstrom -from BioSimSpace.Sandpit.Exscientia.Units.Angle import radian, degree +from BioSimSpace.Sandpit.Exscientia.Units.Angle import degree, radian from BioSimSpace.Sandpit.Exscientia.Units.Energy import kcal_per_mol +from BioSimSpace.Sandpit.Exscientia.Units.Length import angstrom from BioSimSpace.Sandpit.Exscientia.Units.Temperature import kelvin +from tests.conftest import has_gromacs + # Store the tutorial URL. url = BSS.tutorialUrl() @@ -75,6 +77,29 @@ def boresch_restraint_component(): return system, restraint_dict +@pytest.mark.skipif(has_gromacs is False, reason="Requires GROMACS to be installed.") +@pytest.mark.parametrize( + ("protocol", "posres"), + [ + (BSS.Protocol.FreeEnergy, "heavy"), + (BSS.Protocol.FreeEnergy, None), + ], +) +def test_top_write( + boresch_restraint_component, protocol, posres, tmp_path, boresch_restraint +): + system, restraint_dict = boresch_restraint_component + bss_protocol = protocol(restraint=posres) + BSS.Process.Gromacs( + system, + bss_protocol, + work_dir=str(tmp_path), + restraint=boresch_restraint, + ) + with open(str(tmp_path / "gromacs.top"), "r") as f: + assert "intermolecular_interactions" in f.read() + + @pytest.fixture(scope="session") def boresch_restraint(boresch_restraint_component): """Generate the Boresch restraint object.""" diff --git a/tests/Sandpit/Exscientia/Process/test_amber.py b/tests/Sandpit/Exscientia/Process/test_amber.py index 1a6daded5..02ad1239a 100644 --- a/tests/Sandpit/Exscientia/Process/test_amber.py +++ b/tests/Sandpit/Exscientia/Process/test_amber.py @@ -7,7 +7,12 @@ import pytest import BioSimSpace.Sandpit.Exscientia as BSS -from tests.Sandpit.Exscientia.conftest import has_amber, has_pyarrow +from BioSimSpace.Sandpit.Exscientia._Exceptions import IncompatibleError +from tests.Sandpit.Exscientia.conftest import ( + has_amber, + has_pyarrow, + url, +) from tests.conftest import root_fp @@ -373,6 +378,10 @@ def test_parse_fep_output(system, protocol): assert len(records_sc1) != 0 +@pytest.mark.skipif( + has_amber is False or has_pyarrow is False, + reason="Requires AMBER and pyarrow to be installed.", +) class TestsaveMetric: @staticmethod @pytest.fixture() @@ -401,13 +410,29 @@ def setup(alchemical_system): process.saveMetric() return process + def test_incompatible_error(self): + alchemical_system = BSS.IO.readPerturbableSystem( + f"{url}/complex_vac0.prm7.bz2", + f"{url}/complex_vac0.rst7.bz2", + f"{url}/complex_vac1.prm7.bz2", + f"{url}/complex_vac1.rst7.bz2", + ) + with pytest.raises( + IncompatibleError, + match="Perturbable system is not compatible with none free energy protocol.", + ): + BSS.Process.Amber( + alchemical_system, + BSS.Protocol.Minimisation(), + ) + def test_error_alchemlyb_extract(self, alchemical_system): # Create a process using any system and the protocol. process = BSS.Process.Amber( alchemical_system, BSS.Protocol.FreeEnergy(temperature=298 * BSS.Units.Temperature.kelvin), ) - process.wait() + process.saveMetric() with open(process.workDir() + "/amber.err", "r") as f: text = f.read() assert "Exception Information" in text @@ -433,3 +458,37 @@ def test_no_output(self, system): process = BSS.Process.Amber(system, BSS.Protocol.Production()) with pytest.warns(match="Simulation didn't produce any output."): process.saveMetric() + + @staticmethod + @pytest.mark.parametrize( + ("filename", "u_nk", "dHdl"), + [ + ("metric.parquet", None, None), + (None, "u_nk.parquet", None), + (None, None, "dHdl.parquet"), + ], + ) + def test_selective(alchemical_system, filename, u_nk, dHdl): + # Create a process using any system and the protocol. + process = BSS.Process.Amber( + alchemical_system, + BSS.Protocol.FreeEnergy(temperature=298 * BSS.Units.Temperature.kelvin), + ) + shutil.copyfile( + f"{root_fp}/Sandpit/Exscientia/output/amber_fep.out", + process.workDir() + "/amber.out", + ) + process.saveMetric(filename, u_nk, dHdl) + if filename is not None: + assert (Path(process.workDir()) / filename).exists() + else: + assert not (Path(process.workDir()) / "metric.parquet").exists() + if u_nk is not None: + assert (Path(process.workDir()) / u_nk).exists() + else: + assert not (Path(process.workDir()) / "u_nk.parquet").exists() + if dHdl is not None: + # The file doesn't contain dHdl information + pass + else: + assert not (Path(process.workDir()) / "dHdl.parquet").exists() diff --git a/tests/Sandpit/Exscientia/Process/test_gromacs.py b/tests/Sandpit/Exscientia/Process/test_gromacs.py index 9a1886e97..70250a847 100644 --- a/tests/Sandpit/Exscientia/Process/test_gromacs.py +++ b/tests/Sandpit/Exscientia/Process/test_gromacs.py @@ -226,7 +226,27 @@ def test_regular_protocol(self, setup, tmp_path_factory): # Run the process and check that it finishes without error. run_process(system, protocol, restraint=restraint, work_dir=str(tmp_path)) with open(tmp_path / "test.top", "r") as f: - assert "intermolecular_interactions" in f.read() + text = f.read() + assert text.count("intermolecular_interactions") == 1 + + def test_position_restraint_protocol(self, setup, tmp_path_factory): + """Test if the restraint has been written in a way that could be processed + correctly. + """ + tmp_path = tmp_path_factory.mktemp("out") + system, restraint = setup + # Create a short production protocol. + protocol = BSS.Protocol.FreeEnergy( + runtime=BSS.Types.Time(0.0001, "nanoseconds"), + perturbation_type="full", + restraint="heavy", + ) + + # Run the process and check that it finishes without error. + run_process(system, protocol, restraint=restraint, work_dir=str(tmp_path)) + with open(tmp_path / "test.top", "r") as f: + text = f.read() + assert text.count("intermolecular_interactions") == 1 def test_restraint_lambda(self, setup, tmp_path_factory): """Test if the restraint has been written correctly when restraint lambda is evoked.""" @@ -383,11 +403,64 @@ def extract(*args): perturbable_system, BSS.Protocol.FreeEnergy(temperature=298 * BSS.Units.Temperature.kelvin), ) - process.wait() + process.saveMetric() with open(process.workDir() + "/gromacs.err", "r") as f: text = f.read() assert "Exception Information" in text + @staticmethod + @pytest.mark.parametrize( + ("filename", "u_nk", "dHdl"), + [ + ("metric.parquet", None, None), + (None, "u_nk.parquet", None), + (None, None, "dHdl.parquet"), + ], + ) + def test_selective(perturbable_system, filename, u_nk, dHdl): + from alchemtest.gmx import load_ABFE + + protocol = BSS.Protocol.FreeEnergy( + runtime=BSS.Types.Time(60, "picosecond"), + timestep=BSS.Types.Time(4, "femtosecond"), + report_interval=200, + ) + process = BSS.Process.Gromacs(perturbable_system, protocol) + shutil.copyfile( + f"{root_fp}/Sandpit/Exscientia/output/gromacs.edr", + process.workDir() + "/gromacs.edr", + ) + shutil.copyfile( + load_ABFE().data["ligand"][0], + process.workDir() + "/gromacs.xvg", + ) + + process.saveMetric(filename, u_nk, dHdl) + if filename is not None: + assert (Path(process.workDir()) / filename).exists() + else: + assert not (Path(process.workDir()) / "metric.parquet").exists() + if u_nk is not None: + assert (Path(process.workDir()) / u_nk).exists() + else: + assert not (Path(process.workDir()) / "u_nk.parquet").exists() + if dHdl is not None: + assert (Path(process.workDir()) / dHdl).exists() + else: + assert not (Path(process.workDir()) / "dHdl.parquet").exists() + + +@pytest.mark.skipif( + has_gromacs is False or has_pyarrow is False, + reason="Requires GROMACS and pyarrow to be installed.", +) +def test_error_saveMetric(perturbable_system): + protocol = BSS.Protocol.FreeEnergy() + process = BSS.Process.Gromacs(perturbable_system, protocol) + process.saveMetric() + with open(f"{process.workDir()}/{process._name}.err", "r") as f: + assert "Exception Information during saveMetric():" in f.read() + @pytest.mark.skipif( has_amber is False diff --git a/tests/Sandpit/Exscientia/Process/test_position_restraint.py b/tests/Sandpit/Exscientia/Process/test_position_restraint.py index b82ec3c1f..542b71c54 100644 --- a/tests/Sandpit/Exscientia/Process/test_position_restraint.py +++ b/tests/Sandpit/Exscientia/Process/test_position_restraint.py @@ -36,12 +36,20 @@ def alchemical_ion_system(): ) ion = solvated.getMolecule(-1) pert_ion = BSS.Align.merge(ion, ion, mapping={0: 0}) + for lambda_ in [0, 1]: + atomtype = pert_ion._sire_object.property(f"atomtype{lambda_}") + pert_ion._sire_object = ( + pert_ion._sire_object.edit() + .setProperty(f"ambertype{lambda_}", atomtype) + .molecule() + ) pert_ion._sire_object = ( pert_ion.getAtoms()[0] ._sire_object.edit() .setProperty("charge1", 0 * SireUnits.mod_electron) .molecule() ) + alchemcial_ion = _mark_alchemical_ion(pert_ion) solvated.updateMolecule(solvated.getIndex(ion), alchemcial_ion) return solvated @@ -164,6 +172,8 @@ def test_gromacs(protocol, system, ref_system, tmp_path): reason="Requires AMBER and openff to be installed", ) def test_amber(protocol, system, ref_system, tmp_path): + if not isinstance(protocol, BSS.Protocol._FreeEnergyMixin): + pytest.skip("AMBER position restraint only works for free energy protocol") proc = BSS.Process.Amber( system, protocol, reference_system=ref_system, work_dir=str(tmp_path) ) @@ -240,19 +250,27 @@ def test_gromacs_alchemical_ion( reason="Requires AMBER, GROMACS and OpenFF to be installed", ) @pytest.mark.parametrize( - ("restraint", "target"), + ("restraint", "protocol", "target"), [ - ("backbone", "@5-7,9,15-17 | @2148 | @8"), - ("heavy", "@2,5-7,9,11,15-17,19 | @2148 | @8"), - ("all", "@1-22 | @2148 | @8"), - ("none", "@2148 | @8"), + ( + "backbone", + BSS.Protocol.FreeEnergyEquilibration, + "@5-7,9,15-17 | @9,6442,6443", + ), + ( + "heavy", + BSS.Protocol.FreeEnergyEquilibration, + "@2,5-7,9,11,15-17,19 | @9,6442,6443", + ), + ("all", BSS.Protocol.FreeEnergyEquilibration, "@1-22 | @9,6442,6443"), + ("none", BSS.Protocol.FreeEnergyEquilibration, "@9,6442,6443"), ], ) def test_amber_alchemical_ion( - alchemical_ion_system, restraint, target, alchemical_ion_system_psores + alchemical_ion_system, restraint, protocol, target, alchemical_ion_system_psores ): # Create an equilibration protocol with backbone restraints. - protocol = BSS.Protocol.Equilibration(restraint=restraint) + protocol = protocol(restraint=restraint) # Create the process object. process = BSS.Process.Amber( diff --git a/tests/Sandpit/Exscientia/Process/test_restart.py b/tests/Sandpit/Exscientia/Process/test_restart.py index 64c00f58d..c4873f7ce 100644 --- a/tests/Sandpit/Exscientia/Process/test_restart.py +++ b/tests/Sandpit/Exscientia/Process/test_restart.py @@ -116,6 +116,8 @@ def test_gromacs(protocol, system, tmp_path): reason="Requires AMBER and OpenFF to be installed.", ) def test_amber(protocol, system, tmp_path): + if not isinstance(protocol, BSS.Protocol._FreeEnergyMixin): + pytest.skip("AMBER position restraint only works for free energy protocol") BSS.Process.Amber(system, protocol, work_dir=str(tmp_path)) with open(tmp_path / "amber.cfg", "r") as f: cfg = f.read()