6.1. How to use HORTON for numerical integration?¶
6.1.1. Building an integration grid¶
HORTON primarily makes use of BeckeLebedev grids for numerical integrations over molecular volumes. It can automatically construct BeckeLebedev integration grids for a given molecular geometry. This is needed for DFT computations (typically to evaluate the exchangecorrelation functional) or atomsinmolecules analyses. To familiarize yourself with these grids, please refer to [becke1988_multicenter] and [lebedev1999].
In a nutshell, a BeckeLebedev 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 onedimensional radial grids are possible, and the LebedevLaikov 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 BeckeLebedev grid is constructed by specifying three nonoptional 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 builtin 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 objectdm_full
: a spinsummed density matrixgrid
: a BeckeLebedev 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 HartreeFock/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:
compute_grid_density_dm()
compute_grid_gradient_dm()
compute_grid_gga_dm()
compute_grid_kinetic_dm()
compute_grid_hessian_dm()
compute_grid_mgga_dm()
compute_grid_hartree_dm()
compute_grid_esp_dm()
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
onedimensional 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 wavefunction loaded from a file.
from horton import * # pylint: disable=wildcardimport,unusedwildcardimport
# 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 spinsummed 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 onebody operator from a realspace potential¶
This section assumes that the following objects are already available:
obasis
: an orbital basis set objectlf
: an instance ofDenseLinalgFactory
orCholeskyLinalgFactory
dm_full
: a spinsummed density matrixgrid
: a BeckeLebedev 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 HartreeFock/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 spinsummed 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)
Other chain rules are also implemented:

\[\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.
\[\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})\] 
\[\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.