{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# NCI with a Jupyter notebook #\n",
    "_(Peter Reinhardt, Paris 2020)_\n",
    "\n",
    "We approach the NCI with the simplest possible molecule, H$_2$ in a minimal basis. The orbital is a Slater function $N_a e^{-a r}$ with $a=1$, and thus $N=\\sqrt{a^3/\\pi}$.\n",
    "\n",
    "In the notebook, we \n",
    "1. define an atomic orbital in 3D space, centered  at (R,0,0)\n",
    "2. define a molecular orbital for H$_2$. The norm of the orbital can be approximated by a Gaussian function $a\\, e^{-b R^2} + c\\, e^{-d R^2}$\n",
    "3. calculate grad rho / rho^(4/3) on the line going through the two atoms\n",
    "4. upon varying R, we observe the evolution of this index"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**The notebook is not complete - please fill in the missing parts !!**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "def atomic_orb(x,y,z,R):\n",
    "    import math\n",
    "    from math import exp \n",
    "    from math import sqrt\n",
    "    from math import pi\n",
    "    \n",
    "    r=math.sqrt((x-R)*(x-R)+y*y+z*z)\n",
    "    val=math.sqrt(1./math.pi)*math.exp(-r) \n",
    "    return val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def grad_atomic_orb(x,y,z,R):\n",
    "    import math\n",
    "    from math import exp \n",
    "    from math import sqrt\n",
    "    from math import pi\n",
    "    \n",
    "    denom=math.sqrt((x-R)*(x-R)+y*y+z*z)\n",
    "    if (abs(denom) > 1.e-8):\n",
    "        gradx=-atomic_orb(x,y,z,R)*(x-R)/denom\n",
    "        grady=-atomic_orb(x,y,z,R)*y/denom\n",
    "        gradz=-atomic_orb(x,y,z,R)*z/denom\n",
    "    else:\n",
    "        gradx=0.\n",
    "        grady=0.\n",
    "        gradz=0.\n",
    "        \n",
    "    return [gradx,grady,gradz]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def norm_molorb(R):\n",
    "    import math\n",
    "    from math import exp \n",
    "    from math import sqrt\n",
    "    from math import pi\n",
    "\n",
    "# this a fitting formula for the norm of the orbital\n",
    "# to avoid to recalculate the (non-trivial) overlap integral\n",
    "    afit=0.460515\n",
    "    bfit=0.940009\n",
    "    cfit=0.532886\n",
    "    dfit=0.270818\n",
    "    val=afit*math.exp(-bfit*R*R) + cfit*math.exp(-dfit*R*R)\n",
    "    val=2.+2.*val\n",
    "    val=1./math.sqrt(val)\n",
    "\n",
    "    return val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def molecular_orb(x,y,z,R):\n",
    "    val=norm_molorb(R)*(atomic_orbital(x,y,z,-R/2)+atomic_orbital(x,y,z,R/2))\n",
    "    return val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dens(x,y,z,R):\n",
    "    val=molecular_orb(x,y,z,R)*molecular_orb(x,y,z,R)\n",
    "    val=2.*val # we have 2 elecrons inside the orbital\n",
    "    return val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def grad_dens(x,y,z,R):\n",
    "    gradx=4.*molecular_orb(x,y,z,R)*(grad_atomic_orb(x,y,z,R/2)[0] \\\n",
    "                                     +grad_atomic_orb(x,y,z,-R/2)[0])\n",
    "    grady=4.*molecular_orb(x,y,z,R)*(grad_atomic_orb(x,y,z,R/2)[1] \\\n",
    "                                     +grad_atomic_orb(x,y,z,-R/2)[1])\n",
    "    gradz=4.*molecular_orb(x,y,z,R)*(grad_atomic_orb(x,y,z,R/2)[2] \\\n",
    "                                     +grad_atomic_orb(x,y,z,-R/2)[2])\n",
    "    return [gradx,grady,gradz]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def abs_grad_dens(x,y,z,R):\n",
    "    import math\n",
    "    from math import sqrt\n",
    "    # here we define the norm of the gradient |grad(rho(x,y,z)|\n",
    "    gradx,grady,gradz=grad_dens(x,y,z,R)\n",
    "    val=sqrt(gradx**2.+grady**2.+gradz**2.)\n",
    "    return val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def s(x,y,z,R):    \n",
    "    import math\n",
    "    from math import sqrt\n",
    "    # here we define the reduced gradient s(x,y,z) = |grad(rho(x,y,z))|/rho(x,y,z)^(4/3)\n",
    "    val=abs_grad_dens(x,y,z,R)/dens(x,y,z,R)**(4./3.)\n",
    "    return val"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have everything defined for plotting the NCI for 2 hydrogen atoms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import sys\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# we have to create table for a fixed R with triples (x,dens,s) for y,z=0\n",
    "\n",
    "# in the plot we need x_table, s_table, rho_table for plotting \n",
    "# s aginst x, rho against x, and s against rho\n",
    "\n",
    "R=5.5\n",
    "\n",
    "# we plot from -2R to 2R in nmax+1 points\n",
    "# dx = 4 R/(nmax)\n",
    "# x0 = -2 R\n",
    "\n",
    "nmax=500\n",
    "dx=4*R/nmax\n",
    "x0=-2.*R\n",
    "\n",
    "\n",
    "# initialize the tables \n",
    "x_table=[]\n",
    "rho_table=[]\n",
    "s_table=[]\n",
    "x=x0\n",
    "\n",
    "# the loop\n",
    "for i in range(nmax):\n",
    "    x_table.append(x)\n",
    "    s_table.a\n",
    "\n",
    "\n",
    "\n",
    "# The plots    \n",
    "fig, (ax1, ax2, ax3) = plt.subplots(1, 3,figsize=(18,6))\n",
    "fig.suptitle('rho against x   --  s against x   --       s against rho')\n",
    "ax1.plot(x_table, rho_table)\n",
    "ax2.set_ylim(0,50)\n",
    "ax2.plot(x_table,s_table)\n",
    "ax3.set_ylim(0,50)\n",
    "ax3.plot(rho_table,s_table)"
   ]
  }
 ],
 "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}