Jupyter Notebook Tutorials: Quantum Chemistry

cis

Configuration Interaction with Single Excitations¶

The Jupyter notebook for this tutorial can be found here.

In the preceding tutorials, we have explored strategies for finding approximations to the lowest-energy eigenfunctions of the electronic Hamiltonian (using Hartree-Fock theory and many-body perturbation theory). In this project, we explore the simplest possible wave function ansatz for modeling electronically excited states, which is known as configuration interaction with single excitations (CIS). We will develop a CIS code using the Python interface to the Psi4 electronic structure package. Most of Psi4's core libraries (integral generation, SCF solver, etc.) are exposed to Python, which facilitates the development of short, easy to understand codes for complex electronic structure methods.

The CIS Wave Function and Hamiltonian¶

At the Hartree-Fock level of theory, the $N$-electron wave function is approximated as an antisymmetrized product of $N$ one-electron functions called molecular orbitals (MOs) (a Slater determinant). The simplest representation of an excited-state wave function would be a single Slater determinant comprised of a different set of MOs, but we can obtain a slightly better description of the excited state by expressing it as a linear combination of Slater determinants that differ by the ground-state configuration by one electron (i.e., singly-excited confugurations). What we are describing here is a configuration interaction with single excitations (CIS) wave function; the CIS wave function for the $n$th excited state can be expressed as

$$\begin{align} |\psi(n)\rangle = \sum_{ia} c_i^a(n) \hat{a}^\dagger_a\hat{a}_i |\Phi_0\rangle = \sum_{ia} c_i^a(n) |\Phi_i^a\rangle \end{align}$$

Here, the symbol $|\Phi_0\rangle$ represents the Hartree-Fock configuration, and $|\Phi_i^a\rangle$ is a Slater determinant that is singly excited relative $|\Phi_0\rangle$. The indices $i$ and $a$ represent spin orbitals that are occupied and unoccupied (virtual) in the reference, respectively.

How do we determine the expansion coefficients? The variation theorem tells us that we can determine these coefficients variationally, that is, by minimizing the expectation value of the Hamiltonian with respect to variations in $c_i^a$. Recall that, when the Hamiltonian is expanded within a basis of know functions (in this case, the singly-excited configurations), then the variational problem can be expressed as a matrix eigenvalue problem, $i.e.$,

$$\begin{align} {\bf H}{\bf c}(n) = E_n {\bf c}(n) \end{align}$$

where $E_n$ represents the energy of state $n$ ($n > 0$). Here, ${\bf H}$ is the matrix representation of the Hamiltonian expanded within the basis of singly-excited configurations, with elements

$$\begin{align} H_{ia,jb} = \langle \Phi_i^a | \hat{H} | \Phi_j^b \rangle \end{align}$$

We can use the p$^\dagger$q package to evaluate this matrix element using the following code.

In [1]:

import pdaggerq

pq = pdaggerq.pq_helper('fermi')

pq.set_left_operators([['a*(i)', 'a(a)']])
pq.set_right_operators([['a*(b)', 'a(j)']])

pq.add_operator_product(1.0, ['f'])
pq.add_operator_product(1.0, ['v'])

pq.simplify()

terms = pq.strings()
for term in terms:
    print(term)

['+1.00', 'd(a,b)', 'd(i,j)', 'f(k,k)']
['-1.00', 'd(a,b)', 'f(j,i)']
['+1.00', 'd(i,j)', 'f(a,b)']
['-0.50', 'd(a,b)', 'd(i,j)', '<l,k||l,k>']
['+1.00', '<j,a||b,i>']

In the code snippet, 'f' and 'v' refer to the Fock operator ($\hat{f}$) and fluctuation potential ($\hat{v}$), respectively, and $\hat{H} = \hat{f} + \hat{v}$. In the output, 'd(a,b)' is a Kronecker delta function, 'f(a,b)' is an element of the Fock matrix, and '<j,a||b,i>' refers to an antisymmetrized electron repulsion integral (ERI) in physicists' notation. So, we can see that the matrix elements of the CIS Hamiltonian are

$$\begin{align} H_{ia,jb} = f_{ab}\delta_{ij}-f_{ij}\delta_{ab} + \langle aj ||ib\rangle + \delta_{ij}\delta_{ab} \left ( \sum_k f_{kk} - \frac{1}{2} \sum_{kl} \langle kl||kl \rangle \right ) \end{align}$$

We can simplify this expression in two ways. First, we recall that the Fock matrix for Hartree-Fock theory is diagonal in the MO basis, and the diagonal elements are the orbital energies ($f_{pp} = \epsilon_p$). Second, note that

$$\begin{align} E_0 = \sum_k f_{kk} - \frac{1}{2} \sum_{kl} \langle kl||kl \rangle = \sum_k \epsilon_k - \frac{1}{2} \sum_{kl} \langle kl||kl \rangle \end{align}$$

is the Hartree-Fock energy. We now have

$$\begin{align} H_{ia,jb} = (\epsilon_a - \epsilon_i + E_0) \delta_{ij}\delta_{ab} + \langle aj ||ib\rangle \end{align}$$

Singlet Excitation Energies¶

Let us assume that $|\Phi_0\rangle$ represents a closed-shell restricted Hartree-Fock determinant. In this case, it turns out that the eigenvalues of the CIS Hamiltonian represented in the spin-orbital basis correspond to energies for CIS states of two different spin symmetries: singlets and triplets. If we are only interested in the singlet states, then we can spin adapt our basis functions and represent the Hamiltonian matrix in this new spin-adapted basis. The singlet spin-adapted basis functions are given by

$$\begin{align} |\Phi_i^a\rangle = \frac{1}{\sqrt{2}} \left ( |\Phi^{a_\alpha}_{i_\alpha}\rangle + |\Phi^{a_\beta}_{i_\beta}\rangle \right ) \end{align}$$

and the corresponding matrix elements of the Hamiltonian are

$$\begin{align} H_{ia,jb} = \frac{1}{2} \left ( \langle \Phi^{a_\alpha}_{i_\alpha} | \hat{H} | \Phi_{j_\alpha}^{b_\alpha} \rangle + \langle \Phi^{a_\alpha}_{i_\alpha} | \hat{H} | \Phi_{j_\beta}^{b_\beta} \rangle + \langle \Phi^{a_\beta}_{i_\beta} | \hat{H} | \Phi_{j_\alpha}^{b_\alpha} \rangle +\langle \Phi^{a_\beta}_{i_\beta} | \hat{H} | \Phi_{j_\beta}^{b_\beta} \rangle \right ) \end{align}$$

Comparing to the expression for the Hamiltonian in the spin-orbital basis generated by p$^\dagger$q, we note the following:

Kronecker delta functions involving orbitals with different spin-symmetry labels will be zero.
The exchange part of the antisymmetrized ERIs vanish for terms involving different spin-symmetry labels, e.g., $\langle a_\alpha j_\beta || i_\alpha b_\beta \rangle = \langle a_\alpha j_\beta | i_\alpha b_\beta \rangle$
In restricted Hartree-Fock theory, the $\alpha$- and $\beta$-spin orbitals and orbital energies are equivalent.

With these considerations in mind, we find that

$$\begin{align} H_{ia,jb} = (\epsilon_a - \epsilon_i + E_0) \delta_{ij}\delta_{ab} + 2 \langle aj |ib\rangle - \langle aj | bi \rangle \end{align}$$

Note that the orbital labels $i$, $j$, $a,$ and $b$ no longer represent spin orbitals; they represent spatial orbitals. Using chemists' notation, these matrix elements are

$$\begin{align} H_{ia,jb} = (\epsilon_a - \epsilon_i + E_0) \delta_{ij}\delta_{ab} + 2 ( ai |jb ) - ( ab | ji ) \end{align}$$

The benefits of spin-adaptation here are two-fold. First, the dimension of the Hamiltonian is reduced by a factor of two. Second, because matrix diagonalization scales as the third power of matrix size, the cost of the diagonalization procedure is reduced by a factor of eight!

With the theory out of the way, we can begin developing Python code that diagonalizes the CIS Hamiltonian. First, import the Psi4 and numpy libraries, create a molecule object, and set some basic options.

In [2]:

import psi4
import numpy as np

# set molecule
mol = psi4.geometry("""
o
h 1 1.0
h 1 1.0 2 104.5
symmetry c1
""")   

psi4.set_options({'basis': 'sto-3g',
                  'scf_type': 'pk',
                  'e_convergence': 1e-8,
                  'd_convergence': 1e-8})

psi4.core.be_quiet()

Now, we are ready to perform a Hartree-Fock computation to determine our reference electronic configuration:

In [3]:

# compute the Hartree-Fock energy and wavefunction
scf_e, wfn = psi4.energy('SCF', return_wfn=True)

Note that we asked the energy routine to return a wavefunction. This object contains all of the important information from the Hartree-Fock computation, including the number of electrons, the number of orbitals, the AO/MO transformation matrices, and more! We will need some of this information for our CIS routine:

In [4]:

# Grab data from wavfunction

# number of doubly occupied orbitals (for RHF, nalpha = nbeta = ndocc)
ndocc   = wfn.nalpha()

# total number of orbitals
nmo     = wfn.nmo()

# number of virtual orbitals
nvirt   = nmo - ndocc

# orbital energies
eps     = np.asarray(wfn.epsilon_a())

# occupied orbitals:
Co = wfn.Ca_subset("AO", "OCC")

# virtual orbitals:
Cv = wfn.Ca_subset("AO", "VIR")

Now, we can use Psi4's MintsHelper class to generate the two-electron integrals that we need. The CIS Hamiltonian requires two classes of integrals of the type (ov|ov) and (oo|vv), where "o" represents an orbital that is occupied in the reference function, and "v" represents a virtual orbital. The MintsHelper class can construct tensors containing these specific classes of orbitals, provided we provide to it the corresponding definitions of the molecular orbitals (given by the Co and Cv matrices above):

In [5]:

# use Psi4's MintsHelper to generate ERIs
mints = psi4.core.MintsHelper(wfn.basisset())

# build the (ov|ov) integrals:
g_ovov = np.asarray(mints.mo_eri(Co, Cv, Co, Cv))

# build the (oo|vv) integrals:
g_oovv = np.asarray(mints.mo_eri(Co, Co, Cv, Cv))

Given these tensors, you can access the element (ij|ab) in Python as g_oovv[i,j,a,b]. Here, the indices i and j could run from 0 to ndocc-1, and the indices a and b could run from 0 to nvirt-1. For convenience, we can also strip out the occupied and virtual orbital energies from the full list of orbital energies that we obtained from the wave function object above:

In [6]:

# strip out occupied orbital energies
eps_o = eps[:ndocc]

# strip out virtual orbital energies
eps_v = eps[ndocc:]

Now, we have all of the ingredients to build the CIS Hamiltonian! Using the expressions above, build the matrix representation of the Hamiltonian in the basis of spin-adapted singly-excited functions. The dimensions of this Hamiltonian should be ndocc x nvirt by ndocc x nvirt:

In [7]:

Ham = np.zeros((ndocc*nvirt,ndocc*nvirt))

The composite index that represents each basis function (|$\Phi_i^a\rangle$), can be expressed in terms of the orbitals that define that basis function. For example, the index ia = i*nvirt + a, where i=0..ndocc-1, and a=0..nvirt-1.

In [8]:

# build singlet hamiltonian
for i in range(0,ndocc):
    for a in range(0,nvirt):
        ia = i * nvirt + a
        for j in range(0,ndocc):
            for b in range(0,nvirt):
                jb = j * nvirt + b
                Ham[ia][jb] = 2.0 * g_ovov[i][a][j][b] - g_oovv[i][j][a][b]
        Ham[ia][ia] += eps_v[a] - eps_o[i]

# diagonalize Hamiltonian
eig = np.linalg.eigvals(Ham)

# sort excitation energies
eig.sort()

print("")
print("    ==> CIS singlet excitation energies (eV) <==")
print("")
for ia in range(0,ndocc*nvirt):
    print("    %5i %10.5f" % (ia,eig[ia]*27.21138))
print("")

    ==> CIS singlet excitation energies (eV) <==

        0   12.03295
        1   13.89434
        2   15.79662
        3   17.88952
        4   20.69644
        5   27.65951
        6   38.62665
        7   39.46291
        8  546.26101
        9  547.54054

Note that we did not add the reference energy to the diagonal elements of the Hamiltonian. This choice has the effect of shifting the eigenvalues of the Hamiltonian by this constant. In other words, we have modified the original eigenvalue problem

$$\begin{align} {\bf H}{\bf c}(n) = E_n {\bf c}(n) \end{align}$$

such that the eigenvalues we compute correspond to excitation energies rather than total energies, $i.e.$,

$$\begin{align} ({\bf H} - E_0){\bf c}(n) = (E_n-E_0) {\bf c}(n) = \omega_n {\bf c}(n) \end{align}$$

Triplet Excitation Energies¶

We can also define spin-adapted basis functions for triplet spin states. The triplet spin-adapted basis functions are given by

$$ |\Phi_i^a\rangle = \frac{1}{\sqrt{2}} \left ( |\Phi^{a_\alpha}_{i_\alpha}\rangle - |\Phi^{a_\beta}_{i_\beta}\rangle \right ) $$

and the corresponding matrix elements of the Hamiltonian are

$$H_{ia,jb} = (\epsilon_a - \epsilon_i + E_0) \delta_{ij}\delta_{ab} - ( ab | ji )$$

Now, we can build and diagonalize the Hamiltonian expanded in the basis of triplet spin-adapted functions:

In [9]:

# build triplet hamiltonian
for i in range(0,ndocc):
    for a in range(0,nvirt):
        ia = i * nvirt + a
        for j in range(0,ndocc):
            for b in range(0,nvirt):
                jb = j * nvirt + b
                Ham[ia][jb] = - g_oovv[i][j][a][b]
        Ham[ia][ia] += eps_v[a] - eps_o[i]

# diagonalize Hamiltonian
eig = np.linalg.eigvals(Ham)

# sort excitation energies
eig.sort()

print("")
print("    ==> CIS triplet excitation energies (eV) <==")
print("")
for ia in range(0,ndocc*nvirt):
    print("    %5i %10.5f" % (ia,eig[ia]*27.21138))
print("")

    ==> CIS triplet excitation energies (eV) <==

        0   10.00098
        1   12.10703
        2   12.55493
        3   13.80473
        4   16.74649
        5   18.67018
        6   33.30584
        7   36.30692
        8  544.64456
        9  546.40432