{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "source": [
        "#**Química Computacional - Aula Prática 10**\n",
        "\n",
        "---\n",
        "Bibliografia de Suporte:\n",
        "\n",
        "Attila Szabo and Neil S. Ostlund, Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory, Dover Publications Inc., New York, 1996,  **Capítulo 3**\n",
        "\n",
        "Documentação PySCF: https://pyscf.org\n",
        "\n",
        "---"
      ],
      "metadata": {
        "id": "8crRlPgWWMi6"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "**Aplicação de Cálculos QM**\n",
        "\n",
        "\\\\\n",
        "\n",
        "**Exercício 1**\n",
        "\n",
        "Vimos na aula teórica que a equação\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",
        "pode ser escrita em função da expansão de funções base conhecidas\n",
        "\n",
        "\\begin{equation}\n",
        "\\psi_i(\\vec{r}) = \\sum_{u=1}^K C_{ui} \\phi_u(\\vec{r})\n",
        "\\end{equation}\n",
        "\n",
        "originando\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 a mesma equação usando a notação de Dirac\n",
        "\n",
        "\\\\\n",
        "\n",
        "**Exercício 2**\n",
        "\n",
        "O cálculo de Hartree-Fock restricto para a molécula de $H_2$ numa base mínima. com $R_{12}$ = 0.741 Å origina as seguintes matrizes $\\mathbf{C}$ e $\\boldsymbol{\\varepsilon}$, solução das equações de Roothaan, na forma matricial $\\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) Recorde que\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",
        "Escreva agora $\\psi_i$  em função dos coeficientes da matriz $\\mathbf{C}$\n",
        "\n",
        "\\\\\n",
        "\n",
        "b) Recorde que a energia total do sistema é dada por\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",
        "sabendo que $E_{tot} = E(RHF) = −1.1167$ Hartree, calcule o valor da energia de repulsão nuclear e a energia electrónica\n",
        "\n",
        "\\\\\n",
        "c) Os elementos da matriz densidade são dados por\n",
        "\n",
        "\\begin{equation}\n",
        " P_{uv} = 2 \\sum_{a=1}^{N/2} C_{ua} C^*_{va}\n",
        "\\end{equation}\n",
        "\n",
        "Com base nos dados fornecidos, escreva a matriz densidade para a molécula de $H_2$ numa base mínima.\n",
        "\n",
        "\\\n",
        "\n",
        "d) Use agora a relação dos coeficientes com os elementos da matriz de sobreposição para escrever a matriz de sobreposição $\\mathbf{S}$\n",
        "\n",
        "e) Use em alternativa a relação dos coeficiente com os elementos da matriz densidade $\\mathbf{P}$"
      ],
      "metadata": {
        "id": "CXoD_3r1u7QB"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "**Exercício 3**\n",
        "\n",
        "Utilize agora o programa PySCF para obter 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 ao nível HF/STO-3G. Compare com os valores fornecidos e com os que calculou nas alineas anteriores\n",
        "\n",
        "Nota 1: Use a distância indicada nos dados do problema sem optimizar\n",
        "\n",
        "Nota 2: Para extrair as matrizes de um cálculos podemos usar\n",
        "\n",
        "```\n",
        "# Matriz dos Coeficientes(C): cada coluna é uma OM!\n",
        "C = mf.mo_coeff  # shape (nbasis, nbasis)\n",
        "\n",
        "# Energias das orbitais (Matriz diagonal epsilon)\n",
        "# Vector 1D com os valores próprios\n",
        "epsilon = mf.mo_energy  \n",
        "\n",
        "# Matriz Densidade (P)\n",
        "# No caso RHF: 2 electrões por orbital ocupada\n",
        "# A dimensão é (nbasis, nbasis)\n",
        "P = mf.make_rdm1()  \n",
        "\n",
        "# Matriz de 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": [
        "**Exercício 4**\n",
        "\n",
        "Repita o exercício anterior usando 6-31g, 6-311g, 6-311g** e cc-pVDZ, cc-pVTZ e cc-pVQZ (não precisa copiar para o caderno).\n",
        "Verifique a dimensão das matrizes e justifique. Como varia a energia total?"
      ],
      "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": [
        "**Exercício 5**\n",
        "\n",
        "O integral da densidade de carga total é igual ao número de electrões $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",
        "Vimos também que podemos expressar a densidade de carga total usando a matriz densidade $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",
        "Escreva o integral da densidade de carga total escrito sob a forma da equação 2 e verifique o aparecimento do integral de sobreposição $\\mathbf{S}$\n",
        "\n",
        "**Exercício 6**\n",
        "\n",
        "Não existe uma definição única de quantificar o número de electrões associado a um dado átomo ou núcleo. Contudo, poderá ser útil efectuar uma análise de população como *Análise de População Mulliken*, na qual os electrões associados a um determinado átomo numa molécula são obtidos **atribuindo metade da densidade a cada átomo**. A população num átomo A é então dada por\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",
        "e por conseguinte, a carga de Mulliken no átomo A é:\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",
        "com base nos dados do Exercício 2, calcule a carga de Mulliken em cada átomo de hidrogénio. Era o valor que esperava?\n",
        "\n",
        "Nota: $\\sum_{\\mu \\in A}$ indica a soma em todas as orbitais atómicas centradas no átomo A; $\\sum_{\\nu}$ percorre todos as orbitais atómicos da base, i.e., o conjunto completo da molécula.\n",
        "\n",
        "\n",
        "**Exercício 7**\n",
        "\n",
        "Tendo por base o exercício 3, use agora o programa PySCF para calcular as cargas de Mulliken dos átomos de hidrogénio ao nível HF/STO-3G. Espera que as cargas variem com as funções base?\n",
        "\n",
        "Nota: No pySCF, podemos realizar uma análise de poplação de Mulliken usando `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": [
        "**Exercício 9**\n",
        "\n",
        "Tendo por base o exercício 3 da Aula 9 referente à molécula de água, optimize a geometria ao nível HF usando as bases sto-3g, 6-31g, 6-311g, 6-311g** e cc-pVDZ, cc-pVTZ e cc-pVQZ. Para cada base e respectiva geometria optimizada, calcule as cargas de Mulliken em cada átomo. Comente os resultados. Pode usar informação do Exercício 3 da Aula 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"
      ],
      "metadata": {
        "id": "zruDZjROVrSJ"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}