01_introduction/03_molecular_dynamics.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "encouraging-talent",
   "metadata": {},
   "source": [
    "Author: Lester Hedges<br>\n",
    "Email:&nbsp;&nbsp; lester.hedges@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",
    "\n",
    "# Molecular dynamics\n",
    "\n",
    "## Introduction\n",
    "\n",
    "In this section we will learn how to use BioSimSpace to configure and run some basic molecular dynamics simulations.\n",
    "\n",
    "\n",
    "## Protocols\n",
    "\n",
    "One of the key goals of BioSimSpace was to start a conversation regarding _best practice_ within the biomolecular simulation community and to facilitate the codifying of shareable, re-usable, and extensible simulation protocols.\n",
    "\n",
    "The [BioSimSpace.Protocol](https://biosimspace.org/api/index_Protocol.html) package defines protocols for a range of common molecular dynamics simulations. We can query the package to ee what protocols are available:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "objective-introduction",
   "metadata": {},
   "outputs": [],
   "source": [
    "import BioSimSpace as BSS\n",
    "BSS.Protocol.protocols()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "colonial-shark",
   "metadata": {},
   "source": [
    "Since we require protocols to be _interoperable_, the classes listed above are simple objects that allow you to configure a _limited_ set of options that are handled by _all_ of the molecular dynamics engines that we support. This might seem quite restrictive, but we will see later how it is possible to fully customise a simulation for a particular molecular dynamics engine.\n",
    "\n",
    "Each protocol comes with some default options. To see what those are we can instantiate an object using the default constructor. For example, let's explore the [Equilibration](https://biosimspace.org/api/generated/BioSimSpace.Protocol.Equilibration.html#BioSimSpace.Protocol.Equilibration) protocol."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "mighty-easter",
   "metadata": {},
   "outputs": [],
   "source": [
    "protocol = BSS.Protocol.Equilibration()\n",
    "print(protocol)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "opposed-control",
   "metadata": {},
   "source": [
    "Here we can see that the default protocol performs an equlibration at fixed temperature (`temperature_start == temperature_end`) in the NVT ensemble (`pressure=None`) with no restraints (`restraint=None`). The total simulation time is 0.2 nanoseconds with an integration timestep of 2 femtoseconds. The `report_interval` and `restart_interval` govern how frequently information is written to log and restart (and/or trajectory) files respectively.\n",
    "\n",
    "If any of these defaults are unsuitable, then you are free to change them by passing in appropriate values for each of the arguments when instantiating the object. In some cases it might be desirable to _override_ the default protocols and set specific values of the arguments that are suitable for a particular project or team. This can be achieved by defining a set of function wrappers that configure and return the protocols using your own defaults.\n",
    "\n",
    "As an example, the following configuration could be used to provide alternative defaults for NVT and NPT equlibration protocols. (You could simply re-use the existing protocol name, but here we provide two separate protocols for convenience.)\n",
    "\n",
    "```python\n",
    "# myconfig/Protocol.py\n",
    "import BioSimSpace as BSS\n",
    "\n",
    "# Override the equilibration protocol with some custom defaults. Ideally all\n",
    "# arguments to the BioSimSpace function would be mapped, but here we use a\n",
    "# subset for simplicity.\n",
    "\n",
    "# A custom equlibration in the NVT ensemble.\n",
    "def EquilibrationNVT(runtime=5*BSS.Units.Time.nanosecond,\n",
    "                     report_interval=2500,\n",
    "                     restart_interval=250000,\n",
    "                     restraint=\"backbone\"):\n",
    "    return BSS.Protocol.Equilibration(runtime=runtime,\n",
    "                                      report_interval=report_interval,\n",
    "                                      restart_interval=restart_interval,\n",
    "                                      restraint=restraint)\n",
    "\n",
    "# A custom equlibration in the NPT ensemble.\n",
    "def EquilibrationNPT(runtime=5*BSS.Units.Time.nanosecond,\n",
    "                     pressure=BSS.Units.Pressure.atm,\n",
    "                     report_interval=2500,\n",
    "                     restart_interval=25000,\n",
    "                     restraint=\"backbone\"):\n",
    "    return BSS.Protocol.Equilibration(runtime=runtime,\n",
    "                                      pressure=pressure,\n",
    "                                      report_interval=report_interval,\n",
    "                                      restart_interval=restart_interval,\n",
    "                                      restraint=restraint)\n",
    "\n",
    "```\n",
    "\n",
    "We could then import the customised protocols from our local configuraton and use them instead, e.g.:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bridal-pakistan",
   "metadata": {},
   "outputs": [],
   "source": [
    "from myconfig.Protocol import *\n",
    "\n",
    "protocol = EquilibrationNVT()\n",
    "print(protocol)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "delayed-shuttle",
   "metadata": {},
   "source": [
    "## Processes\n",
    "\n",
    "Once you have created a molecular system and chosen a protocol, then it is time to create a simulation _process_. The [BioSimSpace.Process](https://biosimspace.org/api/index_Process.html) package provides functionality for configuring and running processes with several common molecular dynamics engines.\n",
    "\n",
    "Let's query the package to see what engines are available:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "steady-munich",
   "metadata": {},
   "outputs": [],
   "source": [
    "BSS.Process.engines()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "still-reggae",
   "metadata": {},
   "source": [
    "Before creating a process let us once again load our example alanine-dipeptide system from file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "paperback-network",
   "metadata": {},
   "outputs": [],
   "source": [
    "system = BSS.IO.readMolecules(\"inputs/ala*\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "southern-breach",
   "metadata": {},
   "source": [
    "As a simple example, let us use a short minimisation protocol:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "little-charter",
   "metadata": {},
   "outputs": [],
   "source": [
    "protocol = BSS.Protocol.Minimisation(steps=1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "coral-retro",
   "metadata": {},
   "source": [
    "We'll now create a process to apply the `Protocol` to the `System` using the AMBER molecular dynamics engine:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "registered-license",
   "metadata": {},
   "outputs": [],
   "source": [
    "process = BSS.Process.Amber(system, protocol)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "threatened-typing",
   "metadata": {},
   "source": [
    "A lot of complexity is hidden in this line. BioSimSpace has automatically found an AMBER executable on the underlying operating system, has automatically written AMBER format molecular input files, generated an AMBER configuration file for the minimisation protocol, and configured any command-line arguments that are required.\n",
    "\n",
    "By default, processes are run inside of a temporary working directory hidden away from the user. To see where this is, run:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "discrete-table",
   "metadata": {},
   "outputs": [],
   "source": [
    "process.workDir()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "robust-peeing",
   "metadata": {},
   "source": [
    "N.B. If you want to use a different temporary directory, e.g. one with a faster disk, then simply set the `TMPDIR` environment variable. Alternatively, you can pass the `work_dir` argument to the `Process` constructor to explicitly specify the path. This can be useful when you want named directories, or want to examine the intermediate files from the `Process` for debugging purposes.\n",
    "\n",
    "To see what executable was found, run:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "opposed-sleeve",
   "metadata": {},
   "outputs": [],
   "source": [
    "process.exe()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "neural-split",
   "metadata": {},
   "source": [
    "To see the list of autogenerated input files:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "exposed-broadcasting",
   "metadata": {},
   "outputs": [],
   "source": [
    "process.inputFiles()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "isolated-longer",
   "metadata": {},
   "source": [
    "If you like, we could zip up the input files to use on another occasion. When working on a notebook server it's possible to return a file link so that we can download them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "informal-international",
   "metadata": {},
   "outputs": [],
   "source": [
    "process.getInput(file_link=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "empty-geneva",
   "metadata": {},
   "source": [
    "We can query also query the list of configuration file options:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dominican-colombia",
   "metadata": {},
   "outputs": [],
   "source": [
    "process.getConfig()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "seeing-disaster",
   "metadata": {},
   "source": [
    "And also get command-line argument string for the process:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "victorian-gabriel",
   "metadata": {},
   "outputs": [],
   "source": [
    "process.getArgString()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "olympic-hearts",
   "metadata": {},
   "source": [
    "If you're an expert in a particular package then BioSimSpace allows you to fully customise the process by tweaking the configuration options and command-line arguments. Read the help documentation for `process.setConfig` and `process.setArgs` if you are interested. Once again, it's possible to wrap the instantiation of `Process` objects in your own custom functions, allowing you to tweak the default configuration options for your own requirements. For example, if you always want to wrap coordinates to the minimum image when using AMBER, then this could be achieved as follows:\n",
    "\n",
    "```python\n",
    "# myconfig/Process.py\n",
    "import BioSimSpace as BSS\n",
    "\n",
    "# Wrap the instantiation of BSS.Process.Amber objects to configure them\n",
    "# such that coordinates are always wrapped to the minimum image.\n",
    "def Amber(system, protocol, exe=None, name=\"amber\",\n",
    "            work_dir=None, seed=None, property_map={}):\n",
    "\n",
    "    # Create process using the passed parameters.\n",
    "    process = BSS.Process.Amber(system,\n",
    "                                protocol,\n",
    "                                exe=exe,\n",
    "                                name=name,\n",
    "                                work_dir=work_dir,\n",
    "                                seed=seed,\n",
    "                                property_map=property_map)\n",
    "    \n",
    "    # Get the config.\n",
    "    config = process.getConfig()\n",
    "    \n",
    "    # Add coordinate wrapping to the end of the config.\n",
    "    config[-1] = \"  iwrap=1,\"\n",
    "    config.append(\" /\")\n",
    "\n",
    "    # Set the new config.\n",
    "    process.setConfig(config)\n",
    "    \n",
    "    # Return the process.\n",
    "    return process\n",
    "```\n",
    "\n",
    "Let us know create our custom AMBER process and check the configuration:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "unsigned-policy",
   "metadata": {},
   "outputs": [],
   "source": [
    "from myconfig.Process import *\n",
    "\n",
    "process = Amber(system, protocol)\n",
    "process.getConfig()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "considered-trinidad",
   "metadata": {},
   "source": [
    "N.B. You might want to add additional configuration details to your `Process` wrappers, e.g. to ensure that a specific executable is used.\n",
    "\n",
    "Now that we have a process, let's go ahead and start it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "taken-compact",
   "metadata": {},
   "outputs": [],
   "source": [
    "process.start()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "inappropriate-disorder",
   "metadata": {},
   "source": [
    "BioSimSpace has now launched a minimisation process in the background! When in an interactive session you carry on working and periodically check in on the process to see how its doing.\n",
    "\n",
    "To check whether the process is running:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "driven-melissa",
   "metadata": {},
   "outputs": [],
   "source": [
    "process.isRunning()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "assigned-reliance",
   "metadata": {},
   "source": [
    "We can see how many minutes it has been running for:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "understood-vehicle",
   "metadata": {},
   "outputs": [],
   "source": [
    "process.runTime()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "finite-abortion",
   "metadata": {},
   "source": [
    "Since this is a short minimisation it will likely finish pretty quickly. Let's print the final energy of the system and return the minimised molecular configuration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "european-cabinet",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(process.getTotalEnergy(block=True))\n",
    "minimised = process.getSystem()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "answering-oxford",
   "metadata": {},
   "source": [
    "When working interactively, any time we query a running process we get back the _latest_ information that has been written to disk. This means that we can get an update on how things are progressing, then immediately carry on with what we were doing in our notebook. By passing `block=True`, as we do when we call `getTotalEnergy` above, we request that the process finishes running before returning a result. This means we get the _final_ energy, and the minimised system that is returned afterwards represents the _final_ snapshot that was saved.\n",
    "\n",
    "Let's now re-run the simulation, instead using GROMACS as the MD engine."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ongoing-revision",
   "metadata": {},
   "outputs": [],
   "source": [
    "process = BSS.Process.Gromacs(system, protocol)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "optimum-consultancy",
   "metadata": {},
   "source": [
    "When the process is instantiated, BioSimSpace takes the system that was read from AMBER format files and converts it to GROMACS format ready for simulation. Let's take a look at the list of input files that were autogenerated for us:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "recovered-first",
   "metadata": {},
   "outputs": [],
   "source": [
    "process.inputFiles()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "accepted-identification",
   "metadata": {},
   "source": [
    "Let's start the process running and, once again, wait for it to finish before getting the minimised system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dental-forum",
   "metadata": {},
   "outputs": [],
   "source": [
    "process.start()\n",
    "minimised = process.getSystem(block=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "coastal-driving",
   "metadata": {},
   "source": [
    "## Interactive molecular dynamics\n",
    "\n",
    "The example in the previous section was finished almost as soon as it began. Let's run a more complicated equilibration protocol so that we can learn more about how to monitor processes interactively using BioSimSpace."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "secondary-landscape",
   "metadata": {},
   "outputs": [],
   "source": [
    "protocol = BSS.Protocol.Equilibration(runtime=20*BSS.Units.Time.picosecond,\n",
    "                                      temperature_start=0*BSS.Units.Temperature.kelvin,\n",
    "                                      temperature_end=300*BSS.Units.Temperature.kelvin,\n",
    "                                      restraint=\"backbone\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "random-bishop",
   "metadata": {},
   "source": [
    "This protocol will equlibrate a system for 20 picoseconds, while heating it from 0 to 300 Kelvin and restraining any atoms in the backbone of the molecule. Note that some of the parameters passed have units, e.g. the temperatures are in Kelvin. BioSimSpace has a built in type system for handling variables with units. The `BSS.Units` package provides a convenient way of declaring these, for example `10*BSS.Units.Temperature.kelvin` creates an object of type `BSS.Types.Temperature` with a magnitude of 10 and unit of Kelvin. This allows the user to pass parameters with whatever unit they like. BioSimSpace will simply convert it to the correct unit for the chosen MD engine internally.\n",
    "\n",
    "One again, we now need a `Process` in order to run our simulation. Exectute the cell below to initialise an AMBER process and start it immediately. Note that we pass in the minimised system from the last example, along with our new protocol."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "entitled-olive",
   "metadata": {},
   "outputs": [],
   "source": [
    "process = BSS.Process.Amber(minimised, protocol).start()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "characteristic-anniversary",
   "metadata": {},
   "source": [
    "We can monitor the time, temperature, and energy as the process runs. If you run this multiple times using \"CTRL+Return\" you'll see the temperature slowly increasing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "rocky-twenty",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(process.getTime(), process.getTemperature(), process.getTotalEnergy())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "subtle-landscape",
   "metadata": {},
   "source": [
    "Since all of the values returned above are typed we can easily convert them to other units:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "surrounded-arabic",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(process.getTime().nanoseconds(), process.getTemperature().celsius(), process.getTotalEnergy().kj_per_mol())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "confused-valuation",
   "metadata": {},
   "source": [
    "It's possible to query many other thermodynamic records. What's available depends on type of protocol and the MD package that is used to run the protocol. To get more information, run:\n",
    "\n",
    "N.B. Certain functionality is specific to the process in question, i.e. `BSS.Process.Amber` will have different options to `BSS.Process.Gromacs`, but, for the purposes of interoperability, there is a core set of functionality that is consistent across all `Process` classes, e.g. all classes implement a `getSystem` method.)\n",
    "\n",
    "### Plotting time series data\n",
    "\n",
    "As well as querying the most recent records we can also get a time series of results by passing the `time_series` keyword argument to any of the data record getter methods, e.g.\n",
    "\n",
    "```python\n",
    "# Get a time series of pressure records.\n",
    "pressure = process.getPressure(time_series=True)\n",
    "```\n",
    "\n",
    "The `BSS.Notebook` package provides several useful tools that are available when working inside of a Jupyter notebook. One of these is the plot function, that allows us to create simple x/y plots of time-series data.\n",
    "\n",
    "Let's grab the same record data as above and use it to make some graphs of the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "pressed-physics",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate a plot of time vs temperature.\n",
    "plot1 = BSS.Notebook.plot(process.getTime(time_series=True), process.getTemperature(time_series=True))\n",
    "\n",
    "# Generate a plot of time vs energy.\n",
    "plot2 = BSS.Notebook.plot(process.getTime(time_series=True), process.getTotalEnergy(time_series=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "remarkable-freeze",
   "metadata": {},
   "source": [
    "(Note that, by default, the axis labels axis labels are automatically generated from the types and units of the x and y data that are passed to the function.)\n",
    "\n",
    "Re-run the cell using \"CTRL+Return\" to see the graphs update as the simulation progesses. (Occasionally, you might see a warning that the x and y data sets are mismatched in length, this is because the data was extracted before all records were written to disk.)\n",
    "\n",
    "Being able to query a process in real time is an incredibly useful tool. This could enable us to check for convergence, or spot errors in the simulation. If you ever need to kill a running process (perhaps it was configured incorrectly), run:\n",
    "\n",
    "```python\n",
    "process.kill()\n",
    "```\n",
    "\n",
    "### Visualising the molecular system\n",
    "\n",
    "Another useful tool that is available when working inside of a notebook is the `View` class that can be used to visualise the molecular system while a process is running. To create a `View` object we must attach it to a process (or a molecular system), e.g.:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bored-dimension",
   "metadata": {},
   "outputs": [],
   "source": [
    "view = BSS.Notebook.View(process)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "reduced-multimedia",
   "metadata": {},
   "source": [
    "We can now visualise the system:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "metallic-depression",
   "metadata": {},
   "outputs": [],
   "source": [
    "view.system()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "boxed-forestry",
   "metadata": {},
   "source": [
    "(If you see an empty view, try re-executing the cell.)\n",
    "\n",
    "To only view a specific molecule:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "restricted-major",
   "metadata": {},
   "outputs": [],
   "source": [
    "view.molecule(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "mediterranean-workplace",
   "metadata": {},
   "source": [
    "To view a list of molecules:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cellular-sussex",
   "metadata": {},
   "outputs": [],
   "source": [
    "view.molecules([0, 5, 10])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "quality-floor",
   "metadata": {},
   "source": [
    "If a particular view was of interest it can be reloaded as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "boolean-tulsa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reload the original view.\n",
    "view.reload(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "elementary-scene",
   "metadata": {},
   "source": [
    "To save a specific view as a PDB file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "viral-semiconductor",
   "metadata": {},
   "outputs": [],
   "source": [
    "view.savePDB(\"my_view.pdb\", index=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "welcome-dividend",
   "metadata": {},
   "source": [
    "### Reading and analysing trajectory data¶ \n",
    "\n",
    "The `BSS.Trajectory` package comes with a set of tools for reading and analysis trajectory files. Files can be loaded directly, or if supported, can be read from a running process.\n",
    "\n",
    "For example, to get the trajectory from the process, run:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "massive-christmas",
   "metadata": {},
   "outputs": [],
   "source": [
    "traj = process.getTrajectory()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "northern-governor",
   "metadata": {},
   "source": [
    "(If you get an error, then the trajectory file may be in the process of being written. Simply try again.)\n",
    "\n",
    "To get the current number of frames:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aging-affiliation",
   "metadata": {},
   "outputs": [],
   "source": [
    "traj.nFrames()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "judicial-victor",
   "metadata": {},
   "source": [
    "To get all of the frames as a list of `System` objects:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "featured-infrastructure",
   "metadata": {},
   "outputs": [],
   "source": [
    "frames = traj.getFrames()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "alternative-concert",
   "metadata": {},
   "source": [
    "Specific frames can be extracted by passing a list of indices, e.g. the first and last:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "liquid-harvey",
   "metadata": {},
   "outputs": [],
   "source": [
    "frames = traj.getFrames([0, -1])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "conditional-gross",
   "metadata": {},
   "source": [
    "Like most things in BioSimSpace, the `Trajectory` class is simply a wrapper around existing tools. Internally, trajectories are stored as an [MDTraj](http://mdtraj.org) object. This can be obtained, allowing the user direct access to the full power of MDTraj:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "suited-donor",
   "metadata": {},
   "outputs": [],
   "source": [
    "mdtraj = traj.getTrajectory()\n",
    "type(mdtraj)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "nervous-white",
   "metadata": {},
   "source": [
    "Alternatively, a trajectory can be returned in [MDAnalysis](https://www.mdanalysis.org) format:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "assisted-internship",
   "metadata": {},
   "outputs": [],
   "source": [
    "mdanalysis = traj.getTrajectory(format=\"mdanalysis\")\n",
    "type(mdanalysis)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "talented-miniature",
   "metadata": {},
   "source": [
    "The `Trajectory` class also provides wrappers around some basic MDTraj analysis tools, allowing the user to compute quantities such as the root mean squared displacement (RMSD).\n",
    "\n",
    "Let's measure the RMSD of the alanine-dipeptide molecule with a reference to its configuration in the first trajectory frame. To extract the alanine-dipeptide, we search the system for a residue named ALA. We'll also plot the RMSD for each frame of the trajectory.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "saving-laser",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Search the system for a residue named ALA. Since there is a single match,\n",
    "# we take the first result.\n",
    "molecule = system.search(\"mol with resname ALA\")[0]\n",
    "\n",
    "# Get the indices of the atoms in the molecule, relative to the original system.\n",
    "indices = [system.getIndex(x) for x in molecule.getAtoms()]\n",
    "\n",
    "# Compute the RMSD for each frame and plot the result.\n",
    "BSS.Notebook.plot(y=process.getTrajectory().rmsd(frame=0, atoms=indices),\n",
    "                  xlabel=\"Frame\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}