Reduced particle densities#

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 be evaluated as expectation values of one- and two-particle operators, respectively, using the expressions 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.

As an illustration, we will determine the one- and two-particle densities along the internuclear axis of carbon monoxide at the level of Hartree–Fock theory.

import matplotlib.pyplot as plt
import numpy as np
import veloxchem as vlx

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

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

One-particle density#

Probabilistic interpretation#

The probability density of finding any one electron in an infinitesimal volume element at position \(\mathbf{r}\) regardless of the positions of other electrons is equal to

\[ n(\mathbf{r}) = N \int \cdots \int \Psi^\dagger(\mathbf{r}, \mathbf{r}_2, \ldots, \mathbf{r}_N) \Psi(\mathbf{r}, \mathbf{r}_2, \ldots, \mathbf{r}_N) \, \mathrm{d}^3 \mathbf{r}_2 \cdots \mathrm{d}^3\mathbf{r}_N \]

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.

Illustration#

The one-particle density is available in VeloxChem by means of the VisualizationDriver class. As indicated by the string argument “alpha” in the get_density method, the \(\alpha\)- and \(\beta\)-spin densities can be determined individually, referring to the relation

\[ n(\mathbf{r}) = n^\alpha(\mathbf{r}) + n^\beta(\mathbf{r}) \]

For a closed shell system, the two spin densities are equal and thus equal to half the total one-particle density.

vis_drv = vlx.VisualizationDriver()

# list of coordinates in units of Bohr
n = 2000
coords = np.zeros((n, 3))
z = np.linspace(-1, 4, n)
coords[:, 2] = z

one_part_den = vis_drv.get_density(coords, molecule, basis, scf_drv.density, "alpha")
fig = plt.figure(figsize=(12, 4))

plt.subplot(1, 2, 1)
plt.plot(z, one_part_den, lw=2, label=r"$n^\alpha(\mathbf{r})$")

plt.setp(plt.gca(), xlim=(-1, 4))
plt.legend()

plt.xlabel(r"Internuclear axis coordinate (Bohr)")
plt.ylabel(r"One-particle density (a.u.)")

plt.subplot(1, 2, 2)
plt.plot(z, one_part_den, lw=2, label=r"$n^\alpha(\mathbf{r})$")
plt.plot(z, one_part_den, "ko")
plt.grid(True)

plt.setp(plt.gca(), xlim=(-0.05, 0.05), ylim=(50, 62))
plt.legend()

plt.xlabel(r"Internuclear axis coordinate (Bohr)")
plt.ylabel(r"One-particle density (a.u.)")

plt.show()
../../_images/reduced_density_4_0.png

The right figure panel is an enlargement of the one-particle density in the region of the carbon nucleus, illustrating the fact that \(n(\mathbf{r})\) displays a cusp at any given nuclear position.

Two-particle density#

Probabilistic interpretation#

The probability density of finding any two electrons in separate infinitesimal volume elements at positions \(\mathbf{r}_1\) and \(\mathbf{r}_2\) regardless of the positions of other electrons is equal to

\[ n(\mathbf{r}_1, \mathbf{r}_2) = N (N-1) \int \cdots \int \Psi^\dagger(\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3, \ldots, \mathbf{r}_N) \Psi(\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3, \ldots, \mathbf{r}_N) \, \mathrm{d}^3 \mathbf{r}_3 \cdots \mathrm{d}^3\mathbf{r}_N \]

Two-particle density operator#

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(\mathbf{r}) D_{pqrs} \psi_r^\dagger(\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.

Illustration#

As indicated by the string arguments “alpha” and “beta” in the get_two_particle_density method, the spin densities can be determined individually, referring to the relation

\[ n(\mathbf{r}_1, \mathbf{r}_2) = n^{\alpha\alpha}(\mathbf{r}_1, \mathbf{r}_2) + n^{\alpha\beta}(\mathbf{r}_1, \mathbf{r}_2) + n^{\beta\alpha}(\mathbf{r}_1, \mathbf{r}_2) + n^{\beta\beta}(\mathbf{r}_1, \mathbf{r}_2) \]

For a closed shell system, we have

\[\begin{align*} n^{\alpha\alpha}(\mathbf{r}_1, \mathbf{r}_2) & = n^{\beta\beta}(\mathbf{r}_1, \mathbf{r}_2) \\ n^{\alpha\beta}(\mathbf{r}_1, \mathbf{r}_2) & = n^{\beta\alpha}(\mathbf{r}_1, \mathbf{r}_2) \end{align*}\]

We will determine the two-particle density for electron 1 being poistioned at the carbon nucleus and electron 2 being positioned along the internuclear axis.

n = 100

origin = np.zeros((n, 3))

coords = np.zeros((n, 3))
z = np.linspace(-1, 4, n)
coords[:, 2] = z

two_part_den_aa = vis_drv.get_two_particle_density(
    origin, coords, molecule, basis, scf_drv.density, "alpha", "alpha"
)
two_part_den_ab = vis_drv.get_two_particle_density(
    origin, coords, molecule, basis, scf_drv.density, "alpha", "beta"
)
fig = plt.figure(figsize=(12, 8))

plt.subplot(2, 2, 1)
plt.plot(
    z, two_part_den_aa, lw=2, label=r"$n^{\alpha\alpha}(\mathbf{r}_1, \mathbf{r}_2)$"
)

plt.setp(plt.gca(), xlim=(-1, 4))
plt.legend()

plt.figtext(0.2, 0.75, r"Same spin electrons")

plt.xlabel(r"Internuclear axis coordinate (Bohr)")
plt.ylabel(r"Two-particle density (a.u.)")

plt.subplot(2, 2, 3)
plt.plot(
    z, two_part_den_aa, lw=2, label=r"$n^{\alpha\alpha}(\mathbf{r}_1, \mathbf{r}_2)$"
)
plt.plot(
    z, two_part_den_ab, lw=2, label=r"$n^{\alpha\beta}(\mathbf{r}_1, \mathbf{r}_2)$"
)

plt.setp(plt.gca(), xlim=(-1, 2), ylim=(0, 20))
plt.legend()

plt.figtext(0.22, 0.25, r"Same spin")
plt.figtext(0.15, 0.15, r"Fermi hole")
plt.figtext(0.32, 0.4, r"Opposite spin ")

plt.xlabel(r"Internuclear axis coordinate (Bohr)")
plt.ylabel(r"Two-particle density (a.u.)")

plt.subplot(2, 2, 2)
plt.plot(
    z, two_part_den_ab, lw=2, label=r"$n^{\alpha\beta}(\mathbf{r}_1, \mathbf{r}_2)$"
)

plt.setp(plt.gca(), xlim=(-1, 4))
plt.legend()

plt.figtext(0.63, 0.75, r"Oppsosite spin electrons")

plt.xlabel(r"Internuclear axis coordinate (Bohr)")
plt.ylabel(r"Two-particle density (a.u.)")

plt.show()
../../_images/reduced_density_8_0.png