{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Aligning a trajectory to itself\n",
    "\n",
    "We use `align.AlignTraj` to align a trajectory to a reference frame and write it to a file.\n",
    "\n",
    "**Last updated:** December 2022 with MDAnalysis 2.4.0-dev0\n",
    "\n",
    "**Minimum version of MDAnalysis:** 1.0.0\n",
    "\n",
    "**Packages required:**\n",
    "    \n",
    "* MDAnalysis (<a data-cite=\"michaud-agrawal_mdanalysis_2011\" href=\"https://doi.org/10.1002/jcc.21787\">Michaud-Agrawal *et al.*, 2011</a>, <a data-cite=\"gowers_mdanalysis_2016\" href=\"https://doi.org/10.25080/Majora-629e541a-00e\">Gowers *et al.*, 2016</a>)\n",
    "* MDAnalysisTests\n",
    "\n",
    "**See also**\n",
    "\n",
    "* [Aligning a trajectory to a frame from another](aligning_trajectory.ipynb)\n",
    "* [Aligning a structure to another](aligning_structure_to_another.ipynb)\n",
    "* [RMSD](rmsd.ipynb)\n",
    "   \n",
    "\n",
    "<div class=\"alert alert-info\">\n",
    "    \n",
    "**Note**\n",
    "\n",
    "MDAnalysis implements RMSD calculation using the fast QCP algorithm (<a data-cite=\"theobald_rapid_2005\" href=\"https://doi.org/10.1107/S0108767305015266\">Theobald, 2005</a>) and a rotation matrix *R* that minimises the RMSD (<a data-cite=\"liu_fast_2009\" href=\"https://doi.org/10.1002/jcc.21439\">Liu *et al.*, 2009</a>). Please cite (<a data-cite=\"theobald_rapid_2005\" href=\"https://doi.org/10.1107/S0108767305015266\">Theobald, 2005</a>) and (<a data-cite=\"liu_fast_2009\" href=\"https://doi.org/10.1002/jcc.21439\">Liu *et al.*, 2009</a>) when using the ``MDAnalysis.analysis.align`` module in published work.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-19T05:56:34.905776Z",
     "iopub.status.busy": "2021-05-19T05:56:34.904828Z",
     "iopub.status.idle": "2021-05-19T05:56:35.633469Z",
     "shell.execute_reply": "2021-05-19T05:56:35.634004Z"
    },
    "scrolled": false,
    "ExecuteTime": {
     "end_time": "2023-06-09T13:05:23.592176544Z",
     "start_time": "2023-06-09T13:05:23.151189108Z"
    }
   },
   "outputs": [],
   "source": [
    "import MDAnalysis as mda\n",
    "from MDAnalysis.analysis import align, rms\n",
    "from MDAnalysis.tests.datafiles import PSF, DCD"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading files\n",
    "\n",
    "The test files we will be working with here feature adenylate kinase (AdK), a phosophotransferase enzyme. (<a data-cite=\"beckstein_zipping_2009\" href=\"https://doi.org/10.1016/j.jmb.2009.09.009\">Beckstein *et al.*, 2009</a>) The trajectory samples a transition from a closed to an open conformation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-19T05:56:35.641643Z",
     "iopub.status.busy": "2021-05-19T05:56:35.640355Z",
     "iopub.status.idle": "2021-05-19T05:56:35.990432Z",
     "shell.execute_reply": "2021-05-19T05:56:35.991124Z"
    },
    "ExecuteTime": {
     "end_time": "2023-06-09T13:05:23.942673356Z",
     "start_time": "2023-06-09T13:05:23.593672649Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/pbarletta/mambaforge/envs/guide/lib/python3.9/site-packages/MDAnalysis/coordinates/DCD.py:165: DeprecationWarning: DCDReader currently makes independent timesteps by copying self.ts while other readers update self.ts inplace. This behavior will be changed in 3.0 to be the same as other readers. Read more at https://github.com/MDAnalysis/mdanalysis/issues/3889 to learn if this change in behavior might affect you.\n",
      "  warnings.warn(\"DCDReader currently makes independent timesteps\"\n"
     ]
    }
   ],
   "source": [
    "mobile = mda.Universe(PSF, DCD)\n",
    "ref = mda.Universe(PSF, DCD)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Aligning a trajectory to the first frame\n",
    "\n",
    "While [align.alignto](https://docs.mdanalysis.org/stable/documentation_pages/analysis/align.html#MDAnalysis.analysis.align.alignto) aligns structures, or a frame of a trajectory, `align.AlignTraj` ([API docs](https://docs.mdanalysis.org/stable/documentation_pages/analysis/align.html#MDAnalysis.analysis.align.AlignTraj)) efficiently aligns an entire trajectory to a reference. \n",
    "\n",
    "We first check the [root mean square deviation (RMSD) values](rmsd.ipynb) of our unaligned trajectory, so we can compare results (please see the linked notebook for an explanation of RMSD). The code below sets the ``mobile`` trajectory to the last frame by indexing the last timestep, ``ref`` to the first frame by indexing the first timestep, and computes the root mean squared deviation between the alpha-carbon positions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-19T05:56:35.996879Z",
     "iopub.status.busy": "2021-05-19T05:56:35.996217Z",
     "iopub.status.idle": "2021-05-19T05:56:36.004323Z",
     "shell.execute_reply": "2021-05-19T05:56:36.004701Z"
    },
    "ExecuteTime": {
     "end_time": "2023-06-09T13:05:23.943555057Z",
     "start_time": "2023-06-09T13:05:23.928868940Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Unaligned RMSD: 6.84\n"
     ]
    }
   ],
   "source": [
    "mobile.trajectory[-1]  # set mobile trajectory to last frame\n",
    "ref.trajectory[0]  # set reference trajectory to first frame\n",
    "\n",
    "mobile_ca = mobile.select_atoms('name CA')\n",
    "ref_ca = ref.select_atoms('name CA')\n",
    "unaligned_rmsd = rms.rmsd(mobile_ca.positions, ref_ca.positions, superposition=False)\n",
    "print(f\"Unaligned RMSD: {unaligned_rmsd:.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can align the trajectory. We have already set ``ref`` to the first frame. In the cell below, we load the positions of the trajectory into memory so we can modify the trajectory in Python. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-19T05:56:36.009560Z",
     "iopub.status.busy": "2021-05-19T05:56:36.008743Z",
     "iopub.status.idle": "2021-05-19T05:56:36.090226Z",
     "shell.execute_reply": "2021-05-19T05:56:36.090631Z"
    },
    "ExecuteTime": {
     "end_time": "2023-06-09T13:05:24.023863380Z",
     "start_time": "2023-06-09T13:05:23.932962914Z"
    }
   },
   "outputs": [],
   "source": [
    "aligner = align.AlignTraj(mobile, ref, select='name CA', in_memory=True).run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you don't have enough memory to do that, write the trajectory out to a file and reload it into MDAnalysis (uncomment the cell below). Otherwise, you don't have to run it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-19T05:56:36.094367Z",
     "iopub.status.busy": "2021-05-19T05:56:36.093685Z",
     "iopub.status.idle": "2021-05-19T05:56:36.095476Z",
     "shell.execute_reply": "2021-05-19T05:56:36.095881Z"
    },
    "ExecuteTime": {
     "end_time": "2023-06-09T13:05:24.031106057Z",
     "start_time": "2023-06-09T13:05:24.026374041Z"
    }
   },
   "outputs": [],
   "source": [
    "# aligner = align.AlignTraj(mobile, ref, select='backbone', \n",
    "#                           filename='aligned_to_first_frame.dcd').run()\n",
    "# mobile = mda.Universe(PSF, 'aligned_to_first_frame.dcd')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can see that the RMSD has reduced (minorly)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-19T05:56:36.101933Z",
     "iopub.status.busy": "2021-05-19T05:56:36.100806Z",
     "iopub.status.idle": "2021-05-19T05:56:36.104896Z",
     "shell.execute_reply": "2021-05-19T05:56:36.105359Z"
    },
    "ExecuteTime": {
     "end_time": "2023-06-09T13:05:24.085485673Z",
     "start_time": "2023-06-09T13:05:24.030076994Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Aligned RMSD: 6.81\n"
     ]
    }
   ],
   "source": [
    "mobile.trajectory[-1]  # set mobile trajectory to last frame\n",
    "ref.trajectory[0]  # set reference trajectory to first frame\n",
    "\n",
    "mobile_ca = mobile.select_atoms('name CA')\n",
    "ref_ca = ref.select_atoms('name CA')\n",
    "aligned_rmsd = rms.rmsd(mobile_ca.positions, ref_ca.positions, superposition=False)\n",
    "\n",
    "print(f\"Aligned RMSD: {aligned_rmsd:.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Aligning a trajectory to the third frame\n",
    "\n",
    "We can align the trajectory to any frame: for example, the third one. The procedure is much the same, except that we must set ``ref`` to the third frame by indexing the third timestep."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-19T05:56:36.109732Z",
     "iopub.status.busy": "2021-05-19T05:56:36.108816Z",
     "iopub.status.idle": "2021-05-19T05:56:36.112429Z",
     "shell.execute_reply": "2021-05-19T05:56:36.112811Z"
    },
    "ExecuteTime": {
     "end_time": "2023-06-09T13:05:24.085920220Z",
     "start_time": "2023-06-09T13:05:24.071254242Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Aligned RMSD: 6.73\n"
     ]
    }
   ],
   "source": [
    "mobile.trajectory[-1]  # set mobile trajectory to last frame\n",
    "ref.trajectory[2]  # set reference trajectory to third frame\n",
    "\n",
    "aligned_rmsd_3 = rms.rmsd(mobile.atoms.positions, ref.atoms.positions, superposition=False)\n",
    "\n",
    "print(f\"Aligned RMSD: {aligned_rmsd_3:.2f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-19T05:56:36.116539Z",
     "iopub.status.busy": "2021-05-19T05:56:36.115939Z",
     "iopub.status.idle": "2021-05-19T05:56:36.200995Z",
     "shell.execute_reply": "2021-05-19T05:56:36.201497Z"
    },
    "ExecuteTime": {
     "end_time": "2023-06-09T13:05:24.176553080Z",
     "start_time": "2023-06-09T13:05:24.071471277Z"
    }
   },
   "outputs": [],
   "source": [
    "aligner = align.AlignTraj(mobile, ref, select='all', in_memory=True).run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-19T05:56:36.206527Z",
     "iopub.status.busy": "2021-05-19T05:56:36.205346Z",
     "iopub.status.idle": "2021-05-19T05:56:36.209965Z",
     "shell.execute_reply": "2021-05-19T05:56:36.210328Z"
    },
    "ExecuteTime": {
     "end_time": "2023-06-09T13:05:24.227607962Z",
     "start_time": "2023-06-09T13:05:24.179898326Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Aligned RMSD, all-atom: 6.72\n"
     ]
    }
   ],
   "source": [
    "mobile.trajectory[-1]  # set mobile trajectory to last frame\n",
    "ref.trajectory[2]  # set reference trajectory to third frame\n",
    "\n",
    "aligned_rmsd_3 = rms.rmsd(mobile.atoms.positions, ref.atoms.positions, superposition=False)\n",
    "print(f\"Aligned RMSD, all-atom: {aligned_rmsd_3:.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "[1] Oliver Beckstein, Elizabeth&nbsp;J. Denning, Juan&nbsp;R. Perilla, and Thomas&nbsp;B. Woolf.\n",
    "Zipping and <span class=\"bibtex-protected\">Unzipping</span> of <span class=\"bibtex-protected\">Adenylate</span> <span class=\"bibtex-protected\">Kinase</span>: <span class=\"bibtex-protected\">Atomistic</span> <span class=\"bibtex-protected\">Insights</span> into the <span class=\"bibtex-protected\">Ensemble</span> of <span class=\"bibtex-protected\">Open</span>↔<span class=\"bibtex-protected\">Closed</span> <span class=\"bibtex-protected\">Transitions</span>.\n",
    "<em>Journal of Molecular Biology</em>, 394(1):160–176, November 2009.\n",
    "00107.\n",
    "URL: <a href=\"https://linkinghub.elsevier.com/retrieve/pii/S0022283609011164\">https://linkinghub.elsevier.com/retrieve/pii/S0022283609011164</a>, <a href=\"https://doi.org/10.1016/j.jmb.2009.09.009\">doi:10.1016/j.jmb.2009.09.009</a>.\n",
    "\n",
    "[2] Richard&nbsp;J. Gowers, Max Linke, Jonathan Barnoud, Tyler J.&nbsp;E. Reddy, Manuel&nbsp;N. Melo, Sean&nbsp;L. Seyler, Jan Domański, David&nbsp;L. Dotson, Sébastien Buchoux, Ian&nbsp;M. Kenney, and Oliver Beckstein.\n",
    "<span class=\"bibtex-protected\">MDAnalysis</span>: <span class=\"bibtex-protected\">A</span> <span class=\"bibtex-protected\">Python</span> <span class=\"bibtex-protected\">Package</span> for the <span class=\"bibtex-protected\">Rapid</span> <span class=\"bibtex-protected\">Analysis</span> of <span class=\"bibtex-protected\">Molecular</span> <span class=\"bibtex-protected\">Dynamics</span> <span class=\"bibtex-protected\">Simulations</span>.\n",
    "<em>Proceedings of the 15th Python in Science Conference</em>, pages 98–105, 2016.\n",
    "00152.\n",
    "URL: <a href=\"https://conference.scipy.org/proceedings/scipy2016/oliver_beckstein.html\">https://conference.scipy.org/proceedings/scipy2016/oliver_beckstein.html</a>, <a href=\"https://doi.org/10.25080/Majora-629e541a-00e\">doi:10.25080/Majora-629e541a-00e</a>.\n",
    "\n",
    "[3] Pu&nbsp;Liu, Dimitris&nbsp;K. Agrafiotis, and Douglas&nbsp;L. Theobald.\n",
    "Fast determination of the optimal rotational matrix for macromolecular superpositions.\n",
    "<em>Journal of Computational Chemistry</em>, pages n/a–n/a, 2009.\n",
    "URL: <a href=\"http://doi.wiley.com/10.1002/jcc.21439\">http://doi.wiley.com/10.1002/jcc.21439</a>, <a href=\"https://doi.org/10.1002/jcc.21439\">doi:10.1002/jcc.21439</a>.\n",
    "\n",
    "[4] Naveen Michaud-Agrawal, Elizabeth&nbsp;J. Denning, Thomas&nbsp;B. Woolf, and Oliver Beckstein.\n",
    "<span class=\"bibtex-protected\">MDAnalysis</span>: <span class=\"bibtex-protected\">A</span> toolkit for the analysis of molecular dynamics simulations.\n",
    "<em>Journal of Computational Chemistry</em>, 32(10):2319–2327, July 2011.\n",
    "00778.\n",
    "URL: <a href=\"http://doi.wiley.com/10.1002/jcc.21787\">http://doi.wiley.com/10.1002/jcc.21787</a>, <a href=\"https://doi.org/10.1002/jcc.21787\">doi:10.1002/jcc.21787</a>.\n",
    "\n",
    "[5] Douglas&nbsp;L. Theobald.\n",
    "Rapid calculation of <span class=\"bibtex-protected\">RMSDs</span> using a quaternion-based characteristic polynomial.\n",
    "<em>Acta Crystallographica Section A Foundations of Crystallography</em>, 61(4):478–480, July 2005.\n",
    "00127.\n",
    "URL: <a href=\"http://scripts.iucr.org/cgi-bin/paper?S0108767305015266\">http://scripts.iucr.org/cgi-bin/paper?S0108767305015266</a>, <a href=\"https://doi.org/10.1107/S0108767305015266\">doi:10.1107/S0108767305015266</a>."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Raw Cell Format",
  "kernelspec": {
   "display_name": "Python 3.9.15 ('mda-user-guide')",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.15"
  },
  "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": true
  },
  "vscode": {
   "interpreter": {
    "hash": "7b52aa17ef4e9358c0e91ff3f0bf50d10667bb8b55636d4b9fb967a2ff94bd8c"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}