{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "source": [
        "#**Computational Chemistry - Practical Class 7**\n",
        "\n",
        "**Minimal Basis MO-LCAO for the H$_2$ Molecule**\n",
        "\n",
        "\\\\\n",
        "\n",
        "---\n",
        "Supporting Bibliography:\n",
        "\n",
        "Attila Szabo and Neil S. Ostlund, Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory, Dover Publications Inc., New York, 1996, **Chapter 2, Sections 2.2.5 and 2.3.1**\n",
        "\n",
        "---\n",
        "\\\\\n",
        "**Exercise 1**\n",
        "\n",
        "Consider the following basis functions (radial part) of 1s atomic orbitals centered on hydrogen atoms\n",
        "\n",
        "\\\\\n",
        "$$\n",
        "GTO(r) = (\\frac{2\\alpha}{\\pi})^{3/4} e^{-\\alpha_n r^2} \\text{ with $\\alpha = 0.282942$}\n",
        "$$\n",
        "\n",
        "$$\n",
        "STO(r) = \\sqrt{\\frac{\\alpha^{3}}{\\pi}} e^{-\\alpha r} \\text{ with $\\alpha = 1.000000$}\n",
        "$$\n",
        "\n",
        "\\\\\n",
        "Represent from $-5~u.a.$ to $5~u.a.$ the functions GTO(r) and STO(r) centered on a hydrogen atom located at $r =~0~a.u.$\n",
        "\n",
        "Note 1: Given the radial nature of the function, $f(-r) = f(r)$, that is, the function decays exponentially regardless of direction. We can ensure this by using $f(x) = f(|x|)$\n",
        "\n",
        "Note 2: `np.abs()` calculates the absolute value"
      ],
      "metadata": {
        "id": "8crRlPgWWMi6"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "# Set the X values: 1000 points from -10 to 10\n",
        "x = np.linspace(-10, 10, 1000)\n",
        "\n",
        "# Write your code for STO here"
      ],
      "metadata": {
        "id": "jqPNdh-uYzrG"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "# Set the X values: 1000 points from -10 to 10\n",
        "x = np.linspace(-10, 10, 1000)\n",
        "\n",
        "# Write your code for GTO here"
      ],
      "metadata": {
        "id": "zRIvakliZFGD"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "**Exercise 2**\n",
        "\n",
        "Now repeat the exercise but, in addition to the functions centered at $r = 0.0~u.a.$, simultaneously represent the functions GTO(r) and STO(r) centered at another hydrogen atom located at $r =~0.5~a.u.$\n",
        "\n",
        "Note 1: Note that the function from the generic point of view can be represented as\n",
        "\n",
        "$$\n",
        "STO(r) = N e^{-\\alpha |r-R|}\n",
        "$$\n",
        "\n",
        "$$\n",
        "GTO(r) = N e^{-\\alpha_n |r-R|^2}\n",
        "$$\n",
        "\n",
        "with $R$ equal to the location where the function is centered"
      ],
      "metadata": {
        "id": "n5U01KrcaKXO"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "# Set the X values: 1000 points from -10 to 10\n",
        "x = np.linspace(-10, 10, 1000)\n",
        "\n",
        "# Write your code for STO here"
      ],
      "metadata": {
        "id": "Z1NgrJOxLNNJ"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "# Set the X values: 1000 points from -10 to 10\n",
        "x = np.linspace(-10, 10, 1000)\n",
        "\n",
        "# Write your code for GTO here"
      ],
      "metadata": {
        "id": "5Ngsh2aJaulk"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "**Exercise 4**\n",
        "\n",
        "Now represent\n",
        "\n",
        "\\\\\n",
        "$F(r) = \\phi_1^*(r)\\phi_2(r)$\n",
        "\n",
        "\\\\\n",
        "with $\\phi_{1,2} = STO(r)$ and $\\phi_{1,2} = GTO(r)$, where the index 1,2 corresponds to the functions centered at $0.0~u.a.$ and $0.5~u.a.$. Keep the representation of the original functions.\n",
        "\n",
        "\\\\\n",
        "If we integrate the function $F(r)$ with respect to the volume, what do we get?"
      ],
      "metadata": {
        "id": "NCwiBGm3bHjg"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy.integrate import simpson\n",
        "\n",
        "# Set the X values: 1000 points from -10 to 10\n",
        "x = np.linspace(-10, 10, 1000)\n",
        "\n",
        "# Write your code for STO here"
      ],
      "metadata": {
        "id": "nOMWuMlXlWx4"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy.integrate import simpson\n",
        "\n",
        "# Set the X values: 1000 points from -10 to 10\n",
        "x = np.linspace(-10, 10, 1000)\n",
        "\n",
        "# Write your code for GTO here"
      ],
      "metadata": {
        "id": "uTwKNNANb_Ql"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "**Exercise 5**\n",
        "\n",
        "We saw in the previous classes that when the 2 H atoms approach each other, from the two localized atomic orbitals, we form two delocalized spatial molecular orbitals as a linear combination of $\\phi_1$ and $\\phi_2$.\n",
        "\n",
        "\\\\\n",
        "\\begin{equation}\n",
        "\\begin{aligned}\n",
        "    \\psi_1 &= [2(1+S_{12})]^{-1/2}(\\phi_1 + \\phi_2) \\\\\n",
        "    \\psi_2 &= [2(1-S_{12})]^{-1/2}(\\phi_1 - \\phi_2)\n",
        "\\end{aligned}   \n",
        "\\end{equation}\n",
        "\n",
        "\\\\\n",
        "\n",
        "We also saw that $\\langle \\phi_1\\left|\\phi_1\\right\\rangle = 1$ but that $\\langle \\phi_1\\left|\\phi_2\\right\\rangle = S_{12}$ with\n",
        "\n",
        "$$\n",
        "S_{12} = \\int \\phi_1^*(\\vec{r})\\phi_2(\\vec{r})d(\\vec{r})\n",
        "$$\n",
        "\n",
        "Calculate $S_{12}$ approximately for both cases (STOs vs GTOs) **when the hydrogen atoms are separated** by $0.5~u.a.$, $1~u.a.$ and $5~u.a.$. Comment on the variation of $S_{12}$ with the separation distance $R$\n",
        "\n",
        "\\\\\n",
        "Note 1: There are several ways to numerically calculate the approximate integral of a function. Perhaps the best known is the **Trapezoidal Rule**. We can also use **Simpson's Rule** (more info at https://www.freecodecamp.org/portuguese/news/a-regra-de-simpson-a-formula-e-como-ela-funciona/). The two methods differ in the order of the function used to interpolate (linear vs quadratic)\n",
        "\n",
        "Note 2: There is an implementation of **Simpson's Rule** in the `scipy` library. The syntax is as follows:\n",
        "\n",
        "```\n",
        "# Import the library\n",
        "from scipy.integrate import simpson\n",
        "\n",
        "# define the function\n",
        "f = x**2 + 5 + 2*x\n",
        "\n",
        "# calculate the integral of y=f with respect to x=x\n",
        "integral=simpson(y=f, x=x)\n",
        "\n",
        "```"
      ],
      "metadata": {
        "id": "Mihdk7rvxUMR"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy.integrate import simpson\n",
        "\n",
        "# Set the X values: 1000 points from -10 to 10\n",
        "x = np.linspace(-10, 10, 1000)\n",
        "\n",
        "# Write your code for STO here"
      ],
      "metadata": {
        "id": "z0Sccx4mCbnq"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy.integrate import simpson\n",
        "\n",
        "# Set the X values: 1000 points from -10 to 10\n",
        "x = np.linspace(-10, 10, 1000)\n",
        "\n",
        "# Write your code for GTO here"
      ],
      "metadata": {
        "id": "Xy-252D02T3E"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "**Exercise 6**\n",
        "\n",
        "Without doing any calculations, indicate the value you would expect to obtain for $S_{12}$ when $R = 0.0~u.a.$. Justify the result using the definition of $S_{12}$\n",
        "\n",
        "\\\\\n",
        "**Exercise 7**\n",
        "\n",
        "Now calculate $S_{12}$ approximately (e.g. using Simpson's rule) for $R = 0.0~u.a.$. Is the result obtained what was expected? Why?\n",
        "\n",
        "Tip: always check the variables and integration limits of the functions! Also check under what conditions the normalization constants were obtained."
      ],
      "metadata": {
        "id": "5TK0RBXEInCe"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy.integrate import simpson\n",
        "\n",
        "# Set the X values: 1000 points from -10 to 10\n",
        "x = np.linspace(-10, 10, 1000)\n",
        "\n",
        "# Write your code for STO here"
      ],
      "metadata": {
        "id": "51K8b5g5Jt4t"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy.integrate import simpson\n",
        "\n",
        "# Set the X values: 1000 points from -10 to 10\n",
        "x = np.linspace(-10, 10, 1000)\n",
        "\n",
        "# Write your code for GTO hereO"
      ],
      "metadata": {
        "id": "L1DSOevZRRE7"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "**Exercise 8**\n",
        "\n",
        "Now remember that the PySCF program allows you to perform QM calculations using the Python language. Calculate the ground state energy of the H$_2$ molecule when the hydrogen atoms are spaced ($R$) by $0.5~u.a.$, $1~u.a.$, $1.5~u.a.$ and $5~u.a.$. Use sto-3g as the base function set.\n",
        "\n",
        "The relevant syntax modifications are:\n",
        "\n",
        "```\n",
        "# Define the molecule (or atoms)\n",
        "mol = gto.M(\n",
        "    atom='''\n",
        "        H coordinates x y z (in Å)\n",
        "        H coordinates x y z (in Å)\n",
        "    ''',               \n",
        "    basis='sto-3g', # Basis functions\n",
        "    spin=0, # Multiplicity: 2S\n",
        "    charge=0 # Charge\n",
        ")\n",
        "\n",
        "# Perform a restricted Hartree-Fock (RHF) calculation\n",
        "mf = scf.RHF(mol)\n",
        "\n",
        "```\n",
        "\n",
        "What does \"converted SCF energy\" mean? What is the most favorable distance? How could I find the optimal distance?"
      ],
      "metadata": {
        "id": "9oYorGtHMA5f"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Install the program\n",
        "!pip install pyscf\n",
        "\n",
        "# Import libraries\n",
        "import numpy as np\n",
        "from pyscf import gto, scf\n",
        "\n",
        "# your code here"
      ],
      "metadata": {
        "id": "451XS2cdK-JY"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}