08_huckel.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "<figure>\n",
    "  <IMG SRC=\"gfx/Logo_norsk_pos.png\" WIDTH=100 ALIGN=\"right\">\n",
    "</figure>\n",
    "\n",
    "# H&uuml;ckel Method\n",
    "*Roberto Di Remigio*, *Luca Frediani*\n",
    "\n",
    "With the [H&uuml;ckel method] we aim at calculating the energy levels of conjugated $\\pi$ systems.\n",
    "The method is limited to planar, conjugated hydrocarbons and includes only the $\\pi$ electrons in the system in the quantum mechanical treatment. The $\\sigma$ electrons are ignored, as they are at much lower energies and will not affect molecular properties significantly.\n",
    "The molecular orbitals (MOs) are built as a linear combination of atomic orbitals (LCAO),\n",
    "where the atomic orbitals (AOs) in question are the $p_z$ orbitals on each of the carbon atoms.\n",
    "\n",
    "[H&uuml;ckel method]: https://en.wikipedia.org/wiki/H%C3%BCckel_method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "\n",
    "## Variation principle and secular equations\n",
    "\n",
    "If $N$ is the number of carbon atoms in the conjugated hydrocarbon, the trial wavefunction is thus written as:\n",
    "\\begin{equation}\n",
    "\\Psi = \\sum_{i=1}^{N}c_i\\chi_i\n",
    "\\end{equation}\n",
    "where the index $i$ runs over carbon atoms and $\\chi_i$ is the $p_z$ orbital on the $i$-th atom in the molecule.\n",
    "The energy, i.e. the expectation value of the many-electron Hamiltonian, is:\n",
    "\\begin{equation}\n",
    "E(\\mathbf{c}^\\dagger, \\mathbf{c}) = \\frac{\\langle \\Psi | H | \\Psi \\rangle}{\\langle \\Psi | \\Psi \\rangle}\n",
    "= \\frac{\\sum_{i,j=1}^{N} c_i^*H_{ij}c_j }{\\sum_{i, j=1}^{N}c_i^*S_{ij}c_j}\n",
    "= \\frac{\\mathbf{c}^\\dagger\\mathbf{H}\\mathbf{c}}{\\mathbf{c}^\\dagger\\mathbf{S}\\mathbf{c}}\n",
    "\\end{equation}\n",
    "where the Hamiltonian and overlap matrices were introduced:\n",
    "\\begin{alignat}{2}\n",
    "H_{ij} = \\langle \\chi_i | H |\\chi_j\\rangle \\quad& S_{ij} = \\langle \\chi_i | \\chi_j \\rangle\n",
    "\\end{alignat}\n",
    "The variation principle guarantees that minimization of the energy function will lead to the optimal value of the coefficients, in the chosen LCAO basis.\n",
    "This reduces to the generalized eigenvalue problem:\n",
    "\\begin{equation}\n",
    "\\mathbf{H}\\mathbf{c} = \\mathbf{S}\\mathbf{c}E\n",
    "\\end{equation}\n",
    "where the **eigenvector** $\\mathbf{c}$ is the set of coefficients in the linear combination describing the MO with energy **eigenvalue** $E$.\n",
    "Finding the eigenvalues can be achieved by calculating the determinant of $\\mathbf{H} - E\\mathbf{S}$ and setting it to zero:\n",
    "\\begin{equation}\n",
    "\\det(\\mathbf{H} - E\\mathbf{S}) = 0\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## The eigenvalue problem in H&uuml;ckel method\n",
    "\n",
    "As already mentions, the molecular orbitals (MOs) are built as a linear combination of atomic orbitals (LCAO),\n",
    "where the atomic orbitals (AOs) in question are the $p_z$ orbitals on each of the carbon atoms.\n",
    "A schematic view of the model is given in the Figure below.\n",
    "\n",
    "<figure>\n",
    "  <IMG SRC=\"gfx/huckel.png\">\n",
    "</figure>\n",
    "\n",
    "The Hamiltonian and overlap matrices are constructed from a simple set of rules:\n",
    "\n",
    "* all overlap integrals between distinct orbitals are zero:\n",
    "\\begin{equation}\n",
    "S_{ii} = 1 \\quad S_{ij} = 0 \\,\\,\\,\\text{for}\\, i\\neq j\n",
    "\\end{equation}\n",
    "* all diagonal Hamiltonian matrix elements are equal:\n",
    "\\begin{equation}\n",
    "H_{ii} = \\alpha \\,\\,\\,\\forall i\n",
    "\\end{equation}\n",
    "where $\\alpha$ is a parameter of the model. It is usually set to the energy of an isolated electron in a $2p_z$ orbital. $\\alpha$ is thus a **negative number**.\n",
    "* all off-diagonal Hamiltonian matrix elements are equal. They are nonzero for **adjacent** carbon atoms and zero otherwise:\n",
    "\\begin{equation}\n",
    "H_{i, i+1} = \\beta \\quad H_{i, i+n} = 0 \\,\\,\\,\\text{for}\\,n\\geq 2\n",
    "\\end{equation}\n",
    "where $\\beta$ is a **negative** parameter that can be determined from experiment ($\\beta = -75 \\,\\mathrm{kJ}\\cdot\\mathrm{mol}^{-1}$)\n",
    "\n",
    "Once the connectivity and the values of the $\\alpha$ and $\\beta$ parameters are known one can set up the H&uuml;ckel matrix and the correspoding secular problem.\n",
    "\n",
    "### A simple case: ethene\n",
    "\n",
    "Ethene has two carbon atoms. The H&uuml;ckel matrix is:\n",
    "\\begin{equation}\n",
    "\\begin{pmatrix}\n",
    "\\alpha & \\beta \\\\\n",
    "\\beta  & \\alpha\n",
    "\\end{pmatrix}\n",
    "\\end{equation}\n",
    "and its secular problem:\n",
    "\\begin{equation}\n",
    "\\det \\begin{pmatrix}\n",
    "(\\alpha - E) & \\beta \\\\\n",
    "\\beta  & (\\alpha - E)\n",
    "\\end{pmatrix}\n",
    "= (\\alpha - E)^2 - \\beta^2 = 0\n",
    "\\end{equation}\n",
    "The solution of the second order equation gives the energies for the highest occupied molecular orbital $E_-$ (HOMO) and the lowest unoccupied molecular orbital $E_+$ (LUMO):\n",
    "\\begin{equation}\n",
    "E_- = \\alpha -\\beta \\quad E_+ = \\alpha + \\beta\n",
    "\\end{equation}\n",
    "Inserting these energies in the eigenvalue equation will lead to two different systems of equations whose solution determines the coefficient in the MOs linear combinations.\n",
    "Notice that the coefficients are determined up to normalization, meaning that only by imposing a normalization condition the solution to the eigenvalue problem will be unique."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Exercise: a small script for H&uuml;ckel method\n",
    "\n",
    "Write a Python function to calculate eigenvalues and eigenvectors for the H&uuml;ckel matrix of a generic conjugated hydrocarbon. \n",
    "The function should take four parameters:\n",
    "* the number of carbon atoms in the hydrocarbon\n",
    "* the value of the $\\alpha$ parameter\n",
    "* the value of the $\\beta$ parameter\n",
    "* an array specifying the connectivity of the carbon atoms\n",
    "\n",
    "The connectivity array describes which carbon atoms are connected _via_ a double bond. This is essential, since it determines which entries in the H&uuml;ckel matrix are nonzero!\n",
    "As an example, for the ethene molecule $[1, 2]$ would be the connectivity array. For the butadiene molecule one would instead have $[1, 2, 2, 3, 3, 4]$ which means that atom $1$ is connected to atom $2$ by a double bond, atom $2$ and $3$ are connected by a single bond and atom $3$ and $4$ are again connected by a double bond.\n",
    "Why so? Because butadiene has 4 carbon atoms, hence only two double bonds can be present. Since these have to be _alternating_ carbon $2$ and carbon $3$ must be connected by a single bond.\n",
    "To help with determining the connectivity array, draw the molecule and number all carbon atoms.\n",
    "The function will, of course, use NumPy for manipulating matrices and vectors.\n",
    "The module `linalg` includes the functions necessary for linear algebra manipulations. In our case, we're interested in finding eigenvalues and eigenvectors of the H&uuml;ckel matrix with the [`eig` function]\n",
    "\n",
    "```Python\n",
    "def huckel(n, alpha, beta, connectivity):\n",
    "    \"\"\" Huckel method\n",
    "    \n",
    "    n -- number of carbon atoms\n",
    "    alpha -- Coulomb integral\n",
    "    beta --  resonance integral\n",
    "    connectivity -- an array of integers defining the connectivity\n",
    "    \"\"\"\n",
    "    import numpy as np\n",
    "    from numpy import linalg as LA\n",
    "    \n",
    "    H = np.zeros((n, n))\n",
    "    # Fill the matrix!\n",
    "    ## Fill diagonal\n",
    "    \n",
    "    ## Fill off-diagonal\n",
    "    \n",
    "    # Print H\n",
    "    print('Huckel matrix is {}'.format(H))\n",
    "    \n",
    "    # Diagonalize H\n",
    "    energies, orbitals = LA.eig(H)\n",
    "    # Sort eigenvalues and eigenvectors in ascending order\n",
    "    sort_perm = energies.argsort()\n",
    "    energies.sort()     # <-- This sorts the list in place.\n",
    "    orbitals = evecs[:, sort_perm]\n",
    "    \n",
    "    # Return energies and orbitals\n",
    "    return energies, orbitals\n",
    "```\n",
    "\n",
    "[`eig` function]: http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.eig.html#numpy.linalg.eig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Huckel matrix is \n",
      "[[ -5. -75.   0.   0.]\n",
      " [-75.  -5. -75.   0.]\n",
      " [  0. -75.  -5. -75.]\n",
      " [  0.   0. -75.  -5.]]\n",
      "Huckel model energies \n",
      "[-126.35254916  -51.35254916   41.35254916  116.35254916]\n",
      "Huckel model orbitals \n",
      "[[ 0.37174803 -0.60150096  0.60150096 -0.37174803]\n",
      " [ 0.60150096 -0.37174803 -0.37174803  0.60150096]\n",
      " [ 0.60150096  0.37174803 -0.37174803 -0.60150096]\n",
      " [ 0.37174803  0.60150096  0.60150096  0.37174803]]\n"
     ]
    }
   ],
   "source": [
    "def huckel(n, alpha, beta, connectivity):\n",
    "    \"\"\" Huckel method\n",
    "    \n",
    "    n -- number of carbon atoms\n",
    "    alpha -- Coulomb integral\n",
    "    beta --  resonance integral\n",
    "    connectivity -- an array of integers defining the connectivity\n",
    "    \"\"\"\n",
    "    import numpy as np\n",
    "    from numpy import linalg as LA\n",
    "    \n",
    "    H = np.zeros((n, n))\n",
    "    # Fill the matrix!\n",
    "    ## Fill diagonal\n",
    "    \n",
    "    ## Fill off-diagonal\n",
    "    \n",
    "    # Print H\n",
    "    print('Huckel matrix is \\n{}'.format(H))\n",
    "    \n",
    "    # Diagonalize H\n",
    "    energies, evecs = LA.eig(H)\n",
    "    # Sort eigenvalues and eigenvectors in ascending order:q\n",
    "    sort_perm = energies.argsort()\n",
    "    energies.sort()     # <-- This sorts the list in place.\n",
    "    orbitals = evecs[:, sort_perm]\n",
    "    \n",
    "    # Return energies and orbitals\n",
    "    return energies, orbitals\n",
    "    \n",
    "energies, orbitals = huckel(4, -5.0, -75.0, [1, 2, 2, 3, 3, 4])\n",
    "print('Huckel model energies \\n{}'.format(energies))\n",
    "print('Huckel model orbitals \\n{}'.format(orbitals))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Visualizing the molecular orbitals\n",
    "\n",
    "We will use the packages [`chemview`] and [`chemlab`] to visualize the molecular orbitals, as in [this example]. The coefficients are calculated from our small function for the H&uuml;ckel method. The technicalities of the definition of the basis functions are not important in this context. You can find more details [here].\n",
    "\n",
    "[`chemview`]: http://chemview.readthedocs.org/en/latest/\n",
    "[`chemlab`]: http://chemlab.readthedocs.org/en/latest/index.html\n",
    "[this example]: http://gabrielelanaro.github.io/blog/2015/01/14/visualizing-molecular-orbitals-ipython-notebook.html\n",
    "[here]: https://en.wikipedia.org/wiki/Basis_set_(chemistry)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "HDF backend is not available because pytables is not installed\n",
      "cclib not found. Install cclib for more handlers.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/roberto/.virtualenvs/py3/lib/python3.4/site-packages/traitlets/traitlets.py:565: FutureWarning: comparison to `None` will result in an elementwise object comparison in the future.\n",
      "  silent = bool(old_value == new_value)\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "from chemlab.db.cirdb import CirDB\n",
    "from chemlab.qc import molecular_orbital\n",
    "\n",
    "from chemview import enable_notebook, MolecularViewer\n",
    "enable_notebook()\n",
    "\n",
    "# To get molecule by name\n",
    "ethene = CirDB().get(\"molecule\", \"ethene\")\n",
    "# Huckel's method results for ethene\n",
    "orbitals = np.array([[ 0.70710678,  0.70710678],\n",
    " [-0.70710678,  0.70710678]])\n",
    "\n",
    "# Initialize molecular viewer\n",
    "mv = MolecularViewer(ethene.r_array, { 'atom_types': ethene.type_array,\n",
    "                                         'bonds': ethene.bonds })\n",
    "# Define rendering type\n",
    "mv.wireframe()\n",
    "\n",
    "# Define basis set\n",
    "# p-functions for carbon from STO-3G basis set\n",
    "sto3g = [ \n",
    "    # atom 1\n",
    "[\n",
    "('P', [\n",
    "(2.941249, 0.155916), (0.683483, 0.607684), (0.222290, 0.391957), ]),\n",
    "],\n",
    "    # atom 2\n",
    "[\n",
    "('P', [\n",
    "(2.941249, 0.155916), (0.683483, 0.607684), (0.222290, 0.391957), ]),\n",
    "]\n",
    "]\n",
    "\n",
    "# Define molecular orbital to be plotted \n",
    "# Notice that we pass only the coordinates of the carbon atoms!\n",
    "coords = np.array([ethene.r_array[0], ethene.r_array[1]])\n",
    "# We need to fudge the coefficient arrays with zeros (p_x and p_y are absent!)\n",
    "HOMO_coeff = np.array([0.0, 0.0, orbitals[0][0], 0.0, 0.0, orbitals[0][1]])\n",
    "LUMO_coeff = np.array([0.0, 0.0, orbitals[1][0], 0.0, 0.0, orbitals[1][1]])\n",
    "HOMO = molecular_orbital(coords, HOMO_coeff, sto3g)\n",
    "LUMO = molecular_orbital(coords, LUMO_coeff, sto3g)\n",
    "# Add isosurface of MO with positive isovalue\n",
    "mv.add_isosurface(HOMO, isolevel=0.3, color=0xff0000)\n",
    "# Add isosurface of MO with negative isovalue\n",
    "mv.add_isosurface(HOMO, isolevel=-0.3, color=0x0000ff)\n",
    "# This is necessary to actually get a picture of the molecule\n",
    "mv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.4.3+"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}