04_molecular_setup.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Author: Lester Hedges and Christopher Woods<br>\n",
    "Email:&nbsp;&nbsp; lester.hedges@bristol.ac.uk, christopher.woods@bristol.ac.uk\n",
    "\n",
    "___Jupyter Recap___:\n",
    "* Press Shift+Enter to execute a cell and move to the cell below.\n",
    "* Press Ctrl+Enter to execute a cell and remain in that cell.\n",
    "* Run a shell command on the underlying operating system by prefixing the command with an exclamation mark, !\n",
    "* Remember that the flow is in the order that you execute cells, which is not necessarily linear in the notebook. Keep track of the numbers in brackets to the left of the cell!\n",
    "\n",
    "# Molecular setup\n",
    "\n",
    "So far we have learned how to use BioSimSpace to run molecular dynamics simulations. In this notebook you'll work with molecules, and learn how to use BioSimSpace to parameterise and solvate them ready for simulation. As always, BioSimSpace wraps around existing tools: AmberTools [tLEaP](http://ambermd.org/tutorials/pengfei/index.htm) and GROMACS's [pdb2gmx](http://manual.gromacs.org/2019/onlinehelp/gmx-pdb2gmx.html) for parameterisation, and GROMACS's [solvate](http://manual.gromacs.org/documentation/2019/onlinehelp/gmx-solvate.html) for solvation. We don't yet provide tools for preparing proteins, since this still requires too much human input and intuition to be encoded in a robust and reproducible way.\n",
    "\n",
    "To begin with we'll need to import BioSimSpace:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import BioSimSpace as BSS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's load in benzene molecule. As you learned earlier, the `BSS.IO.readMolecules` function loads molecular information from file and returns a `System` object. This object can be indexed to extract the molecules that it contains. Since we know that the system contains a single benzene molecule, we simply take the first molecule."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "molecule = BSS.IO.readMolecules(\"input/molecules/benzene.pdb\")[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(Just like the BioSimSpace `System`, a `Molecule` is simply a thin wrapper around the molecule class from [Sire](https://github.com/michellab/Sire). For those that are familiar with Sire, you can get access to the underlying object directly using `molecule._getSireObject()`. This private method is hidden from the user.)\n",
    "\n",
    "Now let's paramterise the molecule. We do so by calling the `parameterise` function from the `BSS.Parameters` package, passing the molecule and force field name as arguments. Since parameterisation can be slow, the function returns a handle to a process that runs the parameterisation in the background. To get the parameterised molecule from the process we need to call the `getMolecule` method. This is a blocking operation which waits for the process to finish before grabbing the molecule and returning it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "molecule = BSS.Parameters.parameterise(molecule, \"gaff2\").getMolecule()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(Note that there is also a convenience function, `BSS.Parameters.gaff2` which could be used too. However, it's useful to use the general `parameterise` function in a node since the user might have passed the name of the force field as an input parameter.)\n",
    "\n",
    "If you find that parameterisation fails, It is a good idea to create a parameterisation process, then call `getOutput` on the process in order to get all of the output files, which can then be examined for errors. For example:\n",
    "\n",
    "```python\n",
    "process = BSS.Parameters.parameterise(molecule, \"gaff2\")\n",
    "\n",
    "# This could throw an error.\n",
    "molecule = process.getOutput()\n",
    "\n",
    "# Get the output files of the process.\n",
    "process.getOutput()\n",
    "```\n",
    "\n",
    "(In the example in the previous cell the intermediate process goes out of scope when we call `getMolecule`.)\n",
    "\n",
    "We can see that it is parameterised by printing out all of the atoms. Here we use some functionality from the underlying [Sire](https://github.com/michellab/Sire) object that is normally hidden from the user. Note how symmetrical atoms have symmetrical charges."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for atom in molecule._getSireObject().atoms():\n",
    "    print(atom, atom.property(\"charge\"), atom.property(\"LJ\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we will solvate our molecule in a box of water using the `solvate` function from the `BSS.Solvent` package. This will centre the molecule in an orthorhombic box of a specified size (for simplicity we'll use a cubic box) and surround it by water molecules.\n",
    "\n",
    "Note that the molecule is an optional `keyword` argument to the solvate function. This is because its also possible to create a pure water box, i.e. without any molecules in it.\n",
    "\n",
    "(As for the paremeterisation, there is are also convenience function `BSS.Solvent.tip3p` that we could use instead.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "system = BSS.Solvent.solvate(\n",
    "    \"tip3p\", molecule=molecule, box=3 * [4 * BSS.Units.Length.nanometer]\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now write the output as GROMACS format files representing the parameterised and solvated molecular system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "BSS.IO.saveMolecules(\"solvated_benzene\", system, [\"grotop\", \"gro87\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!head -n 50 solvated_benzene.top"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now take a look at the solvated structure:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "BSS.Notebook.View([\"solvated_benzene.top\", \"solvated_benzene.gro\"]).system()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok, that was a small molecule. How about something bigger? Let's now load up a protein:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "molecule = BSS.IO.readMolecules(\"input/molecules/2JJC.pdb\")[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will parameterise this using FF14SB:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "molecule = BSS.Parameters.parameterise(molecule, \"ff14SB\").getMolecule()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for atom in molecule._getSireObject().atoms():\n",
    "    print(atom, atom.property(\"charge\"), atom.property(\"LJ\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's solvate this in a box of TIP4P:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "system = BSS.Solvent.solvate(\n",
    "    \"tip3p\", molecule=molecule, box=3 * [5 * BSS.Units.Length.nanometer]\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will save this to a Amber PRM7, NETCDF and PDB formats:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "BSS.IO.saveMolecules(\"solvated_protein\", system, [\"prm7\", \"rst\", \"pdb\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!head -500 solvated_protein.prm7"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What does this look like?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "BSS.Notebook.View(\"solvated_protein.pdb\").system()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To view the water molecules, click on the `Representation` tab, then `Add` and select `line` in the dropdown. Next type `water` in the selection box to the right where an asterisk is visible, then hit enter.\n",
    "\n",
    "\n",
    "# Excercise\n",
    "\n",
    "Write a node to parameterise and solvate a molecule, saving it to AMBER format files ready for simulation. As in the previous exercises, remember that functions like `BSS.Parameters.forceFields` are useful for automatically generating allowed options in the GUI and command-line argument parser. Start with parameterising the molecule, then add functionality to solvate it if you have time. Make sure to use the help, e.g. `help(BSS.Parameters.parameterise)` and `help(BSS.Solvent.solvate)` to see what options you might want to support. If you get stuck, a solution can be found [here](nodes/molecular_setup.ipynb)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}