First-order properties

First-order properties#

In Kohn–Sham DFT the accessible quantity is the one-particle density that here is expressed in terms of the atomic orbitals and their associated density matrix

\[ n(\mathbf{r}) = \sum_{\alpha, \beta} D_{\alpha\beta} \chi_\alpha^\ast(\mathbf{r}) \chi_\beta(\mathbf{r}) \]

Let us determine this density matrix for water at the level of B3LYP/cc-pVDZ.

import numpy as np
import veloxchem as vlx

Define the system.

h2o_xyz = """3

O    0.000000000000        0.000000000000        0.000000000000
H    0.000000000000        0.740848095288        0.582094932012
H    0.000000000000       -0.740848095288        0.582094932012
"""
molecule = vlx.Molecule.read_xyz_string(h2o_xyz)
basis = vlx.MolecularBasis.read(molecule, "cc-pvdz", ostream=None)

molecule.show()

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

Define the method and perform an SCF optimization.

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

D = scf_results["D_alpha"] + scf_results["D_beta"]

An observable associated with a scalar one-electron operator in the coordinate basis

\[ \hat{\Omega} = \sum_{i=1}^N \omega(\mathbf{r}_i) \]

takes the form of an expectation value

\[\begin{align*} \langle \hat{\Omega} \rangle & = \int \omega(\mathbf{r}) \, n(\mathbf{r}) \, \mathrm{d}^3 \mathbf{r} \end{align*}\]

We will determine the permanent dipole moment along the \(z\)-axis such that

\[ \hat{\Omega} = \hat{\mu}_z ; \qquad \omega(\mathbf{r}) = - e z \]

Numerical integration#

In DFT, integrals involving the exchange–correlation kernel are evaluated numerically. The same technique can, however, be used to evaluate first-order properties

\[ \langle \hat{\Omega} \rangle \approx \sum_{g=1}^P w_g \, \omega(\mathbf{r}_g) \, n(\mathbf{r}_g) \]
grid_level = scf_drv.grid_level
if grid_level is None:
    grid_level = vlx.dftutils.get_default_grid_level(scf_drv.xcfun)

grid_drv = vlx.veloxchemlib.GridDriver()
grid_drv.set_level(grid_level)
molgrid = grid_drv.generate(molecule)

# grid point weights
w_g = molgrid.w_to_numpy()

# grid point z-coordinates
z_g = molgrid.z_to_numpy()

xc_drv = vlx.XCIntegrator()

# generate AOs on the grid points
chi_g = xc_drv.compute_gto_values(molecule, basis, molgrid)

# determine the density on the grid points
G = np.einsum("ab,bg->ag", D, chi_g)
n_g = np.einsum("ag,ag->g", chi_g, G)

print(f"Dipole moment: {-1.0 * np.dot(w_g, z_g * n_g) : 12.8f}")

Dipole moment:  -1.44687595

The numerical accuracy in this integration is limited and determined by the quality of the grid. We here adopted the default grid level in VeloxChem.

Analytic integration#

If we combine the expressions for the density and the expectation value, we get

\[\begin{align*} \langle \hat{\Omega} \rangle &= \sum_{\alpha,\beta} D_{\alpha\beta} \langle \chi_\alpha | \hat{\omega} | \chi_\beta \rangle \end{align*}\]

In form, this Kohn–Sham DFT expression is identical to the result obtained when considering the expectation value of the operator \(\hat{\Omega}\) with respect to the many-electron Hartree–Fock wave function. It would be incorrect, however, to think of the Kohn–Sham reference state as the wave function of the system.

In practice, this expression is used for the evaluation of first-order properties at the DFT level. The scheme is referred to as analytic since the property integrals are determined to double precision on the computer and the only numerical integrations are those performed in the SCF iterations.

dipole_drv = vlx.ElectricDipoleIntegralsDriver()
dipole_mats = dipole_drv.compute(molecule, basis)

mu_z = -1.0 * dipole_mats.z_to_numpy()

print(f"Dipole moment: {np.einsum('ab,ab->', D, mu_z) : 12.8f}")

Dipole moment:  -1.44687507

The dipole moment obtained from numerical and analytic integrations are seen to be in very close agreement. The error in the numerical integration is in the order of \(10^{-6}\) a.u.

Contents