{ "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 9**\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 Calculations**\n", "\n", "\\\\\n", "\n", "**Exercise 1: Geometry Optimizations**\n", "\n", "As we saw in the previous lesson, we can use quantum mechanics calculations to optimize the geometry of molecules. In short, the program used (in this case, PySCF) evaluates the energy for various geometries, trying to find a minimum on the potential energy surface defined by the electronic energy.\n", "\n", "According to the values ​​reported in the Computational Chemistry Comparison and Benchmark DataBase (https://cccbdb.nist.gov/expgeom2x.asp) the experimental values ​​for the internal coordinates of the formaldehyde molecule\n", "\n", "
\n", "\n", "
\n", "\n", "are as follows: r(C=O) = 1.205 Å ; r(C-H) = 1.111 Å; α(H-C-H) = 116.1°; α(H-C-O) = 121.9°\n", "\n", "Optimize the geometry of the formaldehyde molecule at the HF level using the basis functions sto-3g, 6-31g, and cc-pVDZ. Report the values ​​obtained for the energy and for the angles and distances. Compare, when possible, with the experimental values.\n", "\n", "Note 1: Use the program from **Exercise 3** of the previous lesson as a basis and use the following as the starting coordinates (in Angstrom)\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", "Note 2: You can store the x, y, z coordinates of the optimized molecule in a variable, eg.\n", "\n", "```\n", "# store the coordinates of mol_eq in the variable coords\n", "coords = mol_eq.atom_coords()\n", "```\n", "\n", "Note 3: We can perform operations like calculating distances and angles on the coords matrix. For example\n", "\n", "```\n", "from scipy.spatial.distance import euclidean\n", "import numpy as np\n", "\n", "# stores in the variable d the Euclidean distance between element 0(z,y,z) and element 1(x,y,z)\n", "\n", "d = euclidean(coords[0], coords[1])\n", "\n", "# a function that calculates the angle\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)) # returns the angle in degrees\n", "\n", "# stores in the variable angle the angle between the vectors defined by\n", "# points/elements 0(x,y,z) 2(x,y,z) and 0(x,y,z) 1(x,y,z)\n", "angle = angle(coords[2], coords[0], coords[1])\n", "\n", "```\n", "\n", "Don't forget that the program internally works in atomic units and therefore, the values ​​returned are in this unit!" ], "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", "# Write the code here\n", "\n" ], "metadata": { "id": "ZIEevtteuqli" }, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": [ "**Exercise 2**\n", "\n", "Now visualize the last optimized geometry." ], "metadata": { "id": "sF4EzR7_EDui" } }, { "cell_type": "code", "source": [ "import py3Dmol\n", "\n", "# Write the code here" ], "metadata": { "id": "shAkLknRD2Ty" }, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": [ "**Exercise 3**\n", "\n", "In addition to predicting geometries, these calculations can be used to review other properties such as the dipole moment.\n", "\n", "The dipole moment of an isolated water molecule is 1.855 D (Debye) [https://www.science.org/doi/10.1126/science.275.5301.814].\n", "\n", "Optimize a water molecule at the HF level using the sto-3g, 6-31g, cc-pVDZ, cc-pVDZ and cc-pVQZ bases. For each base and its optimized geometry, calculate the energy and dipole moment. Also report the relevant parameters of the optimized geometry and compare with the experimental values ​​(https://doi.org/10.1002/jcc.20157)\n", "\n", "`rOH = 0.9572; aHOH = 104.52°`\n", "\n", "Comment on the results, in particular, whether the trend was as expected.\n", "\n", "Note 1: Use exercise 5 from the last lesson as a starting point.\n", "\n", "Note 2: The dipole moment can be obtained from `mf_eq.dip_moment()` where mf_eq is the name where the entire SCF calculation is stored. Remember that the dipole moment is a vector. To compare it with the experimental value we need to calculate the norm (see Lesson 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", "# Write your code here\n" ], "metadata": { "id": "ksXQ8I2O06Vm" }, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": [ "**Exercise 4**\n", "\n", "Now repeat the calculations to obtain the dipole moment at the HF level with the same bases using the experimental geometry of water without performing the optimization (\"single point\" calculation). Also check the energy value.\n", "\n", "Note 1: PySCF, instead of Cartesian coordinates, accepts a type of internal coordinates called Z-matrix that has the following format:\n", "\n", "```\n", "mol = gto.M(\n", "atom='''\n", "# atom1\n", "H 1 d_21 # atom 2, is bonded to 1 with a distance d_21\n", "H 1 d_31 2 a123 # atom 3, is bonded to 1 with a distance d_21 and makes an angle of a123 with atom 3\n", "''',\n", "unit='Angstrom'\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", "# Write your code here" ], "metadata": { "id": "zruDZjROVrSJ" }, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": [ "**Exercise 5**\n", "\n", "Consider the Coulomb operator acting on a 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", "Show that the expected value with respect to $\\chi_a$ is equal to\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", "**Exercise 6**\n", "\n", "Consider the Exchange operator acting on a 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", "Show that the expected value with respect to $\\chi_a$ is equal to\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", "**Exercise 7**\n", "\n", "Consider $( \\chi_a, \\chi_b )$ as orthonormal spin orbitals. The Fock operator acting on a spin-orbital $\\chi_a$ is given by:\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", "Write the expression for the elements of the Fock matrix $\\langle \\chi_a | \\hat{f} | \\chi_a \\rangle$ as a function of the mono- and bi-electron integrals. What is the sign of the Coulomb and Exchange terms?\n" ], "metadata": { "id": "LGuP2vbe08Uo" } }, { "cell_type": "markdown", "source": [ "**Exercise 8**\n", "\n", "The Hartree-Fock energy of a closed-shell determinant is\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", "Now consider the helium atom. Write the expectation value of the closed-shell Hartree-Fock energy in terms of the single- and two-electron integrals.\n" ], "metadata": { "id": "taU6ehM4jf-q" } }, { "cell_type": "markdown", "source": [ "**Exercise 9**\n", "\n", "Now use the PySCF program to evaluate the values ​​of the contribution of the core, Coulomb and exchange terms. Based on these, calculate the ground state energy for the helium atom at the RHF/sto-3g level. Compare the value with that obtained directly by the program.\n", "\n", "Note: To get the contributions you can use\n", "\n", "```\n", "# Density matrix: total and 1/2\n", "dm = mf.make_rdm1()\n", "dm_half = dm / 2 # To get J and K in the restricted case\n", "\n", "# Hamiltonian core h(1)\n", "hcore = mf.get_hcore()\n", "\n", "# Get J and K from dm_half\n", "J = mf.get_j(dm=dm_half)\n", "K = mf.get_k(dm=dm_half)\n", "\n", "# Evaluate the energies - sum over the elements of the matrices\n", "# Avaliar as energias - soma ao longo dos elementos das matrizes\n", "E_1e = (np.einsum('ij,ij', dm, hcore))/2 # energy core (1e inegral) per electron\n", "E_J = np.einsum('ij,ij', dm_half, J) # Coulomb energy (2e integrals)\n", "E_K = np.einsum('ij,ij', dm_half, K) # Exchange energy (2e integrals)\n", "```" ], "metadata": { "id": "TdS78_g1nqDo" } }, { "cell_type": "code", "source": [ "from pyscf import gto, scf\n", "import numpy as np\n", "\n", "# Write your code here" ], "metadata": { "id": "E8blioKviUpu" }, "execution_count": null, "outputs": [] } ] }