# Operators and matrix elements#

## Single-electron systems#

Operators acting on single-electron wave functions, or spin orbitals (ignoring the freedom of charge conjugation), take the form

$\hat{\omega} = \hat{\omega}^\mathrm{o} \otimes \hat{\omega}^\mathrm{s}$

Pure spatial (or spin) operators refer to the situation when $$\hat{\omega}^\mathrm{s}$$ (or $$\hat{\omega}^\mathrm{o}$$) equals the identity operator $$\hat{I}^\mathrm{s}$$ (or $$\hat{I}^\mathrm{o}$$).

Note

All electronic operators have components in all underlying vector spaces. However, it is a common practice to leave out a reference to identiy operators.

Example: The kinetic energy operator can be written

$\begin{split} \hat{T} = \frac{\hat{p}^2}{2 m_\mathrm{e}} \otimes \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} \equiv \frac{\hat{p}^2}{2 m_\mathrm{e}} \end{split}$

## Many-electron systems#

### One-electron operators#

One-electron operators take the form

$\hat{\Omega} = \hat{\omega} \otimes \hat{I} \otimes \cdots \otimes \hat{I} + \cdots + \hat{I} \otimes \hat{I} \otimes \cdots \otimes \hat{\omega}$

Despite the name, one-electron operators act on many-electron wave functions, but each term in the sum has a non-trivial action only in a single particle space. Furthermore, the non-trivial one-particle operator $$\hat{\omega}$$ is identical in all terms in the sum.

Note

It is a common practice to leave out a reference to identiy operators.

Example: The kinetic energy operator for a $$N$$-electron system can be written

$\hat{T} = \sum_{i=1}^N \frac{\hat{p}^2(i)}{2 m_\mathrm{e}}$

where the particle index, $$i$$, is used to denote in which particle space the action takes place.

#### Action of one-electron operators on Slater determinants#

Consider a one-electron operator for an $$N$$-electron system

$\hat{\Omega} = \sum_{i=1}^N \hat{\omega}(i)$

The action of this operator on a Slater determinant results in a sum of $$N$$ Slater determinants, each with a single modified spin-orbital

$\hat{\Omega} | \psi_1, \ldots, \psi_N \rangle = | \tilde{\psi}_1, \ldots, \psi_N \rangle + \cdots + | \psi_1, \ldots, \tilde{\psi}_N \rangle$

where the modified orbitals are given by

$| \tilde{\psi}_p \rangle = \hat{\omega} | \psi_p \rangle$

#### Commutator of one-electron operators#

Consider two one-electron operators for an $$N$$-electron system

$\hat{\Omega} = \sum_{i=1}^N \hat{\omega}(i) ; \qquad \hat{\Lambda} = \sum_{i=1}^N \hat{\lambda}(i)$

The commutator of these operators results in another one-electron operator

$\big[ \hat{\Omega}, \hat{\Lambda} \big] = \sum_{i=1}^N [\hat{\omega}, \hat{\lambda}](i)$

### Two-electron operators#

Two-electron operators take the form

$\hat{\Omega} = \sum_{j > i}^N \hat{\omega}(i,j)$

Analogous to one-electron operators, two-electron operators act on many-electron wave functions, but each term in the sum has a non-trivial action only in two particle spaces. Furthermore, the non-trivial two-particle operator $$\hat{\omega}(i,j)$$ is identical in all terms in the sum, where the particle indices, $$i$$ and $$j$$, are used to denote in which two particle spaces the action takes place.

Example: The Coulomb repulsion operator for a $$N$$-electron system can be written

$\hat{V}^\mathrm{e-e} = \sum_{j > i}^N \frac{e^2}{4 \pi \varepsilon_0 |\mathbf{r}_i - \mathbf{r}_j|}$

where $$e$$ is the elementary charge.

## Matrix elements for Slater determinants#

Let $$|\Phi \rangle$$ be a Slater determinant for an $$N$$-electron system.

With $$\hat{\Omega}$$ being a one-electron operator, the pertinent matrix elements become

\begin{align*} \langle \Phi | \hat{\Omega} | \Phi \rangle &= \sum_{i}^\mathrm{occ} \langle i| \hat{\omega} | i \rangle \\ \langle \Phi | \hat{\Omega} | \Phi_{i}^{a} \rangle &= \langle i| \hat{\omega} |a \rangle \end{align*}

With $$\hat{\Omega}$$ being a two-electron operator, the pertinent matrix elements become

\begin{align*} \langle \Phi | \hat{\Omega} | \Phi \rangle &= \frac{1}{2} \sum_{i,j}^\mathrm{occ} \left[\rule{0pt}{12pt} \langle ij| \hat{\omega} | ij \rangle - \langle ij| \hat{\omega} |ji \rangle \right] \\ \langle \Phi | \hat{\Omega} | \Phi_{i}^{a} \rangle &= \sum_{j}^\mathrm{occ} \left[\rule{0pt}{12pt} \langle ij| \hat{\omega} |aj \rangle - \langle ij| \hat{\omega} |ja \rangle \right] \\ \langle \Phi | \hat{\Omega} | \Phi_{ij}^{ab} \rangle &= \langle ij| \hat{\omega} |ab \rangle - \langle ij| \hat{\omega} |ba \rangle \end{align*}

The two-electron integrals are here written in physicist’s notation and summation run over occupied spin orbitals.

### Spin orbitals#

Integrals with respect to spin orbitals are in this book expressed in the physicist’s notation:

$\omega_{ab} = \langle \psi_a | \hat{\omega} | \psi_b \rangle = \int \psi_a^\dagger(\mathbf{r}) \, \hat{\omega} \, \psi_b(\mathbf{r}) \, d^3\mathbf{r}$

and

$\omega_{abcd} = \langle \psi_a \psi_b | \hat{\omega} | \psi_c \psi_d \rangle = \int \psi_a^\dagger(\mathbf{r}_1) \psi_b^\dagger(\mathbf{r}_2) \, \hat{\omega} \, \psi_c(\mathbf{r}_1) \psi_d(\mathbf{r}_2) \, d^3\mathbf{r}_1 d^3\mathbf{r}_2$

### Molecular orbitals#

Integrals with respect to spatial (or molecular) orbitals are in this book expressed in the chemist’s notation:

$\omega_{ab} = ( \phi_a | \hat{\omega} | \phi_b ) = \int \phi_a^\ast(\mathbf{r}) \, \hat{\omega} \, \phi_b(\mathbf{r}) \, d^3\mathbf{r}$

and

$\omega_{abcd} = ( \phi_a \phi_b | \hat{\omega} | \phi_c \phi_d ) = \int \phi_a^\ast(\mathbf{r}_1) \phi_b(\mathbf{r}_1) \, \hat{\omega} \, \phi_c^\ast(\mathbf{r}_2) \phi_d(\mathbf{r}_2) \, d^3\mathbf{r}_1 d^3\mathbf{r}_2$

## Reduced particle density operators#

Consider an $$N$$-electron system in a state represented by a general state vector that can be expressed as a linear combination of Slater determinants

$| \Psi \rangle = \sum_n c_n | \Phi_n \rangle$

The one- and two-particle particle reduced densities can then be evaluated as expectation values of one- and two-particle operators, respectively, using the expressions given above for matrix elements.

The resulting formulas will be expressions involving products of molecular orbitals from which we can identify one- and two-particle density matrices.

### One-particle density operator#

The one-particle density can be written as an expectation of a one-electron operator according to

$n(\mathbf{r}) = \langle \Psi | \hat{n}(\mathbf{r}) | \Psi \rangle$

with

$\hat{n}(\mathbf{r}) = \sum_{i=1}^N \delta(\mathbf{r} - \mathbf{r}_i)$

After evaluation of this expectation value, the one-particle density matrix, $$D$$, is identified from the expression

$n(\mathbf{r}) \equiv \sum_{p,q} \psi_p^\dagger(\mathbf{r}) D_{pq} \psi_q(\mathbf{r})$

For a state represented by a single Slater determinant, the one-electron density becomes

$n(\mathbf{r}) = \sum_{i=1}^N |\psi_i(\mathbf{r})|^2$

where the summation run over the $$N$$ occupied spin orbitals. Consequently, the density matrix is in this case equal to the identity matrix in the occupied–occupied block and zero elsewhere.

### Two-particle density#

The two-particle density can be written as an expectation of a two-electron operator according to

$n(\mathbf{r}, \mathbf{r}') = \langle \Psi | \hat{n}(\mathbf{r}, \mathbf{r}') | \Psi \rangle$

with

$\hat{n}(\mathbf{r}, \mathbf{r}') = \sum_{i=1}^N \sum_{j>i}^N \left[ \delta(\mathbf{r} - \mathbf{r}_i) \delta(\mathbf{r}' - \mathbf{r}_j) + \delta(\mathbf{r} - \mathbf{r}_j) \delta(\mathbf{r}' - \mathbf{r}_i) \right]$

After evaluation of this expectation value, the two-particle density matrix, $$D$$, is identified from the expression

$n(\mathbf{r}, \mathbf{r}') \equiv \sum_{p, q, r, s} \psi_p^\dagger(\mathbf{r}) \psi_q^\dagger(\mathbf{r}') D_{pqrs} \psi_r(\mathbf{r}) \psi_s(\mathbf{r}')$

For a state represented by a single Slater determinant, the two-electron density becomes

$n(\mathbf{r}, \mathbf{r}') = \sum_{i=1}^N \sum_{j=1}^N \left[\rule{0pt}{12pt} |\psi_i(\mathbf{r})|^2 |\psi_j(\mathbf{r}')|^2 - \psi_i^\dagger(\mathbf{r}) \psi_j^\dagger(\mathbf{r}') \psi_j(\mathbf{r}) \psi_i(\mathbf{r}') \right]$

where the summations run over the $$N$$ occupied spin orbitals.

## Molecular Hamiltonian#

Molecular quatum mechanics is concerned with the motions of electrons and nuclei, governed in turn by the molecular Hamiltonian

$\hat{H} = \hat{T}_\mathrm{n} + \hat{H}_\mathrm{e}$

where the nuclear kinetic energy operator is given by

$\hat{T}_\mathrm{n} = - \sum_{A=1}^M \frac{\hbar^2}{2 M_A} \nabla^2_A$

The electromagnetic forces acting on nuclei and electrons are of the same order of magnitude, but since nuclear masses are at least 1800 times larger ($$m_\mathrm{p} / m_\mathrm{e} \approx 1836.15$$), we expect nuclear velocities to be much smaller than those of electrons. To decouple their motions, it is a common practice to work in the basis of the eigenstates of the electronic Hamiltonian

$\hat{H}_\mathrm{e} = % - \sum_{i=1}^N \frac{\hbar^2}{2 m_\mathrm{e}} \nabla^2_i % - \sum_{i=1}^N \sum_{A=1}^M \frac{Z_A e^2}{4 \pi \varepsilon_0 |\mathbf{r}_i - \mathbf{R}_A|} % + \sum_{i=1}^N \sum_{j>i}^N \frac{e^2}{4 \pi \varepsilon_0 |\mathbf{r}_i - \mathbf{r}_j|} % + \sum_{A=1}^M \sum_{B>A}^M \frac{Z_A Z_B e^2}{4 \pi \varepsilon_0 |\mathbf{R}_A - \mathbf{R}_B|}$

introducing in order terms associated with the electronic kinetic energy, electron–nuclear attraction, electron–electron repulsion, and nuclear–nuclear repulsion. We use variables $$\mathbf{r}$$ and $$\mathbf{R}$$ to collectively denote the sets of electronic and nuclear coordinates, respectively.

In a brief notation, the electronic Hamiltonian is expressed in terms of the one- and two-electron components in addition to the nuclear repulsion term

$\hat{H}_\mathrm{e} = % \sum_{i} \hat{h}(i) + \sum_{j>i} \hat{g}(i,j) + V^\mathrm{n,rep}$

For a given fixed nuclear configuration the nuclear repulsion amounts a mere scalar number as we solve the time-independent Schrödinger equation

$\hat{H}_\mathrm{e} \Psi_K^\mathrm{e}(\mathbf{r}; \mathbf{R}) = V_K(\mathbf{R}) \Psi_K^\mathrm{e}(\mathbf{r}; \mathbf{R})$

where it is noted that the electronic wave functions (indexed with $$K$$) have a parametric dependence on the nuclear coordinates. We use these states as a basis in which we express a given total rovibronic molecular wave function according to

$\Psi(\mathbf{r}, \mathbf{R}) = \sum_K \Phi_K(\mathbf{R}) \Psi_K^\mathrm{e}(\mathbf{r}; \mathbf{R})$

## Adiabatic Born–Oppenheimer approximation#

In view of the differences in nuclear and electronic velocities, it is reasonable to assume that the variations of the electronic wave functions with respect to nuclear displacements are small. Moreover, the summation over electronic states is often approximated with a single term. This is known as the adiabatic approximation and it is accurate as long as the electronic states are well separated in energy.

We arrive at the adiabatic Born–Oppenheimer approximation where the total molecular wave function takes the form of a product

$\Psi_{K, k}(\mathbf{r}, \mathbf{R}) = \Phi_{K;k}(\mathbf{R}) \Psi_K^\mathrm{e}(\mathbf{r}; \mathbf{R})$

or, in a state vector notation,

$|K, k \rangle = | K \rangle| k \rangle_K$

In cases where it is clear from the context, the electronic state index on the nuclear state vector $$|k\rangle$$ is often omitted.

## Matrix element of electronic Hamiltonian#

import numpy as np
import veloxchem as vlx

mol_str = """
C        0.00000000    0.00000000    0.00000000
O        0.00000000    0.00000000    1.43
"""
molecule = vlx.Molecule.read_str(mol_str, units="angstrom")
basis = vlx.MolecularBasis.read(molecule, "cc-pVDZ")

scf_drv = vlx.ScfRestrictedDriver()
scf_drv.compute(molecule, basis)


Some information after the SCF optimization.

nocc = molecule.number_of_alpha_electrons()
C = scf_drv.scf_tensors["C"]  # MO coefficient matrix

E_nuc = molecule.nuclear_repulsion_energy()
E_tot = scf_drv.get_scf_energy()
E_elec = E_tot - E_nuc

print("Number of occupied orbitals:", nocc)

print(f"Nuclear repulsion energy: {E_nuc : 19.12f}")
print(f"Electronic energy: {E_elec : 26.12f}")
print(f"Total Hartree-Fock energy: {E_tot : 18.12f}")

Number of occupied orbitals: 7
Nuclear repulsion energy:     17.762591694646
Electronic energy:          -130.366128020400
Total Hartree-Fock energy:  -112.603536325754


### Expectation value of one-electron Hamiltonian#

For a closed-shell Slater determinant, the expectation value of a one-electron operator becomes

$\langle \Phi | \hat{\Omega} | \Phi \rangle = 2 \sum_{i}^\mathrm{occ} ( i| \hat{\omega} | i )$

where the sum now runs over doubly occpuied molecular orbitals, referring to spatial integrals.

kin_drv = vlx.KineticEnergyIntegralsDriver()
T_ao = kin_drv.compute(molecule, basis).to_numpy()
T = np.einsum("ai, ab, bj -> ij", C, T_ao, C)

npot_drv = vlx.NuclearPotentialIntegralsDriver()
V_ao = -npot_drv.compute(molecule, basis).to_numpy()
V = np.einsum("ai, ab, bj -> ij", C, V_ao, C)

# core Hamiltonian
h = T + V

h_aver = 2 * np.einsum("ii", h[:nocc, :nocc])
print("Expectation value of one-electron Hamiltonian:", h_aver)

Expectation value of one-electron Hamiltonian: -188.97928532728577


### Expectation value of two-electron Hamiltonian#

For a closed-shell Slater determinant, the expectation value of a two-electron operator becomes

$\langle \Phi | \hat{\Omega} | \Phi \rangle = \sum_{i,j}^\mathrm{occ} \left[\rule{0pt}{12pt} 2 ( ii| \hat{\omega} | jj ) - ( ij| \hat{\omega} |ji ) \right]$

where the sum now runs over doubly occpuied molecular orbitals, referring to spatial two-electron integrals (Coulomb and exchange) written in the chemist’s notation.

# obtain all-occupied block ("OOOO") of ERIs in physicist's notation
moints_drv = vlx.MOIntegralsDriver()
ijij = moints_drv.compute_in_mem(molecule, basis, scf_drv.mol_orbs, "OOOO")

g_aver = 2 * np.einsum("ijij", ijij[:, :, :, :]) - np.einsum("ijji", ijij[:, :, :, :])
print("Expectation value of two-electron Hamiltonian:", g_aver)

Expectation value of two-electron Hamiltonian: 58.61315730688542


### Expectation value of electronic Hamiltonian#

print(f"Expectation value of electronic Hamiltonian: {h_aver + g_aver:.12f}")

Expectation value of electronic Hamiltonian: -130.366128020400