Reduced particle densities#

Reduced particle densities associated with the probalistic interpretation of wave functions are made available in VeloxChem by means of the VisualizationDriver class. The reduced particle densities can be determined as expectation values of the one- and two-particle reduced density operators.

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#

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, 0, "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/dens_vis_6_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#

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, 0, "alpha", "alpha"
)
two_part_den_ab = vis_drv.get_two_particle_density(
    origin, coords, molecule, basis, scf_drv.density, 0, "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/dens_vis_11_0.png