{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Centering a trajectory in the box\n", "\n", "Here we use MDAnalysis transformations to make a protein whole, center it in the box, and then wrap the water back into the box. We then look at how to do this on-the-fly.\n", "\n", "**Last updated:** July 3, 2020 with MDAnalysis 1.0.0\n", "\n", "**Minimum version of MDAnalysis:** 1.0.0\n", "\n", "**Packages required:**\n", " \n", "* MDAnalysis (Michaud-Agrawal *et al.*, 2011, Gowers *et al.*, 2016)\n", "* MDAnalysisTests\n", "* [nglview](http://nglviewer.org/nglview/latest/) (Nguyen *et al.*, 2018)\n", "\n", "**See also:**\n", "\n", "* [On-the-fly transformations](../../trajectories/transformations.rst)\n", "* [On-the-fly transformations (blog post)](https://www.mdanalysis.org/2020/03/09/on-the-fly-transformations/)\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "33719b2e23374559b80febdfb86e437e", "version_major": 2, "version_minor": 0 }, "text/plain": [ "_ColormakerRegistry()" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import MDAnalysis as mda\n", "from MDAnalysis.tests.datafiles import TPR, XTC\n", "import numpy as np\n", "import nglview as nv" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Loading files\n", "\n", "The test files we will be working with here are trajectories of a adenylate kinase (AdK), a phosophotransferase enzyme. (Beckstein *et al.*, 2009)\n", "\n", "For the step-by-step transformations, we need to load the trajectory into memory so that our changes to the coordinates persist. If your trajectory is too large for that, see the [on-the-fly transformation](#Doing-all-this-on-the-fly) section for how to do this out-of-memory." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "u = mda.Universe(TPR, XTC, in_memory=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Before transformation\n", "\n", "If you have NGLView installed, you can view the trajectory as it currently is." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "view = nv.show_mdanalysis(u)\n", "view.add_representation('point', 'resname SOL')\n", "view.center()\n", "view" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Otherwise, we embed a gif of it below.\n", "\n", "![original](original.gif)\n", "\n", "For easier analysis and nicer visualisation, we want to center this protein in the box." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Unwrapping the protein\n", "\n", "The first step is to \"unwrap\" the protein from the border of the box, to make the protein whole. MDAnalysis provides the `AtomGroup.unwrap` function to do this easily. Note that this function requires your universe to have bonds in it.\n", "\n", "We loop over the trajectory to unwrap for each frame." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "protein = u.select_atoms('protein')\n", "\n", "for ts in u.trajectory:\n", " protein.unwrap(compound='fragments')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As you can see, the protein is now whole, but not centered." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "unwrapped = nv.show_mdanalysis(u)\n", "unwrapped.add_representation('point', 'resname SOL')\n", "unwrapped.center()\n", "unwrapped" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Over the course of the trajectory it leaves the box.\n", "\n", "![unwrapped](unwrapped.gif)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Centering in the box\n", "\n", "The next step is to center the protein in the box. We calculate the center-of-mass of the protein and the center of the box for each timestep. We then move *all* the atoms so that the protein center-of-mass is in the center of the box.\n", "\n", "If you don't have masses in your trajectory, try using the `center_of_geometry`." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "for ts in u.trajectory:\n", " protein_center = protein.center_of_mass(pbc=True)\n", " dim = ts.triclinic_dimensions\n", " box_center = np.sum(dim, axis=0) / 2\n", " u.atoms.translate(box_center - protein_center)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The protein is now in the center of the box, but the solvent is likely outside it, as we have just moved all the atoms." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "centered = nv.show_mdanalysis(u)\n", "centered.add_representation('point', 'resname SOL')\n", "centered.center()\n", "centered" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "![centered](centered.gif)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Wrapping the solvent back into the box\n", "\n", "Luckily, MDAnalysis also has `AtomGroup.wrap` to wrap molecules back into the box. Our trajectory has dimensions defined, which the function will find automatically. If your trajectory does not, or you wish to use a differently sized box, you can pass in a ``box`` with dimensions in the form ``[lx, ly, lz, alpha, beta, gamma]``." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "not_protein = u.select_atoms('not protein')\n", "\n", "for ts in u.trajectory:\n", " not_protein.wrap(compound='residues')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And now it is centered!" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "wrapped = nv.show_mdanalysis(u)\n", "wrapped.add_representation('point', 'resname SOL')\n", "wrapped.center()\n", "wrapped" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "![wrapped](wrapped.gif)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Doing all this on-the-fly\n", "\n", "Running all the transformations above can be difficult if your trajectory is large, or your computational resources are limited. Use on-the-fly transformations to keep your data out-of-memory.\n", "\n", "Some common transformations are defined in `MDAnalysis.transformations`." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "import MDAnalysis.transformations as trans" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We re-load our universe." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "u2 = mda.Universe(TPR, XTC)\n", "\n", "protein2 = u2.select_atoms('protein')\n", "not_protein2 = u2.select_atoms('not protein')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "From version 1.0.0 onwards, the `MDAnalysis.transformations` module contains `wrap` and `unwrap` functions that correspond with the `AtomGroup` versions above. You can only use `add_transformations` *once*, so pass them all at the same time." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "transforms = [trans.unwrap(protein2),\n", " trans.center_in_box(protein2, wrap=True),\n", " trans.wrap(not_protein2)]\n", "\n", "u2.trajectory.add_transformations(*transforms)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "otf = nv.show_mdanalysis(u2)\n", "otf.add_representation('point', 'resname SOL')\n", "otf.center()\n", "otf" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "![on the fly](on_the_fly.gif)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## References\n", "\n", "[1] Richard J. Gowers, Max Linke, Jonathan Barnoud, Tyler J. E. Reddy, Manuel N. Melo, Sean L. Seyler, Jan Domański, David L. Dotson, Sébastien Buchoux, Ian M. Kenney, and Oliver Beckstein.\n", "MDAnalysis: A Python Package for the Rapid Analysis of Molecular Dynamics Simulations.\n", "Proceedings of the 15th Python in Science Conference, pages 98–105, 2016.\n", "00152.\n", "URL: https://conference.scipy.org/proceedings/scipy2016/oliver_beckstein.html, doi:10.25080/Majora-629e541a-00e.\n", "\n", "[2] Naveen Michaud-Agrawal, Elizabeth J. Denning, Thomas B. Woolf, and Oliver Beckstein.\n", "MDAnalysis: A toolkit for the analysis of molecular dynamics simulations.\n", "Journal of Computational Chemistry, 32(10):2319–2327, July 2011.\n", "00778.\n", "URL: http://doi.wiley.com/10.1002/jcc.21787, doi:10.1002/jcc.21787.\n", "\n", "[3] Hai Nguyen, David A Case, and Alexander S Rose.\n", "NGLview–interactive molecular graphics for Jupyter notebooks.\n", "Bioinformatics, 34(7):1241–1242, April 2018.\n", "00024.\n", "URL: https://academic.oup.com/bioinformatics/article/34/7/1241/4721781, doi:10.1093/bioinformatics/btx789.\n", "\n", "[4] Oliver Beckstein, Elizabeth J. Denning, Juan R. Perilla, and Thomas B. Woolf.\n", "Zipping and Unzipping of Adenylate Kinase: Atomistic Insights into the Ensemble of OpenClosed Transitions.\n", "Journal of Molecular Biology, 394(1):160–176, November 2009.\n", "00107.\n", "URL: https://linkinghub.elsevier.com/retrieve/pii/S0022283609011164, doi:10.1016/j.jmb.2009.09.009." ] } ], "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.6" }, "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 }