diff --git a/.gitignore b/.gitignore
index 39421ce87..9e4be271b 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,6 +1,5 @@
*~
.vagrant
-notebooks/.ipynb_checkpoints
notebooks/*.png
notebooks/*.svg
notebooks/*.pdf
@@ -8,4 +7,7 @@ notebooks/*.pdb
notebooks/*.ndx
tmp
doc/build
-doc/source/.ipynb_checkpoints
+**/.ipynb_checkpoints
+**/.vscode
+.vscode
+**/.DS_Store
\ No newline at end of file
diff --git a/doc/examples/README.rst b/doc/examples/README.rst
new file mode 100644
index 000000000..2fc2c67eb
--- /dev/null
+++ b/doc/examples/README.rst
@@ -0,0 +1,6 @@
+
+========
+Examples
+========
+
+Here are some examples.
\ No newline at end of file
diff --git a/doc/examples/analysis/pca.ipynb b/doc/examples/analysis/pca.ipynb
new file mode 100644
index 000000000..8ef60060d
--- /dev/null
+++ b/doc/examples/analysis/pca.ipynb
@@ -0,0 +1,344 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Principal Component Analysis in MDAnalysis\n",
+ "\n",
+ "2019\n",
+ "\n",
+ "Author: [Lily Wang](http://minium.com.au) [(@lilyminium)](https://github.com/lilyminium)\n",
+ "\n",
+ "Inspired by the MDAnalysis PCA tutorial by [Kathleen Clark](https://becksteinlab.physics.asu.edu/people/75/kathleen-clark) [(@kaceyreidy)](https://github.com/kaceyreidy)\n",
+ "\n",
+ "In this tutorial we:\n",
+ "\n",
+ "* use PCA to analyse and visualise large macromolecular conformational changes in the enzyme adenylate kinase (AdK)\n",
+ "* use PCA to compare the conformational differences between ???"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Background"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Principal component analysis (PCA) is a statistical technique that decomposes a system of observations into linearly uncorrelated variables called **principal components**. These components are ordered so that the first principal component accounts for the largest variance in the data, and each following component accounts for less and less variance. PCA is often applied to molecular dynamics trajectories to extract the large-scale conformational motions or \"essential dynamics\" of a protein. The frame-by-frame conformational fluctuation can be considered a linear combination of the essential dynamics yielded by the PCA.\n",
+ "\n",
+ "In MDAnalysis, the method is as follows:\n",
+ "\n",
+ "1. Optionally align each frame in your trajectory to the first frame.\n",
+ "2. Construct a 3N x 3N covariance for the N atoms in your trajectory. Optionally, you can provide a mean; otherwise the covariance is to the averaged structure over the trajectory.\n",
+ "3. Diagonalise the covariance matrix. The eigenvectors are the principal components, and their eigenvalues are the associated variance.\n",
+ "4. Sort the eigenvalues so that the principal components are ordered by variance.\n",
+ "\n",
+ "
\n",
+ " \n",
+ "**Note**\n",
+ " \n",
+ "It should be noted that principal component analysis algorithms are deterministic, but the solutions are not unique. For example, you could easily change the sign of an eigenvector without altering the PCA. Different algorithms are likely to produce different answers, due to variations in implementation. `MDAnalysis` is likely to return different solutions to, say, `cpptraj`. \n",
+ "
"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Large conformational changes in adenylate kinase"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "In MDAnalysis, analysis modules usually need to be imported explicitly. The `pca` module contains the `PCA` class that we will use for analysis. We also import the AdK files from the MDAnalysis test suite."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "application/vnd.jupyter.widget-view+json": {
+ "model_id": "3996de84aab74b3aa903e977a8ac80b7",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ "_ColormakerRegistry()"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "import MDAnalysis as mda\n",
+ "import MDAnalysis.analysis.pca as pca\n",
+ "from MDAnalysisTests.datafiles import PSF, DCD\n",
+ "\n",
+ "import nglview as nv"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "As usual, we start off by creating a universe."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "u = mda.Universe(PSF, DCD)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Unlike other analyses, `pc.PCA` can only be applied to `Universe`s. The default `PCA` arguments are:\n",
+ "\n",
+ "```python\n",
+ "my_pca = pca.PCA(u, select='all', align=False, mean=None, n_components=None)\n",
+ "```\n",
+ "\n",
+ "By default (`align=False`), your trajectory will not be aligned to any structure. If you set `align=True`, every frame will be aligned to the first frame of your trajectory, based on the atoms in your `select` string. \n",
+ "\n",
+ "As PCA is usually used to extract large-scale conformational motions, we select only the backbone atoms here."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "pc = pca.PCA(u, select=\"backbone\", align=True)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Once you set up the class, you can run the analysis with `.run(start=None, stop=None, step=None, verbose=None)`. These allow you to specify the frames to compute the analysis over. The default arguments compute over every frame."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ ""
+ ]
+ },
+ "execution_count": 4,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "pc.run(verbose=True)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The principal components are accessible in `.p_components`."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "(2565, 2565)\n"
+ ]
+ },
+ {
+ "data": {
+ "text/plain": [
+ "array([ 0.02725098, 0.00156086, 0.00816821, ..., -0.01783826,\n",
+ " 0.04746114, 0.04257271])"
+ ]
+ },
+ "execution_count": 5,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "print(pc.p_components.shape)\n",
+ "pc.p_components[0]"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The variance of each principal component is in `.variance`. For example, to get the variance explained by the first principal component:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "281443.5086197605"
+ ]
+ },
+ "execution_count": 6,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "pc.variance[0]"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "This variance is somewhat meaningless by itself. It is much more intuitive to consider the variance of a principal component as a percentage of the total variance in the data. MDAnalysis also tracks the percentage cumulative variance in `.cumulated_variance`. As shown below, the first principal component contains 98.7% the total trajectory variance. The first three components combined account for 99.9% of the total variance."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "0.9873464381554058\n",
+ "0.999419901112709\n"
+ ]
+ }
+ ],
+ "source": [
+ "print(pc.cumulated_variance[0])\n",
+ "print(pc.cumulated_variance[3])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The average structure is also saved as an `AtomGroup` in `.mean_atoms`."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 15,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[[16.297781 6.8397956 -7.622989 ]\n",
+ " [14.900139 7.062459 -7.235277 ]\n",
+ " [14.185768 5.8268375 -6.879689 ]\n",
+ " ...\n",
+ " [13.035071 15.354209 -3.8042812]\n",
+ " [13.695147 15.725297 -4.988666 ]\n",
+ " [12.63667 15.566869 -6.1185045]]\n"
+ ]
+ }
+ ],
+ "source": [
+ "print(pc.mean_atoms.positions)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 14,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "application/vnd.jupyter.widget-view+json": {
+ "model_id": "a33209712ed34b32b969c1f8258aed91",
+ "version_major": 2,
+ "version_minor": 0
+ },
+ "text/plain": [
+ "NGLWidget()"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "mean_structure = mda.Merge(pc.mean_atoms)\n",
+ "nv.show_mdanalysis(mean_structure)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python (mdanalysis)",
+ "language": "python",
+ "name": "mdanalysis"
+ },
+ "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.3"
+ },
+ "toc": {
+ "base_numbering": 1,
+ "nav_menu": {},
+ "number_sections": false,
+ "sideBar": true,
+ "skip_h1_title": false,
+ "title_cell": "Table of Contents",
+ "title_sidebar": "Contents",
+ "toc_cell": false,
+ "toc_position": {},
+ "toc_section_display": true,
+ "toc_window_display": false
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/doc/source/_static/.DS_Store b/doc/source/_static/.DS_Store
new file mode 100644
index 000000000..a1d2192ce
Binary files /dev/null and b/doc/source/_static/.DS_Store differ
diff --git a/doc/source/_static/custom.css b/doc/source/_static/custom.css
index cd04a426b..5fbb18515 100644
--- a/doc/source/_static/custom.css
+++ b/doc/source/_static/custom.css
@@ -5,8 +5,7 @@
/* MDAnalysis gray: #808080 */
/* MDAnalysis white: #FFFFFF */
/* MDAnalysis black: #000000 */
-
-@import url('readable.css');
+/* Darker orange: e76900 */
/* -- page layout --------------------------------------------------------- */
@@ -19,27 +18,57 @@ div.body {
color: #000000;
}
+div.sphinxsidebar a {
+ color: #808080;
+}
+
div.sphinxsidebar a:hover {
color: #FF9200;
}
+div.sphinxsidebar a:visited {
+ color: #808080;
+}
+
div.sphinxsidebar p {
color: #808080;
}
+div.wy-menu-vertical span.caption-text {
+ color: #FF9200;
+}
+
+/* li.current ul {
+ background: #fbd29b;
+} */
+
+
+/* div.sphinxsidebar a {
+ color: #808080;
+}
+
+div.sphinxsidebar a:hover {
+ color: #FF9200;
+}
+
+div.wy-menu-vertical a:visited {
+ color: #808080;
+}
+ */
+
/* -- body styles --------------------------------------------------------- */
-a {
+div.rst-content a {
color: #FF9200;
text-decoration: none;
}
-a:visited {
- color: #FF9200;
+div.rst-content a:visited {
+ color: #e76900;
}
-a:hover {
+div.rst-content a:hover {
color: #FF9200;
text-decoration: underline;
}
@@ -49,6 +78,10 @@ pre, tt, code {
font-family: Menlo, Monaco, 'Courier New', monospace
}
+/* For the backtick tags */
+/* code > span.pre {
+ color: #e76900;
+} */
div.body h1 {
font-weight: bolder;
@@ -75,4 +108,36 @@ div.admonition {
div.admonition p.admonition-title {
font-size: 100%;
font-weight: bolder;
+}
+
+div.admonition.warning {
+ background: #ffdeb3;
+}
+
+div.admonition.warning > p.admonition-title {
+ background: #FF9200;
+}
+
+div.admonition.note {
+ background: #e6e6e6; /* RTD menu light grey */
+}
+
+div.admonition.note > p.admonition-title {
+ background: #808080;
+}
+
+/* -- notebook styles --------------------------------------------------------- */
+.rst-content div[class^="highlight"] pre {
+ white-space: pre-wrap;
+}
+
+/* -- gallery styles ---------------------------------------------------------- */
+div.sphx-glr-download a {
+ border: 1px solid #fff;
+ background: #808080;
+}
+
+div.sphx-glr-download a:hover {
+ border: 1px solid #fff;
+ background: #000;
}
\ No newline at end of file
diff --git a/doc/source/_static/logos/user_guide.png b/doc/source/_static/logos/user_guide.png
new file mode 100644
index 000000000..5d82c3cfd
Binary files /dev/null and b/doc/source/_static/logos/user_guide.png differ
diff --git a/doc/source/_static/logos/user_guide.psd b/doc/source/_static/logos/user_guide.psd
new file mode 100644
index 000000000..428781843
Binary files /dev/null and b/doc/source/_static/logos/user_guide.psd differ
diff --git a/doc/source/auto_examples/index.rst b/doc/source/auto_examples/index.rst
new file mode 100644
index 000000000..74cd10b92
--- /dev/null
+++ b/doc/source/auto_examples/index.rst
@@ -0,0 +1,24 @@
+:orphan:
+
+
+
+.. _sphx_glr_auto_examples:
+
+
+========
+Examples
+========
+
+Here are some examples.
+
+.. raw:: html
+
+
+
+
+
+.. only:: html
+
+ .. rst-class:: sphx-glr-signature
+
+ `Gallery generated by Sphinx-Gallery `_
diff --git a/doc/source/conf.py b/doc/source/conf.py
index fe119d756..54d3689d5 100644
--- a/doc/source/conf.py
+++ b/doc/source/conf.py
@@ -13,7 +13,8 @@
# import os
# import sys
# sys.path.insert(0, os.path.abspath('.'))
-
+import sphinx_rtd_theme
+import sphinx_gallery
# -- Project information -----------------------------------------------------
@@ -35,7 +36,9 @@
'sphinx.ext.viewcode',
'sphinx.ext.githubpages',
'sphinx_sitemap',
- 'nbsphinx'
+ 'nbsphinx',
+ 'sphinx_rtd_theme',
+ 'sphinx_gallery.gen_gallery'
]
mathjax_path = 'https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.0/MathJax.js?config=TeX-AMS-MML_HTMLorMML'
@@ -46,15 +49,33 @@
# List of patterns, relative to source directory, that match files and
# directories to ignore when looking for source files.
# This pattern also affects html_static_path and html_extra_path.
-exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store']
+exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store', '.ipynb_checkpoints']
+
+sphinx_gallery_conf = {
+ 'example_dirs': ['../examples'], # source files
+ 'gallery_dirs': ['auto_examples'], # where files are built
+ 'filename_pattern': '/plot_', # filenames to execute, accepts regex
+ 'ignore_pattern': r'__init__\.py', # filenames to not execute/parse/add, accepts regex
+ 'download_all_examples': False,
+
+ # 'binder': {
+ # 'org': 'MDAnalysis',
+ # 'repo': 'UserGuide',
+ # 'branch': 'master'
+ # 'binderhub_url': 'full URL of binderhub deployment',
+ # 'dependencies': './binder/requirements.txt',
+ # 'use_jupyter_lab': False
+ # }
+}
# -- Options for HTML output -------------------------------------------------
# The theme to use for HTML and HTML Help pages. See the documentation for
# a list of builtin themes.
#
-html_theme = 'alabaster'
+html_theme = 'sphinx_rtd_theme'
+# html_theme_path = ['_themes', ]
# styles/fonts to match http://mdanalysis.org (see public/css)
#
@@ -78,41 +99,32 @@
extra_nav_links['@mdanalysis'] = 'https://twitter.com/mdanalysis'
html_theme_options = {
- 'logo' : "logos/mdanalysistutorial-logo-200x150.png",
- 'github_user': 'MDAnalysis',
- 'github_repo': 'mdanalysis',
- 'github_button': False,
- # 'github_type': 'star',
- 'github_banner': True,
- 'show_related': True,
- 'fixed_sidebar': False,
- 'sidebar_includehidden': True,
- 'sidebar_collapse': True,
- # style
- 'link': color['orange'],
- 'link_hover': color['orange'],
- 'gray_1': color['gray'],
- 'narrow_sidebar_bg': color['gray'],
- 'narrow_sidebar_fg': color['white'],
- # typography
- 'font_family': "'PT Sans', Helvetica, Arial, 'sans-serif'",
- 'head_font_family': "",
- 'code_font_family': "Menlo, Monaco, 'Courier New', monospace",
- 'caption_font_size': "smaller",
- # external links
- 'extra_nav_links': extra_nav_links,
+ 'canonical_url': '',
+ 'logo_only': True,
+ 'display_version': True,
+ 'prev_next_buttons_location': 'bottom',
+ 'style_external_links': False,
+ 'style_nav_header_background': 'white', #'#e76900', # dark orange
+ # Toc options
+ 'collapse_navigation': True,
+ 'sticky_navigation': True,
+ 'navigation_depth': 5,
+ 'includehidden': True,
+ 'titles_only': False,
}
# The name of an image file (within the static path) to use as favicon of the
# docs. This file should be a Windows icon file (.ico) being 16x16 or 32x32
# pixels large.
html_favicon = "_static/logos/mdanalysis-logo.ico"
+html_logo = '_static/logos/user_guide.png'
# Add any paths that contain custom static files (such as style sheets) here,
# relative to this directory. They are copied after the builtin static files,
# so a file named "default.css" will overwrite the builtin "default.css".
html_static_path = ['_static']
+html_css_files = ['custom.css']#, 'readable.css']
# Custom sidebar templates, maps document names to template names.
# alabaster sidebars
diff --git a/doc/source/data_structures.rst b/doc/source/data_structures.rst
index cbd96bba8..e89f41fb5 100644
--- a/doc/source/data_structures.rst
+++ b/doc/source/data_structures.rst
@@ -52,6 +52,23 @@ AtomGroup
AtomGroup instances can be easily created
(e.g., from an AtomGroup.select_atoms() selection or simply by slicing).
+--------------------
+Indexing and slicing
+--------------------
+blah blah
+
+-----------
+Set methods
+-----------
+blah
+
+----------------
+Geometry methods
+----------------
+
+
+UpdatingAtomGroup
+====================
ResidueGroup and SegmentGroup
@@ -63,10 +80,3 @@ Residue or a Segment (typically a whole molecule or all of the solvent)
also exist as containers, as well as groups of these units (
ResidueGroup, SegmentGroup).
-
-
-
-Topology
-====================
-
-MDA Topology =/= 'force field' topology. Connectivity only. See ParmEd etc for other stuff.
\ No newline at end of file
diff --git a/doc/source/faq.rst b/doc/source/faq.rst
new file mode 100644
index 000000000..8ac35ab63
--- /dev/null
+++ b/doc/source/faq.rst
@@ -0,0 +1,30 @@
+.. -*- coding: utf-8 -*-
+
+==========================
+Frequently asked questions
+==========================
+
+
+Trajectories
+==========================
+
+Why do the atom positions change over trajectories?
+---------------------------------------------------
+
+A fundamental concept in MDAnalysis is that at any one time,
+only one time frame of the trajectory is being accessed. The
+:code:`trajectory` attribute of a :code:`Universe` is actually (usually) a file reader.
+Think of the trajectory as a function :math:`X(t)` of the frame index :math:`t`
+that makes the data from this specific frame available. This structure is important
+because it allows MDAnalysis to work with trajectory files too large to fit
+into the computer's memory.
+
+The Timestep data structure :code:`my_universe.trajectory.ts` holds frame information.
+You also see the Timestep when you iterate through a trajectory:
+
+.. code-block:: python
+
+ for ts in my_universe.trajectory:
+ print(ts.frame, ts.time)
+
+This prints the frame index (starting with 0) and the time for all frames.
\ No newline at end of file
diff --git a/doc/source/index.rst b/doc/source/index.rst
index 86cd6dcae..2bd271c35 100644
--- a/doc/source/index.rst
+++ b/doc/source/index.rst
@@ -12,9 +12,16 @@ Link to gallery.
.. toctree::
:maxdepth: 1
- :caption: Contents:
+ :caption: Getting started
installation
quickstart
+ faq
+ auto_examples/index
+
+.. toctree::
+ :maxdepth: 1
+ :caption: User guide
+
data_structures
formats/formats
diff --git a/doc/source/quickstart.ipynb b/doc/source/quickstart.ipynb
index fa3447b39..169f92189 100644
--- a/doc/source/quickstart.ipynb
+++ b/doc/source/quickstart.ipynb
@@ -4,27 +4,1616 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "# Quickstart"
+ "# Quick start guide"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "This guide is designed as a basic introduction to MDAnalysis to get you up and running. You can see more complex tasks in our [Example gallery](placeholder). This page outlines how to:\n",
+ "\n",
+ "* [load a molecular dynamics structure or trajectory](#Loading-a-structure-or-trajectory)\n",
+ "* [work with AtomGroups, a central data structure in MDAnalysis](#Working-with-groups-of-atoms)\n",
+ "* [work with a trajectory](#Working-with-trajectories)\n",
+ "* [write out coordinates](#Writing-out-coordinates)\n",
+ "\n",
+ "* [use the analysis algorithms in MDAnalysis](#Analysis)\n",
+ "* [correct and automated citation of MDAnalysis and algorithms](#References)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "\n",
+ "0.19.2\n"
+ ]
+ }
+ ],
+ "source": [
+ "import MDAnalysis as mda\n",
+ "from MDAnalysis.tests.datafiles import PSF, DCD, GRO, XTC\n",
+ "\n",
+ "print(mda.Universe(PSF, DCD))\n",
+ "print(mda.__version__)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "This tutorial assumes that you already have MDAnalysis installed. Running the cell above should give:\n",
+ "\n",
+ "```console\n",
+ "\n",
+ "0.19.2\n",
+ "```\n",
+ "If you get an error message, you need to install MDAnalysis. If your version is under `0.18.0`, you need to upgrade MDAnalysis. [Instructions for both are here.](installation.rst) After installing, restart this notebook."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Overview"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "MDAnalysis is a Python package that provides tools to access and analyse data in molecular dynamics trajectories. Several key data structures form the backbone of MDAnalysis.\n",
+ "\n",
+ "* A molecular system consists of particles. A particle is represented as an `Atom` object, even if it is a coarse-grained bead. \n",
+ "* `Atom`s are grouped into [AtomGroup](data_structures.rst#atomgroups)s. The `AtomGroup` is probably the most important class in MDAnalysis, as almost everything can be accessed through it. See [Working with atoms](#Working-with-groups-of-atoms) below. \n",
+ "* A `Universe` contains all the particles in a molecular system in an `AtomGroup` accessible at the `.atoms` attribute, and combines it with a trajectory at `.trajectory`. \n",
+ "\n",
+ "A fundamental concept in MDAnalysis is that at any one time, \n",
+ "only one time frame of the trajectory is being accessed. The \n",
+ "`trajectory` attribute of a `Universe` is usually a file reader. \n",
+ "Think of the trajectory as a function $X(t)$ of the frame index $t$ \n",
+ "that only makes the data from this specific frame available. This structure is important \n",
+ "because it allows MDAnalysis to work with trajectory files too large to fit \n",
+ "into the computer's memory. \n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Loading a structure or trajectory"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Working with MDAnalysis typically starts with loading data into a `Universe`, the central data structure in MDAnalysis. \n",
+ "\n",
+ "
\n",
+ "\n",
+ "**See also**\n",
+ "\n",
+ "The [user guide](data_structures.rst#universes) has a complete explanation of ways to create and manipulate a `Universe`.\n",
+ "
\n",
+ "\n",
+ "\n",
+ "The first arguments for creating a `Universe` are topology and trajectory files.\n",
+ "\n",
+ "* A **topology file** is always required for loading data into a Universe. A topology file lists atoms, residues, and their connectivity. MDAnalysis accepts the PSF, PDB, CRD, and GRO formats.\n",
+ "\n",
+ "* A topology file can then be followed by **any number of trajectory files**. A trajectory file contains a list of coordinates in the order defined in the topology. If no trajectory files are given, then only a structure is loaded. If multiple trajectory files are given, the trajectories are concatenated in the given order. MDAnalysis accepts single frames (e.g. PDB, CRD, GRO) and timeseries data (e.g. DCD, XTC, TRR, XYZ). "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {
+ "scrolled": true
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "\n",
+ "False\n"
+ ]
+ }
+ ],
+ "source": [
+ "psf = mda.Universe(PSF)\n",
+ "print(psf)\n",
+ "print(hasattr(psf, 'trajectory'))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "As PSF files don't contain any coordinate information and no trajectory file has been loaded, the `structure_only` universe does not contain a trajectory. If the topology file does contain coordinate information, a trajectory of 1 frame is created."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "\n",
+ "1\n"
+ ]
+ }
+ ],
+ "source": [
+ "gro = mda.Universe(GRO)\n",
+ "print(gro)\n",
+ "print(len(gro.trajectory))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "For the remainder of this guide we will work with the universe `u`, created below. This is a simulation where the enzyme adenylate kinase samples a transition from a closed to an open conformation [(Beckstein *et al.*, 2009)](https://doi.org/10.1016/j.jmb.2009.09.009)."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "\n",
+ "98\n"
+ ]
+ }
+ ],
+ "source": [
+ "u = mda.Universe(PSF, DCD)\n",
+ "print(u)\n",
+ "print(len(u.trajectory))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "
\n",
+ "\n",
+ "**Note**\n",
+ "\n",
+ "The MDAnalysis test suite is packaged with a bunch of test files and trajectories, which are named after their file format. We are using these files throughout this guide for convenience. To analyse your own files, simply replace the `PSF` and `DCD` above with paths to your own files. For example:\n",
+ "\n",
+ "```python\n",
+ "structure_only = mda.Universe(\"my_pdb_file.pdb\")\n",
+ "```\n",
+ "
"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Working with groups of atoms"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Most analysis requires creating and working with an `AtomGroup`, a collection of `Atom`s. For convenience, you can also work with chemically meaningful groups of `Atoms` such as a `Residue` or a `Segment`. These come with analogous containers to `AtomGroup`: `ResidueGroup` and `SegmentGroup`. For instance, the `.residues` attribute of a `Universe` returns a `ResidueGroup`."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ ", , , ..., , , ]>\n"
+ ]
+ }
+ ],
+ "source": [
+ "print(u.residues)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Selecting atoms"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The easiest way to access the particles of your `Universe` is with the `atoms` attribute:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ ""
+ ]
+ },
+ "execution_count": 6,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "u.atoms"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "This returns an `AtomGroup`, which can be thought of as a list of `Atom` objects. Most analysis involves working with groups of atoms in `AtomGroup`s. `AtomGroup`s can easily be created by slicing another `AtomGroup`. For example, the below slice returns the last five atoms."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "metadata": {
+ "scrolled": true
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ ", , , , ]>\n"
+ ]
+ }
+ ],
+ "source": [
+ "last_five = u.atoms[-5:]\n",
+ "print(last_five)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "MDAnalysis supports fancy indexing: passing an array or list of indices to get a new AtomGroup with the atoms at those indices in the old AtomGroup."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ ", , , , , ]>\n"
+ ]
+ }
+ ],
+ "source": [
+ "print(last_five[[0, 3, -1, 1, 3, 0]])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "MDAnalysis has also implemented a [powerful atom selection language](selections.rst) that is similar to existing languages in VMD, PyMol, and other packages. This is available with the `.select_atoms()` function of an `AtomGroup` or `Universe` instance:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "metadata": {
+ "scrolled": true
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ ", , , ..., , , ]>\n"
+ ]
+ }
+ ],
+ "source": [
+ "print(u.select_atoms('resname ASP or resname GLU'))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Numerical ranges can be written as `first-last` or `first:last` where the range is **inclusive**. Note that in slicing, the `last` index is not included."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "51\n",
+ "50\n"
+ ]
+ }
+ ],
+ "source": [
+ "print(u.select_atoms('resid 50-100').n_residues)\n",
+ "print(u.residues[50:100].n_residues)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Selections can also be combined with boolean operators, and allow wildcards. \n",
+ "\n",
+ "For example, the command below selects the $C_{\\alpha}$ atoms of glutamic acid and histidine in the first 100 residues of the protein. Glutamic acid is typically named \"GLU\", but histidine can be named \"HIS\", \"HSD\", or \"HSE\" depending on its protonation state and the force field used."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ ""
+ ]
+ },
+ "execution_count": 11,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "u.select_atoms(\"(resname GLU or resname HS*) and name CA and (resid 1:100)\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "
\n",
+ "\n",
+ "**Note**\n",
+ "\n",
+ "An `AtomGroup` created from a selection is sorted and duplicate elements are removed. This is not true for an `AtomGroup` produced by slicing. Thus, slicing can be used when the order of atoms is crucial.\n",
+ "
\n",
+ "\n",
+ "
\n",
+ "\n",
+ "**See also**\n",
+ "\n",
+ "The [user guide](selections.rst) has a complete rundown of creating AtomGroups through indexing, selection language, and set methods.\n",
+ "
"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Getting atom information from AtomGroups"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "An `AtomGroup` can tell you information about the atoms inside it with a number of convenient attributes."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 12,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "['N' 'HT1' 'HT2' 'HT3' 'CA' 'HA' 'CB' 'HB1' 'HB2' 'CG' 'HG1' 'HG2' 'SD'\n",
+ " 'CE' 'HE1' 'HE2' 'HE3' 'C' 'O' 'N']\n"
+ ]
+ }
+ ],
+ "source": [
+ "print(u.atoms[:20].names)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 13,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[ 1.008 1.008 1.008 12.011 1.008 1.008 12.011 1.008 1.008 1.008\n",
+ " 12.011 15.999 14.007 1.008 12.011 1.008 12.011 1.008 12.011 1.008]\n"
+ ]
+ }
+ ],
+ "source": [
+ "print(u.atoms[50:70].masses)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "It also knows which residues and segments the atoms belong to. The `.residues` and `.segments` return a `ResidueGroup` and `SegmentGroup`, respectively."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 14,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ ", ]>\n",
+ "]>\n"
+ ]
+ }
+ ],
+ "source": [
+ "print(u.atoms[:20].residues)\n",
+ "print(u.atoms[-20:].segments)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Note that there are no duplicates in the `ResidueGroup` and `SegmentGroup` above. To get residue attributes atom-wise, you can access them directly through `AtomGroup`."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 15,
+ "metadata": {
+ "scrolled": true
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "['MET' 'MET' 'MET' 'MET' 'MET' 'MET' 'MET' 'MET' 'MET' 'MET' 'MET' 'MET'\n",
+ " 'MET' 'MET' 'MET' 'MET' 'MET' 'MET' 'MET' 'ARG']\n"
+ ]
+ }
+ ],
+ "source": [
+ "print(u.atoms[:20].resnames)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "You can group atoms together by topology attributes.\n",
+ "\n",
+ "For example, to group atoms with the same residue name and mass together:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 16,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "/Users/lily/anaconda3/envs/mdanalysis/lib/python3.7/site-packages/MDAnalysis/lib/util.py:1747: DeprecationWarning: Using or importing the ABCs from 'collections' instead of from 'collections.abc' is deprecated, and in 3.8 it will stop working\n",
+ " if isinstance(v, collections.MutableMapping):\n"
+ ]
+ },
+ {
+ "data": {
+ "text/plain": [
+ "{('ARG', 'N'): ,\n",
+ " ('ASP', 'C'): ,\n",
+ " ('ASP', 'N'): ,\n",
+ " ('LYS', 'N'): ,\n",
+ " ('ALA', 'C'): ,\n",
+ " ('ALA', 'HN'): ,\n",
+ " ('ASN', 'O'): ,\n",
+ " ('GLN', 'C'): ,\n",
+ " ('THR', 'N'): ,\n",
+ " ('ILE', 'C'): ,\n",
+ " ('GLU', 'N'): ,\n",
+ " ('LEU', 'N'): }"
+ ]
+ },
+ "execution_count": 16,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "near_met = u.select_atoms('not resname MET and (around 2 resname MET)')\n",
+ "near_met.groupby(['resnames', 'names'])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "
\n",
+ "\n",
+ "**See also**\n",
+ "\n",
+ "A complete list of topology attributes can be found [in the user guide.](placeholder)\n",
+ "
"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### AtomGroup positions and methods"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The `.positions` attribute is probably the most important information you can get from an `AtomGroup`: a `numpy.ndarray` of coordinates, with the shape (n_atoms, 3)."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 17,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[[11.664622 8.393473 -8.983231 ]\n",
+ " [11.414839 5.4344215 -6.5134845 ]\n",
+ " [ 8.959755 5.612923 -3.6132305 ]\n",
+ " [ 8.290068 3.075991 -0.79665166]\n",
+ " [ 5.011126 3.7638984 1.130355 ]]\n",
+ "(5, 3)\n"
+ ]
+ }
+ ],
+ "source": [
+ "ca = u.select_atoms('resid 1-5 and name CA')\n",
+ "print(ca.positions)\n",
+ "print(ca.positions.shape)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "A number of other quantities have been defined for an AtomGroup, including:\n",
+ "\n",
+ "* `.center_of_mass()`\n",
+ "* `.center_of_geometry()`\n",
+ "* `.total_mass()`\n",
+ "* `.total_charge()`\n",
+ "* `.radius_of_gyration()`\n",
+ "* `.bsphere()` (the bounding sphere of the selection)\n",
+ "\n",
+ "See the [user guide](data_structures.rst#atomgroups) for a complete list and description of AtomGroup methods."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 18,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[ 9.06808195 5.25614133 -3.75524844]\n"
+ ]
+ }
+ ],
+ "source": [
+ "print(ca.center_of_mass())"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "
\n",
+ "\n",
+ "**Note**\n",
+ "\n",
+ "The `.center_of_mass()` function, like many of the analysis modules in MDAnalysis, relies on having accurate mass properties available. **Particle masses may not always be available or accurate!**\n",
+ "\n",
+ "Currently, MDAnalysis assigns masses to particles based on their element or 'atom type', which is guessed from the particle name. If MDAnalysis guesses incorrectly (e.g. a calcium atom called CA is treated as a $C_{\\alpha}$), the mass of that atom will be inaccurate. If MDAnalysis has no idea what the particle is (e.g. coarse-grained beads), it will raise a warning, and give that particle a mass of 0. \n",
+ "\n",
+ "To be certain that MDAnalysis is using the correct masses, you can set them manually.\n",
+ "
"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "MDAnalysis can also create topology geometries such as bonds, angles, dihedral angles, and improper angles from an `AtomGroup`. This `AtomGroup` has a special requirement: only the atoms involved in the geometry can be in the group. For example, an `AtomGroup` used to create a bond can only have 2 atoms in it; an `AtomGroup` used to create a dihedral or improper angle must have 4 atoms."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 19,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "['N' 'HT1' 'HT2']\n"
+ ]
+ }
+ ],
+ "source": [
+ "nhh = u.atoms[:3]\n",
+ "print(nhh.names)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "After a topology object such as an angle is created, the value of the angle (in degrees) can be calculated based on the positions of the atoms."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 20,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "37.99234750892497\n"
+ ]
+ }
+ ],
+ "source": [
+ "angle_nhh = nhh.angle\n",
+ "print(angle_nhh.value())"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Note that the order of the atoms matters for angles, dihedrals, and impropers. The value returned for an angle is the angle between first and third atom, with the apex at the second. Fancy indexing is one way to get an ordered AtomGroup.\n",
+ "```\n",
+ " 3\n",
+ " /\n",
+ " /\n",
+ "2------1\n",
+ "```"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 21,
+ "metadata": {
+ "scrolled": true
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "['HT1' 'N' 'HT2']\n"
+ ]
+ }
+ ],
+ "source": [
+ "hnh = u.atoms[[1, 0, 2]]\n",
+ "print(hnh.names)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 22,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "106.20364651944931\n"
+ ]
+ }
+ ],
+ "source": [
+ "angle_hnh = hnh.angle\n",
+ "print(angle_hnh.value())"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Working with trajectories"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The trajectory of a Universe contains the changing coordinate information. The number of frames in a trajectory is its length:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 23,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "98\n"
+ ]
+ }
+ ],
+ "source": [
+ "print(len(u.trajectory))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The standard way to assess the information of each frame in a trajectory is to iterate over it. When the timestep changes, the universe only contains information associated with that timestep."
]
},
{
"cell_type": "code",
- "execution_count": 1,
+ "execution_count": 24,
+ "metadata": {
+ "scrolled": true
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Frame: 0, Time: 1 ps, Rgyr: 16.669 A\n",
+ "Frame: 1, Time: 2 ps, Rgyr: 16.6732 A\n",
+ "Frame: 2, Time: 3 ps, Rgyr: 16.7315 A\n",
+ "Frame: 3, Time: 4 ps, Rgyr: 16.7223 A\n",
+ "Frame: 4, Time: 5 ps, Rgyr: 16.744 A\n",
+ "Frame: 5, Time: 6 ps, Rgyr: 16.7185 A\n",
+ "Frame: 6, Time: 7 ps, Rgyr: 16.7741 A\n",
+ "Frame: 7, Time: 8 ps, Rgyr: 16.7764 A\n",
+ "Frame: 8, Time: 9 ps, Rgyr: 16.7894 A\n",
+ "Frame: 9, Time: 10 ps, Rgyr: 16.8289 A\n",
+ "Frame: 10, Time: 11 ps, Rgyr: 16.8521 A\n",
+ "Frame: 11, Time: 12 ps, Rgyr: 16.8549 A\n",
+ "Frame: 12, Time: 13 ps, Rgyr: 16.8723 A\n",
+ "Frame: 13, Time: 14 ps, Rgyr: 16.9108 A\n",
+ "Frame: 14, Time: 15 ps, Rgyr: 16.9494 A\n",
+ "Frame: 15, Time: 16 ps, Rgyr: 16.981 A\n",
+ "Frame: 16, Time: 17 ps, Rgyr: 17.0033 A\n",
+ "Frame: 17, Time: 18 ps, Rgyr: 17.0196 A\n",
+ "Frame: 18, Time: 19 ps, Rgyr: 17.0784 A\n",
+ "Frame: 19, Time: 20 ps, Rgyr: 17.1265 A\n"
+ ]
+ }
+ ],
+ "source": [
+ "for ts in u.trajectory[:20]:\n",
+ " time = u.trajectory.time\n",
+ " rgyr = u.atoms.radius_of_gyration()\n",
+ " print(\"Frame: {:3d}, Time: {:4.0f} ps, Rgyr: {:g} A\".format(ts.frame, time, rgyr))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "After iteration, the trajectory 'resets' back to the first frame."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 25,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "0\n"
+ ]
+ }
+ ],
+ "source": [
+ "print(u.trajectory.frame)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "You can set the timestep of the trajectory with the frame index:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 26,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "10\n"
+ ]
+ }
+ ],
+ "source": [
+ "print(u.trajectory[10].frame)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "This persists until the timestep is next changed."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 27,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Frame: 10, Time: 11 ps, Rgyr: 16.8521 A\n"
+ ]
+ }
+ ],
+ "source": [
+ "frame = u.trajectory.frame\n",
+ "time = u.trajectory.time\n",
+ "rgyr = u.atoms.radius_of_gyration()\n",
+ "print(\"Frame: {:3d}, Time: {:4.0f} ps, Rgyr: {:g} A\".format(frame, time, rgyr))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Generally, trajectory analysis first collects frame-wise data in a list."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 28,
"metadata": {},
"outputs": [],
"source": [
- "import MDAnalysis as mda"
+ "rgyr = []\n",
+ "time = []\n",
+ "protein = u.select_atoms(\"protein\")\n",
+ "for ts in u.trajectory:\n",
+ " time.append(u.trajectory.time)\n",
+ " rgyr.append(protein.radius_of_gyration())"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "This can then be converted into other data structures, such as a numpy array or a pandas DataFrame. It can be plotted (as below), or used for further analysis."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 29,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "
"
+ ]
+ },
+ "metadata": {
+ "needs_background": "light"
+ },
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "%matplotlib inline\n",
+ "\n",
+ "rgyr_df.plot(title='Radius of gyration')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Dynamic selection"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "As seen above, coordinates change while iterating over the trajectory. Therefore, properties calculated from the coordinates, such as the radius of gyration, also change. \n",
+ "\n",
+ "Selections are often defined on static properties that do not change when moving through a trajectory. Above, the static selection is all the atoms that are in a protein. You can define the selection once and then recalculate the property of interest for each frame of the trajectory.\n",
+ "\n",
+ "However, some selections contain distance-dependent queries (such as `around` or `point`, [see selection keywords for more details](selections.rst)). In this case, the selection should be updated for each time step using a dynamic selection by setting the keyword `updating=True`. This command gives an `UpdatingAtomGroup` rather than a static `AtomGroup`.\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 31,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "\n"
+ ]
+ },
+ {
+ "data": {
+ "text/plain": [
+ ""
+ ]
+ },
+ "execution_count": 31,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "dynamic = u.select_atoms('around 2 resname ALA', updating=True)\n",
+ "print(type(dynamic))\n",
+ "dynamic"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 32,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "\n"
+ ]
+ },
+ {
+ "data": {
+ "text/plain": [
+ ""
+ ]
+ },
+ "execution_count": 32,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "static = u.select_atoms('around 2 resname ALA')\n",
+ "print(type(static))\n",
+ "static"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "When you call the next frame of the universe, the atoms in the dynamic selection are updated, whereas the atoms in the static selection remain the same."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 33,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ ""
+ ]
+ },
+ "execution_count": 33,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "u.trajectory.next()\n",
+ "dynamic"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 34,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ ""
+ ]
+ },
+ "execution_count": 34,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "static"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Writing out coordinates"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "MDAnalysis supports writing data out into a range of file formats, including both single frame formats (e.g. PDB, GRO) and trajectory writers (e.g. XTC, DCD, and multi-frame PDB files)."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Single frame"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The most straightforward way to write to a file that can only hold a single frame is to use the `write()` method of any `AtomGroup`. MDAnalysis uses the file extension to determine the output file format. For instance, to only write out the $C_{\\alpha}$ atoms to a file in GRO format:\n",
+ "\n",
+ "```python\n",
+ "ca = u.select_atoms('name CA')\n",
+ "ca.write('calphas.gro')\n",
+ "```"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Trajectories"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The standard way to write out trajectories is to:\n",
+ "\n",
+ "1. Open a trajectory `Writer` and specify how many atoms a frame will contain\n",
+ "2. Iterate through the trajectory and write coordinates frame-by-frame with `Writer.write()`\n",
+ "3. If you do not use the context manager and the `with` statement below, you then need to close the trajectory with `.close()`.\n",
+ "\n",
+ "For instance, to write out the $C_{\\alpha}$ atoms to a trajectory in the XTC format:\n",
+ "\n",
+ "```python\n",
+ "ca = u.select_atoms('name CA')\n",
+ "with mda.Writer('calphas.xtc', ca.n_atoms) as w:\n",
+ " for ts in u.trajectory:\n",
+ " w.write(ca)\n",
+ "```"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Analysis"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "MDAnalysis comes with a diverse set of analysis modules, and the building blocks to implement your own.\n",
+ "\n",
+ "The majority of these follow a common interface:\n",
+ "\n",
+ "1. Initialise the analysis with a `Universe` and other required parameters.\n",
+ "2. Run the analysis with `.run()`. Optional arguments are the `start` frame index, `stop` frame index, `step` size, and toggling `verbose`. The default is to run analysis on the whole trajectory.\n",
+ "3. Results are stored within the class.\n",
+ "4. Often, a function is available to operate on single frames.\n",
+ "\n",
+ "**However, not all analysis uses this model. It is important to check the documentation for each analysis.** You can also see examples in the [Example gallery](placeholder).\n",
+ "\n",
+ "Below, simple RMSD analysis is shown. The `rms` module follows the interface above."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### RMSD"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Not all sub-modules of MDAnalysis are imported with `import MDAnalysis`. Most analysis modules have to be imported explicitly. "
]
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 35,
"metadata": {},
"outputs": [],
- "source": []
+ "source": [
+ "from MDAnalysis.analysis import rms"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "MDAnalysis provides a `rmsd()` function for calculating the RMSD between two numpy arrays of coordinates."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 36,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "6.852774844656239"
+ ]
+ },
+ "execution_count": 36,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "bb = u.select_atoms('backbone')\n",
+ "\n",
+ "u.trajectory[0] # first frame\n",
+ "first = bb.positions\n",
+ "\n",
+ "u.trajectory[-1] #last frame\n",
+ "last = bb.positions\n",
+ "\n",
+ "rms.rmsd(first, last)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "An `RMSD` class is also provided to operate on trajectories. \n",
+ "\n",
+ "Below, the `RMSD` class is created. \n",
+ "* The first argument is the `AtomGroup` or `Universe` for which the RMSD is calculated. \n",
+ "* As a reference `AtomGroup` or `Universe` is not given as the second argument, the default is to align to the current frame of the first argument. Here it is set to the first frame. \n",
+ "* We choose to align the trajectory over the backbone atoms, and then compute the RMSD for the backbone atoms (`select`). \n",
+ "* Then, without re-aligning the trajectory, the RMSD is also computed (`groupselections`) for the $C_{\\alpha}$ atoms (`name CA`) and every protein atom (`protein`).\n",
+ "\n",
+ "The RMSD is computed when we call `.run()`."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 37,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ ""
+ ]
+ },
+ "execution_count": 37,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "u.trajectory[0] # set to first frame\n",
+ "rmsd_analysis = rms.RMSD(u, select='backbone', groupselections=['name CA', 'protein'])\n",
+ "rmsd_analysis.run(verbose=True)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The results are stored in the `.rmsd` attribute. This is an array with the shape (n_frames, 2 + n_selections)."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 38,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "(98, 5)\n"
+ ]
+ }
+ ],
+ "source": [
+ "print(rmsd_analysis.rmsd.shape)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We can interpret this as an array with 98 rows and 5 columns. Each row is the RMSD associated with a frame in the trajectory. The columns are as follows:\n",
+ "\n",
+ "1. Frame number\n",
+ "2. Time (ps)\n",
+ "3. RMSD (backbone)\n",
+ "4. RMSD (C-alpha)\n",
+ "5. RMSD (protein)\n",
+ "\n",
+ "We can turn this into a `pandas` DataFrame and plot the results."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 39,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "
\n",
+ "\n",
+ "**See also**\n",
+ "\n",
+ "See the [analysis documentation](analysis.rst) for more information.\n",
+ "
"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## References"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "When using MDAnalysis in published work, please cite both these papers:\n",
+ "\n",
+ "* N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein. MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. *J. Comput. Chem.* **32** (2011), 2319–2327. doi:[10.1002/jcc.21787](http://dx.doi.org/10.1002/jcc.21787)\n",
+ "\n",
+ "* R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N.\n",
+ " Melo, S. L. Seyler, D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney,\n",
+ " and O. Beckstein. [MDAnalysis: A Python package for the rapid analysis of\n",
+ " molecular dynamics simulations](http://conference.scipy.org/proceedings/scipy2016/oliver_beckstein.html). In S. Benthall and S. Rostrup, editors,\n",
+ " *Proceedings of the 15th Python in Science Conference*, pages 98-105,\n",
+ " Austin, TX, 2016. SciPy. doi:[10.25080/Majora-629e541a-00e](https://doi.org/10.25080/Majora-629e541a-00e)\n",
+ " \n",
+ "MDAnalysis includes many algorithms and modules that should also be individually cited. For example, if you use the `MDAnalysis.analysis.rms` or `MDAnalysis.analysis.align` modules, please cite:\n",
+ "\n",
+ "* Douglas L. Theobald. Rapid calculation of RMSD using a\n",
+ " quaternion-based characteristic polynomial. *Acta Crystallographica A*\n",
+ " **61** (2005), 478-480.\n",
+ "* Pu Liu, Dmitris K. Agrafiotis, and Douglas L. Theobald. Fast\n",
+ " determination of the optimal rotational matrix for macromolecular\n",
+ " superpositions. *J. Comput. Chem.* **31** (2010), 1561–1563.\n",
+ " \n",
+ "The primary sources of each module will be in their documentation."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Automatic citations with duecredit"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Citations can be also be automatically generated using `duecredit`. [Complete installation and usage instructions can be found in the user guide](placeholder), but it is simple to generate a list of citations for your python script `my_script.py`.\n",
+ "\n",
+ "```console\n",
+ "$ python -m duecredit my_script.py\n",
+ "```\n",
+ "\n",
+ "This extracts citations into a hidden file, which can then be exported to different formats. For example, to display them as BibTeX, use the command:\n",
+ "\n",
+ "```console\n",
+ "$ duecredit summary --format=bibtex\n",
+ "```"
+ ]
}
],
"metadata": {
+ "celltoolbar": "Raw Cell Format",
"kernelspec": {
"display_name": "Python (mdanalysis)",
"language": "python",
@@ -42,6 +1631,10 @@
"pygments_lexer": "ipython3",
"version": "3.7.3"
},
+ "nbsphinx-toctree": {
+ "maxdepth": 2,
+ "numbered": true
+ },
"toc": {
"base_numbering": 1,
"nav_menu": {},
diff --git a/environment.yml b/environment.yml
index f324ce5c2..e00c49080 100644
--- a/environment.yml
+++ b/environment.yml
@@ -6,4 +6,8 @@ dependencies:
- jupyter_contrib_nbextensions
- nglview
- pip:
- - mdanalysis
\ No newline at end of file
+ - mdanalysis
+ - mdanalysistests
+ - sphinx_rtd_theme
+ - nbsphinx
+ - sphinx-gallery
\ No newline at end of file