diff --git a/.github/workflows/Sandpit_exs.yml b/.github/workflows/Sandpit_exs.yml index 76b9a3242..f5b0fc0d3 100644 --- a/.github/workflows/Sandpit_exs.yml +++ b/.github/workflows/Sandpit_exs.yml @@ -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.4" "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=2023.5" "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/README.rst b/README.rst index 562af1769..b1753915e 100644 --- a/README.rst +++ b/README.rst @@ -63,35 +63,35 @@ Conda package The easiest way to install BioSimSpace is using our `conda channel `__. BioSimSpace is built using dependencies from `conda-forge `__, so please ensure that the channel takes strict priority. We recommend using -`Mambaforge `__. +`Miniforge `__. To create a new environment: .. code-block:: bash - mamba create -n openbiosim -c conda-forge -c openbiosim biosimspace - mamba activate openbiosim + conda create -n openbiosim -c conda-forge -c openbiosim biosimspace + conda activate openbiosim To install the latest development version you can use: .. code-block:: bash - mamba create -n openbiosim-dev -c conda-forge -c openbiosim/label/dev biosimspace - mamba activate openbiosim-dev + conda create -n openbiosim-dev -c conda-forge -c openbiosim/label/dev biosimspace + conda activate openbiosim-dev When updating the development version it is generally advised to update `Sire `_ at the same time: .. code-block:: bash - mamba update -c conda-forge -c openbiosim/label/dev biosimspace sire + conda update -c conda-forge -c openbiosim/label/dev biosimspace sire Unless you add the required channels to your Conda configuration, then you'll need to add them when updating, e.g., for the development package: .. code-block:: bash - mamba update -c conda-forge -c openbiosim/label/dev biosimspace + conda update -c conda-forge -c openbiosim/label/dev biosimspace Installing from source ^^^^^^^^^^^^^^^^^^^^^^ @@ -146,8 +146,8 @@ latest development code into that. .. code-block:: bash - mamba create -n openbiosim-dev -c conda-forge -c openbiosim/label/dev biosimspace --only-deps - mamba activate openbiosim-dev + conda create -n openbiosim-dev -c conda-forge -c openbiosim/label/dev biosimspace --only-deps + conda activate openbiosim-dev git clone https://github.com/openbiosim/biosimspace cd biosimspace/python BSS_SKIP_DEPENDENCIES=1 python setup.py develop diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 5420e438a..0e8765d6c 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -9,6 +9,19 @@ company supporting open-source development of fostering academic/industrial coll within the biomolecular simulation community. Our software is hosted via the `OpenBioSim` `GitHub `__ organisation. +`2023.5.1 `_ - Mar 20 2024 +------------------------------------------------------------------------------------------------- + +* Fixed path to user links file in the :func:`generateNetwork ` function (`#233 `__). +* Fixed redirection of stderr (`#233 `__). +* Switched to using ``AtomCoordMatcher`` to map parameterised molecules back to their original topology. This resolves issues where atoms moved between residues following parameterisation (`#235 `__). +* Make the GROMACS ``_generate_binary_run_file`` function static so that it can be used when initialising free energy simulations in setup-only mode (`#237 `__). +* Improve error handling and message when attempting to extract an all dummy atom selection (`#251 `__). +* Don't set SOMD specific end-state properties when decoupling a molecule (`#253 `__). +* Only convert to a end-state system when not running a free energy protocol with GROMACS so that hybrid topology isn't lost when using position restraints (`#257 `__). +* Exclude standard free ions from the AMBER position restraint mask (`#260 `__). +* Update the ``BioSimSpace.Types._GeneralUnit.__pow__`` operator to support fractional exponents (`#260 `__). + `2023.5.0 `_ - Dec 16 2023 ------------------------------------------------------------------------------------------------- diff --git a/doc/source/contributing/packaging.rst b/doc/source/contributing/packaging.rst index b98bbbf06..5c1425fde 100644 --- a/doc/source/contributing/packaging.rst +++ b/doc/source/contributing/packaging.rst @@ -4,12 +4,11 @@ Development process =================== -:mod:`BioSimSpace` uses a ``main``, ``devel`` and ``future`` development process, +:mod:`BioSimSpace` uses a ``main`` and ``devel`` development process, using feature branches for all code development. * ``main`` - this always contains the latest official release. * ``devel`` - this always contains the latest development release, which will become the next official release. -* ``future`` - this contains pull requests that have been accepted, but which are targetted for a future release (i.e. not the next official release) Code should be developed on a fork or in a feature branch called ``feature_{feature}``. When your feature is ready, please submit a pull request against ``devel``. This @@ -27,29 +26,8 @@ tests, examples and/or tutorial instructions. the tutorials or writing a detailed description for the website. Assuming the CI completes successfully, then one of the release team will -conduct a code review. The outcome of the review will be one of the following; - -1. This feature is ready, and should be part of the next official release. The pull request - will be accepted into ``devel``. This will trigger our CI/CD process, building the new dev - package and uploading it to `anaconda.org `__ - for everyone to use. - -2. This feature is good, but it is not yet ready to be part of the next offical release. This - could be because the feature is part of a series, and all of the series need to be finished - before release. Or because we are in a feature freeze period. Or because you want more time - for people to explore and play with the feature before it is officially released (and would - then need to be supported, and backwards compatibility maintained). If this is the case (or - it is your request) then the pull request will be redirected into the ``future`` branch. - Once it (and features that depend on it) are ready, you can then issue a pull request for - all of the features at once into ``devel``. It will be noted that each of the individual - parts have already been code reviewed, so the process to accept the combination - into ``devel`` should be more straightforward. - -3. This feature is good, but more work is needed before it can be accepted. This could be - because some of the unit tests haven't passed, or the latest version of ``devel`` hasn't - been merged. Or there may be changes that are requested that would make the code easier - to maintain or to preserve backwards compatibility. If this is the case, then we - will engage in conversation with you and will work together to rectify any issues. +conduct a code review, with the code being merged into ``devel`` if it is +approved. Bug fixes or issue fixes are developed on fix branches, called ``fix_{number}`` (again in either the main repository or forks). If no `issue thread `__ diff --git a/doc/source/contributing/roadmap.rst b/doc/source/contributing/roadmap.rst index 808af2c3e..5a13a2b17 100644 --- a/doc/source/contributing/roadmap.rst +++ b/doc/source/contributing/roadmap.rst @@ -76,8 +76,6 @@ You can keep up with what we are working on in several ways; * Keep an eye on the various ``feature_X`` branches as they appear in the repository. Feel free to initiate a conversation on GitHub with the developer who is working on that branch if you want to learn more, or want to make suggestions or offer a helping hand. -* Clone and build your own copy of the ``future`` branch. This is the bleeding edge, and things may change and break. - But it is the earliest way to use the future version of :mod:`BioSimSpace`. Wishlists / suggestions ======================= diff --git a/doc/source/install.rst b/doc/source/install.rst index 87a91a704..728491e61 100644 --- a/doc/source/install.rst +++ b/doc/source/install.rst @@ -46,27 +46,25 @@ The easiest way to install :mod:`BioSimSpace` is in a new `conda environment `__. You can use any conda environment or installation. We recommend using -`mambaforge `__, -as this is pre-configured to use `conda-forge `__, -and bundles `mamba `__, which -is a fast drop-in replacement for `conda `__. +`Miniforge `__, +as this is pre-configured to use `conda-forge `__. -.. _Install_Mambaforge: -Either... Install a new copy of ``mambaforge`` +.. _Install_Miniforge: +Either... Install a new copy of ``Miniforge`` ---------------------------------------------- To install a new copy of -`mambaforge `__, -first download a ``Mambaforge`` from -`this page `__ that +`Miniforge `__, +first download a ``Miniforge`` from +`this page `__ that matches your operating system and processor. -Install ``Mambaforge`` following the +Install ``Miniforge`` following the `instructions here `__. -Once installed, you should be able to run the ``mamba`` command to -install other packages (e.g. ``mamba -h`` will print out help on -how to use the ``mamba`` command). +Once installed, you should be able to run the ``conda`` command to +install other packages (e.g. ``conda -h`` will print out help on +how to use the ``conda`` command). Or... Use an existing anaconda/miniconda install ------------------------------------------------ @@ -80,21 +78,7 @@ the full path to your anaconda or miniconda installation. You should now be able to run the ``conda`` command to install other packages (e.g. ``conda -h`` will print out help on how to use the -``conda`` command). We highly recommend that you use ``mamba`` as a -drop-in replacement for ``conda``, so first install ``mamba``. - -.. code-block:: bash - - $ conda install -c conda-forge mamba - -This should install mamba. If this fails, then your anaconda or miniconda -environment is likely quite full, or else it is outdated. We recommend -going back and following `the instructions <_Install_Mambaforge>` -to install a new copy of ``mambaforge``. - -If this works, then you should now be able to run the ``mamba`` command -to install other packages (e.g. ``mamba -h`` will print out help -on how to use the ``mamba`` command). +``conda`` command). And then... Install BioSimSpace into a new environment ------------------------------------------------------ @@ -107,7 +91,7 @@ by creating a Python 3.9 environment that we will call ``openbiosim``. .. code-block:: bash - $ mamba create -n openbiosim "python<3.10" + $ conda create -n openbiosim "python<3.10" .. note:: @@ -118,27 +102,27 @@ We can now install :mod:`BioSimSpace` into that environment by typing .. code-block:: bash - $ mamba install -n openbiosim -c openbiosim biosimspace + $ conda install -n openbiosim -c openbiosim biosimspace .. note:: - The option ``-n openbiosim`` tells ``mamba`` to install :mod:`BioSimSpace` + The option ``-n openbiosim`` tells ``conda`` to install :mod:`BioSimSpace` into the ``openbiosim`` environment. The option ``-c openbiosim`` - tells ``mamba`` to install :mod:`BioSimSpace` from the ``openbiosim`` + tells ``conda`` to install :mod:`BioSimSpace` from the ``openbiosim`` conda channel. If you want the latest development release, then install by typing .. code-block:: bash - $ mamba install -n openbiosim -c "openbiosim/label/dev" biosimspace + $ conda install -n openbiosim -c "openbiosim/label/dev" biosimspace To install the latest development version you can use: .. code-block:: bash - mamba create -n openbiosim-dev -c conda-forge -c openbiosim/label/dev biosimspace - mamba activate openbiosim-dev + conda create -n openbiosim-dev -c conda-forge -c openbiosim/label/dev biosimspace + conda activate openbiosim-dev To run :mod:`BioSimSpace`, you must now activate the ``openbiosim`` environment. You can do this by typing @@ -268,8 +252,8 @@ latest development code into that. .. code-block:: bash - mamba create -n openbiosim-dev -c conda-forge -c openbiosim/label/dev biosimspace --only-deps - mamba activate openbiosim-dev + conda create -n openbiosim-dev -c conda-forge -c openbiosim/label/dev biosimspace --only-deps + conda activate openbiosim-dev git clone https://github.com/openbiosim/biosimspace cd biosimspace/python BSS_SKIP_DEPENDENCIES=1 python setup.py develop diff --git a/python/BioSimSpace/FreeEnergy/_relative.py b/python/BioSimSpace/FreeEnergy/_relative.py index a54a6d707..4e197f1dc 100644 --- a/python/BioSimSpace/FreeEnergy/_relative.py +++ b/python/BioSimSpace/FreeEnergy/_relative.py @@ -406,7 +406,7 @@ def getData(self, name="data", file_link=False, work_dir=None): ) # Write to the zip file. - with _zipfile.ZipFile(cwd + f"/{zipname}", "w") as zip: + with _zipfile.Zipfile(_os.join(cwd, zipname), "w") as zip: for file in files: zip.write(file) @@ -2074,15 +2074,15 @@ def _initialise_runner(self, system): process._system = first_process._system.copy() process._protocol = self._protocol process._work_dir = new_dir - process._std_out_file = new_dir + "/somd.out" - process._std_err_file = new_dir + "/somd.err" - process._rst_file = new_dir + "/somd.rst7" - process._top_file = new_dir + "/somd.prm7" - process._traj_file = new_dir + "/traj000000001.dcd" - process._restart_file = new_dir + "/latest.rst" - process._config_file = new_dir + "/somd.cfg" - process._pert_file = new_dir + "/somd.pert" - process._gradients_file = new_dir + "/gradients.dat" + process._std_out_file = _os.path.join(new_dir, "somd.out") + process._std_err_file = _os.path.join(new_dir, "somd.err") + process._rst_file = _os.path.join(new_dir, "somd.rst7") + process._top_file = _os.path.join(new_dir, "somd.prm7") + process._traj_file = _os.path.join(new_dir, "traj000000001.dcd") + process._restart_file = _os.path.join(new_dir, "latest.rst") + process._config_file = _os.path.join(new_dir, "somd.cfg") + process._pert_file = _os.path.join(new_dir, "somd.pert") + process._gradients_file = _os.path.join(new_dir, "gradients.dat") process._input_files = [ process._config_file, process._rst_file, @@ -2106,10 +2106,10 @@ def _initialise_runner(self, system): for line in new_config: f.write(line) - mdp = new_dir + "/gromacs.mdp" - gro = new_dir + "/gromacs.gro" - top = new_dir + "/gromacs.top" - tpr = new_dir + "/gromacs.tpr" + mdp = _os.path.join(new_dir, "gromacs.mdp") + gro = _os.path.join(new_dir, "gromacs.gro") + top = _os.path.join(new_dir, "gromacs.top") + tpr = _os.path.join(new_dir, "gromacs.tpr") # Use grompp to generate the portable binary run input file. _Process.Gromacs._generate_binary_run_file( @@ -2129,14 +2129,14 @@ def _initialise_runner(self, system): process._system = first_process._system.copy() process._protocol = self._protocol process._work_dir = new_dir - process._std_out_file = new_dir + "/gromacs.out" - process._std_err_file = new_dir + "/gromacs.err" - process._gro_file = new_dir + "/gromacs.gro" - process._top_file = new_dir + "/gromacs.top" - process._ref_file = new_dir + "/gromacs_ref.gro" - process._traj_file = new_dir + "/gromacs.trr" - process._config_file = new_dir + "/gromacs.mdp" - process._tpr_file = new_dir + "/gromacs.tpr" + process._std_out_file = _os.path.join(new_dir, "gromacs.out") + process._std_err_file = _os.path.join(new_dir, "gromacs.err") + process._gro_file = _os.path.join(new_dir, "gromacs.gro") + process._top_file = _os.path.join(new_dir, "gromacs.top") + process._ref_file = _os.path.join(new_dir, "gromacs_ref.gro") + process._traj_file = _os.path.join(new_dir, "gromacs.trr") + process._config_file = _os.path.join(new_dir, "gromacs.mdp") + process._tpr_file = _os.path.join(new_dir, "gromacs.tpr") process._input_files = [ process._config_file, process._gro_file, @@ -2165,14 +2165,14 @@ def _initialise_runner(self, system): process._system = first_process._system.copy() process._protocol = self._protocol process._work_dir = new_dir - process._std_out_file = new_dir + "/amber.out" - process._std_err_file = new_dir + "/amber.err" - process._rst_file = new_dir + "/amber.rst7" - process._top_file = new_dir + "/amber.prm7" - process._ref_file = new_dir + "/amber_ref.rst7" - process._traj_file = new_dir + "/amber.nc" - process._config_file = new_dir + "/amber.cfg" - process._nrg_file = new_dir + "/amber.nrg" + process._std_out_file = _os.path.join(new_dir, "amber.out") + process._std_err_file = _os.path.join(new_dir, "amber.err") + process._rst_file = _os.path.join(new_dir, "amber.rst7") + process._top_file = _os.path.join(new_dir, "amber.prm7") + process._ref_file = _os.path.join(new_dir, "amber_ref.rst7") + process._traj_file = _os.path.join(new_dir, "amber.nc") + process._config_file = _os.path.join(new_dir, "amber.cfg") + process._nrg_file = _os.path.join(new_dir, "amber.nrg") process._input_files = [ process._config_file, process._rst_file, diff --git a/python/BioSimSpace/Parameters/_Protocol/_amber.py b/python/BioSimSpace/Parameters/_Protocol/_amber.py index 69816c1e4..2ce0b15f9 100644 --- a/python/BioSimSpace/Parameters/_Protocol/_amber.py +++ b/python/BioSimSpace/Parameters/_Protocol/_amber.py @@ -316,9 +316,6 @@ def run(self, molecule, work_dir=None, queue=None): else: is_smiles = False - # Create the file prefix. - prefix = work_dir + "/" - if not is_smiles: # Create a copy of the molecule. new_mol = molecule.copy() @@ -352,7 +349,10 @@ def run(self, molecule, work_dir=None, queue=None): ) # Prepend the working directory to the output file names. - output = [prefix + output[0], prefix + output[1]] + output = [ + _os.path.join(str(work_dir), output[0]), + _os.path.join(str(work_dir), output[1]), + ] try: # Load the parameterised molecule. (This could be a system of molecules.) @@ -443,9 +443,6 @@ def _run_tleap(self, molecule, work_dir): else: _molecule = molecule - # Create the file prefix. - prefix = work_dir + "/" - # Write the system to a PDB file. try: # LEaP expects residue numbering to be ascending and continuous. @@ -454,7 +451,7 @@ def _run_tleap(self, molecule, work_dir): )[0] renumbered_molecule = _Molecule(renumbered_molecule) _IO.saveMolecules( - prefix + "leap", + _os.path.join(str(work_dir), "leap"), renumbered_molecule, "pdb", property_map=self._property_map, @@ -500,7 +497,7 @@ def _run_tleap(self, molecule, work_dir): pruned_bond_records.append(bond) # Write the LEaP input file. - with open(prefix + "leap.txt", "w") as file: + with open(_os.path.join(str(work_dir), "leap.txt"), "w") as file: file.write("source %s\n" % ff) if self._water_model is not None: if self._water_model in ["tip4p", "tip5p"]: @@ -528,14 +525,14 @@ def _run_tleap(self, molecule, work_dir): # Generate the tLEaP command. command = "%s -f leap.txt" % _tleap_exe - with open(prefix + "README.txt", "w") as file: + with open(_os.path.join(str(work_dir), "README.txt"), "w") as file: # Write the command to file. file.write("# tLEaP was run with the following command:\n") file.write("%s\n" % command) # Create files for stdout/stderr. - stdout = open(prefix + "leap.out", "w") - stderr = open(prefix + "leap.err", "w") + stdout = open(_os.path.join(str(work_dir), "leap.out"), "w") + stderr = open(_os.path.join(str(work_dir), "leap.err"), "w") # Run tLEaP as a subprocess. proc = _subprocess.run( @@ -550,12 +547,12 @@ def _run_tleap(self, molecule, work_dir): # tLEaP doesn't return sensible error codes, so we need to check that # the expected output was generated. - if _os.path.isfile(prefix + "leap.top") and _os.path.isfile( - prefix + "leap.crd" - ): + if _os.path.isfile( + _os.path.join(str(work_dir), "leap.top") + ) and _os.path.isfile(_os.path.join(str(work_dir), "leap.crd")): # Check the output of tLEaP for missing atoms. if self._ensure_compatible: - if _has_missing_atoms(prefix + "leap.out"): + if _has_missing_atoms(_os.path.join(str(work_dir), "leap.out")): raise _ParameterisationError( "tLEaP added missing atoms. The topology is now " "inconsistent with the original molecule. Please " @@ -604,13 +601,13 @@ def _run_pdb2gmx(self, molecule, work_dir): else: _molecule = molecule - # Create the file prefix. - prefix = work_dir + "/" - # Write the system to a PDB file. try: _IO.saveMolecules( - prefix + "leap", _molecule, "pdb", property_map=self._property_map + _os.path.join(str(work_dir), "input"), + _molecule, + "pdb", + property_map=self._property_map, ) except Exception as e: msg = "Failed to write system to 'PDB' format." @@ -626,14 +623,14 @@ def _run_pdb2gmx(self, molecule, work_dir): % (_gmx_exe, supported_ff[self._forcefield]) ) - with open(prefix + "README.txt", "w") as file: + with open(_os.path.join(str(work_dir), "README.txt"), "w") as file: # Write the command to file. file.write("# pdb2gmx was run with the following command:\n") file.write("%s\n" % command) # Create files for stdout/stderr. - stdout = open(prefix + "pdb2gmx.out", "w") - stderr = open(prefix + "pdb2gmx.err", "w") + stdout = open(_os.path.join(str(work_dir), "pdb2gmx.out"), "w") + stderr = open(_os.path.join(str(work_dir), "pdb2gmx.err"), "w") # Run pdb2gmx as a subprocess. proc = _subprocess.run( @@ -647,9 +644,9 @@ def _run_pdb2gmx(self, molecule, work_dir): stderr.close() # Check for the expected output. - if _os.path.isfile(prefix + "output.gro") and _os.path.isfile( - prefix + "output.top" - ): + if _os.path.isfile( + _os.path.join(str(work_dir), "output.gro") + ) and _os.path.isfile(_os.path.join(str(work_dir), "output.top")): return ["output.gro", "output.top"] else: raise _ParameterisationError("pdb2gmx failed!") @@ -1010,9 +1007,6 @@ def run(self, molecule, work_dir=None, queue=None): if work_dir is None: work_dir = _os.getcwd() - # Create the file prefix. - prefix = work_dir + "/" - # Convert SMILES to a molecule. if isinstance(molecule, str): is_smiles = True @@ -1092,7 +1086,10 @@ def run(self, molecule, work_dir=None, queue=None): # Write the system to a PDB file. try: _IO.saveMolecules( - prefix + "antechamber", new_mol, "pdb", property_map=self._property_map + _os.path.join(str(work_dir), "antechamber"), + new_mol, + "pdb", + property_map=self._property_map, ) except Exception as e: msg = "Failed to write system to 'PDB' format." @@ -1108,14 +1105,14 @@ def run(self, molecule, work_dir=None, queue=None): + "-o antechamber.mol2 -fo mol2 -c %s -s 2 -nc %d" ) % (_antechamber_exe, self._version, self._charge_method.lower(), charge) - with open(prefix + "README.txt", "w") as file: + 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") file.write("%s\n" % command) # Create files for stdout/stderr. - stdout = open(prefix + "antechamber.out", "w") - stderr = open(prefix + "antechamber.err", "w") + stdout = open(_os.path.join(str(work_dir), "antechamber.out"), "w") + stderr = open(_os.path.join(str(work_dir), "antechamber.err"), "w") # Run Antechamber as a subprocess. proc = _subprocess.run( @@ -1130,20 +1127,20 @@ def run(self, molecule, work_dir=None, queue=None): # Antechamber doesn't return sensible error codes, so we need to check that # the expected output was generated. - if _os.path.isfile(prefix + "antechamber.mol2"): + if _os.path.isfile(_os.path.join(str(work_dir), "antechamber.mol2")): # Run parmchk to check for missing parameters. command = ( "%s -s %d -i antechamber.mol2 -f mol2 " + "-o antechamber.frcmod" ) % (_parmchk_exe, self._version) - with open(prefix + "README.txt", "a") as file: + with open(_os.path.join(str(work_dir), "README.txt"), "a") as file: # Write the command to file. file.write("\n# ParmChk was run with the following command:\n") file.write("%s\n" % command) # Create files for stdout/stderr. - stdout = open(prefix + "parmchk.out", "w") - stderr = open(prefix + "parmchk.err", "w") + stdout = open(_os.path.join(str(work_dir), "parmchk.out"), "w") + stderr = open(_os.path.join(str(work_dir), "parmchk.err"), "w") # Run parmchk as a subprocess. proc = _subprocess.run( @@ -1157,7 +1154,7 @@ def run(self, molecule, work_dir=None, queue=None): stderr.close() # The frcmod file was created. - if _os.path.isfile(prefix + "antechamber.frcmod"): + if _os.path.isfile(_os.path.join(str(work_dir), "antechamber.frcmod")): # Now call tLEaP using the partially parameterised molecule and the frcmod file. # tLEap will run in the same working directory, using the Mol2 file generated by # Antechamber. @@ -1169,7 +1166,7 @@ def run(self, molecule, work_dir=None, queue=None): ff = _find_force_field("gaff2") # Write the LEaP input file. - with open(prefix + "leap.txt", "w") as file: + with open(_os.path.join(str(work_dir), "leap.txt"), "w") as file: file.write("source %s\n" % ff) file.write("mol = loadMol2 antechamber.mol2\n") file.write("loadAmberParams antechamber.frcmod\n") @@ -1179,14 +1176,14 @@ def run(self, molecule, work_dir=None, queue=None): # Generate the tLEaP command. command = "%s -f leap.txt" % _tleap_exe - with open(prefix + "README.txt", "a") as file: + with open(_os.path.join(str(work_dir), "README.txt"), "a") as file: # Write the command to file. file.write("\n# tLEaP was run with the following command:\n") file.write("%s\n" % command) # Create files for stdout/stderr. - stdout = open(prefix + "leap.out", "w") - stderr = open(prefix + "leap.err", "w") + stdout = open(_os.path.join(str(work_dir), "leap.out"), "w") + stderr = open(_os.path.join(str(work_dir), "leap.err"), "w") # Run tLEaP as a subprocess. proc = _subprocess.run( @@ -1201,12 +1198,12 @@ def run(self, molecule, work_dir=None, queue=None): # tLEaP doesn't return sensible error codes, so we need to check that # the expected output was generated. - if _os.path.isfile(prefix + "leap.top") and _os.path.isfile( - prefix + "leap.crd" - ): + if _os.path.isfile( + _os.path.join(str(work_dir), "leap.top") + ) and _os.path.isfile(_os.path.join(str(work_dir), "leap.crd")): # Check the output of tLEaP for missing atoms. if self._ensure_compatible: - if _has_missing_atoms(prefix + "leap.out"): + if _has_missing_atoms(_os.path.join(str(work_dir), "leap.out")): raise _ParameterisationError( "tLEaP added missing atoms. The topology is now " "inconsistent with the original molecule. Please " @@ -1217,7 +1214,10 @@ def run(self, molecule, work_dir=None, queue=None): # Load the parameterised molecule. (This could be a system of molecules.) try: par_mol = _IO.readMolecules( - [prefix + "leap.top", prefix + "leap.crd"] + [ + _os.path.join(str(work_dir), "leap.top"), + _os.path.join(str(work_dir), "leap.crd"), + ], ) # Extract single molecules. if par_mol.nMolecules() == 1: diff --git a/python/BioSimSpace/Parameters/_Protocol/_openforcefield.py b/python/BioSimSpace/Parameters/_Protocol/_openforcefield.py index 018f43e4e..c94e7d349 100644 --- a/python/BioSimSpace/Parameters/_Protocol/_openforcefield.py +++ b/python/BioSimSpace/Parameters/_Protocol/_openforcefield.py @@ -214,9 +214,6 @@ def run(self, molecule, work_dir=None, queue=None): if work_dir is None: work_dir = _os.getcwd() - # Create the file prefix. - prefix = work_dir + "/" - # Flag whether the molecule is a SMILES string. if isinstance(molecule, str): is_smiles = True @@ -256,7 +253,7 @@ def run(self, molecule, work_dir=None, queue=None): # Write the molecule to SDF format. try: _IO.saveMolecules( - prefix + "molecule", + _os.path.join(str(work_dir), "molecule"), molecule, "sdf", property_map=self._property_map, @@ -275,7 +272,7 @@ def run(self, molecule, work_dir=None, queue=None): # Write the molecule to a PDB file. try: _IO.saveMolecules( - prefix + "molecule", + _os.path.join(str(work_dir), "molecule"), molecule, "pdb", property_map=self._property_map, @@ -291,7 +288,7 @@ def run(self, molecule, work_dir=None, queue=None): # Create an RDKit molecule from the PDB file. try: rdmol = _Chem.MolFromPDBFile( - prefix + "molecule.pdb", removeHs=False + _os.path.join(str(work_dir), "molecule.pdb"), removeHs=False ) except Exception as e: msg = "RDKit was unable to read the molecular PDB file!" @@ -303,7 +300,9 @@ def run(self, molecule, work_dir=None, queue=None): # Use RDKit to write back to SDF format. try: - writer = _Chem.SDWriter(prefix + "molecule.sdf") + writer = _Chem.SDWriter( + _os.path.join(str(work_dir), "molecule.sdf") + ) writer.write(rdmol) writer.close() except Exception as e: @@ -317,7 +316,9 @@ def run(self, molecule, work_dir=None, queue=None): # Create the Open Forcefield Molecule from the intermediate SDF file, # as recommended by @j-wags and @mattwthompson. try: - off_molecule = _OpenFFMolecule.from_file(prefix + "molecule.sdf") + off_molecule = _OpenFFMolecule.from_file( + _os.path.join(str(work_dir), "molecule.sdf") + ) except Exception as e: msg = "Unable to create OpenFF Molecule!" if _isVerbose(): @@ -383,8 +384,8 @@ def run(self, molecule, work_dir=None, queue=None): # Export AMBER format files. try: - interchange.to_prmtop(prefix + "interchange.prm7") - interchange.to_inpcrd(prefix + "interchange.rst7") + interchange.to_prmtop(_os.path.join(str(work_dir), "interchange.prm7")) + interchange.to_inpcrd(_os.path.join(str(work_dir), "interchange.rst7")) except Exception as e: msg = "Unable to write Interchange object to AMBER format!" if _isVerbose(): @@ -396,7 +397,10 @@ def run(self, molecule, work_dir=None, queue=None): # Load the parameterised molecule. (This could be a system of molecules.) try: par_mol = _IO.readMolecules( - [prefix + "interchange.prm7", prefix + "interchange.rst7"] + [ + _os.path.join(str(work_dir), "interchange.prm7"), + _os.path.join(str(work_dir), "interchange.rst7"), + ], ) # Extract single molecules. if par_mol.nMolecules() == 1: diff --git a/python/BioSimSpace/Process/_amber.py b/python/BioSimSpace/Process/_amber.py index 466216e87..8140729a4 100644 --- a/python/BioSimSpace/Process/_amber.py +++ b/python/BioSimSpace/Process/_amber.py @@ -235,15 +235,15 @@ def __init__( self._is_header = False # The names of the input files. - self._rst_file = "%s/%s.rst7" % (self._work_dir, name) - self._top_file = "%s/%s.prm7" % (self._work_dir, name) - self._ref_file = "%s/%s_ref.rst7" % (self._work_dir, name) + self._rst_file = _os.path.join(str(self._work_dir), f"{name}.rst7") + self._top_file = _os.path.join(str(self._work_dir), f"{name}.prm7") + self._ref_file = _os.path.join(str(self._work_dir), f"{name}_ref.rst7") # The name of the trajectory file. - self._traj_file = "%s/%s.nc" % (self._work_dir, name) + self._traj_file = _os.path.join(str(self._work_dir), f"{name}.nc") # Set the path for the AMBER configuration file. - self._config_file = "%s/%s.cfg" % (self._work_dir, name) + self._config_file = _os.path.join(str(self._work_dir), f"{name}.cfg") # Create the list of input files. self._input_files = [self._config_file, self._rst_file, self._top_file] @@ -400,7 +400,9 @@ def _generate_config(self): if auxiliary_files is not None: for file in auxiliary_files: file_name = _os.path.basename(file) - _shutil.copyfile(file, self._work_dir + f"/{file_name}") + _shutil.copyfile( + file, _os.path.join(str(self._work_dir), file_name) + ) self._input_files.append(self._plumed_config_file) # Expose the PLUMED specific member functions. @@ -534,7 +536,7 @@ def getSystem(self, block="AUTO"): _warnings.warn("The process exited with an error!") # Create the name of the restart CRD file. - restart = "%s/%s.crd" % (self._work_dir, self._name) + restart = _os.path.join(str(self._work_dir), "%s.crd" % self._name) # Check that the file exists. if _os.path.isfile(restart): diff --git a/python/BioSimSpace/Process/_gromacs.py b/python/BioSimSpace/Process/_gromacs.py index a35fc3531..7de88ac88 100644 --- a/python/BioSimSpace/Process/_gromacs.py +++ b/python/BioSimSpace/Process/_gromacs.py @@ -208,18 +208,18 @@ def __init__( self._energy_file = "%s/%s.edr" % (self._work_dir, name) # The names of the input files. - self._gro_file = "%s/%s.gro" % (self._work_dir, name) - self._top_file = "%s/%s.top" % (self._work_dir, name) - self._ref_file = "%s/%s_ref.gro" % (self._work_dir, name) + self._gro_file = _os.path.join(str(self._work_dir), f"{name}.gro") + self._top_file = _os.path.join(str(self._work_dir), f"{name}.top") + self._ref_file = _os.path.join(str(self._work_dir), f"{name}_ref.gro") # The name of the trajectory file. - self._traj_file = "%s/%s.trr" % (self._work_dir, name) + self._traj_file = _os.path.join(str(self._work_dir), f"{name}.trr") # The name of the output coordinate file. - self._crd_file = "%s/%s_out.gro" % (self._work_dir, name) + self._crd_file = _os.path.join(str(self._work_dir), f"{name}_out.gro") # Set the path for the GROMACS configuration file. - self._config_file = "%s/%s.mdp" % (self._work_dir, name) + self._config_file = _os.path.join(str(self._work_dir), f"{name}.mdp") # Create the list of input files. self._input_files = [self._config_file, self._gro_file, self._top_file] @@ -314,7 +314,7 @@ def _setup(self, **kwargs): ) # Create the binary input file name. - self._tpr_file = "%s/%s.tpr" % (self._work_dir, self._name) + self._tpr_file = _os.path.join(str(self._work_dir), f"{self._name}.tpr") self._input_files.append(self._tpr_file) # Generate the GROMACS configuration file. @@ -397,7 +397,9 @@ def _generate_config(self): if auxiliary_files is not None: for file in auxiliary_files: file_name = _os.path.basename(file) - _shutil.copyfile(file, self._work_dir + f"/{file_name}") + _shutil.copyfile( + file, _os.path.join(str(self._work_dir), file_name) + ) self._input_files.append(self._plumed_config_file) # Expose the PLUMED specific member functions. @@ -423,7 +425,9 @@ def _generate_config(self): if auxiliary_files is not None: for file in auxiliary_files: file_name = _os.path.basename(file) - _shutil.copyfile(file, self._work_dir + f"/{file_name}") + _shutil.copyfile( + file, _os.path.join(str(self._work_dir), file_name) + ) self._input_files.append(self._plumed_config_file) # Expose the PLUMED specific member functions. @@ -2122,7 +2126,9 @@ def _add_position_restraints(self): if len(restrained_atoms) > 0: # Create the file names. include_file = "posre_%04d.itp" % num_restraint - restraint_file = "%s/%s" % (self._work_dir, include_file) + restraint_file = _os.path.join( + str(self._work_dir), include_file + ) with open(restraint_file, "w") as file: # Write the header. @@ -2206,7 +2212,9 @@ def _add_position_restraints(self): if len(atom_idxs) > 0: # Create the file names. include_file = "posre_%04d.itp" % num_restraint - restraint_file = "%s/%s" % (self._work_dir, include_file) + restraint_file = _os.path.join( + str(self._work_dir), include_file + ) with open(restraint_file, "w") as file: # Write the header. @@ -2735,11 +2743,12 @@ def _getFrame(self, time): return old_system except: + raise _warnings.warn( "Failed to extract trajectory frame with trjconv. " "Try running 'getSystem' again." ) - frame = "%s/frame.gro" % self._work_dir + frame = _os.path.join(str(self._work_dir), "frame.gro") if _os.path.isfile(frame): _os.remove(frame) return None @@ -2759,7 +2768,7 @@ def _find_trajectory_file(self): # Check that the current trajectory file is found. if not _os.path.isfile(self._traj_file): # If not, first check for any trr extension. - traj_file = _glob.glob("%s/*.trr" % self._work_dir) + traj_file = _glob.glob(_os.path.join(str(self._work_dir), "*.trr")) # Store the number of trr files. num_trr = len(traj_file) @@ -2769,7 +2778,7 @@ def _find_trajectory_file(self): return traj_file[0] else: # Now check for any xtc files. - traj_file = _glob.glob("%s/*.xtc" % self._work_dir) + traj_file = _glob.glob(_os.path.join(str(self._work_dir), "*.xtc")) if len(traj_file) == 1: return traj_file[0] diff --git a/python/BioSimSpace/Process/_namd.py b/python/BioSimSpace/Process/_namd.py index 9811d26ab..787452905 100644 --- a/python/BioSimSpace/Process/_namd.py +++ b/python/BioSimSpace/Process/_namd.py @@ -148,17 +148,17 @@ def __init__( self._stdout_title = None # The names of the input files. - self._psf_file = "%s/%s.psf" % (self._work_dir, name) - self._top_file = "%s/%s.pdb" % (self._work_dir, name) - self._param_file = "%s/%s.params" % (self._work_dir, name) + self._psf_file = _os.path.join(str(self._work_dir), f"{name}.psf") + self._top_file = _os.path.join(str(self._work_dir), f"{name}.pdb") + self._param_file = _os.path.join(str(self._work_dir), f"{name}.params") self._velocity_file = None self._restraint_file = None # The name of the trajectory file. - self._traj_file = "%s/%s_out.dcd" % (self._work_dir, name) + self._traj_file = _os.path.join(str(self._work_dir), f"{name}_out.dcd") # Set the path for the NAMD configuration file. - self._config_file = "%s/%s.cfg" % (self._work_dir, name) + self._config_file = _os.path.join(str(self._work_dir), f"{name}.cfg") # Create the list of input files. self._input_files = [ @@ -443,9 +443,8 @@ def _generate_config(self): p = _SireIO.PDB2(restrained._sire_object, {prop: "restrained"}) # File name for the restraint file. - self._restraint_file = "%s/%s.restrained" % ( - self._work_dir, - self._name, + self._restraint_file = _os.path.join( + str(self._work_dir), f"{self._name}.restrained" ) # Write the PDB file. @@ -733,13 +732,19 @@ def getSystem(self, block="AUTO"): has_coor = False # First check for final configuration. - if _os.path.isfile("%s/%s_out.coor" % (self._work_dir, self._name)): - coor_file = "%s/%s_out.coor" % (self._work_dir, self._name) + if _os.path.isfile( + _os.path.join(str(self._work_dir), f"{self._name}_out.coor") + ): + coor_file = _os.path.join(str(self._work_dir), f"{self._name}_out.coor") has_coor = True # Otherwise check for a restart file. - elif _os.path.isfile("%s/%s_out.restart.coor" % (self._work_dir, self._name)): - coor_file = "%s/%s_out.restart.coor" % (self._work_dir, self._name) + elif _os.path.isfile( + _os.path.join(str(self._work_dir), f"{self._name}_out.restart.coor") + ): + coor_file = _os.path.join( + str(self._work_dir), f"{self._name}_out.restart.coor" + ) has_coor = True # Try to find an XSC file. @@ -747,13 +752,17 @@ def getSystem(self, block="AUTO"): has_xsc = False # First check for final XSC file. - if _os.path.isfile("%s/%s_out.xsc" % (self._work_dir, self._name)): - xsc_file = "%s/%s_out.xsc" % (self._work_dir, self._name) + if _os.path.isfile(_os.path.join(str(self._work_dir), f"{self._name}_out.xsc")): + xsc_file = _os.path.join(str(self._work_dir), f"{self._name}_out.xsc") has_xsc = True # Otherwise check for a restart XSC file. - elif _os.path.isfile("%s/%s_out.restart.xsc" % (self._work_dir, self._name)): - xsc_file = "%s/%s_out.restart.xsc" % (self._work_dir, self._name) + elif _os.path.isfile( + _os.path.join(str(self._work_dir), f"{self._name}_out.restart.xsc") + ): + xsc_file = _os.path.join( + str(self._work_dir), f"{self._name}_out.restart.xsc" + ) has_xsc = True # We found a coordinate file. diff --git a/python/BioSimSpace/Process/_openmm.py b/python/BioSimSpace/Process/_openmm.py index d11974d22..e19cc1f40 100644 --- a/python/BioSimSpace/Process/_openmm.py +++ b/python/BioSimSpace/Process/_openmm.py @@ -174,7 +174,7 @@ def __init__( self._stdout_dict = _process._MultiDict() # Store the name of the OpenMM log file. - self._log_file = "%s/%s.log" % (self._work_dir, name) + self._log_file = _os.path.join(str(self._work_dir), f"{name}.log") # Initialise the log file separator. self._record_separator = None @@ -184,16 +184,16 @@ def __init__( # The names of the input files. We choose to use AMBER files since they # are self-contained, but could equally work with GROMACS files. - self._rst_file = "%s/%s.rst7" % (self._work_dir, name) - self._top_file = "%s/%s.prm7" % (self._work_dir, name) - self._ref_file = "%s/%s_ref.rst7" % (self._work_dir, name) + self._rst_file = _os.path.join(str(self._work_dir), f"{name}.rst7") + self._top_file = _os.path.join(str(self._work_dir), f"{name}.prm7") + self._ref_file = _os.path.join(str(self._work_dir), f"{name}_ref.rst7") # The name of the trajectory file. - self._traj_file = "%s/%s.dcd" % (self._work_dir, name) + self._traj_file = _os.path.join(str(self._work_dir), f"{name}.dcd") # Set the path for the OpenMM Python script. (We use the concept of a # config file for consistency with other Process classes.) - self._config_file = "%s/%s_script.py" % (self._work_dir, name) + self._config_file = _os.path.join(str(self._work_dir), f"{name}_script.py") # Create the list of input files. self._input_files = [self._config_file, self._rst_file, self._top_file] @@ -772,7 +772,7 @@ def _generate_config(self): ) # Copy the file into the working directory. - _shutil.copyfile(path, self._work_dir + f"/{aux_file}") + _shutil.copyfile(path, _os.path.join(str(self._work_dir), aux_file)) # The following OpenMM native implementation of the funnel metadynamics protocol # is adapted from funnel_maker.py by Dominykas Lukauskis. @@ -970,9 +970,13 @@ def _generate_config(self): # Get the number of steps to date. step = 0 - if _os.path.isfile(f"{self._work_dir}/{self._name}.xml"): - if _os.path.isfile(f"{self._work_dir}/{self._name}.log"): - with open(f"{self._work_dir}/{self._name}.log", "r") as f: + if _os.path.isfile(_os.path.join(str(self._work_dir), f"{self._name}.xml")): + if _os.path.isfile( + _os.path.join(str(self._work_dir), f"{self._name}.log") + ): + with open( + _os.path.join(str(self._work_dir), f"{self._name}.log"), "r" + ) as f: lines = f.readlines() last_line = lines[-1].split() try: @@ -2031,8 +2035,10 @@ def _add_config_restart(self): self.addToConfig("else:") self.addToConfig(" is_restart = False") - if _os.path.isfile(f"{self._work_dir}/{self._name}.xml"): - with open(f"{self._work_dir}/{self._name}.log", "r") as f: + if _os.path.isfile(_os.path.join(str(self._work_dir), f"{self._name}.xml")): + with open( + _os.path.join(str(self._work_dir), f"{self._name}.log"), "r" + ) as f: lines = f.readlines() last_line = lines[-1].split() step = int(last_line[0]) diff --git a/python/BioSimSpace/Process/_plumed.py b/python/BioSimSpace/Process/_plumed.py index aad824087..7ee712696 100644 --- a/python/BioSimSpace/Process/_plumed.py +++ b/python/BioSimSpace/Process/_plumed.py @@ -117,8 +117,8 @@ def __init__(self, work_dir): self._work_dir = work_dir # Set the location of the HILLS and COLVAR files. - self._hills_file = "%s/HILLS" % self._work_dir - self._colvar_file = "%s/COLVAR" % self._work_dir + self._hills_file = _os.path.join(str(self._work_dir), "HILLS") + self._colvar_file = _os.path.join(str(self._work_dir), "COLVAR") # The number of collective variables and total number of components. self._num_colvar = 0 @@ -250,11 +250,11 @@ def _createMetadynamicsConfig(self, system, protocol, property_map={}): # Always remove pygtail offset files. try: - _os.remove("%s/COLVAR.offset" % self._work_dir) + _os.remove(_os.path.join(str(self._work_dir), "COLVAR.offset")) except: pass try: - _os.remove("%s/HILLS.offset" % self._work_dir) + _os.remove(_os.path.join(str(self._work_dir), "HILLS.offset")) except: pass @@ -627,7 +627,9 @@ def _createMetadynamicsConfig(self, system, protocol, property_map={}): colvar_string += " TYPE=%s" % colvar.getAlignmentType().upper() # Write the reference PDB file. - with open("%s/reference.pdb" % self._work_dir, "w") as file: + with open( + _os.path.join(str(self._work_dir), "reference.pdb"), "w" + ) as file: for line in colvar.getReferencePDB(): file.write(line + "\n") @@ -870,7 +872,9 @@ def _createMetadynamicsConfig(self, system, protocol, property_map={}): metad_string += ( " GRID_WFILE=GRID GRID_WSTRIDE=%s" % protocol.getHillFrequency() ) - if is_restart and _os.path.isfile(f"{self._work_dir}/GRID"): + if is_restart and _os.path.isfile( + _os.path.join(str(self._work_dir), "GRID") + ): metad_string += " GRID_RFILE=GRID" metad_string += " CALC_RCT" @@ -940,7 +944,7 @@ def _createSteeringConfig(self, system, protocol, property_map={}): # Always remove pygtail offset files. try: - _os.remove("%s/COLVAR.offset" % self._work_dir) + _os.remove(_os.path.join(str(self._work_dir), "COLVAR.offset")) except: pass @@ -1225,7 +1229,8 @@ def _createSteeringConfig(self, system, protocol, property_map={}): # Write the reference PDB file. with open( - "%s/reference_%i.pdb" % (self._work_dir, num_rmsd), "w" + _os.path.join(str(self._work_dir), "reference_%i.pdb" % num_rmsd), + "w", ) as file: for line in colvar.getReferencePDB(): file.write(line + "\n") @@ -1449,8 +1454,8 @@ def getFreeEnergy(self, index=None, stride=None, kt=_Types.Energy(1.0, "kt")): raise ValueError("'kt' must have value > 0") # Delete any existing FES directotry and create a new one. - _shutil.rmtree(f"{self._work_dir}/fes", ignore_errors=True) - _os.makedirs(f"{self._work_dir}/fes") + _shutil.rmtree(_os.path.join(str(self._work_dir), "fes"), ignore_errors=True) + _os.makedirs(_os.path.join(str(self._work_dir), "fes")) # Create the command string. command = "%s sum_hills --hills ../HILLS --mintozero" % self._exe @@ -1466,7 +1471,7 @@ def getFreeEnergy(self, index=None, stride=None, kt=_Types.Energy(1.0, "kt")): free_energies = [] # Move to the working directory. - with _Utils.cd(self._work_dir + "/fes"): + with _Utils.cd(_os.path.join(str(self._work_dir), "fes")): # Run the sum_hills command as a background process. proc = _subprocess.run( _Utils.command_split(command), @@ -1556,7 +1561,7 @@ def getFreeEnergy(self, index=None, stride=None, kt=_Types.Energy(1.0, "kt")): _os.remove(fes) # Remove the FES output directory. - _shutil.rmtree(f"{self._work_dir}/fes", ignore_errors=True) + _shutil.rmtree(_os.path.join(str(self._work_dir), "fes"), ignore_errors=True) return tuple(free_energies) diff --git a/python/BioSimSpace/Process/_process.py b/python/BioSimSpace/Process/_process.py index a4ef0ebf2..3457e492c 100644 --- a/python/BioSimSpace/Process/_process.py +++ b/python/BioSimSpace/Process/_process.py @@ -281,11 +281,11 @@ def __init__( self._work_dir = _Utils.WorkDir(work_dir) # Files for redirection of stdout and stderr. - self._stdout_file = "%s/%s.out" % (self._work_dir, name) - self._stderr_file = "%s/%s.err" % (self._work_dir, name) + self._stdout_file = _os.path.join(str(self._work_dir), f"{name}.out") + self._stderr_file = _os.path.join(str(self._work_dir), f"{name}.err") # Files for metadynamics simulation with PLUMED. - self._plumed_config_file = "%s/plumed.dat" % self._work_dir + self._plumed_config_file = _os.path.join(str(self._work_dir), "plumed.dat") self._plumed_config = None # Initialise the configuration file string list. @@ -349,13 +349,13 @@ def _clear_output(self): self._stderr = [] # Clean up any existing offset files. - offset_files = _glob.glob("%s/*.offset" % self._work_dir) + offset_files = _glob.glob(_os.path.join(str(self._work_dir), "*.offset")) # Remove any HILLS or COLVAR files from the list. These will be dealt # with by the PLUMED interface. try: - offset_files.remove("%s/COLVAR.offset" % self._work_dir) - offset_files.remove("%s/HILLS.offset" % self._work_dir) + offset_files.remove(_os.path.join(str(self._work_dir), "COLVAR.offset")) + offset_files.remove(_os.path.join(str(self._work_dir), "HILLS.offset")) except: pass @@ -1240,7 +1240,7 @@ def getOutput(self, name=None, block="AUTO", file_link=False): zipname = "%s.zip" % name # Glob all of the output files. - output = _glob.glob("%s/*" % self._work_dir) + output = _glob.glob(_os.path.join(str(self._work_dir), "*")) with _zipfile.ZipFile(zipname, "w") as zip: # Loop over all of the file outputs. diff --git a/python/BioSimSpace/Process/_process_runner.py b/python/BioSimSpace/Process/_process_runner.py index 5c45c4f04..fa1f08cee 100644 --- a/python/BioSimSpace/Process/_process_runner.py +++ b/python/BioSimSpace/Process/_process_runner.py @@ -843,7 +843,9 @@ def _nest_directories(self, processes): # Loop over each process. for process in processes: # Create the new working directory name. - new_dir = "%s/%s" % (self._work_dir, _os.path.basename(process._work_dir)) + new_dir = _os.path.join( + self._work_dir, _os.path.basename(process._work_dir) + ) # Create a new process object using the nested directory. if process._package_name == "SOMD": diff --git a/python/BioSimSpace/Process/_somd.py b/python/BioSimSpace/Process/_somd.py index 1f3094daf..4b66b8646 100644 --- a/python/BioSimSpace/Process/_somd.py +++ b/python/BioSimSpace/Process/_somd.py @@ -227,23 +227,23 @@ def __init__( raise IOError("SOMD executable doesn't exist: '%s'" % exe) # The names of the input files. - self._rst_file = "%s/%s.rst7" % (self._work_dir, name) - self._top_file = "%s/%s.prm7" % (self._work_dir, name) + self._rst_file = _os.path.join(str(self._work_dir), f"{name}.rst7") + self._top_file = _os.path.join(str(self._work_dir), f"{name}.prm7") # The name of the trajectory file. - self._traj_file = "%s/traj000000001.dcd" % self._work_dir + self._traj_file = _os.path.join(str(self._work_dir), "traj000000001.dcd") # The name of the restart file. - self._restart_file = "%s/latest.pdb" % self._work_dir + self._restart_file = _os.path.join(str(self._work_dir), "latest.pdb") # Set the path for the SOMD configuration file. - self._config_file = "%s/%s.cfg" % (self._work_dir, name) + self._config_file = _os.path.join(str(self._work_dir), f"{name}.cfg") # Set the path for the perturbation file. - self._pert_file = "%s/%s.pert" % (self._work_dir, name) + self._pert_file = _os.path.join(str(self._work_dir), f"{name}.pert") # Set the path for the gradient file and create the gradient list. - self._gradient_file = "%s/gradients.dat" % self._work_dir + self._gradient_file = _os.path.join(str(self._work_dir), "gradients.dat") self._gradients = [] # Create the list of input files. @@ -871,26 +871,26 @@ def _clear_output(self): # Delete any restart and trajectory files in the working directory. - file = "%s/sim_restart.s3" % self._work_dir + file = _os.path.join(str(self._work_dir), "sim_restart.s3") if _os.path.isfile(file): _os.remove(file) - file = "%s/SYSTEM.s3" % self._work_dir + file = _os.path.join(str(self._work_dir), "SYSTEM.s3") if _os.path.isfile(file): _os.remove(file) - files = _glob.glob("%s/traj*.dcd" % self._work_dir) + files = _glob.glob(_os.path.join(str(self._work_dir), "traj*.dcd")) for file in files: if _os.path.isfile(file): _os.remove(file) # Additional files for free energy simulations. if isinstance(self._protocol, _Protocol.FreeEnergy): - file = "%s/gradients.dat" % self._work_dir + file = _os.path.join(str(self._work_dir), "gradients.dat") if _os.path.isfile(file): _os.remove(file) - file = "%s/simfile.dat" % self._work_dir + file = _os.path.join(str(self._work_dir), "simfile.dat") if _os.path.isfile(file): _os.remove(file) diff --git a/python/BioSimSpace/Protocol/_metadynamics.py b/python/BioSimSpace/Protocol/_metadynamics.py index 99b49320e..8e6b25e82 100644 --- a/python/BioSimSpace/Protocol/_metadynamics.py +++ b/python/BioSimSpace/Protocol/_metadynamics.py @@ -48,6 +48,7 @@ def __init__( runtime=_Types.Time(1, "nanosecond"), temperature=_Types.Temperature(300, "kelvin"), pressure=_Types.Pressure(1, "atmosphere"), + thermostat_time_constant=_Types.Time(1, "picosecond"), hill_height=_Types.Energy(1, "kj per mol"), hill_frequency=1000, report_interval=1000, @@ -76,6 +77,9 @@ def __init__( pressure : :class:`Pressure ` The pressure. Pass pressure=None to use the NVT ensemble. + thermostat_time_constant : :class:`Time ` + Time constant for thermostat coupling. + hill_height : :class:`Energy ` The height of the Gaussian hills. @@ -123,6 +127,9 @@ def __init__( else: self._pressure = None + # Set the thermostat time constant. + self.setThermostatTimeConstant(thermostat_time_constant) + # Set the hill parameters: height, frequency. self.setHillHeight(hill_height) self.setHillFrequency(hill_frequency) @@ -384,6 +391,42 @@ def setPressure(self, pressure): "'pressure' must be of type 'str' or 'BioSimSpace.Types.Pressure'" ) + def getThermostatTimeConstant(self): + """ + Return the time constant for the thermostat. + + Returns + ------- + + runtime : :class:`Time ` + The time constant for the thermostat. + """ + return self._thermostat_time_constant + + def setThermostatTimeConstant(self, thermostat_time_constant): + """ + Set the time constant for the thermostat. + + Parameters + ---------- + + thermostat_time_constant : str, :class:`Time ` + The time constant for the thermostat. + """ + if isinstance(thermostat_time_constant, str): + try: + self._thermostat_time_constant = _Types.Time(thermostat_time_constant) + except: + raise ValueError( + "Unable to parse 'thermostat_time_constant' string." + ) from None + elif isinstance(thermostat_time_constant, _Types.Time): + self._thermostat_time_constant = thermostat_time_constant + else: + raise TypeError( + "'thermostat_time_constant' must be of type 'str' or 'BioSimSpace.Types.Time'" + ) + def getHillHeight(self): """ Return the height of the Gaussian hills. diff --git a/python/BioSimSpace/Protocol/_position_restraint_mixin.py b/python/BioSimSpace/Protocol/_position_restraint_mixin.py index 8374bf8ad..3211558d1 100644 --- a/python/BioSimSpace/Protocol/_position_restraint_mixin.py +++ b/python/BioSimSpace/Protocol/_position_restraint_mixin.py @@ -110,13 +110,14 @@ def __eq__(self, other): ) def getRestraint(self): - """Return the type of restraint.. + """ + Return the type of restraint. Returns ------- restraint : str, [int] - The type of restraint. + The type of restraint, either a keyword or a list of atom indices. """ return self._restraint diff --git a/python/BioSimSpace/Protocol/_protocol.py b/python/BioSimSpace/Protocol/_protocol.py index 9d37e004f..ad52a55e1 100644 --- a/python/BioSimSpace/Protocol/_protocol.py +++ b/python/BioSimSpace/Protocol/_protocol.py @@ -40,6 +40,23 @@ def __init__(self): # Flag that the protocol hasn't been customised. self._is_customised = False + def getRestraint(self): + """ + Return the type of restraint. + + Returns + ------- + + restraint : str, [int] + The type of restraint, either a keyword or a list of atom indices. + """ + from ._position_restraint_mixin import _PositionRestraintMixin + + if isinstance(self, _PositionRestraintMixin): + return self._restraint + else: + return None + def _setCustomised(self, is_customised): """ Internal function to flag whether a protocol has been customised. diff --git a/python/BioSimSpace/Sandpit/Exscientia/Align/_alch_ion.py b/python/BioSimSpace/Sandpit/Exscientia/Align/_alch_ion.py index 14773adf6..cd75db2a3 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/Align/_alch_ion.py +++ b/python/BioSimSpace/Sandpit/Exscientia/Align/_alch_ion.py @@ -1,6 +1,6 @@ import warnings -from .._SireWrappers import Molecule as _Molecule +from .._SireWrappers import Molecule as _Molecule, System as _System def _mark_alchemical_ion(molecule): @@ -50,3 +50,23 @@ def _mark_alchemical_ion(molecule): mol._sire_object = mol_edit.commit() return mol + + +def _get_protein_com_idx(system: _System) -> int: + """return the index of the atom that is closest to the center of + mass of the biggest molecule in the system. + + Args: + system: The input system. + + Returns: + atom_index + """ + biggest_mol_idx = max(range(system.nMolecules()), key=lambda x: system[x].nAtoms()) + + atom_offset = 0 + for i, mol in enumerate(system): + if biggest_mol_idx == i: + return atom_offset + mol.getCOMIdx() + else: + atom_offset += mol.nAtoms() diff --git a/python/BioSimSpace/Sandpit/Exscientia/FreeEnergy/_restraint.py b/python/BioSimSpace/Sandpit/Exscientia/FreeEnergy/_restraint.py index ab961ddb2..54421eed9 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/FreeEnergy/_restraint.py +++ b/python/BioSimSpace/Sandpit/Exscientia/FreeEnergy/_restraint.py @@ -22,29 +22,54 @@ """ A class for holding restraints. """ -from math import e -from scipy import integrate as _integrate -import numpy as _np +import math as _math import warnings as _warnings +from typing import Literal -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 - -from .._SireWrappers import Atom as _Atom -from .._SireWrappers import System as _System -from ..Types import Angle as _Angle -from ..Types import Length as _Length -from ..Types import Temperature as _Temperature -from ..Units.Angle import degree as _degree -from ..Units.Angle import radian as _radian +import numpy as _np +from scipy import integrate as _integrate +from scipy.special import erf as _erf +from sire.legacy.Units import ( + angstrom3 as _Sire_angstrom3, + k_boltz as _k_boltz, + meter3 as _Sire_meter3, + mole as _Sire_mole, +) +from sire.units import GeneralUnit as _sire_GeneralUnit + +from BioSimSpace.Sandpit.Exscientia.Types._general_unit import ( + GeneralUnit as _GeneralUnit, +) +from ..Types import Angle as _Angle, Length as _Length, Temperature as _Temperature +from ..Units.Angle import degree as _degree, radian as _radian from ..Units.Area import angstrom2 as _angstrom2 -from ..Units.Energy import kcal_per_mol as _kcal_per_mol -from ..Units.Energy import kj_per_mol as _kj_per_mol -from ..Units.Length import angstrom as _angstrom -from ..Units.Length import nanometer as _nanometer +from ..Units.Energy import kcal_per_mol as _kcal_per_mol, kj_per_mol as _kj_per_mol +from ..Units.Length import angstrom as _angstrom, nanometer as _nanometer from ..Units.Temperature import kelvin as _kelvin +from ..Units.Volume import angstrom3 as _angstrom3 +from .._SireWrappers import Atom as _Atom, System as _System + + +def sqrt(u): + dims = u._sire_unit.dimensions() + for dim in dims: + if dim % 2 != 0: + raise ValueError( + "Square root not possible on dimension that is not divisible by 2!" + ) + return _GeneralUnit( + _sire_GeneralUnit(_math.sqrt(u.value()), [int(0.5 * dim) for dim in dims]) + ) + + +def exp(u): + dims = u._sire_unit.dimensions() + return _GeneralUnit(_sire_GeneralUnit(_math.exp(u.value()), dims)) + + +def erf(u): + dims = u._sire_unit.dimensions() + return _GeneralUnit(_sire_GeneralUnit(_erf(u.value()), dims)) class Restraint: @@ -380,7 +405,7 @@ def system(self, system): # Store a copy of solvated system. self._system = system.copy() - def _gromacs_boresch(self, perturbation_type=None): + def _gromacs_boresch(self, perturbation_type=None, restraint_lambda=False): """Format the Gromacs string for boresch restraint.""" # Format the atoms into index list @@ -398,9 +423,18 @@ def format_index(key_list): return " ".join(formated_index) parameters_string = "{eq0:<10} {fc0:<10} {eq1:<10} {fc1:<10}" + # The Gromacs dihedral restraints has the format of + # phi dphi fc and we don't want dphi for the restraint, it is hence zero + dihedral_restraints_parameters_string = ( + "{eq0:<10} 0.00 {fc0:<10} {eq1:<10} 0.00 {fc1:<10}" + ) # Format the parameters for the bonds def format_bond(equilibrium_values, force_constants): + """ + Format the bonds equilibrium values and force constant + in into the Gromacs topology format. + """ converted_equ_val = ( self._restraint_dict["equilibrium_values"][equilibrium_values] / _nanometer @@ -416,14 +450,33 @@ def format_bond(equilibrium_values, force_constants): ) # Format the parameters for the angles and dihedrals - def format_angle(equilibrium_values, force_constants): + def format_angle(equilibrium_values, force_constants, restraint_lambda): + """ + Format the angle equilibrium values and force constant + in into the Gromacs topology format. + + For Boresch restraint, we might want the dihedral to be stored + under the [ dihedral_restraints ] and controlled by restraint-lambdas. + Instead of under the [ dihedrals ] directive and controlled by bonded-lambdas. + + However, for dihedrals, the [ dihedral_restraints ] has a different function type + compared with [ dihedrals ] and more values for the force constant, so we need + to format them differently. + + When restraint_lambda is True, the dihedrals will be stored in the dihedral_restraints. + """ converted_equ_val = ( self._restraint_dict["equilibrium_values"][equilibrium_values] / _degree ) converted_fc = self._restraint_dict["force_constants"][force_constants] / ( _kj_per_mol / (_radian * _radian) ) - return parameters_string.format( + par_string = ( + dihedral_restraints_parameters_string + if restraint_lambda + else parameters_string + ) + return par_string.format( eq0="{:.3f}".format(converted_equ_val), fc0="{:.2f}".format(0), eq1="{:.3f}".format(converted_equ_val), @@ -444,14 +497,28 @@ def write_angle(key_list, equilibrium_values, force_constants): return master_string.format( index=format_index(key_list), func_type=1, - parameters=format_angle(equilibrium_values, force_constants), + parameters=format_angle( + equilibrium_values, force_constants, restraint_lambda=False + ), ) - def write_dihedral(key_list, equilibrium_values, force_constants): + def write_dihedral( + key_list, equilibrium_values, force_constants, restraint_lambda + ): + if restraint_lambda: + # In [ dihedral_restraints ], function type 1 + # means the dihedral is restrained harmonically. + func_type = 1 + else: + # In [ dihedrals ], function type 2 + # means the dihedral is restrained harmonically. + func_type = 2 return master_string.format( index=format_index(key_list), - func_type=2, - parameters=format_angle(equilibrium_values, force_constants), + func_type=func_type, + parameters=format_angle( + equilibrium_values, force_constants, restraint_lambda + ), ) # Writing the string @@ -473,16 +540,43 @@ def write_dihedral(key_list, equilibrium_values, force_constants): # Angles: r1-l1-l2 (thetaB0, kthetaB) output.append(write_angle(("r1", "l1", "l2"), "thetaB0", "kthetaB")) - output.append("[ dihedrals ]") + if restraint_lambda: + output.append("[ dihedral_restraints ]") + output.append( + "; ai aj ak al type phiA dphiA fcA phiB dphiB fcB" + ) + else: + output.append("[ dihedrals ]") + output.append( + "; ai aj ak al type phiA fcA phiB fcB" + ) + # Dihedrals: r3-r2-r1-l1 (phiA0, kphiA) output.append( - "; ai aj ak al type phiA fcA phiB fcB" + write_dihedral( + ("r3", "r2", "r1", "l1"), + "phiA0", + "kphiA", + restraint_lambda=restraint_lambda, + ) ) - # Dihedrals: r3-r2-r1-l1 (phiA0, kphiA) - output.append(write_dihedral(("r3", "r2", "r1", "l1"), "phiA0", "kphiA")) # Dihedrals: r2-r1-l1-l2 (phiB0, kphiB) - output.append(write_dihedral(("r2", "r1", "l1", "l2"), "phiB0", "kphiB")) + output.append( + write_dihedral( + ("r2", "r1", "l1", "l2"), + "phiB0", + "kphiB", + restraint_lambda=restraint_lambda, + ) + ) # Dihedrals: r1-l1-l2-l3 (phiC0, kphiC) - output.append(write_dihedral(("r1", "l1", "l2", "l3"), "phiC0", "kphiC")) + output.append( + write_dihedral( + ("r1", "l1", "l2", "l3"), + "phiC0", + "kphiC", + restraint_lambda=restraint_lambda, + ) + ) return "\n".join(output) @@ -702,7 +796,7 @@ def _add_restr_to_str(restr, restr_string): ) return standard_restr_string + permanent_restr_string[:-2] + "}" - def toString(self, engine, perturbation_type=None): + def toString(self, engine, perturbation_type=None, restraint_lambda=False): """ The method for convert the restraint to a format that could be used by MD Engines. @@ -721,27 +815,34 @@ def toString(self, engine, perturbation_type=None): turned on when the perturbation type is "restraint", but for which the permanent distance restraint is always active if the perturbation type is "release_restraint" (or any other perturbation type). + restraint_lambda : str, optional, default=False + Whether to use restraint_lambda in Gromacs, this would move the dihedral restraints + from [ dihedrals ], which is controlled by the bonded-lambda to + [ dihedral_restraints ], which is controlled by restraint-lambda. """ - to_str_functions = { - "boresch": {"gromacs": self._gromacs_boresch, "somd": self._somd_boresch}, - "multiple_distance": { - "gromacs": self._gromacs_multiple_distance, - "somd": self._somd_multiple_distance, - }, - } - engine = engine.strip().lower() - try: - str_fn = to_str_functions[self._restraint_type][engine] - except KeyError: - raise NotImplementedError( - f"Restraint type {self._restraint_type} not implemented " - f"yet for {engine}." - ) - - return str_fn(perturbation_type) + match (self._restraint_type, engine): + case "boresch", "gromacs": + return self._gromacs_boresch( + perturbation_type, restraint_lambda=restraint_lambda + ) + case "boresch", "somd": + return self._somd_boresch(perturbation_type) + case "multiple_distance", "gromacs": + return self._gromacs_multiple_distance(perturbation_type) + case "multiple_distance", "somd": + return self._somd_multiple_distance(perturbation_type) + case _: + raise NotImplementedError( + f"Restraint type {self._restraint_type} not implemented " + f"yet for {engine}." + ) - def getCorrection(self, method="analytical"): + def getCorrection( + self, + method="analytical", + flavour: Literal["boresch", "schrodinger"] = "boresch", + ): """ Calculate the free energy of releasing the restraint to the standard state volume. @@ -755,6 +856,11 @@ def getCorrection(self, method="analytical"): correction can introduce errors when the restraints are weak, restrained angles are close to 0 or pi radians, or the restrained distance is close to 0. + flavour : str + When analytical correction is used, one could either use + Boresch's derivation or Schrodinger's derivation. Both of + them usually agrees quite well with each other to the extent + of 0.2 kcal/mol. Returns ---------- @@ -903,59 +1009,10 @@ def numerical_dihedral_integrand(phi, phi0, kphi): return dg elif method == "analytical": - # Only need three equilibrium values for the analytical correction - r0 = ( - self._restraint_dict["equilibrium_values"]["r0"] / _angstrom - ) # Distance in A - thetaA0 = ( - self._restraint_dict["equilibrium_values"]["thetaA0"] / _radian - ) # Angle in radians - thetaB0 = ( - self._restraint_dict["equilibrium_values"]["thetaB0"] / _radian - ) # Angle in radians - - force_constants = [] - - # Loop through and correct for force constants of zero, - # which break the analytical correction. To account for this, - # divide the prefactor accordingly. Note that setting - # certain force constants to zero while others are non-zero - # will result in unstable restraints, but this will be checked when - # the restraint object is created - for k, val in self._restraint_dict["force_constants"].items(): - if val.value() == 0: - if k == "kr": - raise ValueError("The force constant kr must not be zero") - if k == "kthetaA": - prefactor /= 2 / _np.sin(thetaA0) - if k == "kthetaB": - prefactor /= 2 / _np.sin(thetaB0) - if k[:4] == "kphi": - prefactor /= 2 * _np.pi - else: - if k == "kr": - force_constants.append(val / (_kcal_per_mol / _angstrom2)) - else: - force_constants.append( - val / (_kcal_per_mol / (_radian * _radian)) - ) - - # Calculation - n_nonzero_k = len(force_constants) - prod_force_constants = _np.prod(force_constants) - numerator = prefactor * _np.sqrt(prod_force_constants) - denominator = ( - (r0**2) - * _np.sin(thetaA0) - * _np.sin(thetaB0) - * (2 * _np.pi * R * T) ** (n_nonzero_k / 2) - ) - - # Compute dg and attach unit - dg = -R * T * _np.log(numerator / denominator) - dg *= _kcal_per_mol - - return dg + if flavour.lower() == "schrodinger": + return self._schrodinger_analytical_correction() + elif flavour.lower() == "boresch": + return self._boresch_analytical_correction() else: raise ValueError( @@ -1055,6 +1112,125 @@ def _get_correction(r0, r_fb, kr): ) return _get_correction(r0, r_fb, kr) + def _schrodinger_analytical_correction(self): + # Adapted from DOI: 10.1021/acs.jcim.3c00013 + k_boltz = _GeneralUnit(_k_boltz) + beta = 1 / (k_boltz * self.T) + V = 1660 * _angstrom3 + + r = self._restraint_dict["equilibrium_values"]["r0"] + # Schrodinger uses k(b-b0)**2 + kr = self._restraint_dict["force_constants"]["kr"] / 2 + + Z_dist = r / (2 * beta * kr) * _np.exp(-beta * kr * r**2) + _np.sqrt( + _np.pi + ) / (4 * beta * kr * sqrt(beta * kr)) * (1 + 2 * beta * kr * r**2) * ( + 1 + _erf(sqrt(beta * kr) * r) + ) + + Z_angles = [] + for angle in ["A", "B"]: + theta = ( + self._restraint_dict["equilibrium_values"][f"theta{angle}0"] / _radian + ) # Angle in radians + # Schrodinger uses k instead of k/2 + ktheta = self._restraint_dict["force_constants"][f"ktheta{angle}"] / 2 + Z_angle = ( + sqrt(_np.pi / (beta * ktheta)) + * exp(-1 / (4 * beta * ktheta)) + * _np.sin(theta) + ) + Z_angle /= _radian**3 + Z_angles.append(Z_angle) + + Z_dihedrals = [] + for dihedral in ["A", "B", "C"]: + # Schrodinger uses k instead of k/2 + kphi = self._restraint_dict["force_constants"][f"kphi{dihedral}"] / 2 + Z_dihedral = sqrt(_np.pi / (beta * kphi)) * erf(_np.pi * sqrt(beta * kphi)) + Z_dihedrals.append(Z_dihedral) + + dG = ( + k_boltz + * self.T + * _np.log( + Z_angles[0] + * Z_angles[1] + * Z_dist + * Z_dihedrals[0] + * Z_dihedrals[1] + * Z_dihedrals[2] + / (8 * _np.pi**2 * V) + ) + ) + return dG + + def _boresch_analytical_correction(self): + R = ( + _k_boltz.value() * _kcal_per_mol / _kelvin + ).value() # molar gas constant in kcal mol-1 K-1 + + # Parameters + T = self.T / _kelvin # Temperature in Kelvin + v0 = ( + ((_Sire_meter3 / 1000) / _Sire_mole) / _Sire_angstrom3 + ).value() # standard state volume in A^3 + prefactor = ( + 8 * (_np.pi**2) * v0 + ) # In A^3. Divide this to account for force constants of 0 in the + # analytical correction + # Only need three equilibrium values for the analytical correction + r0 = ( + self._restraint_dict["equilibrium_values"]["r0"] / _angstrom + ) # Distance in A + thetaA0 = ( + self._restraint_dict["equilibrium_values"]["thetaA0"] / _radian + ) # Angle in radians + thetaB0 = ( + self._restraint_dict["equilibrium_values"]["thetaB0"] / _radian + ) # Angle in radians + + force_constants = [] + + # Loop through and correct for force constants of zero, + # which break the analytical correction. To account for this, + # divide the prefactor accordingly. Note that setting + # certain force constants to zero while others are non-zero + # will result in unstable restraints, but this will be checked when + # the restraint object is created + for k, val in self._restraint_dict["force_constants"].items(): + if val.value() == 0: + if k == "kr": + raise ValueError("The force constant kr must not be zero") + if k == "kthetaA": + prefactor /= 2 / _np.sin(thetaA0) + if k == "kthetaB": + prefactor /= 2 / _np.sin(thetaB0) + if k[:4] == "kphi": + prefactor /= 2 * _np.pi + else: + if k == "kr": + force_constants.append(val / (_kcal_per_mol / _angstrom2)) + else: + force_constants.append(val / (_kcal_per_mol / (_radian * _radian))) + + # Calculation + n_nonzero_k = len(force_constants) + prod_force_constants = _np.prod(force_constants) + numerator = prefactor * _np.sqrt(prod_force_constants) + denominator = ( + (r0**2) + * _np.sin(thetaA0) + * _np.sin(thetaB0) + * (2 * _np.pi * R * T) ** (n_nonzero_k / 2) + ) + + # Compute dg and attach unit + dg = -R * T * _np.log(numerator / denominator) + dg *= _kcal_per_mol + + return dg + @property def correction(self): """Give the free energy of removing the restraint.""" diff --git a/python/BioSimSpace/Sandpit/Exscientia/Parameters/_Protocol/_amber.py b/python/BioSimSpace/Sandpit/Exscientia/Parameters/_Protocol/_amber.py index 69816c1e4..2ce0b15f9 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/Parameters/_Protocol/_amber.py +++ b/python/BioSimSpace/Sandpit/Exscientia/Parameters/_Protocol/_amber.py @@ -316,9 +316,6 @@ def run(self, molecule, work_dir=None, queue=None): else: is_smiles = False - # Create the file prefix. - prefix = work_dir + "/" - if not is_smiles: # Create a copy of the molecule. new_mol = molecule.copy() @@ -352,7 +349,10 @@ def run(self, molecule, work_dir=None, queue=None): ) # Prepend the working directory to the output file names. - output = [prefix + output[0], prefix + output[1]] + output = [ + _os.path.join(str(work_dir), output[0]), + _os.path.join(str(work_dir), output[1]), + ] try: # Load the parameterised molecule. (This could be a system of molecules.) @@ -443,9 +443,6 @@ def _run_tleap(self, molecule, work_dir): else: _molecule = molecule - # Create the file prefix. - prefix = work_dir + "/" - # Write the system to a PDB file. try: # LEaP expects residue numbering to be ascending and continuous. @@ -454,7 +451,7 @@ def _run_tleap(self, molecule, work_dir): )[0] renumbered_molecule = _Molecule(renumbered_molecule) _IO.saveMolecules( - prefix + "leap", + _os.path.join(str(work_dir), "leap"), renumbered_molecule, "pdb", property_map=self._property_map, @@ -500,7 +497,7 @@ def _run_tleap(self, molecule, work_dir): pruned_bond_records.append(bond) # Write the LEaP input file. - with open(prefix + "leap.txt", "w") as file: + with open(_os.path.join(str(work_dir), "leap.txt"), "w") as file: file.write("source %s\n" % ff) if self._water_model is not None: if self._water_model in ["tip4p", "tip5p"]: @@ -528,14 +525,14 @@ def _run_tleap(self, molecule, work_dir): # Generate the tLEaP command. command = "%s -f leap.txt" % _tleap_exe - with open(prefix + "README.txt", "w") as file: + with open(_os.path.join(str(work_dir), "README.txt"), "w") as file: # Write the command to file. file.write("# tLEaP was run with the following command:\n") file.write("%s\n" % command) # Create files for stdout/stderr. - stdout = open(prefix + "leap.out", "w") - stderr = open(prefix + "leap.err", "w") + stdout = open(_os.path.join(str(work_dir), "leap.out"), "w") + stderr = open(_os.path.join(str(work_dir), "leap.err"), "w") # Run tLEaP as a subprocess. proc = _subprocess.run( @@ -550,12 +547,12 @@ def _run_tleap(self, molecule, work_dir): # tLEaP doesn't return sensible error codes, so we need to check that # the expected output was generated. - if _os.path.isfile(prefix + "leap.top") and _os.path.isfile( - prefix + "leap.crd" - ): + if _os.path.isfile( + _os.path.join(str(work_dir), "leap.top") + ) and _os.path.isfile(_os.path.join(str(work_dir), "leap.crd")): # Check the output of tLEaP for missing atoms. if self._ensure_compatible: - if _has_missing_atoms(prefix + "leap.out"): + if _has_missing_atoms(_os.path.join(str(work_dir), "leap.out")): raise _ParameterisationError( "tLEaP added missing atoms. The topology is now " "inconsistent with the original molecule. Please " @@ -604,13 +601,13 @@ def _run_pdb2gmx(self, molecule, work_dir): else: _molecule = molecule - # Create the file prefix. - prefix = work_dir + "/" - # Write the system to a PDB file. try: _IO.saveMolecules( - prefix + "leap", _molecule, "pdb", property_map=self._property_map + _os.path.join(str(work_dir), "input"), + _molecule, + "pdb", + property_map=self._property_map, ) except Exception as e: msg = "Failed to write system to 'PDB' format." @@ -626,14 +623,14 @@ def _run_pdb2gmx(self, molecule, work_dir): % (_gmx_exe, supported_ff[self._forcefield]) ) - with open(prefix + "README.txt", "w") as file: + with open(_os.path.join(str(work_dir), "README.txt"), "w") as file: # Write the command to file. file.write("# pdb2gmx was run with the following command:\n") file.write("%s\n" % command) # Create files for stdout/stderr. - stdout = open(prefix + "pdb2gmx.out", "w") - stderr = open(prefix + "pdb2gmx.err", "w") + stdout = open(_os.path.join(str(work_dir), "pdb2gmx.out"), "w") + stderr = open(_os.path.join(str(work_dir), "pdb2gmx.err"), "w") # Run pdb2gmx as a subprocess. proc = _subprocess.run( @@ -647,9 +644,9 @@ def _run_pdb2gmx(self, molecule, work_dir): stderr.close() # Check for the expected output. - if _os.path.isfile(prefix + "output.gro") and _os.path.isfile( - prefix + "output.top" - ): + if _os.path.isfile( + _os.path.join(str(work_dir), "output.gro") + ) and _os.path.isfile(_os.path.join(str(work_dir), "output.top")): return ["output.gro", "output.top"] else: raise _ParameterisationError("pdb2gmx failed!") @@ -1010,9 +1007,6 @@ def run(self, molecule, work_dir=None, queue=None): if work_dir is None: work_dir = _os.getcwd() - # Create the file prefix. - prefix = work_dir + "/" - # Convert SMILES to a molecule. if isinstance(molecule, str): is_smiles = True @@ -1092,7 +1086,10 @@ def run(self, molecule, work_dir=None, queue=None): # Write the system to a PDB file. try: _IO.saveMolecules( - prefix + "antechamber", new_mol, "pdb", property_map=self._property_map + _os.path.join(str(work_dir), "antechamber"), + new_mol, + "pdb", + property_map=self._property_map, ) except Exception as e: msg = "Failed to write system to 'PDB' format." @@ -1108,14 +1105,14 @@ def run(self, molecule, work_dir=None, queue=None): + "-o antechamber.mol2 -fo mol2 -c %s -s 2 -nc %d" ) % (_antechamber_exe, self._version, self._charge_method.lower(), charge) - with open(prefix + "README.txt", "w") as file: + 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") file.write("%s\n" % command) # Create files for stdout/stderr. - stdout = open(prefix + "antechamber.out", "w") - stderr = open(prefix + "antechamber.err", "w") + stdout = open(_os.path.join(str(work_dir), "antechamber.out"), "w") + stderr = open(_os.path.join(str(work_dir), "antechamber.err"), "w") # Run Antechamber as a subprocess. proc = _subprocess.run( @@ -1130,20 +1127,20 @@ def run(self, molecule, work_dir=None, queue=None): # Antechamber doesn't return sensible error codes, so we need to check that # the expected output was generated. - if _os.path.isfile(prefix + "antechamber.mol2"): + if _os.path.isfile(_os.path.join(str(work_dir), "antechamber.mol2")): # Run parmchk to check for missing parameters. command = ( "%s -s %d -i antechamber.mol2 -f mol2 " + "-o antechamber.frcmod" ) % (_parmchk_exe, self._version) - with open(prefix + "README.txt", "a") as file: + with open(_os.path.join(str(work_dir), "README.txt"), "a") as file: # Write the command to file. file.write("\n# ParmChk was run with the following command:\n") file.write("%s\n" % command) # Create files for stdout/stderr. - stdout = open(prefix + "parmchk.out", "w") - stderr = open(prefix + "parmchk.err", "w") + stdout = open(_os.path.join(str(work_dir), "parmchk.out"), "w") + stderr = open(_os.path.join(str(work_dir), "parmchk.err"), "w") # Run parmchk as a subprocess. proc = _subprocess.run( @@ -1157,7 +1154,7 @@ def run(self, molecule, work_dir=None, queue=None): stderr.close() # The frcmod file was created. - if _os.path.isfile(prefix + "antechamber.frcmod"): + if _os.path.isfile(_os.path.join(str(work_dir), "antechamber.frcmod")): # Now call tLEaP using the partially parameterised molecule and the frcmod file. # tLEap will run in the same working directory, using the Mol2 file generated by # Antechamber. @@ -1169,7 +1166,7 @@ def run(self, molecule, work_dir=None, queue=None): ff = _find_force_field("gaff2") # Write the LEaP input file. - with open(prefix + "leap.txt", "w") as file: + with open(_os.path.join(str(work_dir), "leap.txt"), "w") as file: file.write("source %s\n" % ff) file.write("mol = loadMol2 antechamber.mol2\n") file.write("loadAmberParams antechamber.frcmod\n") @@ -1179,14 +1176,14 @@ def run(self, molecule, work_dir=None, queue=None): # Generate the tLEaP command. command = "%s -f leap.txt" % _tleap_exe - with open(prefix + "README.txt", "a") as file: + with open(_os.path.join(str(work_dir), "README.txt"), "a") as file: # Write the command to file. file.write("\n# tLEaP was run with the following command:\n") file.write("%s\n" % command) # Create files for stdout/stderr. - stdout = open(prefix + "leap.out", "w") - stderr = open(prefix + "leap.err", "w") + stdout = open(_os.path.join(str(work_dir), "leap.out"), "w") + stderr = open(_os.path.join(str(work_dir), "leap.err"), "w") # Run tLEaP as a subprocess. proc = _subprocess.run( @@ -1201,12 +1198,12 @@ def run(self, molecule, work_dir=None, queue=None): # tLEaP doesn't return sensible error codes, so we need to check that # the expected output was generated. - if _os.path.isfile(prefix + "leap.top") and _os.path.isfile( - prefix + "leap.crd" - ): + if _os.path.isfile( + _os.path.join(str(work_dir), "leap.top") + ) and _os.path.isfile(_os.path.join(str(work_dir), "leap.crd")): # Check the output of tLEaP for missing atoms. if self._ensure_compatible: - if _has_missing_atoms(prefix + "leap.out"): + if _has_missing_atoms(_os.path.join(str(work_dir), "leap.out")): raise _ParameterisationError( "tLEaP added missing atoms. The topology is now " "inconsistent with the original molecule. Please " @@ -1217,7 +1214,10 @@ def run(self, molecule, work_dir=None, queue=None): # Load the parameterised molecule. (This could be a system of molecules.) try: par_mol = _IO.readMolecules( - [prefix + "leap.top", prefix + "leap.crd"] + [ + _os.path.join(str(work_dir), "leap.top"), + _os.path.join(str(work_dir), "leap.crd"), + ], ) # Extract single molecules. if par_mol.nMolecules() == 1: diff --git a/python/BioSimSpace/Sandpit/Exscientia/Parameters/_Protocol/_openforcefield.py b/python/BioSimSpace/Sandpit/Exscientia/Parameters/_Protocol/_openforcefield.py index 018f43e4e..c94e7d349 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/Parameters/_Protocol/_openforcefield.py +++ b/python/BioSimSpace/Sandpit/Exscientia/Parameters/_Protocol/_openforcefield.py @@ -214,9 +214,6 @@ def run(self, molecule, work_dir=None, queue=None): if work_dir is None: work_dir = _os.getcwd() - # Create the file prefix. - prefix = work_dir + "/" - # Flag whether the molecule is a SMILES string. if isinstance(molecule, str): is_smiles = True @@ -256,7 +253,7 @@ def run(self, molecule, work_dir=None, queue=None): # Write the molecule to SDF format. try: _IO.saveMolecules( - prefix + "molecule", + _os.path.join(str(work_dir), "molecule"), molecule, "sdf", property_map=self._property_map, @@ -275,7 +272,7 @@ def run(self, molecule, work_dir=None, queue=None): # Write the molecule to a PDB file. try: _IO.saveMolecules( - prefix + "molecule", + _os.path.join(str(work_dir), "molecule"), molecule, "pdb", property_map=self._property_map, @@ -291,7 +288,7 @@ def run(self, molecule, work_dir=None, queue=None): # Create an RDKit molecule from the PDB file. try: rdmol = _Chem.MolFromPDBFile( - prefix + "molecule.pdb", removeHs=False + _os.path.join(str(work_dir), "molecule.pdb"), removeHs=False ) except Exception as e: msg = "RDKit was unable to read the molecular PDB file!" @@ -303,7 +300,9 @@ def run(self, molecule, work_dir=None, queue=None): # Use RDKit to write back to SDF format. try: - writer = _Chem.SDWriter(prefix + "molecule.sdf") + writer = _Chem.SDWriter( + _os.path.join(str(work_dir), "molecule.sdf") + ) writer.write(rdmol) writer.close() except Exception as e: @@ -317,7 +316,9 @@ def run(self, molecule, work_dir=None, queue=None): # Create the Open Forcefield Molecule from the intermediate SDF file, # as recommended by @j-wags and @mattwthompson. try: - off_molecule = _OpenFFMolecule.from_file(prefix + "molecule.sdf") + off_molecule = _OpenFFMolecule.from_file( + _os.path.join(str(work_dir), "molecule.sdf") + ) except Exception as e: msg = "Unable to create OpenFF Molecule!" if _isVerbose(): @@ -383,8 +384,8 @@ def run(self, molecule, work_dir=None, queue=None): # Export AMBER format files. try: - interchange.to_prmtop(prefix + "interchange.prm7") - interchange.to_inpcrd(prefix + "interchange.rst7") + interchange.to_prmtop(_os.path.join(str(work_dir), "interchange.prm7")) + interchange.to_inpcrd(_os.path.join(str(work_dir), "interchange.rst7")) except Exception as e: msg = "Unable to write Interchange object to AMBER format!" if _isVerbose(): @@ -396,7 +397,10 @@ def run(self, molecule, work_dir=None, queue=None): # Load the parameterised molecule. (This could be a system of molecules.) try: par_mol = _IO.readMolecules( - [prefix + "interchange.prm7", prefix + "interchange.rst7"] + [ + _os.path.join(str(work_dir), "interchange.prm7"), + _os.path.join(str(work_dir), "interchange.rst7"), + ], ) # Extract single molecules. if par_mol.nMolecules() == 1: diff --git a/python/BioSimSpace/Sandpit/Exscientia/Process/_amber.py b/python/BioSimSpace/Sandpit/Exscientia/Process/_amber.py index 257deef09..64c7a173f 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/Process/_amber.py +++ b/python/BioSimSpace/Sandpit/Exscientia/Process/_amber.py @@ -44,7 +44,6 @@ from sire.legacy import Base as _SireBase from sire.legacy import IO as _SireIO from sire.legacy import Mol as _SireMol -import pandas as pd from .._Utils import _assert_imported, _have_imported, _try_import @@ -236,8 +235,12 @@ def _setup(self): ) # Create the reference file - if self._ref_system is not None and self._protocol.getRestraint() is not None: - self._write_system(self._ref_system, ref_file=self._ref_file) + if self._ref_system is not None: + if ( + self._system.getAlchemicalIon() + or self._protocol.getRestraint() is not None + ): + self._write_system(self._ref_system, ref_file=self._ref_file) else: _shutil.copy(self._rst_file, self._ref_file) @@ -543,7 +546,10 @@ def _generate_args(self): if not isinstance(self._protocol, _Protocol.Custom): # Append a reference file if this a restrained simulation. if isinstance(self._protocol, _Protocol._PositionRestraintMixin): - if self._protocol.getRestraint() is not None: + if ( + self._protocol.getRestraint() is not None + or self._system.getAlchemicalIon() + ): self.setArg("-ref", "%s_ref.rst7" % self._name) # Append a trajectory file if this anything other than a minimisation. diff --git a/python/BioSimSpace/Sandpit/Exscientia/Process/_gromacs.py b/python/BioSimSpace/Sandpit/Exscientia/Process/_gromacs.py index 6c17685f4..788450c84 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/Process/_gromacs.py +++ b/python/BioSimSpace/Sandpit/Exscientia/Process/_gromacs.py @@ -28,9 +28,6 @@ import glob as _glob import os as _os -import warnings as _warnings - -import pandas as pd from .._Utils import _try_import @@ -261,8 +258,12 @@ def _setup(self): ) # Create the reference file - if self._ref_system is not None and self._protocol.getRestraint() is not None: - self._write_system(self._ref_system, ref_file=self._ref_file) + if self._ref_system is not None: + if ( + self._system.getAlchemicalIon() + or self._protocol.getRestraint() is not None + ): + self._write_system(self._ref_system, ref_file=self._ref_file) else: _shutil.copy(self._gro_file, self._ref_file) @@ -407,6 +408,8 @@ def _write_system(self, system, coord_file=None, topol_file=None, ref_file=None) self._restraint.toString( engine="GROMACS", perturbation_type=self._protocol.getPerturbationType(), + restraint_lambda="restraint" + in self._protocol.getLambda(type="series"), ) ) @@ -2078,7 +2081,7 @@ def _add_position_restraints(self, config_options): # Get the restraint type. restraint = self._protocol.getRestraint() - if restraint is not None: + if restraint is not None or self._system.getAlchemicalIon(): # Get the force constant in units of kJ_per_mol/nanometer**2 force_constant = self._protocol.getForceConstant()._sire_unit force_constant = force_constant.to( @@ -2140,8 +2143,15 @@ def _add_position_restraints(self, config_options): moltypes_sys_idx[mol_type].append(idx) sys_idx_moltypes[idx] = mol_type + if self._system.getAlchemicalIon(): + biggest_mol_idx = max( + range(system.nMolecules()), key=lambda x: system[x].nAtoms() + ) + else: + biggest_mol_idx = -1 + # A keyword restraint. - if isinstance(restraint, str): + if isinstance(restraint, str) or self._system.getAlchemicalIon(): # The number of restraint files. num_restraint = 1 @@ -2158,12 +2168,36 @@ def _add_position_restraints(self, config_options): for idx, mol_idx in enumerate(mol_idxs): # Get the indices of any restrained atoms in this molecule, # making sure that indices are relative to the molecule. - atom_idxs = self._system.getRestraintAtoms( - restraint, - mol_index=mol_idx, - is_absolute=False, - allow_zero_matches=True, - ) + if restraint is not None: + atom_idxs = self._system.getRestraintAtoms( + restraint, + mol_index=mol_idx, + is_absolute=False, + allow_zero_matches=True, + ) + else: + atom_idxs = [] + + if self._system.getMolecule(mol_idx).isAlchemicalIon(): + alch_ion = self._system.getMolecule(mol_idx).getAtoms() + alch_idx = alch_ion[0].index() + if alch_idx != 0 or len(alch_ion) != 1: + # The alchemical ions should only contain 1 atom + # and the relative index should thus be 0. + raise ValueError( + f"{self._system.getMolecule(mol_idx)} is marked as an alchemical ion but has more than 1 atom." + ) + else: + atom_idxs.append(alch_idx) + + if mol_idx == biggest_mol_idx: + # Only triggered when there is alchemical ion present. + # The biggest_mol_idx is -1 when there is no alchemical ion. + protein_com_idx = self._system.getMolecule( + mol_idx + ).getCOMIdx() + if protein_com_idx not in atom_idxs: + atom_idxs.append(protein_com_idx) # Store the atom index if it hasn't already been recorded. for atom_idx in atom_idxs: diff --git a/python/BioSimSpace/Sandpit/Exscientia/Protocol/_config.py b/python/BioSimSpace/Sandpit/Exscientia/Protocol/_config.py index 5b9f53f6f..99446e14b 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/Protocol/_config.py +++ b/python/BioSimSpace/Sandpit/Exscientia/Protocol/_config.py @@ -1,16 +1,15 @@ +import math as _math import warnings as _warnings -import math as _math from sire.legacy import Units as _SireUnits -from ..Units.Time import nanosecond as _nanosecond -from .. import Protocol as _Protocol -from .. import _gmx_version -from .._Exceptions import IncompatibleError as _IncompatibleError +from .. import Protocol as _Protocol, _gmx_version +from ..Align._alch_ion import _get_protein_com_idx from ..Align._squash import _amber_mask_from_indices, _squashed_atom_mapping from ..FreeEnergy._restraint import Restraint as _Restraint from ..Units.Energy import kj_per_mol as _kj_per_mol from ..Units.Length import nanometer as _nanometer +from .._Exceptions import IncompatibleError as _IncompatibleError class ConfigFactory: @@ -248,6 +247,13 @@ def generateAmberConfig(self, extra_options=None, extra_lines=None): # Restrain the backbone. restraint = self.protocol.getRestraint() + 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}" + else: + alchemical_ion_mask = None + if restraint is not None: # Get the indices of the atoms that are restrained. if type(restraint) is str: @@ -304,13 +310,22 @@ def generateAmberConfig(self, extra_options=None, extra_lines=None): "AMBER atom 'restraintmask' exceeds 256 character limit!" ) - protocol_dict["ntr"] = 1 - force_constant = self.protocol.getForceConstant()._sire_unit - force_constant = force_constant.to( - _SireUnits.kcal_per_mol / _SireUnits.angstrom2 - ) - protocol_dict["restraint_wt"] = force_constant - protocol_dict["restraintmask"] = f'"{restraint_mask}"' + else: + restraint_mask = None + + if restraint_mask or alchemical_ion_mask: + if restraint_mask and alchemical_ion_mask: + restraint_mask = f"{restraint_mask} | {alchemical_ion_mask}" + elif alchemical_ion_mask: + restraint_mask = alchemical_ion_mask + + protocol_dict["ntr"] = 1 + force_constant = self.protocol.getForceConstant()._sire_unit + force_constant = force_constant.to( + _SireUnits.kcal_per_mol / _SireUnits.angstrom2 + ) + protocol_dict["restraint_wt"] = force_constant + protocol_dict["restraintmask"] = f'"{restraint_mask}"' # Pressure control. if not isinstance(self.protocol, _Protocol.Minimisation): @@ -516,8 +531,11 @@ def generateGromacsConfig( # Temperature control. if not isinstance(self.protocol, _Protocol.Minimisation): - protocol_dict["integrator"] = "md" # leap-frog dynamics. - protocol_dict["tcoupl"] = "v-rescale" + if isinstance(self.protocol, _Protocol._FreeEnergyMixin): + protocol_dict["integrator"] = "sd" # langevin dynamics. + else: + protocol_dict["integrator"] = "md" # leap-frog dynamics. + protocol_dict["tcoupl"] = "v-rescale" protocol_dict["tc-grps"] = ( "system" # A single temperature group for the entire system. ) @@ -565,34 +583,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." diff --git a/python/BioSimSpace/Sandpit/Exscientia/Trajectory/_trajectory.py b/python/BioSimSpace/Sandpit/Exscientia/Trajectory/_trajectory.py index 7b5f9d45e..1b7ac59ed 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/Trajectory/_trajectory.py +++ b/python/BioSimSpace/Sandpit/Exscientia/Trajectory/_trajectory.py @@ -157,7 +157,7 @@ def getFrame(trajectory, topology, index, system=None, property_map={}): errors = [] is_sire = False is_mdanalysis = False - pdb_file = work_dir + f"/{str(_uuid.uuid4())}.pdb" + pdb_file = _os.path.join(str(work_dir), f"{str(_uuid.uuid4())}.pdb") try: frame = _sire_load( [trajectory, topology], @@ -169,7 +169,7 @@ def getFrame(trajectory, topology, index, system=None, property_map={}): except Exception as e: errors.append(f"Sire: {str(e)}") try: - frame_file = work_dir + f"/{str(_uuid.uuid4())}.rst7" + frame_file = _os.path.join(str(work_dir), f"{str(_uuid.uuid4())}.rst7") frame = _mdtraj.load_frame(trajectory, index, top=topology) frame.save(frame_file, force_overwrite=True) frame.save(pdb_file, force_overwrite=True) @@ -178,7 +178,7 @@ def getFrame(trajectory, topology, index, system=None, property_map={}): errors.append(f"MDTraj: {str(e)}") # Try to load the frame with MDAnalysis. try: - frame_file = work_dir + f"/{str(_uuid.uuid4())}.gro" + frame_file = _os.path.join(str(work_dir), f"{str(_uuid.uuid4())}.gro") universe = _mdanalysis.Universe(topology, trajectory) universe.trajectory.trajectory[index] with _warnings.catch_warnings(): @@ -615,7 +615,9 @@ def getTrajectory(self, format="auto"): # If this is a PRM7 file, copy to PARM7. if extension == ".prm7": # Set the path to the temporary topology file. - top_file = self._work_dir + f"/{str(_uuid.uuid4())}.parm7" + top_file = _os.path.join( + str(self._work_dir), f"{str(_uuid.uuid4())}.parm7" + ) # Copy the topology to a file with the correct extension. _shutil.copyfile(self._top_file, top_file) @@ -761,16 +763,20 @@ def getFrames(self, indices=None): # Write the current frame to file. - pdb_file = self._work_dir + f"/{str(_uuid.uuid4())}.pdb" + pdb_file = _os.path.join(str(self._work_dir), f"{str(_uuid.uuid4())}.pdb") if self._backend == "SIRE": frame = self._trajectory[x] elif self._backend == "MDTRAJ": - frame_file = self._work_dir + f"/{str(_uuid.uuid4())}.rst7" + frame_file = _os.path.join( + str(self._work_dir), f"{str(_uuid.uuid4())}.rst7" + ) self._trajectory[x].save(frame_file, force_overwrite=True) self._trajectory[x].save(pdb_file, force_overwrite=True) elif self._backend == "MDANALYSIS": - frame_file = self._work_dir + f"/{str(_uuid.uuid4())}.gro" + frame_file = _os.path.join( + str(self._work_dir), f"{str(_uuid.uuid4())}.gro" + ) self._trajectory.trajectory[x] with _warnings.catch_warnings(): _warnings.simplefilter("ignore") @@ -1110,8 +1116,8 @@ def _split_molecules(frame, pdb, reference, work_dir, property_map={}): formats = reference.fileFormat() # Write the frame coordinates/velocities to file. - coord_file = work_dir + f"/{str(_uuid.uuid4())}.coords" - top_file = work_dir + f"/{str(_uuid.uuid4())}.top" + coord_file = _os.path.join(str(work_dir), f"{str(_uuid.uuid4())}.coords") + top_file = _os.path.join(str(work_dir), f"{str(_uuid.uuid4())}.top") frame.writeToFile(coord_file) # Whether we've parsed as a PDB file. diff --git a/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_molecule.py b/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_molecule.py index 8a381622b..6c9012c23 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_molecule.py +++ b/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_molecule.py @@ -32,20 +32,22 @@ from math import isclose as _isclose from warnings import warn as _warn -from sire.legacy import Base as _SireBase -from sire.legacy import IO as _SireIO -from sire.legacy import MM as _SireMM -from sire.legacy import Maths as _SireMaths -from sire.legacy import Mol as _SireMol -from sire.legacy import System as _SireSystem -from sire.legacy import Units as _SireUnits +import numpy as _np +from sire.legacy import ( + Base as _SireBase, + IO as _SireIO, + MM as _SireMM, + Maths as _SireMaths, + Mol as _SireMol, + System as _SireSystem, + Units as _SireUnits, +) +from ._sire_wrapper import SireWrapper as _SireWrapper from .. import _isVerbose +from ..Types import Coordinate as _Coordinate, Length as _Length, Vector as _BSSVector +from ..Units.Length import angstrom as _angstrom from .._Exceptions import IncompatibleError as _IncompatibleError -from ..Types import Coordinate as _Coordinate -from ..Types import Length as _Length - -from ._sire_wrapper import SireWrapper as _SireWrapper class Molecule(_SireWrapper): @@ -1887,6 +1889,20 @@ def _getPerturbationIndices(self): return idxs + def getCOMIdx(self): + """Get the index of the atom that closest to the center of mass.""" + if self.isPerturbable(): + property_map = {"coordinates": "coordinates0", "mass": "mass0"} + else: + property_map = {"coordinates": "coordinates", "mass": "mass"} + coords = self.coordinates(property_map=property_map) + com = self._getCenterOfMass(property_map=property_map) + com = _BSSVector(*[e / _angstrom for e in com]) + diffs = [coord.toVector() - com for coord in coords] + sq_distances = [diff.dot(diff) for diff in diffs] + idx = int(_np.argmin(sq_distances)) + return idx + # Import at bottom of module to avoid circular dependency. from ._atom import Atom as _Atom diff --git a/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py b/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py index 7dc882868..64426775e 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py +++ b/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py @@ -88,12 +88,20 @@ def __init__(self, system): sire_object = _SireSystem.System("BioSimSpace_System.") super().__init__(sire_object) self.addMolecules(_Molecule(system)) + if "fileformat" in system.propertyKeys(): + self._sire_object.setProperty( + "fileformat", system.property("fileformat") + ) # A BioSimSpace Molecule object. elif isinstance(system, _Molecule): sire_object = _SireSystem.System("BioSimSpace_System.") super().__init__(sire_object) self.addMolecules(system) + if "fileformat" in system._sire_object.propertyKeys(): + self._sire_object.setProperty( + "fileformat", system._sire_object.property("fileformat") + ) # A BioSimSpace Molecules object. elif isinstance(system, _Molecules): diff --git a/python/BioSimSpace/Trajectory/_trajectory.py b/python/BioSimSpace/Trajectory/_trajectory.py index 7b5f9d45e..1b7ac59ed 100644 --- a/python/BioSimSpace/Trajectory/_trajectory.py +++ b/python/BioSimSpace/Trajectory/_trajectory.py @@ -157,7 +157,7 @@ def getFrame(trajectory, topology, index, system=None, property_map={}): errors = [] is_sire = False is_mdanalysis = False - pdb_file = work_dir + f"/{str(_uuid.uuid4())}.pdb" + pdb_file = _os.path.join(str(work_dir), f"{str(_uuid.uuid4())}.pdb") try: frame = _sire_load( [trajectory, topology], @@ -169,7 +169,7 @@ def getFrame(trajectory, topology, index, system=None, property_map={}): except Exception as e: errors.append(f"Sire: {str(e)}") try: - frame_file = work_dir + f"/{str(_uuid.uuid4())}.rst7" + frame_file = _os.path.join(str(work_dir), f"{str(_uuid.uuid4())}.rst7") frame = _mdtraj.load_frame(trajectory, index, top=topology) frame.save(frame_file, force_overwrite=True) frame.save(pdb_file, force_overwrite=True) @@ -178,7 +178,7 @@ def getFrame(trajectory, topology, index, system=None, property_map={}): errors.append(f"MDTraj: {str(e)}") # Try to load the frame with MDAnalysis. try: - frame_file = work_dir + f"/{str(_uuid.uuid4())}.gro" + frame_file = _os.path.join(str(work_dir), f"{str(_uuid.uuid4())}.gro") universe = _mdanalysis.Universe(topology, trajectory) universe.trajectory.trajectory[index] with _warnings.catch_warnings(): @@ -615,7 +615,9 @@ def getTrajectory(self, format="auto"): # If this is a PRM7 file, copy to PARM7. if extension == ".prm7": # Set the path to the temporary topology file. - top_file = self._work_dir + f"/{str(_uuid.uuid4())}.parm7" + top_file = _os.path.join( + str(self._work_dir), f"{str(_uuid.uuid4())}.parm7" + ) # Copy the topology to a file with the correct extension. _shutil.copyfile(self._top_file, top_file) @@ -761,16 +763,20 @@ def getFrames(self, indices=None): # Write the current frame to file. - pdb_file = self._work_dir + f"/{str(_uuid.uuid4())}.pdb" + pdb_file = _os.path.join(str(self._work_dir), f"{str(_uuid.uuid4())}.pdb") if self._backend == "SIRE": frame = self._trajectory[x] elif self._backend == "MDTRAJ": - frame_file = self._work_dir + f"/{str(_uuid.uuid4())}.rst7" + frame_file = _os.path.join( + str(self._work_dir), f"{str(_uuid.uuid4())}.rst7" + ) self._trajectory[x].save(frame_file, force_overwrite=True) self._trajectory[x].save(pdb_file, force_overwrite=True) elif self._backend == "MDANALYSIS": - frame_file = self._work_dir + f"/{str(_uuid.uuid4())}.gro" + frame_file = _os.path.join( + str(self._work_dir), f"{str(_uuid.uuid4())}.gro" + ) self._trajectory.trajectory[x] with _warnings.catch_warnings(): _warnings.simplefilter("ignore") @@ -1110,8 +1116,8 @@ def _split_molecules(frame, pdb, reference, work_dir, property_map={}): formats = reference.fileFormat() # Write the frame coordinates/velocities to file. - coord_file = work_dir + f"/{str(_uuid.uuid4())}.coords" - top_file = work_dir + f"/{str(_uuid.uuid4())}.top" + coord_file = _os.path.join(str(work_dir), f"{str(_uuid.uuid4())}.coords") + top_file = _os.path.join(str(work_dir), f"{str(_uuid.uuid4())}.top") frame.writeToFile(coord_file) # Whether we've parsed as a PDB file. diff --git a/python/BioSimSpace/_Config/_amber.py b/python/BioSimSpace/_Config/_amber.py index 861fa5d72..34d6ad84d 100644 --- a/python/BioSimSpace/_Config/_amber.py +++ b/python/BioSimSpace/_Config/_amber.py @@ -321,10 +321,7 @@ def createConfig( ) # Temperature control. - if not isinstance( - self._protocol, - (_Protocol.Metadynamics, _Protocol.Steering, _Protocol.Minimisation), - ): + if not isinstance(self._protocol, _Protocol.Minimisation): # Langevin dynamics. protocol_dict["ntt"] = 3 # Collision frequency (1 / ps). diff --git a/python/BioSimSpace/_SireWrappers/_system.py b/python/BioSimSpace/_SireWrappers/_system.py index 2d458890b..dbb8db8fc 100644 --- a/python/BioSimSpace/_SireWrappers/_system.py +++ b/python/BioSimSpace/_SireWrappers/_system.py @@ -88,12 +88,20 @@ def __init__(self, system): sire_object = _SireSystem.System("BioSimSpace_System.") super().__init__(sire_object) self.addMolecules(_Molecule(system)) + if "fileformat" in system.propertyKeys(): + self._sire_object.setProperty( + "fileformat", system.property("fileformat") + ) # A BioSimSpace Molecule object. elif isinstance(system, _Molecule): sire_object = _SireSystem.System("BioSimSpace_System.") super().__init__(sire_object) self.addMolecules(system) + if "fileformat" in system._sire_object.propertyKeys(): + self._sire_object.setProperty( + "fileformat", system._sire_object.property("fileformat") + ) # A BioSimSpace Molecules object. elif isinstance(system, _Molecules): diff --git a/requirements.txt b/requirements.txt index 81dcc0ce2..a5caf0018 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,7 +1,7 @@ # BioSimSpace runtime requirements. # main -#sire~=2023.5.0 +#sire~=2023.5.2 # devel sire==2024.1.0.dev diff --git a/tests/Convert/test_convert.py b/tests/Convert/test_convert.py index f2f4efbb8..8823e4ef9 100644 --- a/tests/Convert/test_convert.py +++ b/tests/Convert/test_convert.py @@ -5,11 +5,6 @@ import pytest -@pytest.fixture(scope="session") -def system(): - return BSS.IO.readMolecules(["tests/input/ala.crd", "tests/input/ala.top"]) - - def test_system(system): """ Check that system conversions work as expected. diff --git a/tests/FreeEnergy/test_relative.py b/tests/FreeEnergy/test_relative.py index 6d01fe81f..6808b0ed7 100644 --- a/tests/FreeEnergy/test_relative.py +++ b/tests/FreeEnergy/test_relative.py @@ -9,17 +9,6 @@ from tests.conftest import url, has_alchemlyb, has_gromacs -@pytest.fixture(scope="module") -def perturbable_system(): - """Re-use the same perturbable system for each test.""" - return BSS.IO.readPerturbableSystem( - f"{url}/perturbable_system0.prm7", - f"{url}/perturbable_system0.rst7", - f"{url}/perturbable_system1.prm7", - f"{url}/perturbable_system1.rst7", - ) - - @pytest.fixture(scope="module") def fep_output(): """Path to a temporary directory containing FEP output.""" diff --git a/tests/Metadynamics/test_metadynamics.py b/tests/Metadynamics/test_metadynamics.py new file mode 100644 index 000000000..5266e1c5d --- /dev/null +++ b/tests/Metadynamics/test_metadynamics.py @@ -0,0 +1,161 @@ +import pytest +import socket + +import BioSimSpace as BSS + + +@pytest.mark.skipif( + socket.gethostname() != "porridge", + reason="Local test requiring PLUMED patched GROMACS.", +) +def test_metadynamics(system): + # Search for the first molecule containing ALA. + molecule = system.search("resname ALA").molecules()[0] + + # Store the torsion indices. + phi_idx = [4, 6, 8, 14] + psi_idx = [6, 8, 14, 16] + + # Create the collective variables. + phi = BSS.Metadynamics.CollectiveVariable.Torsion(atoms=phi_idx) + psi = BSS.Metadynamics.CollectiveVariable.Torsion(atoms=psi_idx) + + # Create the metadynamics protocol. + protocol = BSS.Protocol.Metadynamics( + collective_variable=[phi, psi], runtime=100 * BSS.Units.Time.picosecond + ) + + # Run the metadynamics simulation. + process = BSS.Metadynamics.run(molecule.toSystem(), protocol, gpu_support=True) + + # Wait for the process to finish. + process.wait() + + # Check if the process has finished successfully. + assert not process.isError() + + free_nrg = process.getFreeEnergy(kt=BSS.Units.Energy.kt) + + # Check if the free energy is not None. + assert free_nrg is not None + + +@pytest.mark.skipif( + socket.gethostname() != "porridge", + reason="Local test requiring PLUMED patched GROMACS.", +) +def test_steering(system): + # Find the indices of the atoms in the ALA residue. + rmsd_idx = [x.index() for x in system.search("resname ALA").atoms()] + + # Create the collective variable. + cv = BSS.Metadynamics.CollectiveVariable.RMSD(system, system[0], rmsd_idx) + + # Add some stages. + start = 0 * BSS.Units.Time.nanosecond + apply_force = 4 * BSS.Units.Time.picosecond + steer = 50 * BSS.Units.Time.picosecond + relax = 100 * BSS.Units.Time.picosecond + + # Create some restraints. + nm = BSS.Units.Length.nanometer + restraint_1 = BSS.Metadynamics.Restraint(cv.getInitialValue(), 0) + restraint_2 = BSS.Metadynamics.Restraint(cv.getInitialValue(), 3500) + restraint_3 = BSS.Metadynamics.Restraint(0 * nm, 3500) + restraint_4 = BSS.Metadynamics.Restraint(0 * nm, 0) + + # Create the steering protocol. + protocol = BSS.Protocol.Steering( + cv, + [start, apply_force, steer, relax], + [restraint_1, restraint_2, restraint_3, restraint_4], + runtime=100 * BSS.Units.Time.picosecond, + ) + + # Create the steering process. + process = BSS.Process.Gromacs(system, protocol, extra_args={"-ntmpi": 1}) + + # Start the process and wait for it to finish. + process.start() + process.wait() + + # Check if the process has finished successfully. + assert not process.isError() + + +@pytest.mark.skipif( + socket.gethostname() != "porridge", + reason="Local test requiring PLUMED.", +) +def test_funnel_metadynamics(): + # Load the protein-ligand system. + system = BSS.IO.readMolecules( + BSS.IO.expand(BSS.tutorialUrl(), ["funnel_system.rst7", "funnel_system.prm7"]) + ) + + # Get the p0 and p1 points for defining the funnel. + p0, p1 = BSS.Metadynamics.CollectiveVariable.makeFunnel(system) + + # Expected p0 and p1 points. + expected_p0 = [1017, 1031, 1050, 1186, 1205, 1219, 1238, 2585, 2607, 2623] + expected_p1 = [ + 519, + 534, + 553, + 572, + 583, + 597, + 608, + 619, + 631, + 641, + 1238, + 1254, + 1280, + 1287, + 1306, + 1313, + 1454, + 1473, + 1480, + 1863, + 1879, + 1886, + 1906, + 2081, + 2116, + 2564, + 2571, + 2585, + 2607, + ] + + # Make sure the p0 and p1 points are as expected. + assert p0 == expected_p0 + assert p1 == expected_p1 + + # Set the upper bound for the funnel collective variable. + upper_bound = BSS.Metadynamics.Bound(value=3.5 * BSS.Units.Length.nanometer) + + # Create the funnel collective variable. + cv = BSS.Metadynamics.CollectiveVariable.Funnel(p0, p1, upper_bound=upper_bound) + + # Create the metadynamics protocol. + protocol = BSS.Protocol.Metadynamics( + cv, + runtime=1 * BSS.Units.Time.picosecond, + hill_height=1.5 * BSS.Units.Energy.kj_per_mol, + hill_frequency=500, + restart_interval=1000, + bias_factor=10, + ) + + # Create the metadynamics process. + process = BSS.Process.OpenMM(system, protocol) + + # Start the process and wait for it to finish. + process.start() + process.wait() + + # Check if the process has finished successfully. + assert not process.isError() diff --git a/tests/Process/test_amber.py b/tests/Process/test_amber.py index a6626047e..18a23df2a 100644 --- a/tests/Process/test_amber.py +++ b/tests/Process/test_amber.py @@ -13,13 +13,7 @@ restraints = BSS.Protocol._position_restraint_mixin._PositionRestraintMixin.restraints() -@pytest.fixture(scope="session") -def system(): - """Re-use the same molecuar system for each test.""" - return BSS.IO.readMolecules(["tests/input/ala.top", "tests/input/ala.crd"]) - - -@pytest.fixture(scope="session") +@pytest.fixture(scope="module") def rna_system(): """An RNA system for re-use.""" return BSS.IO.readMolecules( @@ -27,7 +21,7 @@ def rna_system(): ) -@pytest.fixture(scope="session") +@pytest.fixture(scope="module") def large_protein_system(): """A large protein system for re-use.""" return BSS.IO.readMolecules( @@ -35,28 +29,6 @@ def large_protein_system(): ) -@pytest.fixture(scope="module") -def perturbable_system(): - """Re-use the same perturbable system for each test.""" - return BSS.IO.readPerturbableSystem( - f"{url}/perturbable_system0.prm7", - f"{url}/perturbable_system0.rst7", - f"{url}/perturbable_system1.prm7", - f"{url}/perturbable_system1.rst7", - ) - - -@pytest.fixture(scope="module") -def solvated_perturbable_system(): - """Re-use the same solvated perturbable system for each test.""" - return BSS.IO.readPerturbableSystem( - f"{url}/solvated_perturbable_system0.prm7", - f"{url}/solvated_perturbable_system0.rst7", - f"{url}/solvated_perturbable_system1.prm7", - f"{url}/solvated_perturbable_system1.rst7", - ) - - @pytest.mark.skipif(has_amber is False, reason="Requires AMBER to be installed.") @pytest.mark.parametrize("restraint", restraints) def test_minimise(system, restraint): diff --git a/tests/Process/test_gromacs.py b/tests/Process/test_gromacs.py index a321f033d..5de10a95a 100644 --- a/tests/Process/test_gromacs.py +++ b/tests/Process/test_gromacs.py @@ -19,23 +19,6 @@ restraints = BSS.Protocol._position_restraint_mixin._PositionRestraintMixin.restraints() -@pytest.fixture(scope="session") -def system(): - """Re-use the same molecuar system for each test.""" - return BSS.IO.readMolecules(["tests/input/ala.top", "tests/input/ala.crd"]) - - -@pytest.fixture(scope="session") -def perturbable_system(): - """Re-use the same perturbable system for each test.""" - return 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", - ) - - @pytest.mark.skipif(has_gromacs is False, reason="Requires GROMACS to be installed.") @pytest.mark.parametrize("restraint", restraints) def test_minimise(system, restraint): diff --git a/tests/Process/test_namd.py b/tests/Process/test_namd.py index 03c550d46..7c9eb19b7 100644 --- a/tests/Process/test_namd.py +++ b/tests/Process/test_namd.py @@ -8,8 +8,8 @@ restraints = BSS.Protocol._position_restraint_mixin._PositionRestraintMixin.restraints() -@pytest.fixture(scope="session") -def system(): +@pytest.fixture(scope="module") +def namd_system(): """Re-use the same molecuar system for each test.""" return BSS.IO.readMolecules( [ @@ -22,19 +22,19 @@ def system(): @pytest.mark.skipif(has_namd is False, reason="Requires NAMD to be installed.") @pytest.mark.parametrize("restraint", restraints) -def test_minimise(system, restraint): +def test_minimise(namd_system, restraint): """Test a minimisation protocol.""" # Create a short minimisation protocol. protocol = BSS.Protocol.Minimisation(steps=100, restraint=restraint) # Run the process, check that it finished without error, and returns a system. - run_process(system, protocol) + run_process(namd_system, protocol) @pytest.mark.skipif(has_namd is False, reason="Requires NAMD to be installed.") @pytest.mark.parametrize("restraint", restraints) -def test_equilibrate(system, restraint): +def test_equilibrate(namd_system, restraint): """Test an equilibration protocol.""" # Create a short equilibration protocol. @@ -43,11 +43,11 @@ def test_equilibrate(system, restraint): ) # Run the process, check that it finished without error, and returns a system. - run_process(system, protocol) + run_process(namd_system, protocol) @pytest.mark.skipif(has_namd is False, reason="Requires NAMD to be installed.") -def test_heat(system): +def test_heat(namd_system): """Test a heating protocol.""" # Create a short heating protocol. @@ -58,11 +58,11 @@ def test_heat(system): ) # Run the process, check that it finished without error, and returns a system. - run_process(system, protocol) + run_process(namd_system, protocol) @pytest.mark.skipif(has_namd is False, reason="Requires NAMD to be installed.") -def test_cool(system): +def test_cool(namd_system): """Test a cooling protocol.""" # Create a short heating protocol. @@ -73,12 +73,12 @@ def test_cool(system): ) # Run the process, check that it finished without error, and returns a system. - run_process(system, protocol) + run_process(namd_system, protocol) @pytest.mark.skipif(has_namd is False, reason="Requires NAMD to be installed.") @pytest.mark.parametrize("restraint", restraints) -def test_production(system, restraint): +def test_production(namd_system, restraint): """Test a production protocol.""" # Create a short production protocol. @@ -87,14 +87,14 @@ def test_production(system, restraint): ) # Run the process, check that it finished without error, and returns a system. - run_process(system, protocol) + run_process(namd_system, protocol) -def run_process(system, protocol): +def run_process(namd_system, protocol): """Helper function to run various simulation protocols.""" # Initialise the NAMD process. - process = BSS.Process.Namd(system, protocol, name="test") + process = BSS.Process.Namd(namd_system, protocol, name="test") # Start the NAMD simulation. process.start() diff --git a/tests/Process/test_openmm.py b/tests/Process/test_openmm.py index 7925bd91b..a7968173a 100644 --- a/tests/Process/test_openmm.py +++ b/tests/Process/test_openmm.py @@ -8,12 +8,6 @@ restraints = BSS.Protocol._position_restraint_mixin._PositionRestraintMixin.restraints() -@pytest.fixture(scope="session") -def system(): - """Re-use the same molecuar system for each test.""" - return BSS.IO.readMolecules(["tests/input/ala.top", "tests/input/ala.crd"]) - - @pytest.mark.parametrize("restraint", restraints) def test_minimise(system, restraint): """Test a minimisation protocol.""" diff --git a/tests/Process/test_single_point_energy.py b/tests/Process/test_single_point_energy.py index e252cfb99..8fc3a604e 100644 --- a/tests/Process/test_single_point_energy.py +++ b/tests/Process/test_single_point_energy.py @@ -5,8 +5,8 @@ from tests.conftest import url, has_amber, has_gromacs -@pytest.fixture(scope="session") -def system(): +@pytest.fixture(scope="module") +def ubiquitin_system(): """Re-use the same molecuar system for each test.""" return BSS.IO.readMolecules( [f"{url}/ubiquitin.prm7.bz2", f"{url}/ubiquitin.rst7.bz2"] @@ -17,17 +17,19 @@ def system(): has_amber is False or has_gromacs is False, reason="Requires that both AMBER and GROMACS are installed.", ) -def test_amber_gromacs(system): +def test_amber_gromacs(ubiquitin_system): """Single point energy comparison between AMBER and GROMACS.""" # Create a single-step minimisation protocol. protocol = BSS.Protocol.Minimisation(steps=1) # Create a process to run with AMBER. - process_amb = BSS.Process.Amber(system, protocol) + process_amb = BSS.Process.Amber(ubiquitin_system, protocol) # Create a process to run with GROMACS. - process_gmx = BSS.Process.Gromacs(system, protocol, extra_options={"nsteps": 0}) + process_gmx = BSS.Process.Gromacs( + ubiquitin_system, protocol, extra_options={"nsteps": 0} + ) # Run the AMBER process and wait for it to finish. process_amb.start() @@ -57,23 +59,25 @@ def test_amber_gromacs(system): has_amber is False or has_gromacs is False, reason="Requires that both AMBER and GROMACS are installed.", ) -def test_amber_gromacs_triclinic(system): +def test_amber_gromacs_triclinic(ubiquitin_system): """Single point energy comparison between AMBER and GROMACS in a triclinic box.""" # Swap the space for a triclinic cell (truncated octahedron). from sire.legacy.Vol import TriclinicBox triclinic_box = TriclinicBox.truncatedOctahedron(50) - system._sire_object.setProperty("space", triclinic_box) + ubiquitin_system._sire_object.setProperty("space", triclinic_box) # Create a single-step minimisation protocol. protocol = BSS.Protocol.Minimisation(steps=1) # Create a process to run with AMBER. - process_amb = BSS.Process.Amber(system, protocol) + process_amb = BSS.Process.Amber(ubiquitin_system, protocol) # Create a process to run with GROMACS. - process_gmx = BSS.Process.Gromacs(system, protocol, extra_options={"nsteps": 0}) + process_gmx = BSS.Process.Gromacs( + ubiquitin_system, protocol, extra_options={"nsteps": 0} + ) # Run the AMBER process and wait for it to finish. process_amb.start() diff --git a/tests/Process/test_somd.py b/tests/Process/test_somd.py index 2220e589e..f7eedc52c 100644 --- a/tests/Process/test_somd.py +++ b/tests/Process/test_somd.py @@ -9,23 +9,6 @@ url = BSS.tutorialUrl() -@pytest.fixture(scope="session") -def system(): - """Re-use the same molecuar system for each test.""" - return BSS.IO.readMolecules(["tests/input/ala.top", "tests/input/ala.crd"]) - - -@pytest.fixture(scope="session") -def perturbable_system(): - """Re-use the same perturbable system for each test.""" - return BSS.IO.readPerturbableSystem( - f"{url}/perturbable_system0.prm7", - f"{url}/perturbable_system0.rst7", - f"{url}/perturbable_system1.prm7", - f"{url}/perturbable_system1.rst7", - ) - - def test_minimise(system): """Test a minimisation protocol.""" diff --git a/tests/Protocol/test_protocol.py b/tests/Protocol/test_protocol.py index c38e8e553..b0614b6fa 100644 --- a/tests/Protocol/test_protocol.py +++ b/tests/Protocol/test_protocol.py @@ -6,12 +6,6 @@ # using strings of unit-based types. -@pytest.fixture(scope="session") -def system(): - """Re-use the same molecuar system for each test.""" - return BSS.IO.readMolecules(["tests/input/ala.top", "tests/input/ala.crd"]) - - def test_equilibration(): # Instantiate from types. p0 = BSS.Protocol.Equilibration( diff --git a/tests/Sandpit/Exscientia/Align/test_alchemical_ion.py b/tests/Sandpit/Exscientia/Align/test_alchemical_ion.py index fa2a0aeb6..a3320e5da 100644 --- a/tests/Sandpit/Exscientia/Align/test_alchemical_ion.py +++ b/tests/Sandpit/Exscientia/Align/test_alchemical_ion.py @@ -1,10 +1,12 @@ import pytest import BioSimSpace.Sandpit.Exscientia as BSS -from BioSimSpace.Sandpit.Exscientia.Align._alch_ion import _mark_alchemical_ion +from BioSimSpace.Sandpit.Exscientia.Align._alch_ion import ( + _get_protein_com_idx, + _mark_alchemical_ion, +) from BioSimSpace.Sandpit.Exscientia._SireWrappers import Molecule - -from tests.conftest import root_fp, has_gromacs +from tests.conftest import has_gromacs, root_fp @pytest.fixture @@ -50,3 +52,9 @@ def test_getAlchemicalIon(input_system, isalchem, request): def test_getAlchemicalIonIdx(alchemical_ion_system): index = alchemical_ion_system.getAlchemicalIonIdx() assert index == 680 + + +@pytest.mark.skipif(has_gromacs is False, reason="Requires GROMACS to be installed.") +def test_get_protein_com_idx(alchemical_ion_system): + index = _get_protein_com_idx(alchemical_ion_system) + assert index == 8 diff --git a/tests/Sandpit/Exscientia/FreeEnergy/test_restraint.py b/tests/Sandpit/Exscientia/FreeEnergy/test_restraint.py index e5d7ce24d..0149a64e3 100644 --- a/tests/Sandpit/Exscientia/FreeEnergy/test_restraint.py +++ b/tests/Sandpit/Exscientia/FreeEnergy/test_restraint.py @@ -97,11 +97,22 @@ def test_numerical_correction_boresch(boresch_restraint): def test_analytical_correction_boresch(boresch_restraint): - dG = boresch_restraint.getCorrection(method="analytical") / kcal_per_mol + dG = ( + boresch_restraint.getCorrection(method="analytical", flavour="boresch") + / kcal_per_mol + ) assert np.isclose(-7.2, dG, atol=0.1) assert isinstance(boresch_restraint, Restraint) +def test_analytical_schrodinger_correction_boresch(boresch_restraint): + dG = ( + boresch_restraint.getCorrection(method="analytical", flavour="schrodinger") + / kcal_per_mol + ) + assert np.isclose(-7.2, dG, atol=0.1) + + test_force_constants_boresch = [ ({"kr": 0}, ValueError), ({"kthetaA": 0}, ValueError), @@ -171,7 +182,9 @@ def test_dihedral(self, Topology): assert aj == "2" assert ak == "1" assert al == "1496" + assert type == "2" assert phiA == "148.396" + assert kA == "0.00" assert phiB == "148.396" assert kB == "41.84" ai, aj, ak, al, type, phiA, kA, phiB, kB = Topology[11].split() @@ -186,6 +199,30 @@ def test_dihedral(self, Topology): assert al == "1498" +class TestGromacsOutputBoreschRestraintLambda(TestGromacsOutputBoresch): + @staticmethod + @pytest.fixture(scope="class") + def Topology(boresch_restraint): + return boresch_restraint.toString( + engine="Gromacs", restraint_lambda=True + ).split("\n") + + def test_dihedral(self, Topology): + assert "dihedral_restraints" in Topology[8] + ai, aj, ak, al, type, phiA, dphiA, kA, phiB, dphiB, kB = Topology[10].split() + assert ai == "3" + assert aj == "2" + assert ak == "1" + assert al == "1496" + assert type == "1" + assert phiA == "148.396" + assert dphiA == "0.00" + assert dphiB == "0.00" + assert kA == "0.00" + assert phiB == "148.396" + assert kB == "41.84" + + class TestSomdOutputBoresch: @staticmethod @pytest.fixture(scope="class") diff --git a/tests/Sandpit/Exscientia/Process/test_amber.py b/tests/Sandpit/Exscientia/Process/test_amber.py index 73866d01a..1a6daded5 100644 --- a/tests/Sandpit/Exscientia/Process/test_amber.py +++ b/tests/Sandpit/Exscientia/Process/test_amber.py @@ -5,11 +5,9 @@ import numpy as np import pandas as pd import pytest -import shutil import BioSimSpace.Sandpit.Exscientia as BSS - -from tests.Sandpit.Exscientia.conftest import url, has_amber, has_pyarrow +from tests.Sandpit.Exscientia.conftest import has_amber, has_pyarrow from tests.conftest import root_fp diff --git a/tests/Sandpit/Exscientia/Process/test_gromacs.py b/tests/Sandpit/Exscientia/Process/test_gromacs.py index ec3b9b7f7..9a1886e97 100644 --- a/tests/Sandpit/Exscientia/Process/test_gromacs.py +++ b/tests/Sandpit/Exscientia/Process/test_gromacs.py @@ -1,13 +1,12 @@ import math -import numpy as np -import pytest -import pandas as pd import shutil - from pathlib import Path -import BioSimSpace.Sandpit.Exscientia as BSS +import numpy as np +import pandas as pd +import pytest +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.Angle import radian @@ -17,15 +16,13 @@ from BioSimSpace.Sandpit.Exscientia.Units.Temperature import kelvin from BioSimSpace.Sandpit.Exscientia.Units.Time import picosecond from BioSimSpace.Sandpit.Exscientia.Units.Volume import nanometer3 - from tests.Sandpit.Exscientia.conftest import ( - url, - has_alchemlyb, has_alchemtest, has_amber, has_gromacs, has_openff, has_pyarrow, + url, ) from tests.conftest import root_fp @@ -167,57 +164,87 @@ def test_restraints(perturbable_system, restraint): has_gromacs is False or has_pyarrow is False, reason="Requires GROMACS and pyarrow to be installed.", ) -def test_write_restraint(system, tmp_path): - """Test if the restraint has been written in a way that could be processed - correctly. - """ - ligand = ligand = BSS.IO.readMolecules( - [f"{url}/ligand01.prm7.bz2", f"{url}/ligand01.rst7.bz2"] - ).getMolecule(0) - decoupled_ligand = decouple(ligand) - l1 = decoupled_ligand.getAtoms()[0] - l2 = decoupled_ligand.getAtoms()[1] - l3 = decoupled_ligand.getAtoms()[2] - ligand_2 = BSS.IO.readMolecules( - [f"{url}/ligand04.prm7.bz2", f"{url}/ligand04.rst7.bz2"] - ).getMolecule(0) - r1 = ligand_2.getAtoms()[0] - r2 = ligand_2.getAtoms()[1] - r3 = ligand_2.getAtoms()[2] - system = (decoupled_ligand + ligand_2).toSystem() - - restraint_dict = { - "anchor_points": {"r1": r1, "r2": r2, "r3": r3, "l1": l1, "l2": l2, "l3": l3}, - "equilibrium_values": { - "r0": 7.84 * angstrom, - "thetaA0": 0.81 * radian, - "thetaB0": 1.74 * radian, - "phiA0": 2.59 * radian, - "phiB0": -1.20 * radian, - "phiC0": 2.63 * radian, - }, - "force_constants": { - "kr": 10 * kcal_per_mol / angstrom**2, - "kthetaA": 10 * kcal_per_mol / (radian * radian), - "kthetaB": 10 * kcal_per_mol / (radian * radian), - "kphiA": 10 * kcal_per_mol / (radian * radian), - "kphiB": 10 * kcal_per_mol / (radian * radian), - "kphiC": 10 * kcal_per_mol / (radian * radian), - }, - } - restraint = Restraint( - system, restraint_dict, 300 * kelvin, restraint_type="Boresch" - ) +class TestRestraint: + @pytest.fixture(scope="class") + def setup(self): + ligand = BSS.IO.readMolecules( + [f"{url}/ligand01.prm7.bz2", f"{url}/ligand01.rst7.bz2"] + ).getMolecule(0) + decoupled_ligand = decouple(ligand) + l1 = decoupled_ligand.getAtoms()[0] + l2 = decoupled_ligand.getAtoms()[1] + l3 = decoupled_ligand.getAtoms()[2] + ligand_2 = BSS.IO.readMolecules( + [f"{url}/ligand04.prm7.bz2", f"{url}/ligand04.rst7.bz2"] + ).getMolecule(0) + r1 = ligand_2.getAtoms()[0] + r2 = ligand_2.getAtoms()[1] + r3 = ligand_2.getAtoms()[2] + system = (decoupled_ligand + ligand_2).toSystem() + restraint_dict = { + "anchor_points": { + "r1": r1, + "r2": r2, + "r3": r3, + "l1": l1, + "l2": l2, + "l3": l3, + }, + "equilibrium_values": { + "r0": 7.84 * angstrom, + "thetaA0": 0.81 * radian, + "thetaB0": 1.74 * radian, + "phiA0": 2.59 * radian, + "phiB0": -1.20 * radian, + "phiC0": 2.63 * radian, + }, + "force_constants": { + "kr": 10 * kcal_per_mol / angstrom**2, + "kthetaA": 10 * kcal_per_mol / (radian * radian), + "kthetaB": 10 * kcal_per_mol / (radian * radian), + "kphiA": 10 * kcal_per_mol / (radian * radian), + "kphiB": 10 * kcal_per_mol / (radian * radian), + "kphiC": 10 * kcal_per_mol / (radian * radian), + }, + } + restraint = Restraint( + system, restraint_dict, 300 * kelvin, restraint_type="Boresch" + ) + return system, restraint + + def test_regular_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" + ) - # Create a short production protocol. - protocol = BSS.Protocol.FreeEnergy( - runtime=BSS.Types.Time(0.0001, "nanoseconds"), perturbation_type="full" - ) + # 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() + + def test_restraint_lambda(self, setup, tmp_path_factory): + """Test if the restraint has been written correctly when restraint lambda is evoked.""" + 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"), + lam=pd.Series(data={"bonded": 0.0, "restraint": 0.0}), + lam_vals=pd.DataFrame(data={"bonded": [0.0, 1.0], "restraint": [0.0, 1.0]}), + ) - # 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() + # 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 "dihedral_restraints" in f.read() + with open(tmp_path / "test.mdp", "r") as f: + assert "restraint-lambdas" in f.read() def run_process(system, protocol, **kwargs): diff --git a/tests/Sandpit/Exscientia/Process/test_position_restraint.py b/tests/Sandpit/Exscientia/Process/test_position_restraint.py index 7355adc70..b82ec3c1f 100644 --- a/tests/Sandpit/Exscientia/Process/test_position_restraint.py +++ b/tests/Sandpit/Exscientia/Process/test_position_restraint.py @@ -1,17 +1,20 @@ -from difflib import unified_diff - import itertools import os +from difflib import unified_diff import pandas as pd import pytest +import sire as sr +from sire.legacy import Units as SireUnits from sire.legacy.IO import AmberRst import BioSimSpace.Sandpit.Exscientia as BSS +from BioSimSpace.Sandpit.Exscientia.Align._alch_ion import _mark_alchemical_ion from BioSimSpace.Sandpit.Exscientia.Units.Energy import kj_per_mol from BioSimSpace.Sandpit.Exscientia.Units.Length import angstrom - +from BioSimSpace.Sandpit.Exscientia._SireWrappers import Molecule from tests.Sandpit.Exscientia.conftest import has_amber, has_gromacs, has_openff +from tests.conftest import root_fp @pytest.fixture @@ -22,6 +25,41 @@ def system(): return BSS.Align.merge(mol0, mol1).toSystem() +@pytest.fixture(scope="session") +def alchemical_ion_system(): + """A large protein system for re-use.""" + system = BSS.IO.readMolecules( + [f"{root_fp}/input/ala.top", f"{root_fp}/input/ala.crd"] + ) + solvated = BSS.Solvent.tip3p( + system, box=[4 * BSS.Units.Length.nanometer] * 3, ion_conc=0.15 + ) + ion = solvated.getMolecule(-1) + pert_ion = BSS.Align.merge(ion, ion, mapping={0: 0}) + 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 + + +@pytest.fixture(scope="session") +def alchemical_ion_system_psores(alchemical_ion_system): + # Create a reference system with different coordinate + system = alchemical_ion_system.copy() + mol = system.getMolecule(0) + sire_mol = mol._sire_object + atoms = sire_mol.cursor().atoms() + atoms[0]["coordinates"] = sr.maths.Vector(0, 0, 0) + new_mol = atoms.commit() + system.updateMolecule(0, Molecule(new_mol)) + return system + + @pytest.fixture def ref_system(system): mol = system[0] @@ -146,3 +184,99 @@ def test_amber(protocol, system, ref_system, tmp_path): # We are pointing the reference to the correct file assert f"{proc._work_dir}/{proc.getArgs()['-ref']}" == proc._ref_file + + +@pytest.mark.skipif( + has_gromacs is False or has_openff is False, + reason="Requires GROMACS and openff to be installed", +) +@pytest.mark.parametrize( + "restraint", + ["backbone", "heavy", "all", "none"], +) +def test_gromacs_alchemical_ion( + alchemical_ion_system, restraint, alchemical_ion_system_psores +): + protocol = BSS.Protocol.FreeEnergy(restraint=restraint) + process = BSS.Process.Gromacs( + alchemical_ion_system, + protocol, + name="test", + reference_system=alchemical_ion_system_psores, + ignore_warnings=True, + ) + + # Test the position restraint for protein center + with open(f"{process.workDir()}/posre_0001.itp", "r") as f: + posres = f.read().split("\n") + posres = [tuple(line.split()) for line in posres] + + assert ("9", "1", "4184.0", "4184.0", "4184.0") in posres + + # Test the position restraint for alchemical ion + with open(f"{process.workDir()}/test.top", "r") as f: + top = f.read() + lines = top[top.index("Merged_Molecule") :].split("\n") + assert lines[6] == '#include "posre_0002.itp"' + + with open(f"{process.workDir()}/posre_0002.itp", "r") as f: + posres = f.read().split("\n") + + assert posres[2].split() == ["1", "1", "4184.0", "4184.0", "4184.0"] + + # Test if the original coordinate is correct + with open(f"{process.workDir()}/test.gro", "r") as f: + gro = f.read().splitlines() + assert gro[2].split() == ["1ACE", "HH31", "1", "1.791", "1.610", "2.058"] + + # Test if the reference coordinate is passed + with open(f"{process.workDir()}/test_ref.gro", "r") as f: + gro = f.read().splitlines() + assert gro[2].split() == ["1ACE", "HH31", "1", "0.000", "0.000", "0.000"] + + +@pytest.mark.skipif( + has_amber is False or has_gromacs is False or has_openff is False, + reason="Requires AMBER, GROMACS and OpenFF to be installed", +) +@pytest.mark.parametrize( + ("restraint", "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"), + ], +) +def test_amber_alchemical_ion( + alchemical_ion_system, restraint, target, alchemical_ion_system_psores +): + # Create an equilibration protocol with backbone restraints. + protocol = BSS.Protocol.Equilibration(restraint=restraint) + + # Create the process object. + process = BSS.Process.Amber( + alchemical_ion_system, + protocol, + name="test", + reference_system=alchemical_ion_system_psores, + ) + + # Check that the correct restraint mask is in the config. + config = " ".join(process.getConfig()) + assert target in config + # Check is the reference file is passed to the cmd + assert "-ref test_ref.rst7" in process.getArgString() + + # Test if the original coordinate is correct + original = BSS.IO.readMolecules( + [f"{process.workDir()}/test.prm7", f"{process.workDir()}/test.rst7"] + ) + original_crd = original.getMolecule(0).coordinates()[0] + assert str(original_crd) == "(17.9138 A, 16.0981 A, 20.5786 A)" + # Test if the reference coordinate is passed + ref = BSS.IO.readMolecules( + [f"{process.workDir()}/test.prm7", f"{process.workDir()}/test_ref.rst7"] + ) + ref_crd = ref.getMolecule(0).coordinates()[0] + assert str(ref_crd) == "(0.0000e+00 A, 0.0000e+00 A, 0.0000e+00 A)" diff --git a/tests/Sandpit/Exscientia/Protocol/test_config.py b/tests/Sandpit/Exscientia/Protocol/test_config.py index efb52898a..32446dcda 100644 --- a/tests/Sandpit/Exscientia/Protocol/test_config.py +++ b/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 @@ -110,6 +111,16 @@ def test_tau_t(self, system, protocol): expected_res = {"tau-t = 2.00000"} assert expected_res.issubset(res) + @pytest.mark.parametrize("protocol", [Production, FreeEnergy]) + def test_integrator(self, system, protocol): + config = ConfigFactory(system, protocol(tau_t=BSS.Types.Time(2, "picosecond"))) + res = config.generateGromacsConfig() + if isinstance(protocol(), BSS.Protocol._FreeEnergyMixin): + expected_res = {"integrator = sd"} + else: + expected_res = {"integrator = md", "tcoupl = v-rescale"} + assert expected_res.issubset(res) + @pytest.mark.skipif( has_gromacs is False, reason="Requires GROMACS to be installed." ) @@ -301,6 +312,50 @@ def test_decouple_vdw_q(self, system): assert "couple-lambda1 = none" in mdp_text assert "couple-intramol = yes" in mdp_text + @pytest.mark.skipif( + has_gromacs is False, reason="Requires GROMACS to be installed." + ) + 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" not in mdp_text + assert "couple-lambda0" not in mdp_text + assert "couple-lambda1" not in mdp_text + assert "couple-intramol" not in mdp_text + @pytest.mark.skipif( has_gromacs is False, reason="Requires GROMACS to be installed." ) diff --git a/tests/Solvent/test_solvent.py b/tests/Solvent/test_solvent.py index 185bb6373..02afe1f06 100644 --- a/tests/Solvent/test_solvent.py +++ b/tests/Solvent/test_solvent.py @@ -11,7 +11,7 @@ @pytest.fixture(scope="module") -def system(): +def kigaki_system(): return BSS.IO.readMolecules( BSS.IO.expand( BSS.tutorialUrl(), ["kigaki_xtal_water.gro", "kigaki_xtal_water.top"] @@ -25,7 +25,7 @@ def system(): [partial(BSS.Solvent.solvate, "tip3p"), BSS.Solvent.tip3p], ) @pytest.mark.skipif(not has_gromacs, reason="Requires GROMACS to be installed") -def test_crystal_water(system, match_water, function): +def test_crystal_water(kigaki_system, match_water, function): """ Test that user defined crystal waters can be preserved during solvation and on write to GroTop format. @@ -35,13 +35,13 @@ def test_crystal_water(system, match_water, function): if match_water: num_matches = 0 else: - num_matches = len(system.search("resname COF").molecules()) + num_matches = len(kigaki_system.search("resname COF").molecules()) # Create the box parameters. box, angles = BSS.Box.cubic(5.5 * BSS.Units.Length.nanometer) # Create the solvated system. - solvated = function(system, box, angles, match_water=match_water) + solvated = function(kigaki_system, box, angles, match_water=match_water) # Search for the crystal waters in the solvated system. try: diff --git a/tests/Stream/test_stream.py b/tests/Stream/test_stream.py index 12f169dec..73788bd7e 100644 --- a/tests/Stream/test_stream.py +++ b/tests/Stream/test_stream.py @@ -4,11 +4,6 @@ import BioSimSpace as BSS -@pytest.fixture -def system(scope="session"): - return BSS.IO.readMolecules(["tests/input/ala.top", "tests/input/ala.crd"]) - - @pytest.fixture(autouse=True) def run_around_tests(): yield diff --git a/tests/Trajectory/test_trajectory.py b/tests/Trajectory/test_trajectory.py index b03a554ed..08457529d 100644 --- a/tests/Trajectory/test_trajectory.py +++ b/tests/Trajectory/test_trajectory.py @@ -13,12 +13,6 @@ def wrap(arg): from tests.conftest import url, has_mdanalysis, has_mdtraj -@pytest.fixture(scope="session") -def system(): - """A system object with the same topology as the trajectories.""" - return BSS.IO.readMolecules(["tests/input/ala.top", "tests/input/ala.crd"]) - - @pytest.fixture(scope="session") def traj_sire(system): """A trajectory object using the Sire backend.""" diff --git a/tests/_SireWrappers/test_molecule.py b/tests/_SireWrappers/test_molecule.py index 28b6c28df..8152c038a 100644 --- a/tests/_SireWrappers/test_molecule.py +++ b/tests/_SireWrappers/test_molecule.py @@ -5,11 +5,6 @@ from tests.conftest import url, has_amber -@pytest.fixture(scope="session") -def system(): - return BSS.IO.readMolecules(["tests/input/ala.top", "tests/input/ala.crd"]) - - @pytest.mark.skipif(has_amber is False, reason="Requires AMBER to be installed.") def test_makeCompatibleWith(): # Load the original PDB file. In this representation the system contains diff --git a/tests/_SireWrappers/test_search_result.py b/tests/_SireWrappers/test_search_result.py index b38cb40c2..55aeda098 100644 --- a/tests/_SireWrappers/test_search_result.py +++ b/tests/_SireWrappers/test_search_result.py @@ -6,12 +6,6 @@ url = BSS.tutorialUrl() -@pytest.fixture(scope="session") -def system(): - """Re-use the same molecuar system for each test.""" - return BSS.IO.readMolecules(["tests/input/ala.top", "tests/input/ala.crd"]) - - @pytest.fixture(scope="session") def molecule(system): """Re-use the same molecule for each test.""" diff --git a/tests/_SireWrappers/test_system.py b/tests/_SireWrappers/test_system.py index 81ff55a5f..2073f4b26 100644 --- a/tests/_SireWrappers/test_system.py +++ b/tests/_SireWrappers/test_system.py @@ -8,12 +8,6 @@ from tests.conftest import url, has_amber, has_openff -@pytest.fixture(scope="session") -def system(): - """Re-use the same molecuar system for each test.""" - return BSS.IO.readMolecules(["tests/input/ala.top", "tests/input/ala.crd"]) - - @pytest.fixture(scope="session") def rna_system(): """An RNA system for re-use.""" diff --git a/tests/conftest.py b/tests/conftest.py index 2ef298ea0..46be2d10a 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -1,6 +1,7 @@ collect_ignore_glob = ["*/out_test*.py"] import os +import pytest from pathlib import Path @@ -55,3 +56,33 @@ # Allow tests to be run from any directory. root_fp = Path(__file__).parent.resolve() + +# Fixtures for tests. + + +@pytest.fixture(scope="session") +def system(): + """Solvated alanine dipeptide system.""" + return BSS.IO.readMolecules(["tests/input/ala.top", "tests/input/ala.crd"]) + + +@pytest.fixture(scope="module") +def perturbable_system(): + """A vacuum perturbable system.""" + return BSS.IO.readPerturbableSystem( + f"{url}/perturbable_system0.prm7", + f"{url}/perturbable_system0.rst7", + f"{url}/perturbable_system1.prm7", + f"{url}/perturbable_system1.rst7", + ) + + +@pytest.fixture(scope="module") +def solvated_perturbable_system(): + """A solvated perturbable system.""" + return BSS.IO.readPerturbableSystem( + f"{url}/solvated_perturbable_system0.prm7", + f"{url}/solvated_perturbable_system0.rst7", + f"{url}/solvated_perturbable_system1.prm7", + f"{url}/solvated_perturbable_system1.rst7", + )