Configuration Interaction with Single Excitations¶

In this project, you will learn how to write a restricted configuration interaction with single excitations (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. For example, you will soon find that it is possible to write a working CIS code in Python in less than 100 lines!

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 molecular orbitals, 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. 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

$$|\psi(n)\rangle = \sum_{ia} c_i^a(n) |\phi_i^a\rangle$$

Here, the symbol φ_i^a represents a Slater determinant that is singly excited relative to the reference determinant, and the indices i and a represent a spin orbitals that are occupied and unoccupied (virtual) in the reference, respectively.

The expansion coefficients and energy of the excited state can be determined by testing the Schrödinger equation with one of the singly-excited configurations to obtain

$$\sum_{jb} \langle \phi_i^a | \hat{H} | \phi_j^b \rangle c_j^b(n) = E(n) c_i^a(n)$$

We can form similar expressions for all excited states, and, taken together, these expressions become an eigenvalue equation

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

where H is the matrix representation of the Hamiltonian in the basis of singly-excited determinants. The elements of this matrix are given by

$$H_{ia,jb} = F_{ab}\delta_{ij}-F_{ij}\delta_{ab} + \langle aj ||ib\rangle$$

where $F$ represents the Fock matrix, and $\langle aj || ib \rangle$ represents an antisymmetrized electron-repulsion integral in physicists' notation. Recall that the Fock matrix for restricted Hartree-Fock theory is diagonal in the MO basis, so this expression simplifies sligthly:

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

where $\epsilon_i$ and $\epsilon_a$ represent occupied and virtual orbital energies, respectively.

Singlet Excitation Energies¶

It turns out that when the eigenvalues of the CIS Hamiltonian represented in this 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

$$ |\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) \delta_{ij}\delta_{ab} + 2 \langle aj |ib\rangle - \langle aj | bi \rangle$$

Using chemists' notation, these matrix elements are

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

Note that the orbital labels $i$, $j$, $a,$ and $b$ no longer represent spin orbitals; they represent spatial orbitals. 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 [6]:

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 [7]:

# 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 [8]:

# Grab data from wavfunction

# number of doubly occupied orbitals
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 [9]:

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

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

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

Given these tensors, you can access the element (ij|ab) in Python as 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 [10]:

# 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 [11]:

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

The composite index that represents each basis function (φia), 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 [16]:

# 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 * ovov[i][a][j][b] - 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

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) \delta_{ij}\delta_{ab} - ( ab | ji )$$

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

In [19]:

# 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] = - 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

In [ ]: