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}\) (\(\hat{\omega}^\mathrm{o}\)) equals the identity operator \(\hat{I}^\mathrm{s}\) (\(\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 identity 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.

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 \(|\Psi \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 \Psi | \hat{\Omega} | \Psi \rangle &= \sum_{i}^\mathrm{occ} \langle i| \hat{\omega} | i \rangle \\ \langle \Psi | \hat{\Omega} | \Psi_{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 \Psi | \hat{\Omega} | \Psi \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 \Psi | \hat{\Omega} | \Psi_{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 \Psi | \hat{\Omega} | \Psi_{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 summations run over occupied spin orbitals.

Spin orbitals

Integrals with respect to spin orbitals are here 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 here 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 \]

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

\[\begin{align*} \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|}\\ % &\qquad + \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|} \end{align*}\]

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 = """2

C        -0.71500000    0.00000000   0.00000000
O         0.71500000    0.00000000   0.00000000
"""
molecule = vlx.Molecule.read_xyz_string(mol_str)
basis = vlx.MolecularBasis.read(molecule, "cc-pVDZ", ostream=None)

scf_drv = vlx.ScfRestrictedDriver()
scf_drv.ostream.mute()
scf_results = scf_drv.compute(molecule, basis)

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

From the converged SCF optimization we can retrieve information such as

nocc = molecule.number_of_alpha_electrons()
C = scf_results["C_alpha"]  # 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.97928532728582

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 the oooo-block of ERIs
eri_drv = vlx.MOIntegralsDriver()
g_oooo = eri_drv.compute_in_memory(molecule, basis, scf_drv.mol_orbs, "chem_oooo")

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

Expectation value of two-electron Hamiltonian: 58.61315730688549

Expectation value of electronic Hamiltonian

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

Expectation value of electronic Hamiltonian: -130.366128020400