Kernel integration

In DFT, the exchange-correlation (xc) kernel integrations are performed numerically

\[ E_\mathrm{xc}[n(\mathbf{r})] = \int n(\mathbf{r}) \, \varepsilon_\mathrm{xc}(n, \nabla n, \nabla^2 n) \, \mathrm{d}^3\mathbf{r} \approx \sum_{i=1}^P w_i \, n(\mathbf{r}_i) \, \varepsilon_\mathrm{xc}(n, \nabla n, \nabla^2 n) \,\Big|_{\mathbf{r} = \mathbf{r}_i} \]

In VeloxChem, the set of grid points and the associated grid-point weights, \(w_i\) are constructed using modified versions of the Log3-quadrature [TA95] and SSF-partitioning schemes [SSF96]. In practical work, the number of grid points, \(P\), amounts to in the order of \(10^6\).

Integration of the density itself results in the number of electrons in the system. Let us explore this simple case to assess the accuracy in the numerical integration scheme. We will use the electron density of the water molecule in our example calculations.

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

Define the molecular 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)

molecule.show()

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol

Define the method.

scf_drv = vlx.ScfRestrictedDriver()
scf_drv.ostream.mute()
scf_drv.xcfun = "b3lyp"

Molecular grid points

In the Python layer, VeloxChem provides access to the creation of molecular grids with the GridDriver class.

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 coordinates
x = molgrid.x_to_numpy()
y = molgrid.y_to_numpy()
z = molgrid.z_to_numpy()

P = molgrid.number_of_points()
print("Number of grid points:", P)

Number of grid points: 40720

Visualize all grid points outside the sphere of radius 5.0 \(a_0\).

grid_xyz = "\nGrid points\n"

r = np.sqrt(x**2 + y**2 + z**2)
au2ang = 0.52911721

i = 0
for p in range(P):
    if r[p] > 5.0:
        grid_xyz += f"He {x[p] * au2ang : 16.8f} {y[p] * au2ang : 16.8f} {z[p] * au2ang : 16.8f}\n"
        i += 1

grid_xyz = str(i) + grid_xyz

print("Number of grid point outside sphere:", i)
print("Number of grid point inside sphere:", P - i)

Number of grid point outside sphere: 5398
Number of grid point inside sphere: 35322

Not surprisingly, the grid point density is much larger inside the van der Waals volume of the molecule where the electron density is large.

Let us visualize all grid points outside the sphere of radius 5.0 \(a_0\).

import py3Dmol as p3d

v = p3d.view(width=400, height=400)

v.addModel(h2o_xyz, "xyz")
v.setStyle({"stick": {}})
v.addSurface(p3d.VDW, {"opacity": 0.7, "color": "green"})

v.addModel(grid_xyz, "xyz")
v.setStyle({"elem": "He"}, {"sphere": {"radius": 0.05, "color": "red", "opacity": 0.5}})

v.zoomTo()
v.show()

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol

Another important view of the grid points is from the perspective of their atomic radial density. It can be clearly seen that the grid density is much larger in the regions of the density maxima of the atomic shells.

count, rad_dist, patches = plt.hist(r, bins=75)

plt.grid(True)
plt.xlabel("Radial distance (a.u.)")
plt.ylabel("Grid point count")

plt.show()

Density on grid points

Obtain the total density matrix.

basis = vlx.MolecularBasis.read(molecule, "def2-svpd", ostream=None)
scf_results = scf_drv.compute(molecule, basis)

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

We will first study the total ground-state electron density and perform the integration with molecular grids of varying size and thereby quality. The quality of the grid is determined by the grid level and VeloxChem presently offers grid levels ranging from one to eight and for water that amounts to a number of grid points ranging from 4,400 to 1,107,900.

For a given density matrix, the density on grid points can be determined as follows

\[ n(\mathbf{r}_g) = \sum_{\alpha,\beta} \chi^*_\alpha(\mathbf{r}_g) D_{\alpha\beta} \chi_\beta(\mathbf{r}_g) \]

xc_drv = vlx.XCIntegrator()

n_grid_points = []
n_elec = []

for grid_level in range(1, 9):
    grid_drv.set_level(grid_level)  # available grid levels are between 1-8

    # generate grid points and weights for molecule
    molgrid = grid_drv.generate(molecule)
    weights = molgrid.w_to_numpy()
    n_grid_points.append(molgrid.number_of_points())

    # 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)
    n_elec.append(np.dot(weights, n_g))

print("Number of grid points:", n_grid_points)

Number of grid points: [4404, 10536, 20678, 40720, 63844, 147024, 348468, 1107906]

As measure of accuracy in the numerical integration, we will use the absolute value of the difference between the integer number of electron (ten for water) and the value of the integral with varing grid levels. This error measure is illustrated below in the plot.

fig = plt.figure(figsize=(8, 4))
ax1 = fig.add_subplot(111)

error = abs(np.array(n_elec) - 10)
ax1.plot(range(1, 9), error, "o")

plt.annotate(
    "Default grid level for GGA",
    (4.04, 1e-6),
    xytext=(3.3, 2e-5),
    arrowprops={"arrowstyle": "->"},
)
ax1.set_xlim((0.75, 8.25))

plt.grid()
plt.yscale("log")
plt.xlabel("Grid level")
plt.ylabel("Error in number of electrons")

ax2 = ax1.twiny()
ax2.set_xlim((0.75, 8.25))
ax2.set_xticks(range(1, 9))
ax2.set_xticklabels(n_grid_points)
ax2.set_xlabel("Number of grid points")
plt.show()

It is seen that the error in the numerical integration ranges from \(10^{-2}\) to \(10^{-9}\) and the error at the default grid level amounts to \(10^{-6}\) electrons. These results are rather insensitive to the choice of basis set as the location of radial density maxima for the occupied orbitals do not change much with the choice of basis set, and the numerical grids are optimized to cover well these spatial regoins.

It is less clear, however, how well the grid points cover the important spatial regions of the unoccupied (virtual) orbitals that become important in the calculation of excited states and therefore also in simulations of spectroscopies. To analyze this aspect further, we lock the grid level to four and instead perform numerical integrations of the densities of occupied and unoccupied orbitals. The reference integer number of electrons is in this case equal to two.

We will compare the results for two standard basis sets that are frequently used in DFT spectroscopy calculations namely def2-svpd and aug-cc-pvdz.

grid_drv.set_level(4)  # default level is 4 in VeloxChem

# generate grid points and weights for molecule
molgrid = grid_drv.generate(molecule)
weights = molgrid.w_to_numpy()

n_elec = {}
norb = {}

basis_sets = ["def2-svpd", "aug-cc-pvdz"]
for basis_set in basis_sets:
    basis = vlx.MolecularBasis.read(molecule, basis_set, ostream=None)
    norb[basis_set] = basis.get_dimensions_of_basis(molecule)

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

    # get the MOs
    scf_results = scf_drv.compute(molecule, basis)  # perform an SCF optimization
    C = scf_results["C_alpha"]  # MO coefficients
    n_elec[basis_set] = []

    for mo in range(norb[basis_set]):
        D = np.einsum("a, b -> ab", C[:, mo], C[:, mo])  # orbital density

        # determine the density on the grid points
        G = np.einsum("ab,bg->ag", D, chi_g)
        n_g = 2 * np.einsum("ag,ag->g", chi_g, G)
        n_elec[basis_set].append(np.dot(weights, n_g))

fig = plt.figure(figsize=(8, 4))
ax = plt.axes(xlim=(-1, 41), ylim=(1e-9, 1e-1))

plt.grid(True)
plt.yscale("log")

for basis_set in basis_sets:
    error = abs(np.array(n_elec[basis_set]) - 2)
    plt.scatter(range(norb[basis_set]), error, label=basis_set)

plt.legend()

ax.fill_between(
    [4.6, 41], 0, 1, color="green", alpha=0.1, transform=ax.get_xaxis_transform()
)
ax.text(15, 0.85, "Unoccupied MOs", size=14, transform=ax.get_xaxis_transform())

plt.xlabel("Molecular orbital")
plt.ylabel("Error in number of electrons")

plt.show()

The lowest five orbitals are occupied and for those we note that the largest error is seen in the numerical integration of the core 1s-orbital of oxygen and it is noted that the errors are very much alike for the two basis sets. For the unoccupied orbitals, on the other hand, the two basis sets give very different results. While the errors obtained with use of the def2-svpd basis set remain small for all virtual orbitals , they occationally reach above \(10^{-4}\) electrons for the aug-cc-pvdz basis set. Somewhat alarmingly, it is the lower virtual orbitals that suffer the largest errors at the same time as these are the most important ones in spectrum simulations.