05 The Closed-Shell CCSD energy¶

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 (though notations are different, and this project is not discussing extended systems), and PySCF code (rintermediates.py, rccsd.py) (dimension convention is similar to PySCF implementation).

This project will use spacial orbitals (like previous projects, SCF and MP2 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.

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.

# Following os.chdir code is only for thebe (live code), since only in thebe default directory is /home/jovyan
import os
if os.getcwd().split("/")[-1] != "Project_05":
    os.chdir("source/Project_05")
from solution_05 import Molecule as SolMol

from pyscf import gto
import numpy as np
from typing import Tuple
import pickle
np.set_printoptions(precision=7, linewidth=120, suppress=True)

# Solution mol only uses PySCF approach
sol_mole = SolMol()
sol_mole.construct_from_dat_file("input/h2o/STO-3G/geom.dat")
sol_mole.obtain_mol_instance(basis="STO-3G")
sol_mole.obtain_nao()
sol_mole.obtain_nocc()
sol_mole.obtain_eri_ao()

Molecule Object Initialization¶

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 03 and Project 04.

Note that since details of SCF iteration are no longer required, we can make some simplification to SCF code.

class Molecule:
    def __init__(self):
        # Project 03 existed
        self.atom_charges = NotImplemented  # type: np.ndarray
        self.atom_coords = NotImplemented  # type: np.ndarray
        self.natm = NotImplemented  # type: int
        self.mol = NotImplemented  # type: gto.Mole
        self.nao = NotImplemented  # type: int
        self.charge = 0  # type: int
        self.nocc = NotImplemented  # type: int
        # Project 04 added
        self.mo_coeff = NotImplemented  # type: np.ndarray
        self.mo_energy = NotImplemented  # type: np.ndarray
        self.eri_ao = NotImplemented  # type: np.ndarray
        self.eri_mo = NotImplemented  # type: np.ndarray
        self.energy_rhf = NotImplemented  # type: np.ndarray
        self.energy_corr = NotImplemented  # type: np.ndarray
        # Project 05 added
        self.w = NotImplemented  # type: np.ndarray
        self.w_ovoo = NotImplemented  # type: np.ndarray
        self.w_ovvo = NotImplemented  # type: np.ndarray
        self.w_ovov = NotImplemented  # type: np.ndarray
        self.w_ovvv = NotImplemented  # type: np.ndarray
        self.v_oooo = NotImplemented  # type: np.ndarray
        self.v_ovoo = NotImplemented  # type: np.ndarray
        self.v_ovvo = NotImplemented  # type: np.ndarray
        self.v_ovov = NotImplemented  # type: np.ndarray
        self.v_oovv = NotImplemented  # type: np.ndarray
        self.v_ovvv = NotImplemented  # type: np.ndarray
        self.v_vvvv = NotImplemented  # type: np.ndarray

    def construct_from_dat_file(self, file_path: str):
        # Same to Project 01
        with open(file_path, "r") as f:
            dat = np.array([line.split() for line in f.readlines()][1:])
            self.atom_charges = np.array(dat[:, 0], dtype=float).astype(int)
            self.atom_coords = np.array(dat[:, 1:4], dtype=float)
            self.natm = self.atom_charges.shape[0]

    def obtain_mol_instance(self, basis: str, verbose=0):
        # Same to Project 03
        mol = gto.Mole()
        mol.unit = "Bohr"
        mol.atom = "\n".join([("{:3d} " + " ".join(["{:25.18f}"] * 3)).format(chg, *coord) for chg, coord in zip(self.atom_charges, self.atom_coords)])
        mol.basis = basis
        mol.charge = self.charge
        mol.spin = 0
        mol.verbose = verbose
        self.mol = mol.build()

    def obtain_nao(self):
        self.nao = self.mol.nao_nr()

    def obtain_nocc(self):
        assert (self.atom_charges.sum() - self.charge) % 2 == 0
        self.nocc = (self.atom_charges.sum() - self.charge) // 2

    def obtain_eri_ao(self):
        self.eri_ao = self.mol.intor("int2e")

    def get_hcore(self) -> np.ndarray:
        return self.mol.intor("int1e_kin") + self.mol.intor("int1e_nuc")

    def get_fock(self, dm: np.ndarray) -> np.ndarray:
        return self.get_hcore() + (self.eri_ao * dm).sum(axis=(-1, -2)) - 0.5 * (self.eri_ao * dm[:, None, :]).sum(axis=(-1, -3))

    def make_rdm1(self, coeff: np.ndarray) -> np.ndarray:
        return 2 * coeff[:, :self.nocc] @ coeff[:, :self.nocc].T

    def eng_total(self, dm: np.ndarray) -> float:
        return (0.5 * (self.get_hcore() + self.get_fock(dm)) * dm).sum() + self.mol.energy_nuc()

    def scf_process(self, dm_guess: np.ndarray=None) -> Tuple[float, np.ndarray]:
        eng, dm = 0., np.zeros((self.nao, self.nao)) if dm_guess is None else np.copy(dm_guess)
        max_iter, thresh_eng, thresh_dm = 64, 1e-10, 1e-8
        for epoch in range(max_iter):
            eng_next, dm_next = self.eng_total(dm), self.make_rdm1(scipy.linalg.eigh(self.get_fock(dm), self.mol.intor("int1e_ovlp"))[1])
            if np.abs(eng_next - eng) < thresh_eng and np.linalg.norm(dm_next - dm) < thresh_dm:
                eng, dm = eng_next, dm_next
                break
            eng, dm = eng_next, dm_next
        return eng, dm

    def obtain_scf_intermediates(self, dm_guess: np.ndarray=None):
        eng, dm = self.scf_process(dm_guess)
        self.energy_rhf = eng
        self.mo_energy, self.mo_coeff = scipy.linalg.eigh(self.get_fock(dm), self.mol.intor("int1e_ovlp"))

    def get_eri_mo_einsum(self):
        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)

    def obtain_eri_mo(self):
        self.eri_mo = self.get_eri_mo_einsum()

Lots of attributes like w_ovov is introduced in this project. We will illustrate meaning of these attributes later.

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.

Thus, in this project, we may use pre-computed water/STO-3G results to check whether calculation intermediates is correctly implemented.

mo_coeff: pre-computed \(C_{\mu p}\) molecular orbital coefficients
mo_energy: pre-computed \(\varepsilon_p\) molecular orbital energies
t1_precomput: pre-computed \(t_i^a\) CCSD single excitation amplitude
t2_precomput: pre-computed \(t_{ij}^{ab}\) CCSD double excitation amplitude

You may wish to run the following Molecule initialization code after running the code-cell below:

mole = Molecule()
mole.construct_from_dat_file("input/h2o/STO-3G/geom.dat")
mole.obtain_mol_instance(basis="STO-3G")
mole.obtain_nao()
mole.obtain_nocc()
mole.obtain_eri_ao()
mole.mo_coeff = mo_coeff
mole.mo_energy = mo_energy
mole.obtain_eri_mo()

with open("demo_data_h2o_sto3g.dat", "rb") as f:
    d = pickle.load(f)
mo_coeff = d["mo_coeff"]
mo_energy = d["mo_energy"]
t1_precomput, t2_precomput = d["t1"], d["t2"]

sol_mole.mo_coeff = mo_coeff
sol_mole.mo_energy = mo_energy
sol_mole.obtain_eri_mo()

Step 1: ERI and Its Biorthogonal Form¶

In closed-shell CCSD prototype programming, electron repulsion integral (ERI) can be ambiguous. We usually use

\[ v^{pq}_{rs} = (pr|qs) \]

to denote electron repulsion integrals (ERI), as previous projects did before. Here we define biorthogonal electron repulsion integral (biorthogonal ERI) as (Hirata, eq 23)

\[ w^{pq}_{rs} = 2 v^{pq}_{rs} - v^{pq}_{sr} \]

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.

Dimension convention of ERI

Dimension convention to these tensors is the same to Molecule.eri_mo, i.e.

Molecule.eri_mo \(v^{pq}_{rs}\), dim: \((p, r, q, s)\)
Molecule.w \(w^{pq}_{rs}\), dim: \((p, r, q, s)\)

Take caution when transposing those tensors.

Implementation¶

def obtain_w(mole: Molecule):
    # Attribute Modification: `w` biorthogonal electron repulsion integral
    raise NotImplementedError("Exactly 1 line of code")

Molecule.obtain_w = obtain_w

Solution¶

sol_mole.obtain_w()
sol_mole.w[3, 5]

array([[ 0.5311873,  0.0277178,  0.       ,  0.0183313, -0.       , -0.0017709,  0.       ],
       [ 0.0253682,  0.218575 , -0.       ,  0.1602001,  0.       , -0.060252 , -0.       ],
       [ 0.       , -0.       ,  0.1339514,  0.       ,  0.       , -0.       ,  0.1456495],
       [ 0.0012042,  0.0477327, -0.       ,  0.1319441, -0.       ,  0.0955699,  0.       ],
       [-0.       ,  0.       ,  0.       , -0.       ,  0.2900729,  0.       , -0.       ],
       [-0.0215503, -0.0719668, -0.       , -0.352678 ,  0.       ,  0.0524259, -0.       ],
       [ 0.       , -0.       ,  0.1382612,  0.       , -0.       ,  0.       ,  0.1620096]])

Step 2: ERI Slices¶

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.

Variable naming convention

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.

mole.v_ovov = mole.eri_mo[:mole.nocc, mole.nocc:, :mole.nocc, mole.nocc:]

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.

mole.w_ovoo = mole.w[:mole.nocc, mole.nocc:, :mole.nocc, :mole.nocc]

Slice object in Python

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.

so, sv = slice(0, mole.nocc), slice(mole.nocc, mole.nao)

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

mole.v_ovov = mole.eri_mo[so, sv, so, sv]

Implementation¶

def obtain_wv_slices(mole: Molecule):
    # Attribute Modification: Various (biorthogonal) ERI slices
    nocc, nmo = mole.nocc, mole.nao
    so, sv = slice(0, nocc), slice(nocc, nmo)
    mole.w_ovoo = self.w[so, sv, so, so]
    mole.w_ovvo = ____  # Fill this line
    mole.w_ovov = ____  # Fill this line
    mole.w_ovvv = ____  # Fill this line
    mole.v_oooo = ____  # Fill this line
    mole.v_ovoo = ____  # Fill this line
    mole.v_ovvo = ____  # Fill this line
    mole.v_ovov = self.eri_mo[so, sv, so, sv]
    mole.v_oovv = ____  # Fill this line
    mole.v_ovvv = ____  # Fill this line
    mole.v_vvvv = ____  # Fill this line

Molecule.obtain_wv_slices = obtain_wv_slices

Solution¶

sol_mole.obtain_wv_slices()
print(sol_mole.v_ovvo.shape)
print(sol_mole.v_ovvo[3, 1])

(5, 2, 2, 5)
[[ 0.        -0.         0.0541371  0.        -0.       ]
 [-0.0112477  0.0171886  0.         0.0709527  0.       ]]

Step 3: \(\boldsymbol{\mathscr{F}}\) Matrices¶

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\)).

\[\begin{split} \begin{align} \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 \\ \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 \\ \mathscr{F}^k_c &= \sum_{ld} w^{kl}_{cd} t_l^d \end{align} \end{split}\]

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\).

Dimension convention of CCSD amplitudes

Dimension convention of CCSD single excitation amplitude \(t_i^a\) is \((i, a)\):

>>> t1_precomput.shape
(5, 2)

Dimension of CCSD double excitation amplitude \(t_{ij}^{ab}\) is \((i, j, a, b)\):

>>> t2_precomput.shape
(5, 5, 2, 2)

This is quite different to convention of ERIs. So take caution when dealing those tensors!

Implementation¶

def cc_Foo(mole: Molecule, t1: np.ndarray, t2: np.ndarray) -> np.ndarray:
    # Ref: Hitara, eq 37, generate F^k_i
    # This is a reference implementation. There is no need to change this code.
    Fki = np.diag(self.mo_energy[:self.nocc])
    Fki += np.einsum("kcld, ilcd   -> ki", self.w_ovov, t2,     optimize=True)
    Fki += np.einsum("kcld, ic, ld -> ki", self.w_ovov, t1, t1, optimize=True)
    return Fki

Molecule.cc_Foo = cc_Foo

def cc_Fvv(mole: Molecule, t1: np.ndarray, t2: np.ndarray) -> np.ndarray:
    # Ref: Hitara, eq 38
    raise NotImplementedError("No more than 4 lines of code")

Molecule.cc_Fvv = cc_Fvv

def cc_Fov(mole: Molecule, t1: np.ndarray, t2: np.ndarray) -> np.ndarray:
    # Ref: Hitara, eq 39
    # Note that amplitude t2 is actually not taken into account,
    # but for signature consistency, we still include this useless amplitude
    raise NotImplementedError("Exactly 1 line of code")

Molecule.cc_Fov = cc_Fov

Solution¶

sol_mole.cc_Foo(t1_precomput, t2_precomput)

array([[-20.2629403,   0.0004753,  -0.       ,   0.0008499,   0.       ],
       [  0.0000935,  -1.2196514,  -0.       ,   0.014577 ,   0.       ],
       [ -0.       ,  -0.       ,  -0.5825172,   0.       ,   0.       ],
       [  0.0000376,   0.0098341,   0.       ,  -0.4614169,  -0.       ],
       [  0.       ,   0.       ,   0.       ,  -0.       ,  -0.3888219]])

sol_mole.cc_Fvv(t1_precomput, t2_precomput)

array([[ 0.5124191, -0.       ],
       [-0.       ,  0.624019 ]])

sol_mole.cc_Fov(t1_precomput, t2_precomput)

array([[ 0.0001061, -0.       ],
       [ 0.0012685,  0.       ],
       [ 0.       , -0.0032164],
       [-0.0022473, -0.       ],
       [-0.       ,  0.       ]])

Step 4: \(\boldsymbol{\mathscr{L}}\) Matrices¶

In this step, we will calculate \(\boldsymbol{\mathscr{L}}\) matrices (Hirata, eq 40-41).

\[\begin{split} \begin{align} \mathscr{L}^k_i &= \mathscr{F}^k_i + \sum_{lc} w^{lk}_{ci} t_l^c \\ \mathscr{L}^a_c &= \mathscr{F}^a_c + \sum_{kd} w^{ka}_{dc} t_k^d \end{align} \end{split}\]

Dimension convention of \(\mathscr{L}^k_i\) is \((k, i)\). Same to \(\mathscr{L}^a_c\).

Implementation¶

def Loo(mole: Molecule, t1: np.ndarray, t2: np.ndarray) -> np.ndarray:
    # Ref: Hitara, eq 40
    # This is a reference implementation. There is no need to change this code.
    L_ki = self.cc_Foo(t1, t2)
    L_ki += np.einsum("lcki, lc -> ki", self.w_ovoo, t1, optimize=True)
    return L_ki

Molecule.Loo = Loo

def Lvv(mole: Molecule, t1: np.ndarray, t2: np.ndarray) -> np.ndarray:
    # Ref: Hitara, eq 41
    raise NotImplementedError("About 3 lines of code")

Molecule.Lvv = Lvv

Solution¶

sol_mole.Loo(t1_precomput, t2_precomput)

array([[-20.2746448,  -0.0001825,  -0.       ,   0.000496 ,   0.       ],
       [ -0.0005209,  -1.2246573,  -0.       ,   0.0106743,   0.       ],
       [ -0.       ,  -0.       ,  -0.5850956,   0.       ,   0.       ],
       [  0.0001375,   0.0087952,   0.       ,  -0.4644571,  -0.       ],
       [  0.       ,   0.       ,   0.       ,  -0.       ,  -0.3953228]])

sol_mole.Lvv(t1_precomput, t2_precomput)

array([[ 0.5111931, -0.       ],
       [-0.       ,  0.6206725]])

Step 5: \(\boldsymbol{\mathscr{W}}\) Tensors¶

In this step, we will calculate \(\boldsymbol{\mathscr{W}}\) Tensors (Hirata, eq 42-45).

\[\begin{split} \begin{align} \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 \\ \mathscr{W}^{ab}_{cd} &= v^{ab}_{cd} - \sum_k v^{ka}_{dc} t_k^b - \sum_k v^{kb}_{cd} t_k^a\\ \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} \\ \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 \end{align} \end{split}\]

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}\).

Einstein summation convention

We have learned using numpy.einsum. Actually, the original purpose of this function is to make Einstein summation by program. In Einstein summation convention, 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,

\[ \mathscr{W}^{ak}_{ic} \leftarrow v_{{\color{orange}{d}} c}^{{\color{blue}{l}} k} t_{i {\color{blue}{l}}}^{{\color{orange}{d}} a} \]

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.

So, the \(\mathscr{W}^{ak}_{ic}\) could be recasted by Einstein summation convention as

\[ \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} \]

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.

Hint 1: Use deep copy when necessary

A very common FAULTY implementation to \(\mathscr{W}^{ab}_{cd}\) could be

def cc_Wvvvv(mole: Molecule, t1: np.ndarray, t2: np.ndarray) -> np.ndarray:
    # Ref: Hitara, eq 43
    W_abcd = mole.v_vvvv.transpose((0, 2, 1, 3))  # This line of code is actually WRONG
    W_abcd -= np.einsum("kdac, kb -> abcd", mole.v_ovvv, t1, optimize=True)
    W_abcd -= np.einsum("kcbd, ka -> abcd", mole.v_ovvv, t1, optimize=True)
    return W_abcd

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.

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.

For demonstration, we can use the following code to illustrate this point:

>>> arr_base = np.zeros((3, 4), dtype=int)
>>> arr_view = arr_base[:2, 1:4].transpose((0, 1))
>>> arr_view -= 1
>>> arr_base
array([[ 0, -1, -1, -1],
       [ 0, -1, -1, -1],
       [ 0,  0,  0,  0]])

So, use numpy.copy (NumPy) 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.

Hint 2: Complicated tensor transpose

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.

However, in this step, multiple indice swaps may occur. What tuple shall we pass into transpose? That could be puzzling.

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

\[ 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) \]

So the first line to generate \(\mathscr{W}^{ak}_{ic}\) could be

W_akic = np.copy(mole.v_ovvo.transpose((2, 0, 3, 1)))

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:

W_akic = np.copy(np.einsum("iack -> akic", mole.v_ovvo))

Note that when doing transposing, numpy.einsum actually does not deep copy the underlying array data, so numpy.copy is required.

Implementation¶

def cc_Woooo(mole: Molecule, t1: np.ndarray, t2: np.ndarray) -> np.ndarray:
    # Ref: Hitara, eq 42
    raise NotImplementedError("About 6 lines of code")

Molecule.cc_Woooo = cc_Woooo

def cc_Wvvvv(self, t1: np.ndarray, t2: np.ndarray) -> np.ndarray:
    # Ref: Hitara, eq 43
    raise NotImplementedError("About 4 lines of code")

Molecule.cc_Wvvvv = cc_Wvvvv

def cc_Wvoov(mole: Molecule, t1: np.ndarray, t2: np.ndarray) -> np.ndarray:
    # Ref: Hitara, eq 44
    # This is a reference implementation. There is no need to change this code.
    W_akic = np.copy(mole.v_ovvo.transpose((2, 0, 3, 1)))
    W_akic -=       np.einsum("kcli, la     -> akic", mole.v_ovoo, t1,     optimize=True)
    W_akic +=       np.einsum("kcad, id     -> akic", mole.v_ovvv, t1,     optimize=True)
    W_akic -= 0.5 * np.einsum("ldkc, ilda   -> akic", mole.v_ovov, t2,     optimize=True)
    W_akic -=       np.einsum("ldkc, id, la -> akic", mole.v_ovov, t1, t1, optimize=True)
    W_akic += 0.5 * np.einsum("ldkc, ilad   -> akic", mole.w_ovov, t2,     optimize=True)
    return W_akic

Molecule.cc_Wvoov = cc_Wvoov

def cc_Wvovo(mole: Molecule, t1: np.ndarray, t2: np.ndarray) -> np.ndarray:
    # Ref: Hitara, eq 45
    raise NotImplementedError("About 6 lines of code")

Molecule.cc_Wvovo = cc_Wvovo

Solution¶

sol_mole.cc_Woooo(t1_precomput, t2_precomput)[1, 3]

array([[ 0.0218411,  0.0081018, -0.       ,  0.0113409, -0.       ],
       [ 0.0126575,  0.023576 , -0.       ,  0.6571315, -0.       ],
       [-0.       ,  0.       , -0.0260943, -0.       ,  0.       ],
       [-0.0180776,  0.1264385, -0.       ,  0.0991944, -0.       ],
       [ 0.       , -0.       ,  0.       ,  0.       ,  0.045617 ]])

sol_mole.cc_Wvvvv(t1_precomput, t2_precomput)[0, 1]

array([[-0.       ,  0.55115  ],
       [ 0.1111799,  0.       ]])

sol_mole.cc_Wvoov(t1_precomput, t2_precomput)[0, 1]

array([[-0.0102725,  0.       ],
       [ 0.0947736,  0.       ],
       [ 0.       , -0.0307268],
       [-0.0575545, -0.       ],
       [-0.       ,  0.       ]])

sol_mole.cc_Wvovo(t1_precomput, t2_precomput)[0, 1]

array([[ 0.0074454,  0.5997484,  0.       , -0.0504503, -0.       ],
       [-0.       , -0.       , -0.0874274, -0.       ,  0.       ]])

Step 6: CCSD Single Excitation Amplitude \(t_i^a\)¶

In this step, we will calculate CCSD single excitation amplitude \(t_i^a\) (Hirata, eq 35).

\[\begin{split} \begin{align} 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 \\ &\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 \end{align} \end{split}\]

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

\[ D_i^a = \varepsilon_i - \varepsilon_a \]

Hint 1: Generate \(\tilde{\mathscr{F}}^a_c = \mathscr{F}^a_c - f^a_c\)

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

Fvv = mole.cc_Fvv(t1, t2)
Fvvt = Fvv - np.diag(mole.mo_energy[mole.nocc:])

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 ;->

Hint 2: Generate \(D_i^a\)

Using numpy boardcasting, we can generate D_ov \(D_i^a\) as (dim: \((i, a)\))

D_ov = mole.mo_energy[:mole.nocc, None] - mole.mo_energy[None, mole.nocc:]

Hint 3: Update \(t_i^a\)

We may simply write the formula above as

\[ t_i^a D_i^a = \mathtt{RHS}_i^a \]

So \(t_i^a\) is updated in form of

\[ t_i^a = (D_i^a)^{-1} \mathtt{RHS}_i^a \]

You can first generate RHS \(\mathtt{RHS}_i^a\) (dim: \((i, a)\)) and D_ov \(D_i^a\) (dim: \((i, a)\)), then

return RHS / D_ov

Hint 4: Tensor transpose without explicit statement

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

RHS += np.einsum("kcai, kc ->ia", self.w_ovvo, t1, optimize=True)

The transpose is already imbeded in string "kcai".

Hint 5: Einstein summation expression

Einstein summation expression can be more concise and clear.

\[\begin{split} \begin{align} 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 \\ &\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 \end{align} \end{split}\]

Implementation¶

def update_t1(mole: Molecule, t1: np.ndarray, t2: np.ndarray) -> np.ndarray:
    # Input: `t1`, `t2` CCSD amplitudes guess
    # Output: Updated `t1` CCSD single amplitude
    nocc = mole.nocc
    Foo, Fvv, Fov = mole.cc_Foo(t1, t2), mole.cc_Fvv(t1, t2), mole.cc_Fov(t1, t2)
    Foot = ____  #  Fill this line
    Fvvt = ____  #  Fill this line
    # Formula Line 1
    RHS  =     np.einsum("ac, ic     -> ia", Fvvt, t1   , optimize=True)
    RHS -=     ____  # Fill this line
    RHS += 2 * np.einsum("kc, kica   -> ia", Fov, t2    , optimize=True)
    RHS -=     ____  # Fill this line
    RHS +=     ____  # Fill this line
    # Formula Line 2
    RHS +=     np.einsum("kcai, kc     ->ia", mole.w_ovvo, t1,     optimize=True)
    RHS +=     ____  # Fill this line
    RHS +=     ____  # Fill this line
    RHS -=     ____  # Fill this line
    RHS -=     ____  # Fill this line

    D_ov = ____  # Fill this line
    return ____  # Fill this line

Molecule.update_t1 = update_t1

Solution¶

sol_mole.update_t1(t1_precomput, t2_precomput)

array([[-0.0000489, -0.       ],
       [-0.0032898, -0.       ],
       [ 0.       , -0.0025011],
       [-0.0217781, -0.       ],
       [-0.       ,  0.       ]])

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) to check whether two arrays are very close to each other.

np.allclose(sol_mole.update_t1(t1_precomput, t2_precomput), t1_precomput)

True

Step 7: CCSD Double Excitation Amplitude \(t_{ij}^{ab}\)¶

In this step, we will calculate CCSD double excitation amplitude \(t_{ij}^{ab}\) (Hirata, eq 36).

\[\begin{split} \begin{align} t_{ij}^{ab} D_{ij}^{ab} = \mathscr{P} (ia, jb) \bigg[ & \frac{1}{2} v^{ij}_{ab} + \frac{1}{2} \sum_{kl} \mathscr{W}^{kl}_{ij} t_{kl}^{ab} + \frac{1}{2} \sum_{kl} \mathscr{W}^{kl}_{ij} t_k^a t_l^b + \frac{1}{2} \sum_{cd} \mathscr{W}^{ab}_{cd} t_{ij}^{cd} + \frac{1}{2} \sum_{cd} \mathscr{W}^{ab}_{cd} t_i^c t_j^d \\ & + \sum_c (\mathscr{L}^a_c - f^a_c) t_{ij}^{cb} - \sum_k (\mathscr{L}^k_i - f^k_i) t_{kj}^{ab} \\ & + \sum_c v^{ab}_{ic} t_j^c - \sum_{kc} v^{kb}_{ic} t_k^a t_j^c - \sum_k v^{ak}_{ij} t_k^b - \sum_{kc} v^{ak}_{ic} t_j^c t_k^b \\ & + 2 \sum_{kc} \mathscr{W}^{ak}_{ic} t_{kj}^{cb} - \sum_{kc} \mathscr{W}^{ak}_{ci} t_{kj}^{cb} - \sum_{kc} \mathscr{W}^{ak}_{ic} t_{kj}^{bc} - \sum_{kc} \mathscr{W}^{bk}_{ci} t_{kj}^{ac} \bigg] \end{align} \end{split}\]

Where

\[ D_{ij}^{ab} = \varepsilon_i + \varepsilon_j - \varepsilon_a - \varepsilon_b \]

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

\[ \big( \mathscr{P} (ia, jb) \circ f \big) (i, a, j, b) = f(i, a, j, b) + f(j, b, i, a) \]

Hint 1: Generate \(\tilde{\mathscr{L}}^a_c = \mathscr{L}^a_c - f^a_c\)

Generate \(\tilde{\mathscr{L}}^a_c\) in the same way of \(\tilde{\mathscr{F}}^a_c\) generation

Lvv = mole.Lvv(t1, t2)
Lvvt = Lvv - np.diag(mole.mo_energy[mole.nocc:])

Hint 2: Generate \(D_{ij}^{ab}\)

Generate D_oovv \(D_{ij}^{ab}\) in the same way of D_ov \(D_i^a\) or what you’ve learned from MP2 calculation.

Hint 3: Update \(t_{ij}^{ab}\)

We may write the formula above as

\[ t_{ij}^{ab} D_{ij}^{ab} = \mathscr{P} (ia, jb) \big[ \mathtt{RHS}_{ij}^{ab} \big] = \mathtt{RHS}_{ij}^{ab} + \mathtt{RHS}_{ji}^{ba} \]

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

\[ t_{ij}^{ab} = (D_{ij}^{ab})^{-1} \big( \mathtt{RHS}_{ij}^{ab} + \mathtt{RHS}_{ji}^{ba} \big) \]

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

return (RHS + RHS.transpose((1, 0, 3, 2))) / D_oovv  # permutation of P(ia, jb) applies here

Einstein summation expression

\[\begin{split} \begin{align} t_{ij}^{ab} D_{ij}^{ab} = \mathscr{P} (ia, jb) \bigg[ & \frac{1}{2} v^{ij}_{ab} + \frac{1}{2} \mathscr{W}^{kl}_{ij} t_{kl}^{ab} + \frac{1}{2} \mathscr{W}^{kl}_{ij} t_k^a t_l^b + \frac{1}{2} \mathscr{W}^{ab}_{cd} t_{ij}^{cd} + \frac{1}{2} \mathscr{W}^{ab}_{cd} t_i^c t_j^d \\ & + \tilde{\mathscr{L}}^a_c t_{ij}^{cb} - \tilde{\mathscr{L}}^k_i t_{kj}^{ab} \\ & + v^{ab}_{ic} t_j^c - v^{kb}_{ic} t_k^a t_j^c - v^{ak}_{ij} t_k^b - v^{ak}_{ic} t_j^c t_k^b \\ & + 2 \mathscr{W}^{ak}_{ic} t_{kj}^{cb} - \mathscr{W}^{ak}_{ci} t_{kj}^{cb} - \mathscr{W}^{ak}_{ic} t_{kj}^{bc} - \mathscr{W}^{bk}_{ci} t_{kj}^{ac} \bigg] \end{align} \end{split}\]

Implementation¶

def update_t2(mole: Molecule, t1: np.ndarray, t2: np.ndarray) -> np.ndarray:
    # Input: `t1`, `t2` CCSD amplitudes guess
    # Output: Updated `t2` CCSD double amplitude
    nocc, e = mole.nocc, mole.mo_energy
    Loot = ____  #  Fill this line
    Lvvt = ____  #  Fill this line
    Woooo = mole.cc_Woooo(t1, t2)
    Wvoov = ____  #  Fill this line
    Wvovo = ____  #  Fill this line
    Wvvvv = ____  #  Fill this line
    # Formula Line 1
    RHS = ____  #  Fill this line
    RHS += 0.5 * np.einsum("klij, klab   -> ijab", Woooo, t2,     optimize=True)
    RHS += ____  #  Fill this line
    RHS += ____  #  Fill this line
    RHS += ____  #  Fill this line
    # Formula Line 2
    RHS += ____  #  Fill this line
    RHS -= ____  #  Fill this line
    # Formula Line 3
    RHS += np.einsum("iacb, jc     -> ijab", mole.v_ovvv, t1,     optimize=True)
    RHS -= ____  #  Fill this line
    RHS -= ____  #  Fill this line
    RHS -= ____  #  Fill this line
    # Formula Line 4
    RHS += 2 * np.einsum("akic, kjcb -> ijab", Wvoov, t2, optimize=True)
    RHS -= ____  #  Fill this line
    RHS -= ____  #  Fill this line
    RHS -= ____  #  Fill this line

    D_oovv = ____  #  Fill this line
    return ____  #  Fill this line

Solution¶

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.

np.allclose(sol_mole.update_t2(t1_precomput, t2_precomput), t2_precomput)

True

Step 8: CCSD Correlation Energy¶

Closed-shell CCSD correlation energy is (Hirata, eq 32)

\[ E_\mathrm{corr}^\mathsf{CCSD} = \sum_{ijab} w^{ij}_{ab} (t_{ij}^{ab} + t_i^a t_j^b) \]

Implementation¶

def eng_ccsd_corr(mol: Molecule, t1: np.ndarray, t2: np.ndarray) -> float:
    # Input: `t1`, `t2` CCSD amplitudes
    # Output: (Closed-shell) CCSD correlation energy for given amplitudes (not converged value)
    raise NotImplementedError("About 3 lines of code")

Molecule.eng_ccsd_corr = eng_ccsd_corr

Solution¶

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.

sol_mole.eng_ccsd_corr(t1_precomput, t2_precomput)

-0.07068008832822766

Step 9: CCSD Amplitudes Initial Guess¶

A common initial guess of CCSD amplitudes could be

\[ t_i^a (\mathtt{Guess}) = 0, \quad t_{ij}^{ab} (\mathtt{Guess}) = \frac{v_{ij}^{ab}}{D_{ij}^{ab}} \]

Feed those amplitudes to CCSD correlation energy evaluation function eng_ccsd_corr should recast MP2 correlation energy of water/STO-3G.

Implementation¶

def get_t1_t2_initial_guess(mole: Molecule) -> Tuple[np.ndarray, np.ndarray]:
    # Output: `t1`, `t2` Initial guess of CCSD amplitudes
    raise NotImplementedError("About 3~10 lines of code")

Molecule.get_t1_t2_initial_guess = get_t1_t2_initial_guess

Solution¶

sol_mole.eng_ccsd_corr(*sol_mole.get_t1_t2_initial_guess())

-0.049149636147146854

Step 10: CCSD Loop and Convergence¶

Just like SCF Loop in Project 03, CCSD amplitudes are updated by similar loop.

Within this function, you may

Initialize \(t_i^a\) and \(t_{ij}^{ab}\) guesses;
Determine maximum iteration number and converge threshold;
While not exceeding the maximum iteration number,
- Update \(t_i^a\) and \(t_{ij}^{ab}\);
- Calculate current \(E_\mathrm{corr}^\mathsf{CCSD}\);
- Print debug information;
- Check whether energy converged;
Return converged \(E_\mathrm{corr}^\mathsf{CCSD}\), \(t_i^a\) and \(t_{ij}^{ab}\).

Implementation¶

def ccsd_process(mole: Molecule) -> Tuple[float, np.ndarray, np.ndarray]:
    # Output: Converged CCSD correlation energy, and density matrix
    raise NotImplementedError("About 15 lines of code")

Molecule.ccsd_process = ccsd_process

Solution¶

sol_mole.ccsd_process()[0]

Epoch     CCSD Corr Energy
    -0.062758205988
    -0.067396582633
    -0.069224536447
    -0.070007757593
    -0.070360041940
    -0.070523820256
    -0.070602032655
    -0.070640293065
    -0.070659428867
    -0.070669194464
    -0.070674268086
    -0.070676945033
    -0.070678375897
    -0.070679148925
    -0.070679570177
    -0.070679801317
    -0.070679928834
    -0.070679999483
    -0.070680038755
    -0.070680060641
    -0.070680072863
    -0.070680079698
    -0.070680083526
    -0.070680085671
    -0.070680086874
    -0.070680087549
    -0.070680087929
    -0.070680088141
    -0.070680088261
    -0.070680088328

-0.07068008832822766

Test Cases¶

The input structures, integrals, etc. for these examples are found in the input directory. You can also use PySCF approach and simply ignore the integral files.

Water/STO-3G (Directory)

sol_mole = SolMol()
sol_mole.construct_from_dat_file("input/h2o/STO-3G/geom.dat")
sol_mole.obtain_mol_instance(basis="STO-3G")
sol_mole.print_solution_05()

=== CCSD Iterations ===
Epoch     CCSD Corr Energy
    -0.062758205988
    -0.067396582633
    -0.069224536447
    -0.070007757593
    -0.070360041940
    -0.070523820256
    -0.070602032655
    -0.070640293065
    -0.070659428867
    -0.070669194464
    -0.070674268086
    -0.070676945033
    -0.070678375897
    -0.070679148925
    -0.070679570177
    -0.070679801317
    -0.070679928834
    -0.070679999483
    -0.070680038755
    -0.070680060641
    -0.070680072863
    -0.070680079698
    -0.070680083526
    -0.070680085671
    -0.070680086874
    -0.070680087549
    -0.070680087929
    -0.070680088141
    -0.070680088261
    -0.070680088328
=== Final Results ===
 MP2 Correlation energy:      -0.049149636147
CCSD Correlation energy:      -0.070680088328
 SCF       Total energy:     -74.942079928192
 MP2       Total energy:     -74.991229564340
CCSD       Total energy:     -75.012760016521

Water/DZ (Directory)

sol_mole = SolMol()
sol_mole.construct_from_dat_file("input/h2o/DZ/geom.dat")
sol_mole.obtain_mol_instance(basis="DZ")
sol_mole.print_solution_05()

=== CCSD Iterations ===
Epoch     CCSD Corr Energy
    -0.153219621576
    -0.157583607647
    -0.158780675284
    -0.159372391280
    -0.159624912085
    -0.159744546130
    -0.159800894374
    -0.159828278196
    -0.159841745959
    -0.159848484681
    -0.159851901946
    -0.159853658737
    -0.159854573259
    -0.159855055055
    -0.159855311727
    -0.159855449898
    -0.159855524996
    -0.159855566175
    -0.159855588937
    -0.159855601610
    -0.159855608712
    -0.159855612715
    -0.159855614984
    -0.159855616275
    -0.159855617013
    -0.159855617437
    -0.159855617681
    -0.159855617821
    -0.159855617903
=== Final Results ===
 MP2 Correlation energy:      -0.152709879014
CCSD Correlation energy:      -0.159855617903
 SCF       Total energy:     -75.977878975377
 MP2       Total energy:     -76.130588854391
CCSD       Total energy:     -76.137734593279

Water/DZP (Directory) For this molecule, some additional code is included. See Project 03 for details.

sol_mole = SolMol()
sol_mole.construct_from_dat_file("input/h2o/DZP/geom.dat")
sol_mole.obtain_mol_instance(basis="DZP_Dunning")
sol_mole.mol._basis["H"][2][1][0] = 0.75
sol_mole.mol.cart = True
sol_mole.mol.build()
sol_mole.print_solution_05()

=== CCSD Iterations ===
Epoch     CCSD Corr Energy
    -0.224897568632
    -0.229765989748
    -0.230690952120
    -0.231206877685
    -0.231397181235
    -0.231490655839
    -0.231532173882
    -0.231552454210
    -0.231562184919
    -0.231567038519
    -0.231569475283
    -0.231570725976
    -0.231571376382
    -0.231571720337
    -0.231571904782
    -0.231572005095
    -0.231572060344
    -0.231572091132
    -0.231572108469
    -0.231572118323
    -0.231572123969
    -0.231572127228
    -0.231572129120
    -0.231572130225
    -0.231572130872
    -0.231572131253
    -0.231572131478
    -0.231572131611
    -0.231572131690
=== Final Results ===
 MP2 Correlation energy:      -0.222519233751
CCSD Correlation energy:      -0.231572131690
 SCF       Total energy:     -76.008821792901
 MP2       Total energy:     -76.231341026652
CCSD       Total energy:     -76.240393924591

Methane/STO-3G (Directory)

sol_mole = SolMol()
sol_mole.construct_from_dat_file("input/ch4/STO-3G/geom.dat")
sol_mole.obtain_mol_instance(basis="STO-3G")
sol_mole.print_solution_05()

=== CCSD Iterations ===
Epoch     CCSD Corr Energy
    -0.070745262119
    -0.075483795485
    -0.077200327129
    -0.077865986642
    -0.078135368845
    -0.078247913793
    -0.078296204232
    -0.078317410766
    -0.078326912287
    -0.078331242310
    -0.078333243407
    -0.078334178691
    -0.078334619738
    -0.078334829162
    -0.078334929133
    -0.078334977048
    -0.078335000083
    -0.078335011182
    -0.078335016539
    -0.078335019128
    -0.078335020380
    -0.078335020987
    -0.078335021280
    -0.078335021423
    -0.078335021492
=== Final Results ===
 MP2 Correlation energy:      -0.056046674662
CCSD Correlation energy:      -0.078335021492
 SCF       Total energy:     -39.726850316359
 MP2       Total energy:     -39.782896991021
CCSD       Total energy:     -39.805185337850

References¶

Hirata, S.; Podeszwa, R.; Tobita, M.; Bartlett, R. J. J. Chem. Phys. 2004, 120 (6), 2581

Coupled-cluster singles and doubles for extended systems

doi: 10.1063/1.1637577

Programming Tutorial in Chemistry by Python documentation

05 The Closed-Shell CCSD energy¶

Molecule Object Initialization¶

Step 1: ERI and Its Biorthogonal Form¶

Implementation¶

Solution¶

Step 2: ERI Slices¶

Implementation¶

Solution¶

Step 3: \(\boldsymbol{\mathscr{F}}\) Matrices¶

Implementation¶

Solution¶

Step 4: \(\boldsymbol{\mathscr{L}}\) Matrices¶

Implementation¶

Solution¶

Step 5: \(\boldsymbol{\mathscr{W}}\) Tensors¶

Implementation¶

Solution¶

Step 6: CCSD Single Excitation Amplitude \(t_i^a\)¶

Implementation¶

Solution¶

Step 7: CCSD Double Excitation Amplitude \(t_{ij}^{ab}\)¶

Implementation¶

Solution¶

Step 8: CCSD Correlation Energy¶

Implementation¶

Solution¶

Step 9: CCSD Amplitudes Initial Guess¶

Implementation¶

Solution¶

Step 10: CCSD Loop and Convergence¶

Implementation¶

Solution¶

Test Cases¶

References¶