notebooks/Lecture1_MDAnalysisBasics_TheUniverse.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# <center>Day 1, Session 1: The Universe</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center>MDAnalysis & MOLSSI Workshop - June 2024 - ASU</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Getting started with MDAnalysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Installing MDAnalysis and related dependencies\n",
    "\n",
    "For this workshop, we have provided specific [install instructions](https://github.com/MDAnalysis/MDAnalysisMolSSIWorkshop-Intermediate2Day/blob/jun24-ws/INSTALL.md) for you to follow.\n",
    "\n",
    "These will lead you through setting up a installation of MDAnalysis v2.7.0, and tooling necessary for this workshop including; the MDAKits cookiecutter, sphinx documentation, and jupyter.\n",
    "\n",
    "Should you wish to install a clean version of MDAnalysis, please look at our [userguide instructions](https://userguide.mdanalysis.org/stable/installation.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_style": "split",
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### The general structure of MDAnalysis\n",
    "\n",
    "The two fundamental classes of MDAnalysis are the `Universe` and the `AtomGroup`.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_style": "center"
   },
   "source": [
    "\n",
    "<center><img src=\"imgs/mdaclasses.png\" alt=\"mda\" style=\"width: 1000px;\"/></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**The `Universe`** contains everything about a molecular dynamics system\n",
    " - Static information: atoms and their connectivities (topology)\n",
    " - Dynamic information: the trajectory (coordinates, frames)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<center><img src=\"imgs/mdaclasses.png\" alt=\"mda\" style=\"width: 1000px;\"/></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "- The atoms in a `Universe` can be accessed through a hierarchy of containers:\n",
    " - *Atoms* can be grouped together into **an `AtomGroup`** \n",
    "    - *Residues* are made up of *atoms*. They can be grouped into `ResidueGroups`\n",
    "      - *Segments* are made up of *residues*. They can be grouped into `SegmentGroups`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<center><img src=\"imgs/mdaclasses.png\" alt=\"mda\" style=\"width: 1000px;\"/></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_style": "center",
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### A (very) basic workflow for an analysis in MDAnalysis:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "1. import MDAnalysis\n",
    "2. create a `Universe`\n",
    "3. define an `AtomGroup`\n",
    "4. collect position data\n",
    "5. analyse!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# The `Universe`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**The basic command for loading a universe is:**\n",
    "\n",
    " `u = mda.Universe(topology, trajectory)`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_style": "split",
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "- The *topology* file must contain the atom information \n",
    "- The (optional) *trajectory* file(s) contains the positions of atoms with time - more on this next Session. \n",
    "\n",
    "Note that some files can double as both a *topology* and a *trajectory*.  \n",
    "\n",
    "MDanalysis supports [over 40 input file types](https://userguide.mdanalysis.org/stable/formats/index.html#formats)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "cell_style": "split",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<Universe with 3341 atoms>\n"
     ]
    }
   ],
   "source": [
    "# Avoid any warning messages\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\") \n",
    "\n",
    "# First we import MDAnalysis\n",
    "import MDAnalysis as mda\n",
    "\n",
    "# Let's get some example data\n",
    "from MDAnalysis.tests.datafiles import PSF, DCD\n",
    "\n",
    "# and now load our universe!\n",
    "u = mda.Universe(PSF, DCD)\n",
    "print(u)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Key properties of a `Universe`:**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_style": "split",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "- `atoms`: an `AtomGroup` containing all of the system's atoms\n",
    "    - similarly, `segments` and `residues`; a `SegmentGroup` and a `ResidueGroup`, respectively\n",
    "    \n",
    "- Various bond and angle information, as `TopologyGroups`: `bonds`, `angles`, `dihedrals`, `impropers` (if found in the topology file)\n",
    "\n",
    "- `trajectory`: (section 2): accessing time-dependent data structures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "cell_style": "split",
    "scrolled": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<TopologyGroup containing 3365 bonds>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u.bonds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# AtomGroups"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**An `AtomGroup` is an \"array\" of atoms.**\n",
    "\n",
    "We can get various properties of each atom contained in an `AtomGroup` through attributes, e.g.:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_style": "split",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "   - `names`\n",
    "   - `resnames`\n",
    "   - `resids` \n",
    "   - `charges`\n",
    "   - `masses`\n",
    "\n",
    "Exactly which properties you can get depend on what is read from the topology (see the [documentation](https://userguide.mdanalysis.org/stable/formats/index.html#formats))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "cell_style": "split"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'MDAnalysis.core.groups.AtomGroup'>\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array(['N', 'HT1', 'HT2', ..., 'C', 'OT1', 'OT2'], dtype=object)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ag = u.atoms\n",
    "print(type(ag))\n",
    "ag.names"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_style": "center",
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**`ResidueGroup`s and `SegmentGroup`s work similarly.**\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_style": "split",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "These can be accessed through the `atoms`, `residues` and `segments` properties."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "cell_style": "split",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'MDAnalysis.core.groups.ResidueGroup'>\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([  1,   2,   3,   4,   5,   6,   7,   8,   9,  10,  11,  12,  13,\n",
       "        14,  15,  16,  17,  18,  19,  20,  21,  22,  23,  24,  25,  26,\n",
       "        27,  28,  29,  30,  31,  32,  33,  34,  35,  36,  37,  38,  39,\n",
       "        40,  41,  42,  43,  44,  45,  46,  47,  48,  49,  50,  51,  52,\n",
       "        53,  54,  55,  56,  57,  58,  59,  60,  61,  62,  63,  64,  65,\n",
       "        66,  67,  68,  69,  70,  71,  72,  73,  74,  75,  76,  77,  78,\n",
       "        79,  80,  81,  82,  83,  84,  85,  86,  87,  88,  89,  90,  91,\n",
       "        92,  93,  94,  95,  96,  97,  98,  99, 100, 101, 102, 103, 104,\n",
       "       105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117,\n",
       "       118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130,\n",
       "       131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143,\n",
       "       144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156,\n",
       "       157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169,\n",
       "       170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182,\n",
       "       183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195,\n",
       "       196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208,\n",
       "       209, 210, 211, 212, 213, 214])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rg = u.residues\n",
    "print(type(rg))\n",
    "rg.resids"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Atom selections"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_style": "center",
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We don't ususally want to work with the whole set of atoms in a trajectory. We need a way to create `AtomGroups` containing selected atoms..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### But first... visualising parts of your system!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_style": "split",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "[nglview](https://github.com/nglviewer/nglview#usage) will allow us to view MDAnalysis Universes and AtomGroups inside Jupyter notebooks.\n",
    "\n",
    "We've intalled if for you on the workshop Jupyter hub, but you can also install it on your own computer (*Beware: it's not always straightforward getting it to actually work. Try making sure to install nglview and launch the notebook from the same conda environment*)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "cell_style": "split",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "772ba8d8f6d64b5591dd3d858c9e78cc",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "NGLWidget(max_frame=97)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# first, import nglview\n",
    "import nglview as nv\n",
    "  \n",
    "# add a universe (or atomgroup)\n",
    "view_u = nv.show_mdanalysis(u)\n",
    "\n",
    "# launch the viewer\n",
    "view_u"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_style": "center",
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Selecting atoms to create AtomGroups"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_style": "split"
   },
   "source": [
    "### Atom indexing\n",
    "\n",
    "`AtomGroup`s can be created by indexing `atoms`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "cell_style": "split"
   },
   "outputs": [],
   "source": [
    "ag = u.atoms[0:2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "cell_style": "split"
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8035168e9e9547c5b23dea0a260e8605",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "NGLWidget(max_frame=97)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "view = nv.show_mdanalysis(ag)\n",
    "view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Selection strings and `select_atoms`\n",
    "\n",
    "We can use the `select_atoms()` method of an `AtomGroup` or `Universe` to return an `AtomGroup` based on a selection string."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_style": "split",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "There's a lot of options for selection strings (see the [UserGuide]( https://userguide.mdanalysis.org/stable/selections.html)); including:\n",
    "\n",
    " - selection by attribute (e.g. residue name (`resname`)), including presets like `protein`\n",
    " - wildcard matching (`*`)\n",
    " - boolean operators (`and`, `or`, `not`)\n",
    " - geometric (e.g. `around`, `sphzone`, ...)\n",
    " - and more!\n",
    " \n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "cell_style": "split",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "dc35ced5de85477287f075342a72d0cf",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "NGLWidget(max_frame=97)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ag = u.select_atoms('protein')\n",
    "view_ag = nv.show_mdanalysis(ag)\n",
    "view_ag.add_licorice()\n",
    "view_ag"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_style": "center",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Working with coordinates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### The most useful attribute of out atoms are their coordinates, available in the `positions` attribute of an `AtomGroup`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_style": "split",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "The positions are returned as a numpy array, which we can then readily manipulate.\n",
    "\n",
    "There are some built-in functions based on position data, e.g. `center_of_mass()`, `center_of_geometry()`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "cell_style": "split",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 11.736044    8.500797  -10.445281 ]\n",
      " [ 12.365119    7.839936  -10.834842 ]\n",
      " [ 12.0919485   9.441535  -10.724611 ]\n",
      " ...\n",
      " [  6.512604   18.447018   -7.134053 ]\n",
      " [  6.300186   19.363485   -7.935916 ]\n",
      " [  5.5854015  17.589624   -6.9656615]]\n"
     ]
    }
   ],
   "source": [
    "pos = u.atoms.positions\n",
    "print(pos)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "This is just data from one frame - in the next session, you'll learn how to work with trajectories to get data across a whole simulation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# A note on units"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_style": "split"
   },
   "source": [
    "The units of MDAnalysis trajectories are the Å (ångström) for length and ps (picosecond) for time. Regardless of how the original MD format stored the trajectory data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "cell_style": "split"
   },
   "source": [
    "<center><img src=\"imgs/mdanalysis_units.png\" alt=\"mda\" style=\"width: 1000px;\"/></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## A summary of Lecture 1\n",
    "\n",
    "- A `Universe` contains all information about a simulation system\n",
    "\n",
    "- An `AtomGroup` contains information about a group of atoms\n",
    "\n",
    "- We can use `Universe.select_atoms()` to create an `AtomGroup` containing specific atoms from a `Universe`\n",
    "\n",
    "- Positions of atoms in an AtomGroup are accessed through `AtomGroup.positions`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Now - on to the first tutorial!\n",
    "\n",
    "Find the tutorial notebook at https://github.com/MDAnalysis/MDAnalysisMolSSIWorkshop-Intermediate2Day/tree/jun24-ws/notebooks - look for `Tutorial1_System_Manipulation.ipynb`\n",
    "\n",
    "**Remember:**\n",
    "- Go at your own pace!\n",
    "- Ask questions!\n",
    "- Be respectful!\n",
    "- Take breaks!"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}