{
  "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 10**\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 3**\n",
        "\n",
        "PySCF Documentation: https://pyscf.org\n",
        "\n",
        "---"
      ],
      "metadata": {
        "id": "8crRlPgWWMi6"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "**Application of QM Calculus**\n",
        "\n",
        "\\\\\n",
        "\n",
        "**Exercise 1**\n",
        "\n",
        "We saw in the theoretical class that the equation\n",
        "\n",
        "\\begin{equation}\n",
        "f(\\vec{r}_1) \\psi_i(\\vec{r}_1) = \\varepsilon_i \\psi_i(\\vec{r}_1)\n",
        "\\end{equation}\n",
        "\n",
        "can be written in terms of the expansion of known basis functions\n",
        "\n",
        "\\begin{equation}\n",
        "\\psi_i(\\vec{r}) = \\sum_{u=1}^K C_{ui} \\phi_u(\\vec{r})\n",
        "\\end{equation}\n",
        "\n",
        "yielding\n",
        "\n",
        "\\begin{equation}\n",
        "\\sum_{v=1}^K C_{vi} \\underbrace{\\int d(\\vec{r}_1) \\phi_u^{*}(\\vec{r}_1) f(\\vec{r}_1) \\phi_v(\\vec{r}_1)}_{F_{uv}} = \\varepsilon_i \\sum_{v=1}^K C_{vi} \\underbrace{\\int d(\\vec{r}_1) \\phi_u^{*}(\\vec{r}_1) \\phi_v(\\vec{r}_1)}_{S_{uv}}\n",
        "\\end{equation}\n",
        "\n",
        "Derive the same equation using Dirac notation\n",
        "\n",
        "\\\\\n",
        "\n",
        "**Exercise 2**\n",
        "\n",
        "The Hartree-Fock calculation restricted to the $H_2$ molecule in a minimal basis. with $R_{12}$ = 0.741 Å gives rise to the following matrices $\\mathbf{C}$ and $\\boldsymbol{\\varepsilon}$, solution of the Roothaan equations, in matrix form $\\mathbf{F} \\mathbf{C} = \\mathbf{S} \\mathbf{C} \\boldsymbol{\\varepsilon}$\n",
        "\n",
        "\\begin{equation}\n",
        "    \\textbf{C} =\n",
        "    \\begin{bmatrix}\n",
        "     0.54895 & 1.12113 \\\\\n",
        "     0.54895 & -1.12113 \\\\\n",
        "\\end{bmatrix}\\\\\n",
        "\\end{equation}\n",
        "\n",
        "\\begin{equation}\n",
        "   \\boldsymbol{\\varepsilon} =\n",
        "    \\begin{bmatrix}\n",
        "     −0.578 & 0 \\\\\n",
        "0 & +0.670 \\\\\n",
        "\\end{bmatrix}\\\\\n",
        "\\end{equation}\n",
        "\n",
        "\\\\\n",
        "a) Recall that\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",
        "Now write $\\psi_i$ in terms of the coefficients of the matrix $\\mathbf{C}$\n",
        "\n",
        "\\\\\n",
        "\n",
        "b) Recall that the total energy of the system is given by\n",
        "\n",
        "\\\\\n",
        "\\begin{equation}\n",
        "\\begin{aligned}\n",
        "E_{tot} = E(RHF) = E_{0} + U^{rep}_{nuc-nuc} = E_0 + \\sum_{A=1}^M \\sum_{B \\neq A}^M \\frac{Z_A Z_B}{R_{AB}}\n",
        "\\end{aligned}\n",
        "\\end{equation}\n",
        "\n",
        "\\\\\n",
        "knowing that $E_{tot} = E(RHF) = −1.1167$ Hartree, calculate the value of the nuclear repulsion energy and the electronic energy\n",
        "\n",
        "\\\\\n",
        "c) The elements of the density matrix are given by\n",
        "\n",
        "\\begin{equation}\n",
        "P_{uv} = 2 \\sum_{a=1}^{N/2} C_{ua} C^*_{va}\n",
        "\\end{equation}\n",
        "\n",
        "Based on the data provided, write the density matrix for the $H_2$ molecule in a minimal basis.\n",
        "\n",
        "\\\n",
        "d) Now use the relationship of the coefficients with the elements of the overlap matrix to write the overlap matrix $\\mathbf{S}$\n",
        "\n",
        "e) Alternatively, use the relationship of the coefficients with the elements of the density matrix $\\mathbf{P}$"
      ],
      "metadata": {
        "id": "CXoD_3r1u7QB"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "**Exercise 3**\n",
        "\n",
        "Now use the PySCF program to obtain the matrices $\\mathbf{P}$, $\\mathbf{C}$, $\\mathbf{S}$ as well as the energies of the orbitals (the eigenvalues ​​or matrix $\\boldsymbol{\\varepsilon}$) and the total energy at the HF/STO-3G level. Compare with the values ​​provided and with those calculated in the previous lines\n",
        "\n",
        "Note 1: Use the distance indicated in the problem data without optimizing\n",
        "\n",
        "Note 2: To extract the matrices from a calculation we can use\n",
        "\n",
        "```\n",
        "# Coefficient Matrix(C): each column is an MO!\n",
        "C = mf.mo_coeff # shape (nbasis, nbasis)\n",
        "\n",
        "# Orbital energies (diagonal epsilon matrix)\n",
        "# 1D vector with eigenvalues\n",
        "epsilon = mf.mo_energy\n",
        "\n",
        "# Density matrix (P)\n",
        "# In the RHF case: 2 electrons per occupied orbital\n",
        "# The dimension is (nbasis, nbasis)\n",
        "P = mf.make_rdm1()\n",
        "\n",
        "# Overlap matrix\n",
        "S = mol.intor('int1e_ovlp')\n",
        "```"
      ],
      "metadata": {
        "id": "DjbSliM5lTbG"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "!pip install pyscf\n",
        "!pip install geometric\n",
        "!pip install py3Dmol\n",
        "\n",
        "from pyscf import gto, scf\n",
        "from pyscf.geomopt.geometric_solver import optimize\n",
        "from scipy.spatial.distance import euclidean\n",
        "import numpy as np\n",
        "\n",
        "# Escreva aqui o seu código\n",
        "\n"
      ],
      "metadata": {
        "id": "ZIEevtteuqli"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "**Exercise 4**\n",
        "\n",
        "Repeat the previous exercise using 6-31g, 6-311g, 6-311g** and cc-pVDZ, cc-pVTZ and cc-pVQZ (no need to copy the obtained values manualy).\n",
        "Check the dimensions of the matrices and justify them. How does the total energy vary?"
      ],
      "metadata": {
        "id": "n5U01KrcaKXO"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "!pip install pyscf\n",
        "!pip install geometric\n",
        "!pip install py3Dmol\n",
        "\n",
        "from pyscf import gto, scf\n",
        "from pyscf.geomopt.geometric_solver import optimize\n",
        "from scipy.spatial.distance import euclidean\n",
        "import numpy as np\n",
        "\n",
        "# Escreva aqui o seu código\n",
        "\n"
      ],
      "metadata": {
        "id": "tJbAQuZJb8LX"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "**Exercise 5**\n",
        "\n",
        "The integral of the total charge density is equal to the number of electrons $N$\n",
        "\n",
        "\\begin{equation}\n",
        "2 \\sum_{a=1}^{N/2} \\int d(\\vec{r}) |\\psi_a(\\vec{r})|^2 = N\n",
        "\\tag{1}\n",
        "\\end{equation}\n",
        "\n",
        "We also saw that we can express the total charge density using the density matrix $P_{\\mu\\nu}$\n",
        "\n",
        "\\begin{equation}\n",
        "\\rho(\\vec{r}) = 2 \\sum_{a=1}^{N/2} |\\psi_a(\\vec{r})|^2 = 2 \\sum_{\\mu\\nu} P_{\\mu\\nu} \\phi_\\mu(\\vec{r}) \\phi_\\nu(\\vec{r})\n",
        "\\tag{2}\n",
        "\\end{equation}\n",
        "\n",
        "Write the integral of the total charge density written in the form of equation 2 and verify the appearance of the overlap integral $\\mathbf{S}$\n",
        "\n",
        "**Exercise 6**\n",
        "\n",
        "There is no single definition for quantifying the number of electrons associated with a given atom or nucleus. However, it may be useful to perform a population analysis such as *Mulliken Population Analysis*, in which the electrons associated with a given atom in a molecule are obtained by **assigning half the density to each atom**. The population on atom A is then given by\n",
        "\n",
        "\\\\\n",
        "\\begin{equation}\n",
        "N_A = \\sum_{\\mu \\in A} \\sum_{\\nu} P_{\\mu\\nu} S_{\\mu\\nu}\n",
        "\\end{equation}\n",
        "\n",
        "and therefore, the Mulliken charge on atom A is:\n",
        "\n",
        "\\begin{equation}\n",
        "q_A^{\\text{Mulliken}} = Z_A - \\sum_{\\mu \\in A} \\sum_{\\nu} P_{\\mu\\nu} S_{\\mu\\nu} = Z_A - N_A\n",
        "\\end{equation}\n",
        "\n",
        "Based on the data from Exercise 2, calculate the Mulliken charge on each hydrogen atom. Is this the value you expected?\n",
        "\n",
        "Note: $\\sum_{\\mu \\in A}$ indicates the sum over all atomic orbitals centered on atom A; $\\sum_{\\nu}$ runs through all the atomic orbitals of the base, i.e., the complete set of the molecule.\n",
        "\n",
        "\n",
        "**Exercise 7**\n",
        "\n",
        "Building on exercise 3, now use the PySCF program to calculate the Mulliken charges of hydrogen atoms at the HF/STO-3G level. Do you expect the charges to vary with the basis functions?\n",
        "\n",
        "Note: In pySCF, we can perform a Mulliken population analysis using `mf.mulliken_pop()`"
      ],
      "metadata": {
        "id": "AtqAZhgVb0Ou"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "!pip install pyscf\n",
        "!pip install geometric\n",
        "!pip install py3Dmol\n",
        "\n",
        "from pyscf import gto, scf\n",
        "from pyscf.geomopt.geometric_solver import optimize\n",
        "\n",
        "from scipy.spatial.distance import euclidean\n",
        "import numpy as np\n",
        "\n",
        "# Escreva aqui o seu código"
      ],
      "metadata": {
        "id": "ksXQ8I2O06Vm"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "**Exercício 8**\n",
        "\n",
        "Usando o PySCF, optimize agora a geometria da molécula de ácido fluorídrico (HF) ao nível HF :-) usando uma base mínima (sto-3g). Obtenha as matrizes $\\mathbf{P}$, $\\mathbf{C}$, $\\mathbf{S}$ assim como as energias das orbitais (os valores próprios ou matriz $\\boldsymbol{\\varepsilon}$) e a energia total. Verifique a dimensão das matrizes na base mínima. Calcule as cargas em ambos os átomos usando uma análise de população de Mulliken. Os valores estão de acordo com o esperado? Comente o resultado à luz do pressuposto da análise de população de Mulliken onde a densidade é atribuída metade a cada átomo.  "
      ],
      "metadata": {
        "id": "QoZi97CMuyBp"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "!pip install pyscf\n",
        "!pip install geometric\n",
        "!pip install py3Dmol\n",
        "\n",
        "from pyscf import gto, scf\n",
        "from pyscf.geomopt.geometric_solver import optimize\n",
        "\n",
        "from scipy.spatial.distance import euclidean\n",
        "import numpy as np\n",
        "\n",
        "# Escreva aqui o seu código\n",
        "\n"
      ],
      "metadata": {
        "id": "hwUpUCKPv0cn"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "**Exercise 9**\n",
        "\n",
        "Based on Exercise 3 of Lecture 9 regarding the water molecule, optimize the geometry at the HF level using the bases sto-3g, 6-31g, 6-311g, 6-311g** and cc-pVDZ, cc-pVTZ and cc-pVQZ. For each base and its optimized geometry, calculate the Mulliken charges on each atom. Comment on the results. You can use information from Exercise 3 of Lecture 9."
      ],
      "metadata": {
        "id": "VWWpPBkGSZl1"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "!pip install pyscf\n",
        "!pip install geometric\n",
        "!pip install py3Dmol\n",
        "\n",
        "from pyscf import gto, scf\n",
        "from pyscf.geomopt.geometric_solver import optimize\n",
        "\n",
        "from scipy.spatial.distance import euclidean\n",
        "import numpy as np\n",
        "\n",
        "# Escreva aqui o seu código\n",
        "\n",
        "\n"
      ],
      "metadata": {
        "id": "zruDZjROVrSJ"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}