6.1. How to use HORTON for numerical integration?

6.1.1. Building an integration grid

HORTON primarily makes use of Becke-Lebedev grids for numerical integrations over molecular volumes. It can automatically construct Becke-Lebedev integration grids for a given molecular geometry. This is needed for DFT computations (typically to evaluate the exchange-correlation functional) or atoms-in-molecules analyses. To familiarize yourself with these grids, please refer to [becke1988_multicenter] and [lebedev1999].

In a nutshell, a Becke-Lebedev integration grid works as follows. Say, you are interested in computing the integral of \(f(\mathbf{r})\) over a molecular volume. In practice, the integrand is often a functional of the electron density, so similar to the electron density, it contains sharp spikes close to the atomic nuclei. In order to properly integrate all these unsmooth spikes, the molecular integral is first split into atomic contributions:

\[\int f(\mathbf{r}) d\mathbf{r} = \sum_A \int w_A(\mathbf{r}) f(\mathbf{r}) d\mathbf{r}\]

where \(w_A(\mathbf{r})\) is the atomic weight function for atom A. This function is 1 close to the nucleus of atom A, and goes to zero upon leaving the domain of atom A. Every atomic integral is then computed on a grid in a spherical coordinate system centered on the nucleus of the corresponding atom. This atomic grid is typically a product grid, where different one-dimensional radial grids are possible, and the Lebedev-Laikov grids are always used for the angular part. Putting all these together, you can always approximate the numerical integration as follows:

\[\int f(\mathbf{r}) d\mathbf{r} \approx \sum_{i=1}^{N_\text{grid}} w_i f(\mathbf{r}_i)\]

where \(N_\text{grid}\) is the number of grid points, \(w_i\) are the integration grid weights and \(\mathbf{r}_i\) are the integration grid points.

The default Becke-Lebedev grid is constructed by specifying three non-optional arguments, as follows:

grid = BeckeMolGrid(coordinates, numbers, pseudo_numbers)

where coordinates is an array containing the Cartesian coordinates of the atoms, numbers is an array containing the (integer) atomic numbers of the atoms, and pseudo_numbers is an array containing the (float) effective core charges of the atoms. These arrays are always available when you load a molecule from a file. For example:

mol = IOData.from_file('water.xyz')
grid = BeckeMolGrid(mol.coordinates, mol.numbers, mol.pseudo_numbers)

The fourth optional argument can be provided to control the accuracy of the integration grid, e.g.:

grid = BeckeMolGrid(coordinates, numbers, pseudo_numbers, 'fine')

The available levels of accuracy for the built-in numerical integration grids are documented in Atomic integration grids. You can also control in more detail the radial and angular components of the integration grids. For more details, please refer to the API documentation of horton.grid.molgrid.BeckeMolGrid and horton.grid.atgrid.AtomicGridSpec.

6.1.2. Computing an integral involving the electron density

This section assumes that the following objects are already available:

  • obasis: an orbital basis set object
  • dm_full: a spin-summed density matrix
  • grid: a Becke-Lebedev integration grid as introduced above.

If you are not familiar with the obasis and dm_full objects, please refer to Specifying the basis set and How to use HORTON as a Hartree-Fock/DFT program, respectively. Note that the density matrix can also be loaded from a file instead of being computed by HORTON. This is demonstrated in an example at the end of this section.

First, you must evaluate the integrand on all the points of the integration grid. In case of the electron density, this can be done as follows:

rho = obasis.compute_grid_density_dm(dm_full, grid.points)

It is assumed that dm_full is symmetric. Several other quantities can also be evaluated on the grid. For more details, please refer to:

Integrating the electron density by itself results in the total number of electrons. This is a simple way to verify the accuracy of the integration grid.

print grid.integrate(rho)

Since rho is simply a Numpy array, it can be manipulated easily to compute functionals of the electron density, e.g.

print grid.integrate(rho**(4.0/3.0))

You can also use the grid.points array to evaluate other expectation values numerically, e.g. the following snippet evaluates the expectation value of \(\vert\mathbf{r}\vert=(x^{2}+y^{2}+z^{2})^{0.5}\):

r = (grid.points[:,0]**2 + grid.points[:,1]**2 + grid.points[:,2]**2)**0.5
grid.integrate(rho, r)

As shown in the above snippet, the integrate method can take multiple one-dimensional arrays that are all multiplied before integration.

The following script is a complete example for computing the expectation value of \(\vert\mathbf{r}\vert=(x^{2}+y^{2}+z^{2})^{0.5}\) for a molecular wave-function loaded from a file.

data/examples/grid/expectation_r.py
from horton import *  # pylint: disable=wildcard-import,unused-wildcard-import


# Load the Gaussian output from file from HORTON's test data directory.
fn_fchk = context.get_fn('test/water_sto3g_hf_g03.fchk')
# Replace the previous line with any other fchk file, e.g. fn_fchk = 'yourfile.fchk'.
mol = IOData.from_file(fn_fchk)

# Specify the integration grid
grid = BeckeMolGrid(mol.coordinates, mol.numbers, mol.pseudo_numbers)

# Get the spin-summed density matrix
dm_full = mol.get_dm_full()

# Compute the density on the grid
rho = mol.obasis.compute_grid_density_dm(dm_full, grid.points)

# Compute the expectation value of |r|.
r = (grid.points[:, 0]**2 + grid.points[:, 1]**2 + grid.points[:, 2]**2)**0.5
expt_r = grid.integrate(rho, r)
if log.do_medium:
    log('EXPECTATION VALUE OF |R|: {0}'.format(expt_r))

6.1.3. Constructing a one-body operator from a real-space potential

This section assumes that the following objects are already available:

  • obasis: an orbital basis set object
  • lf: an instance of DenseLinalgFactory or CholeskyLinalgFactory
  • dm_full: a spin-summed density matrix
  • grid: a Becke-Lebedev integration grid as introduced above.

If you are not familiar with the obasis or lf object, please refer to Specifying the basis set. The density matrix can either be read from a file or computed with HORTON. For more information, please refer to How to use HORTON as a Hartree-Fock/DFT program.

Given a multiplicative potential, its expectation value is written as:

\[\langle V \rangle = \int \rho(\mathbf{r}) V(\mathbf{r}) d\mathbf{r}.\]

Expanding the orbitals in a local basis set results in:

\[\langle V \rangle = \sum_{\mu\nu} D_{\mu\nu} \mathcal{V}_{\nu\mu}\]

where \(D_{\mu\nu}\) is the spin-summed density matrix. The matrix \(\mathcal{V}_{\nu\mu}\) is defined as

\[\mathcal{V}_{\nu\mu} = \int V(\mathbf{r}) b_\nu^*(\mathbf{r}) b_\mu(\mathbf{r}) d\mathbf{r}\]

where \(b_\mu(\mathbf{r})\) are the orbital basis functions. Such matrices can be constructed with the compute_grid_density_fock() method. This method is also useful when applying the chain rule to construct the contribution of a density functional to a Fock matrix:

\[\frac{\partial E[\rho]}{\partial D_{\nu\mu}} = \int \frac{\delta E[\rho]}{\delta \rho(\mathbf{r})} b_\nu^*(\mathbf{r}) b_\mu(\mathbf{r}) d\mathbf{r}\]

The usage pattern is as follows:

# Construct some potential, e.g. a hyperbolic well
rsq = grid.points[:,0]**2 + grid.points[:,1]**2 + grid.points[:,2]**2
pot = np.sqrt(1 + rsq)

# Allocate an output array for the operator
fock = lf.create_two_index()

# Actual computation
obasis.compute_grid_density_fock(grid.points, grid.weights, pot, fock)

Other chain rules are also implemented:

  • compute_grid_gradient_fock()

    \[\frac{\partial E[\nabla\rho]}{\partial D_{\nu\mu}} = \int \frac{\delta E[\nabla\rho]}{\delta \nabla\rho(\mathbf{r})} \cdot \left( \nabla b_\nu^*(\mathbf{r}) b_\mu(\mathbf{r}) + b_\nu^*(\mathbf{r}) \nabla b_\mu(\mathbf{r}) \right) d\mathbf{r}\]

  • compute_grid_gradient_fock() just combines the density and gradient chain rules. This is more efficient than computing the separately.

  • compute_grid_kinetic_fock()

    \[\frac{\partial E[\tau]}{\partial D_{\nu\mu}} = \frac{1}{2}\int \frac{\delta E[\tau]}{\delta \tau(\mathbf{r})} \nabla b_\nu^*(\mathbf{r}) \nabla b_\mu(\mathbf{r}) d\mathbf{r}\]

    where \(\tau(\mathbf{r})\) is the positive kinetic energy density:

    \[\tau(\mathbf{r}) = \frac{1}{2} \sum_{\mu\nu} D_{\mu\nu} \nabla b_\nu^*(\mathbf{r}) \nabla b_\mu(\mathbf{r})\]

  • compute_grid_hessian_fock()

    \[\frac{\partial E[\nabla\nabla\rho]}{\partial D_{\nu\mu}} = \int \frac{\delta E[\nabla\nabla\rho]}{\delta \nabla\nabla\rho(\mathbf{r})} \colon \left( \nabla \nabla b_\nu^*(\mathbf{r}) b_\mu(\mathbf{r}) + 2 \nabla b_\nu^*(\mathbf{r}) \nabla b_\mu(\mathbf{r}) + b_\nu^*(\mathbf{r}) \nabla \nabla b_\mu(\mathbf{r}) \right) d\mathbf{r}\]

  • compute_grid_mgga_fock() just combines several of the previous chain rules: density, gradient, laplacian (trace of the Hessian) and kinetic energy density.