{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculating hydrogen bonds: advanced selections\n",
    "\n",
    "We will find all intramolcular hydrogen bonds of a protein without passing any atom selections. We will also look at how to use more advanced atom selections than we saw in [Calculating hydrogen bonds: the basiscs](hbonds.ipynb).\n",
    "\n",
    "**Last updated:** December 2022\n",
    "\n",
    "**Minimum version of MDAnalysis:** 2.0.0-dev0\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",
    "* [numpy](https://numpy.org)\n",
    "\n",
    "**See also**\n",
    "\n",
    "* [Calculating hydrogen bonds: the basics](hbonds.ipynb)\n",
    "* [Calculating hydrogen bond lifetimes](hbonds-lifetimes.ipynb)\n",
    "\n",
    "<div class=\"alert alert-info\">\n",
    "    \n",
    "**Note**\n",
    "\n",
    "Please cite [Smith et al. (2018)](http://dx.doi.org/10.1039/C9CP01532A) when using HydrogenBondAnaysis in published work.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "import MDAnalysis as mda\n",
    "from MDAnalysis.tests.datafiles import PSF, DCD\n",
    "from MDAnalysis.analysis.hydrogenbonds import HydrogenBondAnalysis\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/pbarletta/mambaforge/envs/mda-user-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 behaviour will be changed in 3.0 to be the same as other readers\n",
      "  warnings.warn(\"DCDReader currently makes independent timesteps\"\n"
     ]
    }
   ],
   "source": [
    "u = mda.Universe(PSF, DCD)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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>)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Find all hydrogen bonds\n",
    "The simplest use case is to allow `HydrogenBondAnalysis` to guess the acceptor and hydrogen atoms, then to identify donor-hydrogen pairs via the bond information in the topology.\n",
    "\n",
    "Accetptor and hydrogen atoms are guessed based on the mass and partial charge of the atoms.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "hbonds = HydrogenBondAnalysis(universe=u)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8c64f05e03d34d7d9f20aca0deb4e3d8",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/98 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<MDAnalysis.analysis.hydrogenbonds.hbond_analysis.HydrogenBondAnalysis at 0x7fce51309190>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hbonds.run(verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Use `guess_acceptors` and `guess_hydrogens` to create atom selections\n",
    "\n",
    "It is also possible to use generated `hydrogens_sel` and `acceptors_sel` via the `guess_hydrogens` and `guess_acceptors` class methods.\n",
    "\n",
    "This selection strings may then be modified prior to calling `run`, or a subset of the universe may be used to guess the atoms. For example, below we will find hydrogens and acceptors belonging to specific residues.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "hbonds = HydrogenBondAnalysis(universe=u)\n",
    "hbonds.hydrogens_sel = hbonds.guess_hydrogens(\"resname ARG HIS LYS\")\n",
    "hbonds.acceptors_sel = hbonds.guess_acceptors(\"resname ARG HIS LYS\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7584a0c1a6c4410cae65b29eaba86553",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/98 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/pbarletta/mambaforge/envs/mda-user-guide/lib/python3.9/site-packages/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py:751: UserWarning: No hydrogen bonds were found given angle of 150 between Donor, None, and Acceptor, (resname ARG and name NE) or (resname ARG and name NH1) or (resname ARG and name NH2) or (resname ARG and name O) or (resname LYS and name O).\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<MDAnalysis.analysis.hydrogenbonds.hbond_analysis.HydrogenBondAnalysis at 0x7fce50983790>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hbonds.run(verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can check which atoms were used in the analysis:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "hydrogen_sel = (resname ARG and name HE) or (resname ARG and name HH11) or (resname ARG and name HH12) or (resname ARG and name HH21) or (resname ARG and name HH22) or (resname ARG and name HN) or (resname LYS and name HN) or (resname LYS and name HZ1) or (resname LYS and name HZ2) or (resname LYS and name HZ3)\n",
      "acceptors_sel = (resname ARG and name NE) or (resname ARG and name NH1) or (resname ARG and name NH2) or (resname ARG and name O) or (resname LYS and name O)\n"
     ]
    }
   ],
   "source": [
    "print(f\"hydrogen_sel = {hbonds.hydrogens_sel}\")\n",
    "print(f\"acceptors_sel = {hbonds.acceptors_sel}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## More advanced selections\n",
    "\n",
    "Slightly more complex selection strings are also possible. For example, to find hydrogen bonds involving the protein and water within 10 Å of the protein:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "hbonds = HydrogenBondAnalysis(universe=u)\n",
    "\n",
    "protein_hydrogens_sel = hbonds.guess_hydrogens(\"protein\")\n",
    "protein_acceptors_sel = hbonds.guess_acceptors(\"protein\")\n",
    "\n",
    "water_hydrogens_sel = \"resname TIP3 and name H1 H2\"\n",
    "water_acceptors_sel = \"resname TIP3 and name OH2\"\n",
    "\n",
    "hbonds.hydrogens_sel = f\"({protein_hydrogens_sel}) or ({water_hydrogens_sel} and around 10 not resname TIP3)\"\n",
    "hbonds.acceptors_sel = f\"({protein_acceptors_sel}) or ({water_acceptors_sel} and around 10 not resname TIP3)\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "84362a6b5026457ca6811f69b26ef9e5",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/98 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<MDAnalysis.analysis.hydrogenbonds.hbond_analysis.HydrogenBondAnalysis at 0x7fce5098a250>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hbonds.run(verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "    \n",
    "**Note**\n",
    "    \n",
    "The Universe we are analysing has the water removed. The above example is for illustrative purposes only.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hydrogen bonds between specific groups\n",
    "\n",
    "To calculate the hydrogen bonds between different groups, for example protein and water, one can use the `between` keyword. Below we will find protein-water and protein-protein hydrogen bonds, but not water-water hydrogen bonds. We do this by passing atom selection for the groups we wish to find hydrogen bonds between."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "hbonds = HydrogenBondAnalysis(\n",
    "    universe=u,\n",
    "    between=[\n",
    "        [\"protein\", \"resname TIP3\"], # for protein-water hbonds\n",
    "        [\"protein\", \"protein\"]       # for protein-protein hbonds\n",
    "    ]\n",
    ")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d910e79f17bc49959f49e2d031ed6b85",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/98 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<MDAnalysis.analysis.hydrogenbonds.hbond_analysis.HydrogenBondAnalysis at 0x7fce5096fd60>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hbonds.run(verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "    \n",
    "**Note**\n",
    "    \n",
    "The Universe we are analysing has the water removed. The above example is for illustrative purposes only.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "[1] 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",
    "[2] 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",
    "[3] Paul Smith, Robert&nbsp;M. Ziolek, Elena Gazzarrini, Dylan&nbsp;M. Owen, and Christian&nbsp;D. Lorenz.\n",
    "On the interaction of hyaluronic acid with synovial fluid lipid membranes.\n",
    "<em>Phys. Chem. Chem. Phys.</em>, 21(19):9845-9857, 2018.\n",
    "URL: <a href=\"http://dx.doi.org/10.1039/C9CP01532A\">http://dx.doi.org/10.1039/C9CP01532A\n",
    "\n",
    "[4] 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"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Raw Cell Format",
  "kernelspec": {
   "display_name": "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 | packaged by conda-forge | (main, Nov 22 2022, 15:55:03) \n[GCC 10.4.0]"
  },
  "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": "5edc5d8d8cbc0935a054a8e44024f729bc376180aae27775d15f2ff38c68f892"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}