notebooks/Bonding and J-coupling network output example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting bonding and NMR J-coupling networks in IPython\n",
    "======================================================="
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's import the `magres.atoms` library and use it to load a [CASTEP](http://www.castep.org) calculation output. This calculation was performed on an ethanol molecule between every nucleus with the following input files (`jcoupling_site` set to `C 1` as an example):\n",
    "\n",
    "### ethanol.cell\n",
    "\n",
    "```\n",
    "%block POSITIONS_ABS\n",
    "ang\n",
    "H -0.869341 -0.175289 -0.649057\n",
    "H 0.183454 -0.923201 0.560623\n",
    "H 0.869130 0.198626 -0.628783\n",
    "H -1.129705 0.975874 1.565773\n",
    "H -0.437769 2.079941 0.372865\n",
    "H 1.061671 0.678653 2.298066\n",
    "C -0.003000 -0.003000 -0.003000\n",
    "C -0.246915 1.165107 0.938396\n",
    "O 0.896314 1.459074 1.742964\n",
    "%endblock POSITIONS_ABS\n",
    "\n",
    "%block SPECIES_POT\n",
    "H H_POT.ASC.DAT\n",
    "C C_POT.ASC.DAT\n",
    "O O_POT.ASC.DAT\n",
    "%endblock SPECIES_POT\n",
    "\n",
    "%block LATTICE_CART\n",
    "ang\n",
    "6.000000 0.000000 0.000000\n",
    "0.000000 6.000000 0.000000\n",
    "0.000000 0.000000 6.000000\n",
    "%endblock LATTICE_CART\n",
    "\n",
    "jcoupling_site: C 1\n",
    "kpoint_mp_grid 1 1 1\n",
    "kpoint_mp_offset 0.25 0.25 0.25\n",
    "```\n",
    "\n",
    "### ethanol.param\n",
    "\n",
    "```\n",
    "opt_strategy   : speed\n",
    "task           : magres\n",
    "magres_task    : jcoupling\n",
    "cut_off_energy : 80 ry \n",
    "xc_functional  : PBE\n",
    "popn_calculate : true\n",
    "spin_polarized : false\n",
    "max_scf_cycles : 60\n",
    "nextra_bands   : 20\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "\n",
    "from magres.atoms import MagresAtoms, MagresAtomsView\n",
    "ethanol_atoms = MagresAtoms.load_magres('../samples/ethanol-all.magres')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "IPython's rich representations means that if you just put a set of atoms to the IPython notebook console, the `magres.atoms` library will automatically produce an image of the structure. This is generated using `pydot`, which is a Python interface to the graph (as in graph theory) plotting program `graphviz`.\n",
    "\n",
    "If present, it will also plot the (unreferenced) magnetic shielding below the atom label, and plot the largest indirect nuclear spin-spin interactions (J-coupling) as lines between nuclei. Red lines mean a positive sign on the reduced isotropic coupling and blue lines a negative one."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's also load in bonding information from a .castep file. This is a bit of a hack right now."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "castep_file = open('../samples/ethanol.castep').read()\n",
    "\n",
    "ethanol_atoms.calculate_bonds()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is done by the `magres.atoms.MagresAtomsView` class implementing a `_repr_png_` method, returning raw png data.\n",
    "\n",
    "You can also show just a subset of atoms, e.g. all atoms bonded to `C1` plus `C1` itself."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "unhashable type: 'MagresAtom'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-3-7364da3b33a7>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0;34m(\u001b[0m\u001b[0methanol_atoms\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mC1\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0methanol_atoms\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mC1\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbonded\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m/Users/christian/Dropbox/PycharmProjects/magres-format/magres/atoms.py\u001b[0m in \u001b[0;36m__radd__\u001b[0;34m(self, b)\u001b[0m\n\u001b[1;32m    327\u001b[0m   \u001b[0;32mdef\u001b[0m \u001b[0m__radd__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    328\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mtype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mb\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0mMagresAtom\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 329\u001b[0;31m       \u001b[0mnew_atoms\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mset\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0matoms\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mb\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    330\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    331\u001b[0m       \u001b[0;32mreturn\u001b[0m \u001b[0mMagresAtomsView\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnew_atoms\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlattice\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: unhashable type: 'MagresAtom'"
     ]
    }
   ],
   "source": [
    "(ethanol_atoms.C1 + ethanol_atoms.C1.bonded)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You might use this to investigate a sub-section of a system.\n",
    "\n",
    "Here's the J-coupling network just between hydrogen atoms."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alanine system\n",
    "==============\n",
    "\n",
    "Let's try another, larger system. From [the Wikipedia page on Alanine](http://en.wikipedia.org/wiki/Alanine):\n",
    "\n",
    "> Alanine (abbreviated as Ala or A) is an α-amino acid with the chemical formula CH3CH(NH2)COOH. The L-isomer is one of the 20 amino acids encoded by the genetic code. Its codons are GCU, GCC, GCA, and GCG. It is classified as a nonpolar amino acid. L-Alanine is second only to leucine in rate of occurrence, accounting for 7.8% of the primary structure in a sample of 1,150 proteins. D-Alanine occurs in bacterial cell walls and in some peptide antibiotics.\n",
    "\n",
    "Here we look at its solid-state crystalline form.\n",
    "\n",
    "Let's load our calculation. This only has J-coupling calculations between some of the nuclei in the system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "alanine_atoms = MagresAtoms.load_magres('../samples/alanine-jc-all.magres')\n",
    "\n",
    "#castep_file = open('../samples/alanine.castep').read()\n",
    "\n",
    "alanine_atoms.calculate_bonds()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, all we need to do print the bonding network is output the atoms object to the IPython console."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<magres.atom.MagresAtoms - 52 atoms>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "alanine_atoms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can already see a few interesting features, such as the H-bonding mediated J-coupling between nitrogen and oxygen molecules, e.g. N2 to O1.\n",
    "\n",
    "Let's look at that more closely. We can select all atoms within 2 Angstroms of N1, N2, O1 and O3 (the two H-bridges) by doing a `.within` query and adding together the sets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "unhashable type: 'MagresAtom'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-6-e6894e7d24fa>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0malanine_atoms\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mN1\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbonded\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0malanine_atoms\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mN2\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbonded\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0malanine_atoms\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mO1\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbonded\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0malanine_atoms\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mO3\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbonded\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0malanine_atoms\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mN1\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0malanine_atoms\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mN2\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0malanine_atoms\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mO1\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0malanine_atoms\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mO3\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m/Users/christian/Dropbox/PycharmProjects/magres-format/magres/atoms.py\u001b[0m in \u001b[0;36m__add__\u001b[0;34m(self, b)\u001b[0m\n\u001b[1;32m    321\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    322\u001b[0m     \u001b[0;32melif\u001b[0m \u001b[0mtype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mb\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0mMagresAtom\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 323\u001b[0;31m       \u001b[0mnew_atoms\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mset\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0matoms\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mb\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    324\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    325\u001b[0m       \u001b[0;32mreturn\u001b[0m \u001b[0mMagresAtomsView\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnew_atoms\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlattice\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: unhashable type: 'MagresAtom'"
     ]
    }
   ],
   "source": [
    "alanine_atoms.N1.bonded + alanine_atoms.N2.bonded + \\\n",
    "alanine_atoms.O1.bonded + alanine_atoms.O3.bonded + \\\n",
    "alanine_atoms.N1 + alanine_atoms.N2 + alanine_atoms.O1 + alanine_atoms.O3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's recurse along the the bonding network to pluck out an individual alanine molecule."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_molecule(atom):\n",
    "    atoms = set()\n",
    "\n",
    "    def _get_molecule(atom1):\n",
    "        for atom2 in atom1.bonded:\n",
    "            if atom2 not in atoms:\n",
    "                atoms.add(atom2)\n",
    "                _get_molecule(atom2)\n",
    "                \n",
    "    _get_molecule(atom)\n",
    "    \n",
    "    return MagresAtomsView(list(atoms), atom.bonded.lattice)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<magres.atom.MagresAtomsView - 164 atoms>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "get_molecule(alanine_atoms.N1)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}