OpenBioSim · lohedges · Oct 13, 2023 · Jun 29, 2023 · Jun 30, 2023 · Jul 3, 2023
diff --git a/.github/ISSUE_TEMPLATE/bug_report.md b/.github/ISSUE_TEMPLATE/bug_report.md
@@ -16,6 +16,9 @@ Steps to reproduce the behavior:
 2. Run the code '...'
 3. This is the exception that was raised / this is what went wrong.
 
+(If possible, please quote code snippets using the markdown formatting
+described [here](https://docs.github.com/en/get-started/writing-on-github/getting-started-with-writing-and-formatting-on-github/basic-writing-and-formatting-syntax#quoting-code). This makes it easy to copy and paste the example, avoiding any translation errors.)
+
 **Expected behavior**
 A clear and concise description of what you expected to happen
 (or, it should not raise an exception if an exception was raised)

diff --git a/.github/workflows/Sandpit_exs.yml b/.github/workflows/Sandpit_exs.yml
@@ -19,7 +19,7 @@ jobs:
       fail-fast: false
       matrix:
         os: ["ubuntu-latest", "macOS-latest",]
-        python-version: ["3.8", ]
+        python-version: ["3.10",]
 
     steps:
       - uses: actions/checkout@v2
@@ -35,12 +35,10 @@ jobs:
 
       - name: Install dependency
         run: |
-          mamba install -c conda-forge -c openbiosim/label/main biosimspace python=3.8 ambertools gromacs "sire=2023.2.2" alchemlyb pytest pyarrow "typing-extensions!=4.6.0" openff-interchange pint=0.21 "rdkit<2023" "jaxlib>0.3.7" tqdm
+          mamba install -c conda-forge -c openbiosim/label/main biosimspace python=3.10 ambertools gromacs "sire=2023.3" "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
-          # Before the new alchemlyb is released
-          python -m pip install git+https://github.com/alchemistry/alchemlyb.git
 
       - name: Install the dev version
         run: |

diff --git a/.github/workflows/main.yaml b/.github/workflows/main.yaml
@@ -16,7 +16,7 @@ jobs:
       max-parallel: 9
       fail-fast: false
       matrix:
-        python-version: ["3.8", "3.9", "3.10"]
+        python-version: ["3.9", "3.10", "3.11"]
         platform:
           - { name: "windows", os: "windows-latest", shell: "pwsh" }
           - { name: "linux", os: "ubuntu-latest", shell: "bash -l {0}" }

diff --git a/doc/source/api/index_Box.rst b/doc/source/api/index_Box.rst
@@ -11,11 +11,7 @@ To support triclinic boxes that work across the range of molecular simulation
 engines that support, we represent the triclinic space in reduced form, using
 the approach documented in Appendix A of Chapter 3 from
 `"Molecular dynamics of sense and sensibility in processing and analysis of data" <https://research.rug.nl/files/2839528/01_c1.pdf>`__
-by Tsjerk A. Wassenaar. Due to the fixed-width format used to represent box
-dimensions and angles in the various molecular input files, repeated reading
-and writing can lead to oscillation of the box angles on reduction due to
-rounding precision errors. To account for this, we add a small bias so that
-we always round in a consistent direction.
+by Tsjerk A. Wassenaar.
 
 .. automodule:: BioSimSpace.Box
 

diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst
@@ -9,6 +9,22 @@ company supporting open-source development of fostering academic/industrial coll
 within the biomolecular simulation community. Our software is hosted via the `OpenBioSim`
 `GitHub <https://github.com/OpenBioSim/biosimspace>`__ organisation.
 
+`2023.4.0 <https://github.com/openbiosim/biosimspace/compare/2023.3.1...2023.4.0>`_ - Oct 13 2023
+-------------------------------------------------------------------------------------------------
+
+* Add support for computing trajectory RMSDs using Sire backend (`#152 <https://github.com/OpenBioSim/biosimspace/pull/152>`__).
+* Add support for setting up systems containing crystal waters (`#154 <https://github.com/OpenBioSim/biosimspace/pull/154>`__).
+* Add unified free-energy perturbation analysis using ``alchemlyb`` (`@annamherz <https://github.com/annamherz>`_) (`#155 <https://github.com/OpenBioSim/biosimspace/pull/155>`__).
+* Fix handling of connectivity changes during molecular perturbations (`#157 <https://github.com/OpenBioSim/biosimspace/pull/157>`__).
+* Fix issues related to new shared properties in Sire (`#160 <https://github.com/OpenBioSim/biosimspace/pull/160>`__).
+* Fix issues in SOMD perturbation files for absolute binding free-energy simulations (`@fkclark <https://github.com/fjclark>`_) (`#164 <https://github.com/OpenBioSim/biosimspace/pull/164>`__).
+* Don't generate velocities when performing a continuation with GROMACS (`#169 <https://github.com/OpenBioSim/biosimspace/pull/169>`__).
+* Decouple custom parameters and additional commands in ``LEaP`` input (`#170 <https://github.com/OpenBioSim/biosimspace/pull/170>`__).
+* Check for periodic space when updating box information from restart file or trajectory (`#173 <https://github.com/OpenBioSim/biosimspace/pull/173>`__).
+* Add functionality to allow manual rotation and reduction of triclinic boxes, rather than performing automatically on read (`#175 <https://github.com/OpenBioSim/biosimspace/pull/175>`__).
+* Allow unit-based protocol options to be passed as strings (`#179 <https://github.com/OpenBioSim/biosimspace/pull/179>`__).
+* Fix assignment of ``gpu`` configuration option for SOMD (`#181 <https://github.com/OpenBioSim/biosimspace/pull/181>`__).
+
 `2023.3.1 <https://github.com/openbiosim/biosimspace/compare/2023.3.0...2023.3.1>`_ - Aug 14 2023
 -------------------------------------------------------------------------------------------------
 

diff --git a/doc/source/contributors.rst b/doc/source/contributors.rst
@@ -21,4 +21,3 @@ can be recognised.
 * `@msuruzon <https://github.com/msuruzhon>`__
 * `@xiki-tempula <https://github.com/xiki-tempula>`__
 * `@kexul <https://github.com/kexul>`__
-* `@mb2055 <https://github.com/mb2055>`__
diff --git a/doc/source/tutorials/crystal_water.rst b/doc/source/tutorials/crystal_water.rst
@@ -0,0 +1,135 @@
+.. _ref_crystal_water:
+=============
+Crystal water
+=============
+
+.. toctree::
+   :maxdepth: 2
+
+Water molecules are important for the structure, function, and dynamics of
+proteins. When setting up a simulation it may be important to preserve
+specific water molecules, often referred to as *crystal waters*. For example,
+many structures from the `Protein Data Bank <https://www.rcsb.org/>`__
+contain the coordinates of water molecules that are resolvable via X-ray
+crystallography. In this example, we will use BioSimSpace to parameterise
+and solvate a `tyrosine kinase 2 <https://en.wikipedia.org/wiki/Tyrosine_kinase_2>`__
+protein structure containing two crystal waters within its binding site.
+
+Firstly, let's load the PDB structure:
+
+>>> import BioSimSpace as BSS
+>>> tyk2_xtal = BSS.IO.readMolecules(
+...    BSS.IO.expand(
+...        BSS.IO.tutorialUrl(), "tyk2_xtal.pdb"
+...    )
+...)[0]
+
+This is a single molecule that contains the oxgyen atoms of two crystal waters
+at the end of the structure:
+
+>>> for residue in tyk2_xtal.getResidues()[-2:]:
+...     for atom in residue:
+...         print(atom, atom.coordinates())
+...
+<BioSimSpace.Atom: name='O', molecule=2, index=4652> (-4.5180 A, -0.0520 A, -15.6070 A)
+<BioSimSpace.Atom: name='O', molecule=2, index=4653> (-7.2040 A, -8.5540 A, -20.7520 A)
+
+We can use BioSimSpace to search for and extract the protein and water
+components of the system. In this case, the waters are part of a residue
+named ``WAT``, so we can use this as our search term:
+
+>>> tyk2 = tyk2_xtal.extract(
+...    [atom.index() for atom in tyk2_xtal.search("not resname WAT").atoms()]
+... )
+>>> xtal_water = tyk2_xtal.extract(
+...    [atom.index() for atom in tyk2_xtal.search("resname WAT").atoms()]
+... )
+
+Now we have the components, we can parameterise them individually. First the
+protein:
+
+>>> tyk2 = BSS.Parameters.ff14SB(tyk2).getMolecule()
+
+And then the water:
+
+>>> xtal_water = BSS.Parameters.tip3p(
+...     xtal_water,
+...     water_model="tip3p",
+...     ensure_compatible=False
+... ).getMolecule()
+
+In this case we need to specify the desired water model, and also that we don't
+want to ensure that the parameterised molecule is compatible with the topology
+of the molecule that we passed in. This is because we are only passing in
+oxygen atoms, so ``tLEaP`` will add the missing hydrogens.
+
+.. note ::
+
+   By default, BioSimSpace ensures that the topology of the parameterised molecule(s)
+   matches that of the input molecule and will raise an exception if this is not
+   the case. This is because the input system is often required as a reference by
+   the user, e.g. they might want to preserve a specific naming convention, chain
+   identifiers, etc.
+
+Let's check one of the parameterised crystal waters:
+
+>>> for atom in xtal_water[0].getAtoms():
+...     print(atom, atom.coordinates())
+...
+<BioSimSpace.Atom: name='O', molecule=9, index=0> (-4.5180 A, -0.0520 A, -15.6070 A)
+<BioSimSpace.Atom: name='H1', molecule=9, index=1> (-3.5608 A, -0.0520 A, -15.6070 A)
+<BioSimSpace.Atom: name='H2', molecule=9, index=2> (-4.7580 A, 0.8746 A, -15.6070 A)
+
+Note that the coordinates of the oxygen atom are preserved.
+
+Once both components are parameterised, we can combine them together to create
+a new system:
+
+>>> system = tyk2.toSystem() + xtal_water
+
+The next step in our setup procedure is to solvate the system in a water box.
+Here we will use a truncated octahedral box with a base length of 7 nanometers:
+
+>>> box, angles = BSS.Box.truncatedOctahedron(7 * BSS.Units.Length.nanometer)
+>>> solvated = BSS.Solvent.tip3p(system, box=box, angles=angles)
+
+Since the protein is charged, during the solvation process ``gmx genion`` will
+have been used to neutralise the solvated system by adding counter ions. We can
+check this as follows:
+
+>>> print(system.charge(), solvated.charge())
+... -2.0000 |e| -1.5091e-07 |e|
+>>> print(len(solvated.search("element Na")))
+... 2
+
+In this case, two sodium ions have been added to neutralise the system. In
+doing so, ``gmx genion`` will have removed two *random* water molecules. In
+order to ensure that crystal waters are not removed, we temporarily tag them
+with a unique residue and molecule name during solvation.
+
+To check that they are preserved we can re-solvate the system, asking to
+preserve the water naming for the existing waters in the system. (By default,
+they will be updated to the default ``GROMACS`` naming convention.)
+
+>>> solvated_no_match = BSS.Solvent.tip3p(
+...     system,
+...     box=box,
+...     angles=angles,
+...     match_water=False
+... )
+
+We can check that the crystal waters are still present as follows:
+
+>>> for residue in solvated_no_match.getWaterMolecules().toSystem().search("resname WAT").residues():
+...     print(residue)
+<BioSimSpace.Residue: name='WAT', molecule=8, index=0, nAtoms=3>
+<BioSimSpace.Residue: name='WAT', molecule=9, index=0, nAtoms=3>
+
+.. note ::
+
+   When solvating, molecules in the original system will be centered within the
+   solvation box, hence the coordinates of the crystal waters after solvation
+   won't necessarily match those from before. By setting ``match_water=False``,
+   it is possible to preserve the naming convention used for any existing
+   waters in the system. However, note that some simulation engines rely on a
+   specific naming convention for water molecules and a single solvent group.
diff --git a/doc/source/tutorials/index.rst b/doc/source/tutorials/index.rst
@@ -20,5 +20,6 @@ please :doc:`ask for support. <../support>`
 .. toctree::
    :maxdepth: 1
 
+   crystal_water
    hydration_freenrg
    metadynamics
diff --git a/nodes/playground/prepareFEP.ipynb b/nodes/playground/prepareFEP.ipynb
@@ -159,6 +159,13 @@
     "    BSS.Gateway.String(\n",
     "        help=\"The root name for the files describing the perturbation input1->input2.\"\n",
     "    ),\n",
+    ")\n",
+    "node.addInput(\n",
+    "    \"somd2\",\n",
+    "    BSS.Gateway.Boolean(\n",
+    "        help=\"Whether to generate input for SOMD2.\",\n",
+    "        default=False,\n",
+    "    ),\n",
     ")"
    ]
   },
@@ -321,18 +328,28 @@
    "source": [
     "# Log the mapping used\n",
     "writeLog(lig1, lig2, mapping)\n",
+    "# Are we saving output for SOMD2?\n",
+    "is_somd2 = node.getInput(\"somd2\")\n",
+    "# File root for all output.\n",
+    "root = node.getInput(\"output\")\n",
     "BSS.IO.saveMolecules(\n",
     "    \"merged_at_lam0.pdb\",\n",
     "    merged,\n",
     "    \"PDB\",\n",
     "    {\"coordinates\": \"coordinates0\", \"bond\": \"bond0\", \"element\": \"element0\"},\n",
     ")\n",
-    "# Generate package specific input\n",
-    "protocol = BSS.Protocol.FreeEnergy(runtime=2 * BSS.Units.Time.femtosecond, num_lam=3)\n",
-    "process = BSS.Process.Somd(system1, protocol)\n",
-    "process.getOutput()\n",
-    "with zipfile.ZipFile(\"somd_output.zip\", \"r\") as zip_hnd:\n",
-    "    zip_hnd.extractall(\".\")"
+    "if is_somd2:\n",
+    "    BSS.Stream.save(system1, root)\n",
+    "    stream_file = \"%s.bss\" % root\n",
+    "else:\n",
+    "    # Generate package specific input\n",
+    "    protocol = BSS.Protocol.FreeEnergy(\n",
+    "        runtime=2 * BSS.Units.Time.femtosecond, num_lam=3\n",
+    "    )\n",
+    "    process = BSS.Process.Somd(system1, protocol)\n",
+    "    process.getOutput()\n",
+    "    with zipfile.ZipFile(\"somd_output.zip\", \"r\") as zip_hnd:\n",
+    "        zip_hnd.extractall(\".\")"
    ]
   },
   {
@@ -341,12 +358,12 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "root = node.getInput(\"output\")\n",
-    "mergedpdb = \"%s.mergeat0.pdb\" % root\n",
-    "pert = \"%s.pert\" % root\n",
-    "prm7 = \"%s.prm7\" % root\n",
-    "rst7 = \"%s.rst7\" % root\n",
-    "mapping_str = \"%s.mapping\" % root"
+    "if not is_somd2:\n",
+    "    mergedpdb = \"%s.mergeat0.pdb\" % root\n",
+    "    pert = \"%s.pert\" % root\n",
+    "    prm7 = \"%s.prm7\" % root\n",
+    "    rst7 = \"%s.rst7\" % root\n",
+    "    mapping_str = \"%s.mapping\" % root"
    ]
   },
   {
@@ -355,18 +372,19 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "os.replace(\"merged_at_lam0.pdb\", mergedpdb)\n",
-    "os.replace(\"somd.pert\", pert)\n",
-    "os.replace(\"somd.prm7\", prm7)\n",
-    "os.replace(\"somd.rst7\", rst7)\n",
-    "os.replace(\"somd.mapping\", mapping_str)\n",
-    "try:\n",
-    "    os.remove(\"somd_output.zip\")\n",
-    "    os.remove(\"somd.cfg\")\n",
-    "    os.remove(\"somd.err\")\n",
-    "    os.remove(\"somd.out\")\n",
-    "except Exception:\n",
-    "    pass"
+    "if not is_somd2:\n",
+    "    os.replace(\"merged_at_lam0.pdb\", mergedpdb)\n",
+    "    os.replace(\"somd.pert\", pert)\n",
+    "    os.replace(\"somd.prm7\", prm7)\n",
+    "    os.replace(\"somd.rst7\", rst7)\n",
+    "    os.replace(\"somd.mapping\", mapping_str)\n",
+    "    try:\n",
+    "        os.remove(\"somd_output.zip\")\n",
+    "        os.remove(\"somd.cfg\")\n",
+    "        os.remove(\"somd.err\")\n",
+    "        os.remove(\"somd.out\")\n",
+    "    except Exception:\n",
+    "        pass"
    ]
   },
   {
@@ -375,7 +393,11 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "node.setOutput(\"nodeoutput\", [mergedpdb, pert, prm7, rst7, mapping_str])"
+    "if is_somd2:\n",
+    "    output = [stream_file]\n",
+    "else:\n",
+    "    output = [mergedpdb, pert, prm7, rst7, mapping_str]\n",
+    "node.setOutput(\"nodeoutput\", output)"
    ]
   },
   {
@@ -390,7 +412,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
    "name": "python3"
   },
@@ -404,7 +426,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.3"
+   "version": "3.10.11"
   }
  },
  "nbformat": 4,