{ "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 9**\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: optimizações de geometria**\n", "\n", "Como verificámos na aula anterior, podemos usar cálculos de mecânica quântica para optimizar a geometria de moléculas. De forma muito sucinta, o programa utilizado (neste caso, PySCF) avalia a energia para várias geometrias, tentando encontrar um mínimo na superfície de energia potencial definida pela energia electrónica.\n", "\n", "De acordo com os valores reportados na Computational Chemistry Comparison and Benchmark DataBase (https://cccbdb.nist.gov/expgeom2x.asp) os valores experimentais para as coordenada internas da molécula de formaldeído\n", "\n", "
\n", " \n", "
\n", "\n", "são as seguintes: r(C=O) = 1.205 Å ; r(C-H) = 1.111\tÅ; α(H-C-H) = 116.1°; α(H-C-O) = 121.9° \n", "\n", "Faça a optimização da geometria da molécula de formaldeído ao nível H usando as funções base sto-3g, 6-31g, and cc-pVDZ. Reporte os valores obtidos para a energia e para os ângulos e distâncias. Compare, quando possível, com os valores experimentais.\n", "\n", "Nota 1: Use o programa do **Exercício 3** da aula anterior como base e utilize o seguinte como coordenadas (em Angstrom) de partida\n", "\n", "```\n", " H 1.0686 -0.1411 1.0408\n", " C 0.5979 0.0151 0.0688\n", " H 1.2687 0.2002 -0.7717\n", " O -0.5960 -0.0151 -0.0686\n", "```\n", "\n", "Nota 2: Pode armazenar as coordenadas x, y, z da molecula otimizada numa variável, eg.\n", "\n", "```\n", "# armazena as coordenadas de mol_eq na variável coords\n", "coords = mol_eq.atom_coords()\n", "```\n", "\n", "Nota 3: Podemos fazer operações como cálcular distâncias e ângulos na matriz coords. Por exemplo\n", "\n", "```\n", "from scipy.spatial.distance import euclidean\n", "import numpy as np\n", "\n", "# armazena na variável d a distância euclidiana entre o elemento 0(z,y,z) e o elemento 1(x,y,z)\n", "\n", "d = euclidean(coords[0], coords[1])\n", "\n", "# uma função que calcula o ângulo\n", "def angle(a, b, c):\n", " ba = a - b # vector from B to A\n", " bc = c - b # vector from B to C\n", " cosine = np.dot(ba, bc) / (np.linalg.norm(ba) * np.linalg.norm(bc))\n", " return np.degrees(np.arccos(cosine)) # devolve o ângulo em graus\n", "\n", "# armazena na variável angulo o angulo entre os vectores definidos pelos\n", "# pontos/elementos 0(x,y,z) 2(x,y,z) e 0(x,y,z) 1(x,y,z)\n", "angulo = angle(coords[2], coords[0], coords[1])\n", "\n", "```\n", "\n", "Não se esqueça que o programa internamente trabalha em unidades atómicas e por conseguinte, os valores devolvidos estão nesta unidade!" ], "metadata": { "id": "CXoD_3r1u7QB" } }, { "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 o seu código aqui\n", "\n" ], "metadata": { "id": "ZIEevtteuqli" }, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": [ "**Exercício 2**\n", "\n", "Visualize agora a última geometria optimizada." ], "metadata": { "id": "sF4EzR7_EDui" } }, { "cell_type": "code", "source": [ "import py3Dmol\n", "\n", "# Escreva o seu código aqui" ], "metadata": { "id": "shAkLknRD2Ty" }, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": [ "**Exercício 3**\n", "\n", "Para além da previsão de geometrias, estes cálculos podem ser usados para rever outras propriedades tais como, o momento dipolar.\n", "\n", "O momento dipolar de uma molécula de água isolada é de 1.855 D (Debye) [https://www.science.org/doi/10.1126/science.275.5301.814].\n", "\n", "Optimize uma molécula de água ao nível HF usando as bases sto-3g, 6-31g, cc-pVDZ, cc-pVDZ e cc-pVQZ. Para cada base e respectiva geometria optimizada, calcule o a energia e o momento dipolar. Reporte também os parâmetros relevantes da geometria optimizada e compare com os valores experimentais (https://doi.org/10.1002/jcc.20157)\n", "\n", "`rOH = 0.9572; aHOH = 104.52°`\n", "\n", "Comente os resultados, em particular, se a tendência era a esperada.\n", "\n", "Nota 1: Use o exercício 5 da aula passada como ponto de partida.\n", "\n", "Nota 2: o momento dipolar pode ser obtido a partir de `mf_eq.dip_moment()` onde mf_eq é o nome onde se encontra armazenado todo o cálculo SCF. Recorde que o momento dipolar é um vector. Para comparar com o valor experimental precisamos de calcular a norma (ver Aula 2).\n" ], "metadata": { "id": "n5U01KrcaKXO" } }, { "cell_type": "code", "source": [ "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 o seu código aqui\n" ], "metadata": { "id": "ksXQ8I2O06Vm" }, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": [ "**Exercício 4**\n", "\n", "Repita agora os cálculos para obter o momento dipolar ao nível HF com as mesmas bases usando a geometria experimental da água sem efectuar a optimização (cálculo \"single point\"). Verifique também o valor da energia.\n", "\n", "Nota 1: O PySCF em vez de coordenadas cartesianas, aceita uma espécie de coordenadas internas chamada de matriz-Z que tem o seguinte formato:\n", "\n", "```\n", "mol = gto.M(\n", " atom='''\n", " O # átomo1\n", " H 1 d_21 # átomo 2, está ligado ao 1 com uma distância d_21\n", " H 1 d_31 2 a123 # átomo 3, está ligado ao 1 com uma distância d_21 e faz com o átomo 3 um ângulo de a123\n", " ''',\n", " unit='Angstrom'\n", " )\n", "```\n" ], "metadata": { "id": "VWWpPBkGSZl1" } }, { "cell_type": "code", "source": [ "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 o seu código aqui" ], "metadata": { "id": "zruDZjROVrSJ" }, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": [ "**Exercício 5**\n", "\n", "Considere o operador de Coulomb a actuar numa spin-orbital $\\chi_a(\\mathbf{x}_1)$:\n", "\n", "\\begin{equation}\n", " \\mathcal{J}_b(\\mathbf{x}_1) \\chi_a(\\mathbf{x}_1) = \\left[ \\int d\\mathbf{x}_2 \\ \\chi_b^*(\\mathbf{x}_2) \\frac{1}{r_{12}} \\chi_b(\\mathbf{x}_2) \\right] \\chi_a(\\mathbf{x}_1)\n", "\\end{equation}\n", "\n", "Mostre que o valor esperado com respeito a $\\chi_a$ é igual a\n", "\n", "\\begin{equation}\n", " \\langle \\chi_a(\\mathbf{x}_1) | \\mathcal{J}_b(\\mathbf{x}_1) | \\chi_a(\\mathbf{x}_1) \\rangle = \\langle ab | ab \\rangle\n", "\\end{equation}\n", "\n", "**Exercício 6**\n", "\n", "Considere o operador de Troca a actuar numa spin-orbital $\\chi_a(\\mathbf{x}_1)$:\n", "\n", "\\begin{equation}\n", " \\mathcal{K}_b(\\mathbf{x}_1) \\chi_a(\\mathbf{x}_1) = \\left[ \\int d\\mathbf{x}_2 \\ \\chi_b^*(\\mathbf{x}_2) \\frac{1}{r_{12}} \\chi_a(\\mathbf{x}_2) \\right] \\chi_b(\\mathbf{x}_1)\n", "\\end{equation}\n", "\n", "Mostre que o valor esperado com respeito a $\\chi_a$ é igual a\n", "\n", "\\begin{equation}\n", " \\langle \\chi_a(\\mathbf{x}_1) | \\mathcal{K}_b(\\mathbf{x}_1) | \\chi_a(\\mathbf{x}_1) \\rangle = \\langle ab | ba \\rangle\n", "\\end{equation}\n", "\n", "**Exercício 7**\n", "\n", "Considere $( \\chi_a, \\chi_b )$ como spin orbitais ortonormais. O operador de Fock a actuar numa spin-orbital $\\chi_a$ é dado por:\n", "\n", "$$\n", "\\hat{f} \\chi_a(\\mathbf{x}_1) = \\hat{h} \\chi_a(\\mathbf{x}_1) + \\sum_{b}^{\\text{occ}} \\left[ \\mathcal{J}_b(\\mathbf{x}_1) - \\mathcal{K}_b(\\mathbf{x}_1) \\right] \\chi_a(\\mathbf{x}_1)\n", "$$\n", "\n", "Escreva a expressão para os elementos da matriz de Fock $\\langle \\chi_a | \\hat{f} | \\chi_a \\rangle$ em função dos integrais mono e bi-electrónicos. Qual é o sinal dos termos de Coulomb e de Troca?\n" ], "metadata": { "id": "LGuP2vbe08Uo" } }, { "cell_type": "markdown", "source": [ "**Exercício 8**\n", "\n", "A energia Hartree-Fock de um determinante de camada fechada é\n", "\n", "\\begin{equation}\n", " E_0 = 2\\sum\\limits_{a}^{N/2} h_{aa} + \\sum\\limits_{a}^{N/2} \\sum\\limits_{b}^{N/2} 2 J_{ab} - K_{ab}\n", "\\end{equation}\n", "\n", "Considere agora o átomo de Hélio. Escreva o valor esperado da energia de Hartree-Fock em camada fechada em função dos integrais mono e bi-eletrónicos.\n" ], "metadata": { "id": "taU6ehM4jf-q" } }, { "cell_type": "markdown", "source": [ "**Exercício 9**\n", "\n", "Utilize agora o programa PySCF para avaliar os valores da contribuição dos termos core, Coulomb e Troca. Com base nestes, calcule a energia do estado fundamental para o átomo de hélio ao nível RHF/sto-3g. Compare o valor com o obtido directamente pelo programa.\n", "\n", "Nota: Para obter as contribuições pode usar\n", "\n", "```\n", "# Matriz densidade: total e 1/2\n", "dm = mf.make_rdm1()\n", "dm_half = dm / 2 # Para obter J e K no caso restrito\n", "\n", "# Hamiltoneano core h(1)\n", "hcore = mf.get_hcore()\n", "\n", "# Obter J e K a partir de dm_half\n", "J = mf.get_j(dm=dm_half)\n", "K = mf.get_k(dm=dm_half)\n", "\n", "# Avaliar as energias - soma ao longo dos elementos das matrizes\n", "# Avaliar as energias - soma ao longo dos elementos das matrizes\n", "E_1e = (np.einsum('ij,ij', dm, hcore))/2 # energia core (integrais 1e) por electrão\n", "E_J = np.einsum('ij,ij', dm_half, J) # energia Coulomb (integrais 2e)\n", "E_K = np.einsum('ij,ij', dm_half, K) # energia de Troca (integrais 2e)\n", "```" ], "metadata": { "id": "TdS78_g1nqDo" } }, { "cell_type": "code", "source": [ "from pyscf import gto, scf\n", "import numpy as np\n", "\n", "# escreva o seu código aqui" ], "metadata": { "id": "E8blioKviUpu" }, "execution_count": null, "outputs": [] } ] }