Configuration Interaction Singles
The configuration interaction singles (CIS) method applies a truncated CI wavefunction up to single excitations to a reference determinant, usually the Hartree-Fock determinant: $$ |\Psi^n_{\mathrm{CIS}}\rangle = c_0^n |\Phi_0\rangle + \sum_{ia} c_{ia}^n |\Phi_i^a\rangle $$ The CIS matrix thus takes the form
+----------------+---------------------------+
| "<𝚽₀|Ĥ|𝚽₀>" | "<𝚽₀|Ĥ|𝚽ᵢᵃ>" |
+----------------+---------------------------+
| | |
| | |
| | |
| | |
| | |
| "<𝚽ᵢᵃ|Ĥ|𝚽₀>" | "<𝚽ᵢᵃ|Ĥ|𝚽ⱼᵇ>" |
| | |
| | |
| | |
| | |
| | |
| | |
+----------------+---------------------------+
We shall assume for simplicity that the reference determinant is the Hartree-Fock determinant. This makes our calculation easier, because the Brillouin condition is satisfied and the reference determinant does not interact with singly excited determinants, i.e. \(\langle \Phi_0 | \hat{H} | \Phi_i^a \rangle = 0\). This makes our CIS matrix block diagonal. Since the first block is just the Hartree-Fock energy, we only have to calculate the second block: \(\langle \Phi_i^a | \hat{H} | \Phi_j^b \rangle\) explicitly.
After some algebraical magic, we arrive at $$ \langle \Phi_i^a | \hat{H} | \Phi_j^b \rangle = (E_0 + \epsilon_a - \epsilon_i) \delta_{ij} \delta_{ab} + (jb|ai) - (ji|ab) $$
It should be noted that the two-electron integrals are required in the MO basis, so we have to transform them from the AO basis: $$ (pq|rs) = \sum_{\mu \nu \sigma \lambda} C_{\mu p}^* C_{\nu q} C_{\sigma r}^* C_{\lambda s} (\mu \nu | \sigma \lambda) $$ where \(C_{\mu p}\) is the \(\mu\)-th AO coefficient of the \(p\)-th MO.
We can now implement the CIS method. Since we need several quantities from
the Hartree-Fock calculation, it is convenient to let the CIS class inherit
from the HartreeFock class.
import numpy as np
from hartree_fock import HartreeFock
HARTREE_TO_EV = 27.211_386_245_988
class CIS(HartreeFock):
def initialize(self):
super().initialize()
self.run_hf(max_iter=500, threshold=1e-8, verbose=0)
print(self.energy)
def get_cis_hamiltonian(self):
# Transform ERIs from AO basis to MO basis
eri_mo = np.einsum('pQRS, pP -> PQRS',
np.einsum('pqRS, qQ -> pQRS',
np.einsum('pqrS, rR -> pqRS',
np.einsum('pqrs, sS -> pqrS', self.twoel,
self.mo_coeff, optimize=True),
self.mo_coeff, optimize=True),
self.mo_coeff, optimize=True),
self.mo_coeff, optimize=True)
# Transform to spin-orbit basis
norb = len(self.mo_energy) * 2
nocc = self.nocc * 2
nvirt = norb - nocc
eps = np.repeat(self.mo_energy, 2)
delta = np.zeros((2, 2, 2, 2))
delta[(0, 0, 0, 0)] = 1
delta[(0, 0, 1, 1)] = 1
delta[(1, 1, 0, 0)] = 1
delta[(1, 1, 1, 1)] = 1
eri_mo = np.kron(eri_mo, delta)
# Obtain orbital labels
self.orb_labels = []
for i in range(0, len(self.mo_energy)):
self.orb_labels.extend([f'{i}a', f'{i}b'])
# Obtain excitations in spin-orbit basis
self.excitations = []
for i in range(0, nocc):
for a in range(nocc, norb):
self.excitations.append((i, a))
# Build the Hamiltonian
hamiltonian = np.zeros((nocc * nvirt, nocc * nvirt))
for p, left_excitation in enumerate(self.excitations):
i, a = left_excitation
for q, right_excitation in enumerate(self.excitations):
j, b = right_excitation
hamiltonian[p, q] = (eps[a] - eps[i]) * (i == j) * (a == b) \
+ eri_mo[j, b, a, i] - eri_mo[j, i, a, b]
return hamiltonian
def run_cis(self, nprint=None):
h = self.get_cis_hamiltonian()
eigval, eigvect = np.linalg.eigh(h)
self.cis_energies = eigval
self.cis_states = eigvect
e_ev = eigval * HARTREE_TO_EV
# Print detailed information on significant excitations
print('CIS:')
if nprint is None:
nstate = len(self.excitations)
elif nprint < 0:
nstate = 0
else:
nstate = min(nprint, len(self.excitations))
for state in range(0, nstate):
print(f'Excited State {state + 1:3d}: '
f'E_exc = {e_ev[state]:10.4f} eV')
for idx, exc in enumerate(self.excitations):
coeff = eigvect[idx, state]
contribution = np.abs(coeff)**2
if contribution > 0.1:
i, a = exc
il, al = (self.orb_labels[x] for x in exc)
print(f'{il:4s} -> {al:4s} '
f'{coeff:12.6f} ({100.0 * contribution:.1f} %)')
print()
return eigval
The method get_cis_hamiltonian is the centerpiece of the CIS class. It
ar first transforms the ERIs to the MO basis, spin-blocks them, and the
enumerates all singly excited determinants. The CIS Hamiltonian is then
constructed from the transformed ERIs and the MO energies.
The spin-blocking is useful, since we can obtain excited states with all possible spin multiplicities. Although it is not necessary for CIS, since only singlets and triplets are obtainable from a singlet reference, it will become very useful when higher excitations are included.
Note: Although spin-blocking will ease the implementation of higher excitations, it makes our algorithm less efficient. First of all, each dimension of the spin-blocked ERIs is twice as large as the original one, which makes it 16 times as big. Secondly, since states with different spin multiplicities are orthogonal in the absence of spin-orbit coupling, we only have to include some determinants (actually some carefully chosen linear combinations of determinants called configuration state functions) in the CI expansion. This will reduced the size of the CI matrix greatly when higher excitations are present.
The method run_cis takes the CIS hamiltonian and diagonalizes it. The
eigenvalues are the excitation energies, and the eigenvectors are the
coefficients of the excited determinants in the CIS wavefunction. This method
also prints details about the loewest excited states.
Now we can test our implementation on the water molecule.
from atom import Atom
# Coordinates are in the unit of Angstrom.
o1 = Atom('O', [ 0.000, 0.000, 0.000], unit='A')
h1 = Atom('H', [ 0.758, 0.587, 0.000], unit='A')
h2 = Atom('H', [-0.758, 0.587, 0.000], unit='A')
water = CIS()
water.set_atomlist([o1, h1, h2])
water.get_basis('sto-3g')
water.initialize()
water.run_cis(nprint=20)
We could again plot the density here for the excited states. Since the Slater determinants are orthogonal, the total (one-electron) density is just the sum of the densities of the individual determinants.
But with excited states, we can plot something more interesting: the transition density. Its definition is strongly motivated by the spin-free one-electron density: $$ \rho_{fi}(r) = \int \Psi_f^{*}(x, x_2, \cdots, x_N) \Psi_i(x, x_2, \cdots, x_N) \ \mathrm{d}\sigma \mathrm{d}x_2 \cdots \mathrm{d}x_N $$ So we have two different wavefunctions in the integral, one for the initial state \(\Psi_i\) and one for the final state \(\Psi_f\), instead of just one in the case of the density. This makes the transition density technically not a density, since it is not positive definite.
For a transition density from the ground state to a singly excited state, we can simplify the expression to $$ \rho_{fi}(r) = \sum_{i,a} c_{ia} \int \phi_a^{*}(x) \phi_i(x)\ \mathrm{d}\sigma $$
Using the functions we have written in section 5.1.1, we can calculate and visualize the transition density of water:
XLIM, YLIM, ZLIM = (-5.0, 5.0), (-5.0, 5.0), (-5.0, 5.0)
NX, NY, NZ = 80, 80, 80
grid = build_grid(XLIM, YLIM, ZLIM, NX, NY, NZ)
mo_grid = evaluate_mo_grid(mol, grid)
ISOSURFACE_COLORS = [(1, 0, 0), (0, 0, 1)]
DENSITY_ISOVALUES = [0.01, -0.01]
ISTATE = 3
td_grid = np.zeros((NX, NY, NZ))
for (exc, c_ia) in zip(mol.excitations, mol.cis_states[:, ISTATE]):
i, a = exc
if (i % 2) == (a % 2):
td_grid += c_ia * mo_grid[i] * mo_grid[a]
fig_td = mlab.figure()
visualize_cube(mol, grid, td_grid, DENSITY_ISOVALUES,
ISOSURFACE_COLORS, fig_td)
# Remove this line if you are using Jupyter notebook
mlab.show()
Since our \(\alpha\) orbitals have even indices and \(\beta\) orbitals have odd indices, the if statement ensures that we only include orbital pairs with the same spin.