# 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:

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:

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)
```

Several other quantities can also be evaluated on the grid. For more
details, please refer to: `compute_grid_density_dm()`

,
`compute_grid_gradient_dm()`

,
`compute_grid_kinetic_dm()`

,
`compute_grid_hartree_dm()`

,
`compute_grid_esp_dm()`

,
and
`compute_grid_orbitals_exp()`

.

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.

```
from horton 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
print 'EXPECTATION VALUE OF |R|:', grid.integrate(rho, 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:

Expanding the orbitals in a local basis set results in:

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

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:

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)
```

Similar methods, like
`compute_grid_gradient_fock()`

and
`compute_grid_kinetic_fock()`

,
are available for the following two chain rules:

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