{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# 05 The Closed-Shell CCSD energy" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The coupled cluster model provides a higher level of accuracy beyond the MP2 approach. The purpose of this project is to understand the fundamental aspects of the calculation of the CCSD (coupled cluster singles and doubles) energy. Reference to this project is [Hirata, ..., Bartlett, JCP 2004](https://dx.doi.org/10.1063/1.1637577) (though notations are different, and this project is not discussing extended systems), and PySCF code ([rintermediates.py](https://github.com/pyscf/pyscf/blob/master/pyscf/cc/rintermediates.py), [rccsd.py](https://github.com/pyscf/pyscf/blob/master/pyscf/cc/rccsd.py)) (dimension convention is similar to PySCF implementation).\n", "\n", "This project will use spacial orbitals (like previous projects, [SCF](../Project_03/Project_03.ipynb) and [MP2](../Project_04/Project_04.ipynb) energy calculation), instead of more computation-costly spin orbitals. However, the latter approach could be also applicable in situation of unrestricted or restricted open-shell. We will discuss spin orbital approach in latter projects.\n", "\n", "This project could be a challenging one, and the coding is quite intensive in this project. Transforming tensor contraction formula to codes is the gist of this project. We will make extensive use of `numpy.einsum` here." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "tags": [ "thebe-init", "hide-cell" ] }, "outputs": [], "source": [ "# Following os.chdir code is only for thebe (live code), since only in thebe default directory is /home/jovyan\n", "import os\n", "if os.getcwd().split(\"/\")[-1] != \"Project_05\":\n", " os.chdir(\"source/Project_05\")\n", "from solution_05 import Molecule as SolMol\n", "\n", "from pyscf import gto\n", "import numpy as np\n", "from typing import Tuple\n", "import pickle\n", "np.set_printoptions(precision=7, linewidth=120, suppress=True)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "tags": [ "thebe-init", "hide-cell" ] }, "outputs": [], "source": [ "# Solution mol only uses PySCF approach\n", "sol_mole = SolMol()\n", "sol_mole.construct_from_dat_file(\"input/h2o/STO-3G/geom.dat\")\n", "sol_mole.obtain_mol_instance(basis=\"STO-3G\")\n", "sol_mole.obtain_nao()\n", "sol_mole.obtain_nocc()\n", "sol_mole.obtain_eri_ao()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Molecule Object Initialization" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this project, we may use the even updated Molecule initialization. Since these are technical details, we toggle the code and explanation below. However, take a look at those instructions could be essential. Most of the method functions are illustrated in [Project 01](../Project_01/Project_01.ipynb), [Project 03](../Project_03/Project_03.ipynb) and [Project 04](../Project_04/Project_04.ipynb).\n", "\n", "Note that since details of SCF iteration are no longer required, we can make some simplification to SCF code." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "tags": [ "thebe-init", "hide-cell" ] }, "outputs": [], "source": [ "class Molecule:\n", " def __init__(self):\n", " # Project 03 existed\n", " self.atom_charges = NotImplemented # type: np.ndarray\n", " self.atom_coords = NotImplemented # type: np.ndarray\n", " self.natm = NotImplemented # type: int\n", " self.mol = NotImplemented # type: gto.Mole\n", " self.nao = NotImplemented # type: int\n", " self.charge = 0 # type: int\n", " self.nocc = NotImplemented # type: int\n", " # Project 04 added\n", " self.mo_coeff = NotImplemented # type: np.ndarray\n", " self.mo_energy = NotImplemented # type: np.ndarray\n", " self.eri_ao = NotImplemented # type: np.ndarray\n", " self.eri_mo = NotImplemented # type: np.ndarray\n", " self.energy_rhf = NotImplemented # type: np.ndarray\n", " self.energy_corr = NotImplemented # type: np.ndarray\n", " # Project 05 added\n", " self.w = NotImplemented # type: np.ndarray\n", " self.w_ovoo = NotImplemented # type: np.ndarray\n", " self.w_ovvo = NotImplemented # type: np.ndarray\n", " self.w_ovov = NotImplemented # type: np.ndarray\n", " self.w_ovvv = NotImplemented # type: np.ndarray\n", " self.v_oooo = NotImplemented # type: np.ndarray\n", " self.v_ovoo = NotImplemented # type: np.ndarray\n", " self.v_ovvo = NotImplemented # type: np.ndarray\n", " self.v_ovov = NotImplemented # type: np.ndarray\n", " self.v_oovv = NotImplemented # type: np.ndarray\n", " self.v_ovvv = NotImplemented # type: np.ndarray\n", " self.v_vvvv = NotImplemented # type: np.ndarray\n", "\n", " def construct_from_dat_file(self, file_path: str):\n", " # Same to Project 01\n", " with open(file_path, \"r\") as f:\n", " dat = np.array([line.split() for line in f.readlines()][1:])\n", " self.atom_charges = np.array(dat[:, 0], dtype=float).astype(int)\n", " self.atom_coords = np.array(dat[:, 1:4], dtype=float)\n", " self.natm = self.atom_charges.shape[0]\n", "\n", " def obtain_mol_instance(self, basis: str, verbose=0):\n", " # Same to Project 03\n", " mol = gto.Mole()\n", " mol.unit = \"Bohr\"\n", " mol.atom = \"\\n\".join([(\"{:3d} \" + \" \".join([\"{:25.18f}\"] * 3)).format(chg, *coord) for chg, coord in zip(self.atom_charges, self.atom_coords)])\n", " mol.basis = basis\n", " mol.charge = self.charge\n", " mol.spin = 0\n", " mol.verbose = verbose\n", " self.mol = mol.build()\n", "\n", " def obtain_nao(self):\n", " self.nao = self.mol.nao_nr()\n", "\n", " def obtain_nocc(self):\n", " assert (self.atom_charges.sum() - self.charge) % 2 == 0\n", " self.nocc = (self.atom_charges.sum() - self.charge) // 2\n", "\n", " def obtain_eri_ao(self):\n", " self.eri_ao = self.mol.intor(\"int2e\")\n", "\n", " def get_hcore(self) -> np.ndarray:\n", " return self.mol.intor(\"int1e_kin\") + self.mol.intor(\"int1e_nuc\")\n", "\n", " def get_fock(self, dm: np.ndarray) -> np.ndarray:\n", " return self.get_hcore() + (self.eri_ao * dm).sum(axis=(-1, -2)) - 0.5 * (self.eri_ao * dm[:, None, :]).sum(axis=(-1, -3))\n", "\n", " def make_rdm1(self, coeff: np.ndarray) -> np.ndarray:\n", " return 2 * coeff[:, :self.nocc] @ coeff[:, :self.nocc].T\n", "\n", " def eng_total(self, dm: np.ndarray) -> float:\n", " return (0.5 * (self.get_hcore() + self.get_fock(dm)) * dm).sum() + self.mol.energy_nuc()\n", "\n", " def scf_process(self, dm_guess: np.ndarray=None) -> Tuple[float, np.ndarray]:\n", " eng, dm = 0., np.zeros((self.nao, self.nao)) if dm_guess is None else np.copy(dm_guess)\n", " max_iter, thresh_eng, thresh_dm = 64, 1e-10, 1e-8\n", " for epoch in range(max_iter):\n", " eng_next, dm_next = self.eng_total(dm), self.make_rdm1(scipy.linalg.eigh(self.get_fock(dm), self.mol.intor(\"int1e_ovlp\"))[1])\n", " if np.abs(eng_next - eng) < thresh_eng and np.linalg.norm(dm_next - dm) < thresh_dm:\n", " eng, dm = eng_next, dm_next\n", " break\n", " eng, dm = eng_next, dm_next\n", " return eng, dm\n", "\n", " def obtain_scf_intermediates(self, dm_guess: np.ndarray=None):\n", " eng, dm = self.scf_process(dm_guess)\n", " self.energy_rhf = eng\n", " self.mo_energy, self.mo_coeff = scipy.linalg.eigh(self.get_fock(dm), self.mol.intor(\"int1e_ovlp\"))\n", "\n", " def get_eri_mo_einsum(self):\n", " return np.einsum(\"uvkl, up, vq, kr, ls -> pqrs\", self.eri_ao, self.mo_coeff, self.mo_coeff, self.mo_coeff, self.mo_coeff, optimize=True)\n", "\n", " def obtain_eri_mo(self):\n", " self.eri_mo = self.get_eri_mo_einsum()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "````{toggle}\n", "\n", "Lots of attributes like `w_ovov` is introduced in this project. We will illustrate meaning of these attributes later.\n", "\n", "Since this project is hard to achieve, and numerical problems could occur. Numerical problem means for molecular orbital energy degenerate or non-degenerate systems, $\\mathscr{C}_{\\mu \\mathscr{p}} = - C_{\\mu p}$ is also valid molecular orbital coefficients, but that will cause $(\\mathscr{p} q | rs) \\neq (pq|rs)$. Although final results (SCF, MP2, CCSD) should left unchanged, calculation intermediates could be different and hard to reproduce.\n", "\n", "Thus, in this project, we may use pre-computed water/STO-3G results to check whether calculation intermediates is correctly implemented.\n", "- `mo_coeff`: pre-computed $C_{\\mu p}$ molecular orbital coefficients\n", "- `mo_energy`: pre-computed $\\varepsilon_p$ molecular orbital energies\n", "- `t1_precomput`: pre-computed $t_i^a$ CCSD single excitation amplitude\n", "- `t2_precomput`: pre-computed $t_{ij}^{ab}$ CCSD double excitation amplitude\n", "\n", "You may wish to run the following `Molecule` initialization code after running the code-cell below:\n", "\n", "```python\n", "mole = Molecule()\n", "mole.construct_from_dat_file(\"input/h2o/STO-3G/geom.dat\")\n", "mole.obtain_mol_instance(basis=\"STO-3G\")\n", "mole.obtain_nao()\n", "mole.obtain_nocc()\n", "mole.obtain_eri_ao()\n", "mole.mo_coeff = mo_coeff\n", "mole.mo_energy = mo_energy\n", "mole.obtain_eri_mo()\n", "```\n", "\n", "````" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "tags": [ "thebe-init", "hide-cell" ] }, "outputs": [], "source": [ "with open(\"demo_data_h2o_sto3g.dat\", \"rb\") as f:\n", " d = pickle.load(f)\n", "mo_coeff = d[\"mo_coeff\"]\n", "mo_energy = d[\"mo_energy\"]\n", "t1_precomput, t2_precomput = d[\"t1\"], d[\"t2\"]\n", "\n", "sol_mole.mo_coeff = mo_coeff\n", "sol_mole.mo_energy = mo_energy\n", "sol_mole.obtain_eri_mo()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 1: ERI and Its Biorthogonal Form" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In closed-shell CCSD prototype programming, *electron repulsion integral* (ERI) can be ambiguous. We usually use\n", "\n", "$$\n", "v^{pq}_{rs} = (pr|qs)\n", "$$\n", "\n", "to denote electron repulsion integrals (ERI), as previous projects did before. Here we define *biorthogonal* electron repulsion integral (biorthogonal ERI) as (Hirata, eq 23)\n", "\n", "$$\n", "w^{pq}_{rs} = 2 v^{pq}_{rs} - v^{pq}_{sr}\n", "$$\n", "\n", "Biorthogonal basis is a conception from closed-shell CCSD derivation, which makes formula expression as well as program implementation more concise. We do not discuss this concept in detail." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "````{admonition} Dimension convention of ERI\n", ":class: dropdown\n", "\n", "Dimension convention to these tensors is the same to `Molecule.eri_mo`, i.e.\n", "- `Molecule.eri_mo` $v^{pq}_{rs}$, dim: $(p, r, q, s)$\n", "- `Molecule.w` $w^{pq}_{rs}$, dim: $(p, r, q, s)$\n", "\n", "Take caution when transposing those tensors.\n", "\n", "````" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Implementation" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "def obtain_w(mole: Molecule):\n", " # Attribute Modification: `w` biorthogonal electron repulsion integral\n", " raise NotImplementedError(\"Exactly 1 line of code\")\n", "\n", "Molecule.obtain_w = obtain_w" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solution" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "data": { "text/plain": [ "array([[ 0.5311873, 0.0277178, 0. , 0.0183313, -0. , -0.0017709, 0. ],\n", " [ 0.0253682, 0.218575 , -0. , 0.1602001, 0. , -0.060252 , -0. ],\n", " [ 0. , -0. , 0.1339514, 0. , 0. , -0. , 0.1456495],\n", " [ 0.0012042, 0.0477327, -0. , 0.1319441, -0. , 0.0955699, 0. ],\n", " [-0. , 0. , 0. , -0. , 0.2900729, 0. , -0. ],\n", " [-0.0215503, -0.0719668, -0. , -0.352678 , 0. , 0.0524259, -0. ],\n", " [ 0. , -0. , 0.1382612, 0. , -0. , 0. , 0.1620096]])" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sol_mole.obtain_w()\n", "sol_mole.w[3, 5]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 2: ERI Slices" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will use various kinds of ERIs slices when programming closed-shell CCSD. It could be convenient to pre-store those slices. To storing these slices one may use class attributes, class properties, dictionary, or generate slices on-the-fly by convenient function. In this project, we may use class attributes to pre-store slices, which is the most intutive one." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "````{admonition} Variable naming convention\n", ":class: dropdown\n", "\n", "The most commonly used ERI slice is $v^{ij}_{ab}$. This slice is named to `v_ovov`, since it's dimension is $(i, a, j, b)$, with $i, j$ occupied orbitals and $a, b$ virtual orbitals.\n", "\n", "```python\n", "mole.v_ovov = mole.eri_mo[:mole.nocc, mole.nocc:, :mole.nocc, mole.nocc:]\n", "```\n", "\n", "For $w^{ij}_{ak}$, this slice is named to `w_ovoo`, since it's dimension is $(i, a, j, k)$, with $i, j, k$ occupied orbitals and $a$ virtual orbitals.\n", "\n", "```python\n", "mole.w_ovoo = mole.w[:mole.nocc, mole.nocc:, :mole.nocc, :mole.nocc]\n", "```\n", "\n", "````" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "````{admonition} Slice object in Python\n", ":class: dropdown\n", "\n", "Code above could be somehow unclear, for that `mole.nocc:` looks and codes very similar to `:mole.nocc`. To mitigate this issue, we could use slice object in python.\n", "\n", "```python\n", "so, sv = slice(0, mole.nocc), slice(mole.nocc, mole.nao)\n", "```\n", "\n", "Then `so` works similar to `0:mole.nocc`, i.e. occupied orbital slice; and `sv` works similar to `mole.nocc:mole.nao`, i.e. virtual orbital slice. Note that $n_\\mathrm{AO}$ `mole.nao` is the same to $n_\\mathrm{MO}$ in cases of this project. So code of `mole.v_ovov` could be\n", "\n", "```python\n", "mole.v_ovov = mole.eri_mo[so, sv, so, sv]\n", "```\n", "\n", "````" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Implementation" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "def obtain_wv_slices(mole: Molecule):\n", " # Attribute Modification: Various (biorthogonal) ERI slices\n", " nocc, nmo = mole.nocc, mole.nao\n", " so, sv = slice(0, nocc), slice(nocc, nmo)\n", " mole.w_ovoo = self.w[so, sv, so, so]\n", " mole.w_ovvo = ____ # Fill this line\n", " mole.w_ovov = ____ # Fill this line\n", " mole.w_ovvv = ____ # Fill this line\n", " mole.v_oooo = ____ # Fill this line\n", " mole.v_ovoo = ____ # Fill this line\n", " mole.v_ovvo = ____ # Fill this line\n", " mole.v_ovov = self.eri_mo[so, sv, so, sv]\n", " mole.v_oovv = ____ # Fill this line\n", " mole.v_ovvv = ____ # Fill this line\n", " mole.v_vvvv = ____ # Fill this line\n", "\n", "Molecule.obtain_wv_slices = obtain_wv_slices" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solution" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(5, 2, 2, 5)\n", "[[ 0. -0. 0.0541371 0. -0. ]\n", " [-0.0112477 0.0171886 0. 0.0709527 0. ]]\n" ] } ], "source": [ "sol_mole.obtain_wv_slices()\n", "print(sol_mole.v_ovvo.shape)\n", "print(sol_mole.v_ovvo[3, 1])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 3: $\\boldsymbol{\\mathscr{F}}$ Matrices" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 3 to Step 5 is to implement CCSD calculation intermediates. In this step, we will calculate $\\boldsymbol{\\mathscr{F}}$ matrices (Hirata, eq 37-39, note that our reference state is canonical RHF, so off-diagonal occupied-virtual block Fock matrix $F_{ia} = 0$).\n", "\n", "$$\n", "\\begin{align}\n", "\\mathscr{F}^k_i &= f^k_i + \\sum_{lcd} w^{kl}_{cd} t_{il}^{cd} + \\sum_{lcd} w^{kl}_{cd} t_i^c t_l^d \\\\\n", "\\mathscr{F}^a_c &= f^a_c - \\sum_{kld} w^{kl}_{cd} t_{kl}^{ad} - \\sum_{kld} w^{kl}_{cd} t_k^a t_l^d \\\\\n", "\\mathscr{F}^k_c &= \\sum_{ld} w^{kl}_{cd} t_l^d\n", "\\end{align}\n", "$$\n", "\n", "Where $f^p_q$ is the same to molecular orbital basis Fock matrix $f^p_q = F_{pq} = \\delta_{pq} \\varepsilon_p$. Dimension convention of $\\mathscr{F}^k_i$ is $(k, i)$. Same to $\\mathscr{F}^a_c$ and $\\mathscr{F}^k_c$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "````{admonition} Dimension convention of CCSD amplitudes\n", ":class: dropdown\n", "\n", "Dimension convention of CCSD single excitation amplitude $t_i^a$ is $(i, a)$:\n", "```python\n", ">>> t1_precomput.shape\n", "(5, 2)\n", "```\n", "Dimension of CCSD double excitation amplitude $t_{ij}^{ab}$ is $(i, j, a, b)$:\n", "```python\n", ">>> t2_precomput.shape\n", "(5, 5, 2, 2)\n", "```\n", "This is quite different to convention of ERIs. So take caution when dealing those tensors!\n", "\n", "````" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Implementation" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "def cc_Foo(mole: Molecule, t1: np.ndarray, t2: np.ndarray) -> np.ndarray:\n", " # Ref: Hitara, eq 37, generate F^k_i\n", " # This is a reference implementation. There is no need to change this code.\n", " Fki = np.diag(self.mo_energy[:self.nocc])\n", " Fki += np.einsum(\"kcld, ilcd -> ki\", self.w_ovov, t2, optimize=True)\n", " Fki += np.einsum(\"kcld, ic, ld -> ki\", self.w_ovov, t1, t1, optimize=True)\n", " return Fki\n", "\n", "Molecule.cc_Foo = cc_Foo" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "def cc_Fvv(mole: Molecule, t1: np.ndarray, t2: np.ndarray) -> np.ndarray:\n", " # Ref: Hitara, eq 38\n", " raise NotImplementedError(\"No more than 4 lines of code\")\n", "\n", "Molecule.cc_Fvv = cc_Fvv" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "def cc_Fov(mole: Molecule, t1: np.ndarray, t2: np.ndarray) -> np.ndarray:\n", " # Ref: Hitara, eq 39\n", " # Note that amplitude t2 is actually not taken into account,\n", " # but for signature consistency, we still include this useless amplitude\n", " raise NotImplementedError(\"Exactly 1 line of code\")\n", "\n", "Molecule.cc_Fov = cc_Fov" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solution" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "data": { "text/plain": [ "array([[-20.2629403, 0.0004753, -0. , 0.0008499, 0. ],\n", " [ 0.0000935, -1.2196514, -0. , 0.014577 , 0. ],\n", " [ -0. , -0. , -0.5825172, 0. , 0. ],\n", " [ 0.0000376, 0.0098341, 0. , -0.4614169, -0. ],\n", " [ 0. , 0. , 0. , -0. , -0.3888219]])" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sol_mole.cc_Foo(t1_precomput, t2_precomput)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "data": { "text/plain": [ "array([[ 0.5124191, -0. ],\n", " [-0. , 0.624019 ]])" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sol_mole.cc_Fvv(t1_precomput, t2_precomput)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "data": { "text/plain": [ "array([[ 0.0001061, -0. ],\n", " [ 0.0012685, 0. ],\n", " [ 0. , -0.0032164],\n", " [-0.0022473, -0. ],\n", " [-0. , 0. ]])" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sol_mole.cc_Fov(t1_precomput, t2_precomput)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 4: $\\boldsymbol{\\mathscr{L}}$ Matrices" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this step, we will calculate $\\boldsymbol{\\mathscr{L}}$ matrices (Hirata, eq 40-41)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align}\n", "\\mathscr{L}^k_i &= \\mathscr{F}^k_i + \\sum_{lc} w^{lk}_{ci} t_l^c \\\\\n", "\\mathscr{L}^a_c &= \\mathscr{F}^a_c + \\sum_{kd} w^{ka}_{dc} t_k^d\n", "\\end{align}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Dimension convention of $\\mathscr{L}^k_i$ is $(k, i)$. Same to $\\mathscr{L}^a_c$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Implementation" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "def Loo(mole: Molecule, t1: np.ndarray, t2: np.ndarray) -> np.ndarray:\n", " # Ref: Hitara, eq 40\n", " # This is a reference implementation. There is no need to change this code.\n", " L_ki = self.cc_Foo(t1, t2)\n", " L_ki += np.einsum(\"lcki, lc -> ki\", self.w_ovoo, t1, optimize=True)\n", " return L_ki\n", "\n", "Molecule.Loo = Loo" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "def Lvv(mole: Molecule, t1: np.ndarray, t2: np.ndarray) -> np.ndarray:\n", " # Ref: Hitara, eq 41\n", " raise NotImplementedError(\"About 3 lines of code\")\n", "\n", "Molecule.Lvv = Lvv" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solution" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "data": { "text/plain": [ "array([[-20.2746448, -0.0001825, -0. , 0.000496 , 0. ],\n", " [ -0.0005209, -1.2246573, -0. , 0.0106743, 0. ],\n", " [ -0. , -0. , -0.5850956, 0. , 0. ],\n", " [ 0.0001375, 0.0087952, 0. , -0.4644571, -0. ],\n", " [ 0. , 0. , 0. , -0. , -0.3953228]])" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sol_mole.Loo(t1_precomput, t2_precomput)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "data": { "text/plain": [ "array([[ 0.5111931, -0. ],\n", " [-0. , 0.6206725]])" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sol_mole.Lvv(t1_precomput, t2_precomput)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 5: $\\boldsymbol{\\mathscr{W}}$ Tensors" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this step, we will calculate $\\boldsymbol{\\mathscr{W}}$ Tensors (Hirata, eq 42-45).\n", "\n", "$$\n", "\\begin{align}\n", "\\mathscr{W}^{kl}_{ij} &= v^{kl}_{ij} + \\sum_c v^{lk}_{ci} t_j^c + \\sum_c v^{kl}_{cj} t_i^c + \\sum_{cd} v^{kl}_{cd} t_{ij}^{cd} + \\sum_{cd} v^{kl}_{cd} t_i^c t_j^d \\\\\n", "\\mathscr{W}^{ab}_{cd} &= v^{ab}_{cd} - \\sum_k v^{ka}_{dc} t_k^b - \\sum_k v^{kb}_{cd} t_k^a\\\\\n", "\\mathscr{W}^{ak}_{ic} &= v^{ak}_{ic} - \\sum_l v^{kl}_{ci} t_l^a + \\sum_d v^{ka}_{cd} t_i^d - \\frac{1}{2} \\sum_{ld} v^{lk}_{dc} t_{il}^{da} - \\sum_{ld} v^{lk}_{dc} t_i^d t_l^a + \\frac{1}{2} \\sum_{ld} w^{lk}_{dc} t_{il}^{ad} \\\\\n", "\\mathscr{W}^{ak}_{ci} &= v^{ak}_{ci} - \\sum_l v^{lk}_{ci} t_l^a + \\sum_d v^{ka}_{dc} t_i^d - \\frac{1}{2} \\sum_{ld} v^{lk}_{cd} t_{il}^{da} - \\sum_{ld} v^{lk}_{cd} t_i^d t_l^a\n", "\\end{align}\n", "$$\n", "\n", "Dimension convention of $\\mathscr{W}^{kl}_{ij}$ is $(k, l, i, j)$. Same to $\\mathscr{W}^{ab}_{cd}$, $\\mathscr{W}^{ak}_{ic}$, $\\mathscr{W}^{ak}_{ci}$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "````{admonition} Einstein summation convention\n", ":class: dropdown\n", "\n", "We have learned using `numpy.einsum`. Actually, the original purpose of this function is to make *Einstein summation* by program. In [Einstein summation convention](http://en.wikipedia.org/wiki/Einstein_notation), summation notations are dropped out. One can know that the indices of summation should be taken, if indices occurs both on subscript and superscript in one term. For example,\n", "\n", "$$\n", "\\mathscr{W}^{ak}_{ic} \\leftarrow v_{{\\color{orange}{d}} c}^{{\\color{blue}{l}} k} t_{i {\\color{blue}{l}}}^{{\\color{orange}{d}} a}\n", "$$\n", "\n", "both ${\\color{orange}{d}}$ and ${\\color{blue}{l}}$ are in subscripts and super scripts, so summation should be taken on $d, l$. Further more, superscripts $k, a$ and subscripts $c, i$ are still left not summed, so these indices would be on the term at LHS.\n", "\n", "So, the $\\mathscr{W}^{ak}_{ic}$ could be recasted by Einstein summation convention as\n", "\n", "$$\n", "\\mathscr{W}^{ak}_{ic} = v^{ak}_{ic} - v^{kl}_{ci} t_l^a + v^{ka}_{cd} t_i^d - \\frac{1}{2} v^{lk}_{dc} t_{il}^{da} - v^{lk}_{dc} t_i^d t_l^a + \\frac{1}{2} w^{lk}_{dc} t_{il}^{ad}\n", "$$\n", "\n", "In this documentation, we may not use Einstein summation convention in most cases. However, notations are extremely important in deriving and implementing formulas, and this convention can make equations very concise and clear. You may want to try out this convention in your future research or courses.\n", "\n", "````" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "````{admonition} Hint 1: Use deep copy when necessary\n", ":class: dropdown\n", "\n", "A very common **FAULTY** implementation to $\\mathscr{W}^{ab}_{cd}$ could be\n", "\n", "```python\n", "def cc_Wvvvv(mole: Molecule, t1: np.ndarray, t2: np.ndarray) -> np.ndarray:\n", " # Ref: Hitara, eq 43\n", " W_abcd = mole.v_vvvv.transpose((0, 2, 1, 3)) # This line of code is actually WRONG\n", " W_abcd -= np.einsum(\"kdac, kb -> abcd\", mole.v_ovvv, t1, optimize=True)\n", " W_abcd -= np.einsum(\"kcbd, ka -> abcd\", mole.v_ovvv, t1, optimize=True)\n", " return W_abcd\n", "```\n", "\n", "This is due to NumPy's boardcasting nature. When calling NumPy array view (e.g. slicing array, kind of operations we've discussed in Step 2) or transposing, the data of array is actually not copied, saving memory space and CPU time.\n", "\n", "However, it could be risky to use these arrays without double check. In the code above, tensor `mole.eri_mo` will actually be modified, since transpose of `mole.v_vvvv` is view of array, so the underlying data of `W_abcd` comes from shared-data of `mole.eri_mo`. Then substraction operation `-=` will actually directly modify values on `mole.eri_mo`.\n", "\n", "For demonstration, we can use the following code to illustrate this point:\n", "```python\n", ">>> arr_base = np.zeros((3, 4), dtype=int)\n", ">>> arr_view = arr_base[:2, 1:4].transpose((0, 1))\n", ">>> arr_view -= 1\n", ">>> arr_base\n", "array([[ 0, -1, -1, -1],\n", " [ 0, -1, -1, -1],\n", " [ 0, 0, 0, 0]])\n", "```\n", "\n", "So, use `numpy.copy` ([NumPy](https://numpy.org/doc/stable/reference/generated/numpy.copy.html)) to make *deep copy* of array (i.e. make copy of underlying data, different to operation `=`) when necessary, alhough this will cost memory and CPU time.\n", "\n", "````" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "````{admonition} Hint 2: Complicated tensor transpose\n", ":class: dropdown\n", "\n", "We've learned how to use `numpy.transpose` or `numpy.swapaxes` to tensors for simple cases, such as matrix transpose or $w^{pq}_{rs}$ generation. For those transposes, only two indices are swaped.\n", "\n", "However, in this step, multiple indice swaps may occur. What tuple shall we pass into transpose? That could be puzzling.\n", "\n", "For example, we need $\\mathscr{W}^{ak}_{ic} \\leftarrow v^{ak}_{ic}$, but we don't have slice of `v_voov` $v^{ak}_{ic}$. Recall that ERI $v^{pq}_{rs}$ is 8-fold symmetric, so `v_voov` could be obtained from transpose of `v_ovvo` $v^{ic}_{ak}$. So the transpose of `v_ovvo` to `Wvoov` is\n", "\n", "$$\n", "v^{ic}_{ak} \\mapsto \\mathscr{W}^{ak}_{ic}, \\quad (i, a, c, k) \\mapsto (a, k, i, c), \\quad (2, 0, 3, 1) \\mapsto (0, 1, 2, 3)\n", "$$\n", "\n", "So the first line to generate $\\mathscr{W}^{ak}_{ic}$ could be\n", "```python\n", "W_akic = np.copy(mole.v_ovvo.transpose((2, 0, 3, 1)))\n", "```\n", "\n", "However, the code above is putting a tuple $(2, 0, 3, 1)$ into `numpy.transpose`, which could be a little puzzling. Another intutive way could be utilizing `numpy.einsum`:\n", "```python\n", "W_akic = np.copy(np.einsum(\"iack -> akic\", mole.v_ovvo))\n", "```\n", "\n", "Note that when doing transposing, `numpy.einsum` actually does not deep copy the underlying array data, so `numpy.copy` is required.\n", "\n", "````" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Implementation" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "def cc_Woooo(mole: Molecule, t1: np.ndarray, t2: np.ndarray) -> np.ndarray:\n", " # Ref: Hitara, eq 42\n", " raise NotImplementedError(\"About 6 lines of code\")\n", "\n", "Molecule.cc_Woooo = cc_Woooo" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "def cc_Wvvvv(self, t1: np.ndarray, t2: np.ndarray) -> np.ndarray:\n", " # Ref: Hitara, eq 43\n", " raise NotImplementedError(\"About 4 lines of code\")\n", "\n", "Molecule.cc_Wvvvv = cc_Wvvvv" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "def cc_Wvoov(mole: Molecule, t1: np.ndarray, t2: np.ndarray) -> np.ndarray:\n", " # Ref: Hitara, eq 44\n", " # This is a reference implementation. There is no need to change this code.\n", " W_akic = np.copy(mole.v_ovvo.transpose((2, 0, 3, 1)))\n", " W_akic -= np.einsum(\"kcli, la -> akic\", mole.v_ovoo, t1, optimize=True)\n", " W_akic += np.einsum(\"kcad, id -> akic\", mole.v_ovvv, t1, optimize=True)\n", " W_akic -= 0.5 * np.einsum(\"ldkc, ilda -> akic\", mole.v_ovov, t2, optimize=True)\n", " W_akic -= np.einsum(\"ldkc, id, la -> akic\", mole.v_ovov, t1, t1, optimize=True)\n", " W_akic += 0.5 * np.einsum(\"ldkc, ilad -> akic\", mole.w_ovov, t2, optimize=True)\n", " return W_akic\n", "\n", "Molecule.cc_Wvoov = cc_Wvoov" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "def cc_Wvovo(mole: Molecule, t1: np.ndarray, t2: np.ndarray) -> np.ndarray:\n", " # Ref: Hitara, eq 45\n", " raise NotImplementedError(\"About 6 lines of code\")\n", "\n", "Molecule.cc_Wvovo = cc_Wvovo" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solution" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "data": { "text/plain": [ "array([[ 0.0218411, 0.0081018, -0. , 0.0113409, -0. ],\n", " [ 0.0126575, 0.023576 , -0. , 0.6571315, -0. ],\n", " [-0. , 0. , -0.0260943, -0. , 0. ],\n", " [-0.0180776, 0.1264385, -0. , 0.0991944, -0. ],\n", " [ 0. , -0. , 0. , 0. , 0.045617 ]])" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sol_mole.cc_Woooo(t1_precomput, t2_precomput)[1, 3]" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "data": { "text/plain": [ "array([[-0. , 0.55115 ],\n", " [ 0.1111799, 0. ]])" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sol_mole.cc_Wvvvv(t1_precomput, t2_precomput)[0, 1]" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "data": { "text/plain": [ "array([[-0.0102725, 0. ],\n", " [ 0.0947736, 0. ],\n", " [ 0. , -0.0307268],\n", " [-0.0575545, -0. ],\n", " [-0. , 0. ]])" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sol_mole.cc_Wvoov(t1_precomput, t2_precomput)[0, 1]" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "data": { "text/plain": [ "array([[ 0.0074454, 0.5997484, 0. , -0.0504503, -0. ],\n", " [-0. , -0. , -0.0874274, -0. , 0. ]])" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sol_mole.cc_Wvovo(t1_precomput, t2_precomput)[0, 1]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 6: CCSD Single Excitation Amplitude $t_i^a$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this step, we will calculate CCSD single excitation amplitude $t_i^a$ (Hirata, eq 35).\n", "\n", "$$\n", "\\begin{align}\n", "t_i^a D_i^a &= \\sum_c (\\mathscr{F}^a_c - f^a_c) t_i^c - \\sum_k (\\mathscr{F}^k_i - f^k_i) t_k^a + \\sum_{kc} \\mathscr{F}^k_c (2 t_{ki}^{ca} - t_{ik}^{ca}) + \\sum_{kc} \\mathscr{F}^k_c t_i^c t_k^a \\\\\n", "&\\quad + \\sum_{kc} w^{ak}_{ic} t_k^c + \\sum_{kcd} w^{ak}_{cd} t_{ik}^{cd} + \\sum_{kcd} w^{ak}_{cd} t_i^c t_k^d - \\sum_{klc} w^{kl}_{ic} t_{kl}^{ac} - \\sum_{klc} w^{kl}_{ic} t_k^a t_l^c\n", "\\end{align}\n", "$$\n", "\n", "We use notation $\\delta^p_q$ as kronecker delta ($\\delta^p_q = \\delta_{pq}$), and $f^p_q$ Fock matrix element $f^p_q = F_{pq} = \\delta_{pq} \\varepsilon_p$. We define\n", "\n", "$$\n", "D_i^a = \\varepsilon_i - \\varepsilon_a\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "````{admonition} Hint 1: Generate $\\tilde{\\mathscr{F}}^a_c = \\mathscr{F}^a_c - f^a_c$\n", ":class: dropdown\n", "\n", "Since Fock matrix is diagonal, $f^a_c = \\delta_{ac} \\varepsilon_a$. So, in the program, we can generate `Fvvt` $\\tilde{\\mathscr{F}}^a_c$ as\n", "```python\n", "Fvv = mole.cc_Fvv(t1, t2)\n", "Fvvt = Fvv - np.diag(mole.mo_energy[mole.nocc:])\n", "```\n", "\n", "The same process applies to `Foo` $\\tilde{\\mathscr{F}}^k_i = \\mathscr{F}^k_i - f^k_i$. Well, this variable name is not something like `foo`, `bar` or `ket` ;->\n", "\n", "````" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "````{admonition} Hint 2: Generate $D_i^a$\n", ":class: dropdown\n", "\n", "Using numpy boardcasting, we can generate `D_ov` $D_i^a$ as (dim: $(i, a)$)\n", "```python\n", "D_ov = mole.mo_energy[:mole.nocc, None] - mole.mo_energy[None, mole.nocc:]\n", "```\n", "\n", "````" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "````{admonition} Hint 3: Update $t_i^a$\n", ":class: dropdown\n", "\n", "We may simply write the formula above as\n", "\n", "$$\n", "t_i^a D_i^a = \\mathtt{RHS}_i^a\n", "$$\n", "\n", "So $t_i^a$ is updated in form of\n", "\n", "$$\n", "t_i^a = (D_i^a)^{-1} \\mathtt{RHS}_i^a\n", "$$\n", "\n", "You can first generate `RHS` $\\mathtt{RHS}_i^a$ (dim: $(i, a)$) and `D_ov` $D_i^a$ (dim: $(i, a)$), then\n", "```python\n", "return RHS / D_ov\n", "```\n", "\n", "````" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "````{admonition} Hint 4: Tensor transpose without explicit statement\n", ":class: dropdown\n", "\n", "You may remember that there is no `w_voov` attribute for $w_{ic}^{ak}$. However, $w_{ic}^{ak} = w_{ci}^{ka}$ is 2-fold symmetric, so $w_{ic}^{ak}$ could be seen as transpose of `w_ovvo`. However, there is no need explicitly stating transpose, and $\\sum_{kc} w^{ak}_{ic} t_k^c$ could be written as\n", "```python\n", "RHS += np.einsum(\"kcai, kc ->ia\", self.w_ovvo, t1, optimize=True)\n", "```\n", "\n", "The transpose is already imbeded in string `\"kcai\"`.\n", "\n", "````" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "````{admonition} Hint 5: Einstein summation expression\n", ":class: dropdown\n", "\n", "Einstein summation expression can be more concise and clear.\n", "\n", "$$\n", "\\begin{align}\n", "t_i^a D_i^a &= \\tilde{\\mathscr{F}}^a_c t_i^c - \\tilde{\\mathscr{F}}^k_i t_k^a + 2 \\mathscr{F}^k_c t_{ki}^{ca} - \\mathscr{F}^k_c t_{ik}^{ca} + \\mathscr{F}^k_c t_i^c t_k^a \\\\\n", "&\\quad + w^{ak}_{ic} t_k^c + w^{ak}_{cd} t_{ik}^{cd} + w^{ak}_{cd} t_i^c t_k^d - w^{kl}_{ic} t_{kl}^{ac} - w^{kl}_{ic} t_k^a t_l^c\n", "\\end{align}\n", "$$\n", "\n", "````" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Implementation" ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "def update_t1(mole: Molecule, t1: np.ndarray, t2: np.ndarray) -> np.ndarray:\n", " # Input: `t1`, `t2` CCSD amplitudes guess\n", " # Output: Updated `t1` CCSD single amplitude\n", " nocc = mole.nocc\n", " Foo, Fvv, Fov = mole.cc_Foo(t1, t2), mole.cc_Fvv(t1, t2), mole.cc_Fov(t1, t2)\n", " Foot = ____ # Fill this line\n", " Fvvt = ____ # Fill this line\n", " # Formula Line 1\n", " RHS = np.einsum(\"ac, ic -> ia\", Fvvt, t1 , optimize=True)\n", " RHS -= ____ # Fill this line\n", " RHS += 2 * np.einsum(\"kc, kica -> ia\", Fov, t2 , optimize=True)\n", " RHS -= ____ # Fill this line\n", " RHS += ____ # Fill this line\n", " # Formula Line 2\n", " RHS += np.einsum(\"kcai, kc ->ia\", mole.w_ovvo, t1, optimize=True)\n", " RHS += ____ # Fill this line\n", " RHS += ____ # Fill this line\n", " RHS -= ____ # Fill this line\n", " RHS -= ____ # Fill this line\n", "\n", " D_ov = ____ # Fill this line\n", " return ____ # Fill this line\n", "\n", "Molecule.update_t1 = update_t1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solution" ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "data": { "text/plain": [ "array([[-0.0000489, -0. ],\n", " [-0.0032898, -0. ],\n", " [ 0. , -0.0025011],\n", " [-0.0217781, -0. ],\n", " [-0. , 0. ]])" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sol_mole.update_t1(t1_precomput, t2_precomput)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As said before, `t1_precomput` and `t2_precomput` are converged CCSD amplitudes, so actually result of the code-cell above should be extremely close to `t1_precomput`. One can use `numpy.allclose` ([NumPy API](https://numpy.org/doc/stable/reference/generated/numpy.allclose.html)) to check whether two arrays are very close to each other." ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.allclose(sol_mole.update_t1(t1_precomput, t2_precomput), t1_precomput)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 7: CCSD Double Excitation Amplitude $t_{ij}^{ab}$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this step, we will calculate CCSD double excitation amplitude $t_{ij}^{ab}$ (Hirata, eq 36)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{align}\n", "t_{ij}^{ab} D_{ij}^{ab} = \\mathscr{P} (ia, jb) \\bigg[\n", "& \\frac{1}{2} v^{ij}_{ab}\n", " + \\frac{1}{2} \\sum_{kl} \\mathscr{W}^{kl}_{ij} t_{kl}^{ab}\n", " + \\frac{1}{2} \\sum_{kl} \\mathscr{W}^{kl}_{ij} t_k^a t_l^b\n", " + \\frac{1}{2} \\sum_{cd} \\mathscr{W}^{ab}_{cd} t_{ij}^{cd}\n", " + \\frac{1}{2} \\sum_{cd} \\mathscr{W}^{ab}_{cd} t_i^c t_j^d\n", "\\\\\n", "& + \\sum_c (\\mathscr{L}^a_c - f^a_c) t_{ij}^{cb} - \\sum_k (\\mathscr{L}^k_i - f^k_i) t_{kj}^{ab}\n", "\\\\\n", "& + \\sum_c v^{ab}_{ic} t_j^c - \\sum_{kc} v^{kb}_{ic} t_k^a t_j^c\n", " - \\sum_k v^{ak}_{ij} t_k^b - \\sum_{kc} v^{ak}_{ic} t_j^c t_k^b\n", "\\\\\n", "& + 2 \\sum_{kc} \\mathscr{W}^{ak}_{ic} t_{kj}^{cb}\n", " - \\sum_{kc} \\mathscr{W}^{ak}_{ci} t_{kj}^{cb}\n", " - \\sum_{kc} \\mathscr{W}^{ak}_{ic} t_{kj}^{bc}\n", " - \\sum_{kc} \\mathscr{W}^{bk}_{ci} t_{kj}^{ac}\n", "\\bigg]\n", "\\end{align}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Where\n", "\n", "$$\n", "D_{ij}^{ab} = \\varepsilon_i + \\varepsilon_j - \\varepsilon_a - \\varepsilon_b\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and **operator** $\\mathscr{P} (ia, jb)$ is a permutation operator. When apply $\\mathscr{P} (ia, jb)$ to a function (or tensor) $f(i, a, j, b)$, then\n", "\n", "$$\n", "\\big( \\mathscr{P} (ia, jb) \\circ f \\big) (i, a, j, b) = f(i, a, j, b) + f(j, b, i, a)\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "````{admonition} Hint 1: Generate $\\tilde{\\mathscr{L}}^a_c = \\mathscr{L}^a_c - f^a_c$\n", ":class: dropdown\n", "\n", "Generate $\\tilde{\\mathscr{L}}^a_c$ in the same way of $\\tilde{\\mathscr{F}}^a_c$ generation\n", "```python\n", "Lvv = mole.Lvv(t1, t2)\n", "Lvvt = Lvv - np.diag(mole.mo_energy[mole.nocc:])\n", "```\n", "\n", "````" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "````{admonition} Hint 2: Generate $D_{ij}^{ab}$\n", ":class: dropdown\n", "\n", "Generate `D_oovv` $D_{ij}^{ab}$ in the same way of `D_ov` $D_i^a$ or what you've learned from MP2 calculation.\n", "\n", "````" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "````{admonition} Hint 3: Update $t_{ij}^{ab}$\n", ":class: dropdown\n", "\n", "We may write the formula above as\n", "\n", "$$\n", "t_{ij}^{ab} D_{ij}^{ab} = \\mathscr{P} (ia, jb) \\big[ \\mathtt{RHS}_{ij}^{ab} \\big] = \\mathtt{RHS}_{ij}^{ab} + \\mathtt{RHS}_{ji}^{ba}\n", "$$\n", "\n", "Where $\\mathtt{RHS}_{ij}^{ab}$ is not exactly right-hand-side of double excitation amplitude equation, but the sum of terms in square brackets. So $t_{ij}^{ab}$ is updated in form of\n", "\n", "$$\n", "t_{ij}^{ab} = (D_{ij}^{ab})^{-1} \\big( \\mathtt{RHS}_{ij}^{ab} + \\mathtt{RHS}_{ji}^{ba} \\big)\n", "$$\n", "\n", "You can first generate `RHS` $\\mathtt{RHS}_{ij}^{ab}$ (dim: $(i, j, a, b)$) and `D_ov` $D_{i, j}^{a, b}$ (dim: $(i, j, a, b)$), then\n", "```python\n", "return (RHS + RHS.transpose((1, 0, 3, 2))) / D_oovv # permutation of P(ia, jb) applies here\n", "```\n", "\n", "````" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "````{admonition} Einstein summation expression\n", ":class: dropdown\n", "\n", "$$\n", "\\begin{align}\n", "t_{ij}^{ab} D_{ij}^{ab} = \\mathscr{P} (ia, jb) \\bigg[\n", "& \\frac{1}{2} v^{ij}_{ab}\n", " + \\frac{1}{2} \\mathscr{W}^{kl}_{ij} t_{kl}^{ab}\n", " + \\frac{1}{2} \\mathscr{W}^{kl}_{ij} t_k^a t_l^b\n", " + \\frac{1}{2} \\mathscr{W}^{ab}_{cd} t_{ij}^{cd}\n", " + \\frac{1}{2} \\mathscr{W}^{ab}_{cd} t_i^c t_j^d\n", "\\\\\n", "& + \\tilde{\\mathscr{L}}^a_c t_{ij}^{cb} - \\tilde{\\mathscr{L}}^k_i t_{kj}^{ab}\n", "\\\\\n", "& + v^{ab}_{ic} t_j^c - v^{kb}_{ic} t_k^a t_j^c\n", " - v^{ak}_{ij} t_k^b - v^{ak}_{ic} t_j^c t_k^b\n", "\\\\\n", "& + 2 \\mathscr{W}^{ak}_{ic} t_{kj}^{cb}\n", " - \\mathscr{W}^{ak}_{ci} t_{kj}^{cb}\n", " - \\mathscr{W}^{ak}_{ic} t_{kj}^{bc}\n", " - \\mathscr{W}^{bk}_{ci} t_{kj}^{ac}\n", "\\bigg]\n", "\\end{align}\n", "$$\n", "\n", "````" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Implementation" ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "def update_t2(mole: Molecule, t1: np.ndarray, t2: np.ndarray) -> np.ndarray:\n", " # Input: `t1`, `t2` CCSD amplitudes guess\n", " # Output: Updated `t2` CCSD double amplitude\n", " nocc, e = mole.nocc, mole.mo_energy\n", " Loot = ____ # Fill this line\n", " Lvvt = ____ # Fill this line\n", " Woooo = mole.cc_Woooo(t1, t2)\n", " Wvoov = ____ # Fill this line\n", " Wvovo = ____ # Fill this line\n", " Wvvvv = ____ # Fill this line\n", " # Formula Line 1\n", " RHS = ____ # Fill this line\n", " RHS += 0.5 * np.einsum(\"klij, klab -> ijab\", Woooo, t2, optimize=True)\n", " RHS += ____ # Fill this line\n", " RHS += ____ # Fill this line\n", " RHS += ____ # Fill this line\n", " # Formula Line 2\n", " RHS += ____ # Fill this line\n", " RHS -= ____ # Fill this line\n", " # Formula Line 3\n", " RHS += np.einsum(\"iacb, jc -> ijab\", mole.v_ovvv, t1, optimize=True)\n", " RHS -= ____ # Fill this line\n", " RHS -= ____ # Fill this line\n", " RHS -= ____ # Fill this line\n", " # Formula Line 4\n", " RHS += 2 * np.einsum(\"akic, kjcb -> ijab\", Wvoov, t2, optimize=True)\n", " RHS -= ____ # Fill this line\n", " RHS -= ____ # Fill this line\n", " RHS -= ____ # Fill this line\n", "\n", " D_oovv = ____ # Fill this line\n", " return ____ # Fill this line" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solution" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Again, `t1_precomput` and `t2_precomput` are converged CCSD amplitudes, so when feed those tensors to `update_t2`, it should return a tensor very close to `t2_precomput`." ] }, { "cell_type": "code", "execution_count": 31, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.allclose(sol_mole.update_t2(t1_precomput, t2_precomput), t2_precomput)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 8: CCSD Correlation Energy" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Closed-shell CCSD correlation energy is (Hirata, eq 32)\n", "\n", "$$\n", "E_\\mathrm{corr}^\\mathsf{CCSD} = \\sum_{ijab} w^{ij}_{ab} (t_{ij}^{ab} + t_i^a t_j^b)\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Implementation" ] }, { "cell_type": "code", "execution_count": 32, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "def eng_ccsd_corr(mol: Molecule, t1: np.ndarray, t2: np.ndarray) -> float:\n", " # Input: `t1`, `t2` CCSD amplitudes\n", " # Output: (Closed-shell) CCSD correlation energy for given amplitudes (not converged value)\n", " raise NotImplementedError(\"About 3 lines of code\")\n", "\n", "Molecule.eng_ccsd_corr = eng_ccsd_corr" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solution" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since we use the converged `t1_precomput` and `t2_precomput` amplitudes, the result of solution should recast the final CCSD energy result of water/STO-3G." ] }, { "cell_type": "code", "execution_count": 33, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "data": { "text/plain": [ "-0.07068008832822766" ] }, "execution_count": 33, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sol_mole.eng_ccsd_corr(t1_precomput, t2_precomput)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 9: CCSD Amplitudes Initial Guess" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A common initial guess of CCSD amplitudes could be\n", "\n", "$$\n", "t_i^a (\\mathtt{Guess}) = 0, \\quad t_{ij}^{ab} (\\mathtt{Guess}) = \\frac{v_{ij}^{ab}}{D_{ij}^{ab}}\n", "$$\n", "\n", "Feed those amplitudes to CCSD correlation energy evaluation function `eng_ccsd_corr` should recast MP2 correlation energy of water/STO-3G." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Implementation" ] }, { "cell_type": "code", "execution_count": 34, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "def get_t1_t2_initial_guess(mole: Molecule) -> Tuple[np.ndarray, np.ndarray]:\n", " # Output: `t1`, `t2` Initial guess of CCSD amplitudes\n", " raise NotImplementedError(\"About 3~10 lines of code\")\n", "\n", "Molecule.get_t1_t2_initial_guess = get_t1_t2_initial_guess" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solution" ] }, { "cell_type": "code", "execution_count": 35, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "data": { "text/plain": [ "-0.049149636147146854" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sol_mole.eng_ccsd_corr(*sol_mole.get_t1_t2_initial_guess())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 10: CCSD Loop and Convergence" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Just like SCF Loop in [Project 03](../Project_03/Project_03.ipynb#Step-10:-SCF-Loop-and-Convergence), CCSD amplitudes are updated by similar loop.\n", "\n", "Within this function, you may\n", "- Initialize $t_i^a$ and $t_{ij}^{ab}$ guesses;\n", "- Determine maximum iteration number and converge threshold;\n", "- While not exceeding the maximum iteration number,\n", " - Update $t_i^a$ and $t_{ij}^{ab}$;\n", " - Calculate current $E_\\mathrm{corr}^\\mathsf{CCSD}$;\n", " - Print debug information;\n", " - Check whether energy converged;\n", "- Return converged $E_\\mathrm{corr}^\\mathsf{CCSD}$, $t_i^a$ and $t_{ij}^{ab}$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Implementation" ] }, { "cell_type": "code", "execution_count": 36, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "def ccsd_process(mole: Molecule) -> Tuple[float, np.ndarray, np.ndarray]:\n", " # Output: Converged CCSD correlation energy, and density matrix\n", " raise NotImplementedError(\"About 15 lines of code\")\n", "\n", "Molecule.ccsd_process = ccsd_process" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solution" ] }, { "cell_type": "code", "execution_count": 37, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Epoch CCSD Corr Energy\n", " 0 -0.062758205988\n", " 1 -0.067396582633\n", " 2 -0.069224536447\n", " 3 -0.070007757593\n", " 4 -0.070360041940\n", " 5 -0.070523820256\n", " 6 -0.070602032655\n", " 7 -0.070640293065\n", " 8 -0.070659428867\n", " 9 -0.070669194464\n", " 10 -0.070674268086\n", " 11 -0.070676945033\n", " 12 -0.070678375897\n", " 13 -0.070679148925\n", " 14 -0.070679570177\n", " 15 -0.070679801317\n", " 16 -0.070679928834\n", " 17 -0.070679999483\n", " 18 -0.070680038755\n", " 19 -0.070680060641\n", " 20 -0.070680072863\n", " 21 -0.070680079698\n", " 22 -0.070680083526\n", " 23 -0.070680085671\n", " 24 -0.070680086874\n", " 25 -0.070680087549\n", " 26 -0.070680087929\n", " 27 -0.070680088141\n", " 28 -0.070680088261\n", " 29 -0.070680088328\n" ] }, { "data": { "text/plain": [ "-0.07068008832822766" ] }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sol_mole.ccsd_process()[0]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Test Cases" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The input structures, integrals, etc. for these examples are found in the [input directory](https://github.com/ajz34/PyCrawfordProgProj/tree/master/source/Project_05/). You can also use PySCF approach and simply ignore the integral files." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Water/STO-3G** ([Directory](https://github.com/ajz34/PyCrawfordProgProj/tree/master/source/Project_05/input/h2o/STO-3G))" ] }, { "cell_type": "code", "execution_count": 38, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "=== CCSD Iterations ===\n", "Epoch CCSD Corr Energy\n", " 0 -0.062758205988\n", " 1 -0.067396582633\n", " 2 -0.069224536447\n", " 3 -0.070007757593\n", " 4 -0.070360041940\n", " 5 -0.070523820256\n", " 6 -0.070602032655\n", " 7 -0.070640293065\n", " 8 -0.070659428867\n", " 9 -0.070669194464\n", " 10 -0.070674268086\n", " 11 -0.070676945033\n", " 12 -0.070678375897\n", " 13 -0.070679148925\n", " 14 -0.070679570177\n", " 15 -0.070679801317\n", " 16 -0.070679928834\n", " 17 -0.070679999483\n", " 18 -0.070680038755\n", " 19 -0.070680060641\n", " 20 -0.070680072863\n", " 21 -0.070680079698\n", " 22 -0.070680083526\n", " 23 -0.070680085671\n", " 24 -0.070680086874\n", " 25 -0.070680087549\n", " 26 -0.070680087929\n", " 27 -0.070680088141\n", " 28 -0.070680088261\n", " 29 -0.070680088328\n", "=== Final Results ===\n", " MP2 Correlation energy: -0.049149636147\n", "CCSD Correlation energy: -0.070680088328\n", " SCF Total energy: -74.942079928192\n", " MP2 Total energy: -74.991229564340\n", "CCSD Total energy: -75.012760016521\n" ] } ], "source": [ "sol_mole = SolMol()\n", "sol_mole.construct_from_dat_file(\"input/h2o/STO-3G/geom.dat\")\n", "sol_mole.obtain_mol_instance(basis=\"STO-3G\")\n", "sol_mole.print_solution_05()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Water/DZ** ([Directory](https://github.com/ajz34/PyCrawfordProgProj/tree/master/source/Project_05/input/h2o/DZ))" ] }, { "cell_type": "code", "execution_count": 39, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "=== CCSD Iterations ===\n", "Epoch CCSD Corr Energy\n", " 0 -0.153219621576\n", " 1 -0.157583607647\n", " 2 -0.158780675284\n", " 3 -0.159372391280\n", " 4 -0.159624912085\n", " 5 -0.159744546130\n", " 6 -0.159800894374\n", " 7 -0.159828278196\n", " 8 -0.159841745959\n", " 9 -0.159848484681\n", " 10 -0.159851901946\n", " 11 -0.159853658737\n", " 12 -0.159854573259\n", " 13 -0.159855055055\n", " 14 -0.159855311727\n", " 15 -0.159855449898\n", " 16 -0.159855524996\n", " 17 -0.159855566175\n", " 18 -0.159855588937\n", " 19 -0.159855601610\n", " 20 -0.159855608712\n", " 21 -0.159855612715\n", " 22 -0.159855614984\n", " 23 -0.159855616275\n", " 24 -0.159855617013\n", " 25 -0.159855617437\n", " 26 -0.159855617681\n", " 27 -0.159855617821\n", " 28 -0.159855617903\n", "=== Final Results ===\n", " MP2 Correlation energy: -0.152709879014\n", "CCSD Correlation energy: -0.159855617903\n", " SCF Total energy: -75.977878975377\n", " MP2 Total energy: -76.130588854391\n", "CCSD Total energy: -76.137734593279\n" ] } ], "source": [ "sol_mole = SolMol()\n", "sol_mole.construct_from_dat_file(\"input/h2o/DZ/geom.dat\")\n", "sol_mole.obtain_mol_instance(basis=\"DZ\")\n", "sol_mole.print_solution_05()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Water/DZP** ([Directory](https://github.com/ajz34/PyCrawfordProgProj/tree/master/source/Project_05/input/h2o/DZP)) For this molecule, some additional code is included. See [Project 03](../Project_03/Project_03.ipynb#Test-Cases) for details." ] }, { "cell_type": "code", "execution_count": 40, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "=== CCSD Iterations ===\n", "Epoch CCSD Corr Energy\n", " 0 -0.224897568632\n", " 1 -0.229765989748\n", " 2 -0.230690952120\n", " 3 -0.231206877685\n", " 4 -0.231397181235\n", " 5 -0.231490655839\n", " 6 -0.231532173882\n", " 7 -0.231552454210\n", " 8 -0.231562184919\n", " 9 -0.231567038519\n", " 10 -0.231569475283\n", " 11 -0.231570725976\n", " 12 -0.231571376382\n", " 13 -0.231571720337\n", " 14 -0.231571904782\n", " 15 -0.231572005095\n", " 16 -0.231572060344\n", " 17 -0.231572091132\n", " 18 -0.231572108469\n", " 19 -0.231572118323\n", " 20 -0.231572123969\n", " 21 -0.231572127228\n", " 22 -0.231572129120\n", " 23 -0.231572130225\n", " 24 -0.231572130872\n", " 25 -0.231572131253\n", " 26 -0.231572131478\n", " 27 -0.231572131611\n", " 28 -0.231572131690\n", "=== Final Results ===\n", " MP2 Correlation energy: -0.222519233751\n", "CCSD Correlation energy: -0.231572131690\n", " SCF Total energy: -76.008821792901\n", " MP2 Total energy: -76.231341026652\n", "CCSD Total energy: -76.240393924591\n" ] } ], "source": [ "sol_mole = SolMol()\n", "sol_mole.construct_from_dat_file(\"input/h2o/DZP/geom.dat\")\n", "sol_mole.obtain_mol_instance(basis=\"DZP_Dunning\")\n", "sol_mole.mol._basis[\"H\"][2][1][0] = 0.75\n", "sol_mole.mol.cart = True\n", "sol_mole.mol.build()\n", "sol_mole.print_solution_05()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Methane/STO-3G** ([Directory](https://github.com/ajz34/PyCrawfordProgProj/tree/master/source/Project_05/input/ch4/STO-3G))" ] }, { "cell_type": "code", "execution_count": 41, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "=== CCSD Iterations ===\n", "Epoch CCSD Corr Energy\n", " 0 -0.070745262119\n", " 1 -0.075483795485\n", " 2 -0.077200327129\n", " 3 -0.077865986642\n", " 4 -0.078135368845\n", " 5 -0.078247913793\n", " 6 -0.078296204232\n", " 7 -0.078317410766\n", " 8 -0.078326912287\n", " 9 -0.078331242310\n", " 10 -0.078333243407\n", " 11 -0.078334178691\n", " 12 -0.078334619738\n", " 13 -0.078334829162\n", " 14 -0.078334929133\n", " 15 -0.078334977048\n", " 16 -0.078335000083\n", " 17 -0.078335011182\n", " 18 -0.078335016539\n", " 19 -0.078335019128\n", " 20 -0.078335020380\n", " 21 -0.078335020987\n", " 22 -0.078335021280\n", " 23 -0.078335021423\n", " 24 -0.078335021492\n", "=== Final Results ===\n", " MP2 Correlation energy: -0.056046674662\n", "CCSD Correlation energy: -0.078335021492\n", " SCF Total energy: -39.726850316359\n", " MP2 Total energy: -39.782896991021\n", "CCSD Total energy: -39.805185337850\n" ] } ], "source": [ "sol_mole = SolMol()\n", "sol_mole.construct_from_dat_file(\"input/ch4/STO-3G/geom.dat\")\n", "sol_mole.obtain_mol_instance(basis=\"STO-3G\")\n", "sol_mole.print_solution_05()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## References" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- Hirata, S.; Podeszwa, R.; Tobita, M.; Bartlett, R. J. *J. Chem. Phys.* **2004**, *120* (6), 2581\n", "\n", " Coupled-cluster singles and doubles for extended systems\n", "\n", " doi: [10.1063/1.1637577](https://dx.doi.org/10.1063/1.1637577)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.1" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 4 }