Contents

1. Riemann sums
2. Polynomial interpolation
3. References
4. External links

1. - Riemann sums

Let f : [a,b] ⇒ ℝ be a function defined on a closed domain [a,b] of the real numbers, ℝ. The basic problem in numerical integration is to compute an approximate solution to a definite integral of the function f as a weighted sum of function values at specified points within the domain of integration [a,b]:

provided that f(q) is continuous on the closed interval [a,b] of the real axis q. The real numbers w_i, i = 1, 2, …, n, are called weights or coefficients. The points t_i, i = 1, 2, …, n, denote selected values of the variable q and usually belong to the interval [a,b]. They are called nodes.

The solution proposed in the right hand side of eq. (1.1) can be interpreted as a Riemann sum.[1Stewart, J. Calculus Early Transcendentals; Sixth Edition; Thomson Brooks/Cole: 2008. Section 5.2, p. 366., L1-3(L1) Weisstein, Eric W. "Riemann Sum." From MathWorld − A Wolfram Web Resource.
(L2) Wikipedia contributors, "Riemann sum," Wikipedia, The Free Encyclopedia.
(L3) Wikipedia contributors, "Riemann integral," Wikipedia, The Free Encyclopedia.] Let's partition the interval [a,b] into n sub-intervals [q_i,q_i+1] by means of a finite sequence of numbers of the form

The mesh size of the partition is defined to be the length of the longest sub-interval, that is,

For each interval [q_i,q_i+1] of the partition arbitrarily choose an element t_i ∈ [q_i,q_i+1] and define the Riemann sum as:

Each term in the sum is the product of the value of the function at a given point and the length of an interval. Consequently, each term represents the (signed) area of a rectangle with height f(t_i) and width w_i = q_i+1 − q_i. The Riemann sum is the (signed) area of all the rectangles. If f happens to be positive, then the Riemann sum can be interpreted as a sum of areas of approximating rectangles. If f takes on both positive and negative values, then the Riemann sum is the sum of the areas of the rectangles that lie above the q-axis and the negatives of the areas of the rectangles that lie below the q-axis.

The definite integral of f from a to b is

provided that this limit exists. If it does exist, we say that f is integrable (or more specifically Riemann-integrable) on [a,b]. The definite integral of f(q) from a to b can be interpreted as a net area, that is, the difference of the area of the region above the q-axis and below the graph of f, and the area of the region below the q-axis and above the graph of f.

2. - Polynomial interpolation

A natural way to approximate the integral (1.1) is to find a polynomial L_n−1(f; q₁,q₂,...,q_n; q)  of degree at most (n − 1) that interpolates[2Apostol, T. M. Calculus; Volume II; Second Edition; John Wiley: New York: 1969. Section 15.6, p. 579., L4Wikipedia contributors, "Polynomial interpolation," Wikipedia, The Free Encyclopedia.] f(q) in n distinct points  q₁, q₂, ..., q_n  in [a, b], i.e. such that

and then write

A polynomial L_n−1(q) of degree n − 1 in a single indeterminate (or variable) q can always be written (or rewritten) in the form:

Using eq. (2.1), we write the linear system

The square matrix in eq (2.4) is a Vandermonde matrix of order n.[L5Wikipedia contributors, "Vandermonde matrix," Wikipedia, The Free Encyclopedia.] Its determinant has a particularly simple form and is equal to the product   ∏1≤j<i≤n (q_i − q_j). Since the q_i's are all distinct, this implies that the Vandermonde determinant is always different from zero and that Vandermonde matrices are nonsingular. Therefore the system of equations (2.4) always has a solution for the coefficients  w₁, w₂, …, w_n.

More generally, given a basis l₁, l₂, …, l_n of the space of polynomials of degree less or equal to n − 1, we write the interpolating polynomial L_n−1(q) as a linear combination of the n basis polynomials

and find the coefficients  w₁, w₂, …, w_n  by solving the linear system

The monomial basis l_i = qⁱ⁻¹ (i = 1, 2, …, n), was used in eqs (2.3) and (2.4). Here, we choose, as a basis, the Lagrange polynomials of degree (n − 1)

which satisfy the condition

where δ_ij is the Kronecker delta. Because of this condition, the matrix in eq. (2.6) is equal to the identity matrix and w_i = f(q_i)  ∀ i.  Therefore, the interpolating polynomial (2.5) is given by

This makes L_n−1(q) a polynomial of degree (n − 1) with L_n−1(q_j) = n ∑i=1 δ_ij f(q_i) = f(q_j). Inserting the expression (2.9) for L_n−1(q) into eq. (2.2), gives the quadrature rule (1.1), with

References

1. Stewart, J. Calculus Early Transcendentals; Sixth Edition; Thomson Brooks/Cole: 2008. Section 5.2, p. 366. ISBN-13: 978-0-495-01166-8; ISBN-10: 0-495-01166-5.
2. Apostol, T. M. Calculus; Volume II; Second Edition; John Wiley: New York: 1969. Section 15.6, p. 579. ISBN 0-471-00007-8

External links

1. Weisstein, Eric W. "Riemann Sum." From MathWorld − A Wolfram Web Resource. https://mathworld.wolfram.com/RiemannSum.html (accessed April 11, 2021).
2. Wikipedia contributors, "Riemann sum," Wikipedia, The Free Encyclopedia, https://en.wikipedia.org/wiki/Riemann_sum (accessed April 11, 2021).
3. Wikipedia contributors, "Riemann integral," Wikipedia, The Free Encyclopedia, https://en.wikipedia.org/wiki/Riemann_integral (accessed April 11, 2021).
4. Wikipedia contributors, "Polynomial interpolation," Wikipedia, The Free Encyclopedia, https://en.wikipedia.org/wiki/Polynomial_interpolation (accessed April 11, 2021).
5. Wikipedia contributors, "Vandermonde matrix," Wikipedia, The Free Encyclopedia, https://en.wikipedia.org/wiki/Vandermonde_matrix (accessed April 11, 2021).

The following material will be moved to other sections (web pages)

1.1. - Newton-Cotes quadrature formulas. The interpolatory procedure outlined above leads to an extremely useful and straightforward family of numerical integration techniques, which are known as Newton-Cotes formulas and apply to n equispaced quadrature points or nodes in the range [a, b], such that q_k = a + (k − 1) h  for  k = 1, 2, ..., n  and step size  h ≡ (b − a)/(n − 1). Newton-Cotes formulas may be "closed" if all points in the interval [a = q₁, b = q_n] are used, "open" if only the points [q₂, q_n−1] are used, or a variation of these two. A few closed Newton-Cotes quadrature formulas for small values of n are reported in Table 1.1.

Using polynomial interpolation with polynomials of high degree over a set of equispaced points may result in an oscillating pattern that magnifies near the ends of the interpolation points and is known as Runge's phenomenon.[1Runge, C.
Zeitschrift für Mathematik und Physik 1901, 46, 224-243.] This implies that going to higher degrees does not always improve accuracy. Runge's phenomenon can be avoided by using a composite or extended rule.

Closed "extended" rules use multiple copies of lower order closed rules to build up higher order rules. For n tabulated points, using the 2-point trapezoidal rule (n − 1) times and adding the results gives

where h ≡ (b − a) / (n − 1) and q_i = a + h (i − 1).   It can be shown that the error associated with the n-point closed extended trapezoidal rule can be expressed in terms of derivatives at the endpoints, according to the Euler-Maclaurin formula

The numbers B_j that appear in eq. (1.11) are Bernoulli numbers (B₀ = 1, B₁ = −¹⁄₂, B₂ = ¹⁄₆, B₄ = −¹⁄₃₀, B₆ = ¹⁄₄₂, B₈ = −¹⁄₃₀, ... ; B₃ = B₅ = B₇ = B₉ = ... = 0), and the remainder E_m(f) is often small.

The extended trapezoidal rule is discussed in section 10 with reference to its use in PAMoC.

1.2. - Gaussian quadrature. An n-point Gaussian quadrature rule is developped by requiring that the evaluation points are the roots of a polynomial belonging to a class of orthogonal polynomials   p₀(q), p₁(q), …, p_n−1(q),   defined in some interval [a, b] and such that

1. p_i(q) is a polynomial of degree i;
2. the inner product of any two polynomials p_i(q) and p_j(q) with respect to a positive function ω(q), called a weighting function, is given by

   <pi|ω|pj> ≔ b ∫ a ω(q) pi(q) pj(q) dq = δij b ∫ a ω(q) pi2(q) dq ≕ δij <pi|ω|pi>

   (1.12)

In addition, if   <p_i|ω|p_i> = 1,   the polynomials are orthonormal.   In particular, eq. (1.12) states that

Since p₀(q) is a polynomial of degree zero, it must be a constant and does not affect the integral, so that we may write

Consider the set of roots q₁, q₂, ..., q_n of the polynomial p_n(q). It is trivial to show that n ∑ i=1 v_i p_n(q_i) = 0   for any set of v_i's, because   p_n(q_i) = 0   for all q_i. However it is possible to show that there exists a unique set of v_i's, for which all polynomials of order   1 ≤ k ≤ n   simultaneously obey the relation

This is done by solving the simultaneous set of equations

It can be shown that the determinant of this square matrix is non zero, therefore, it is invertible, i.e., a unique solution exists.

Let's now project f(q) into the set of orthogonal polynomials   p₀(q), p₁(q), …, p_n−1(q)

where

is an expansion coefficient. Integrating both sides of eq. (1.16) over the interval [a, b] yields:

where

Eq. (1.18) shows that the projection of the integrand f(q) into a basis of orthogonal polynomials occurs implicitly, with no need to explicitly evaluate the expansion coefficients. It states that, given an integer n, we can find a set of weigths w_i and abscissas q_i such that the approximation (1.1) or (1.18) is exact if the integrand f(q) is a polynomial. Then formula (1.18) is suitable for integrands that are smooth and well-behaved, and the best results will be achieved if f(q) is approximately a low order polynomial.

Notice that the weight function ω(q) is not overtly visible in the integration formula (1.18). However, using eq. (1.19) and defining

eq. (1.18) becomes

Eq. (1.21) states that, given ω(q) and given an integer n, we can arrange the choice of weights and abscissas to make the integral exact for a class of integrands "polynomials times some known function ω(q)" rather than for the usual class of integrands "polynomials." Again, the best results will be achieved if g(q) is approximately a low order polynomial. The function ω(q) can then be chosen to remove singularities from the desired integral. Eq. (1.21) represents the well known Gaussian quadrature formulas, that in the particular case of ω(q) = 1 assume the form of eq. (1.18). A few Gauss quadrature formulas for different choices of a, b, and ω(q) are reported in Table 1.2.

PAMoC uses John Burkardt's FORTRAN90 version of IQPACK routines[27Elhay, S.; Kautsky, J.
ACM Transactions on Mathematical Software 1987, 13, 399-415.] to estimate roots and weights of Gaussian quadrature rules.

1.3. - Degree of exactness of a quadrature rule. In the following, we will use I(f) to denote the integral over a finite interval [a, b] that appears in the left hand side of eq. (1.1) or (1.18), and Q_n(f) to denote a generic interpolatory quadrature rule over n points in the same interval [a, b], i.e. the summation in the right hand side of eq. (1.1) or (1.18). A quadrature rule for a function f is exact if I(f) = Q_n(f).

Definition. A quadrature rule has degree of exactness m if it renders exact results when the integrand is any polynomial of degree ≤ m but not exact for at least one polynomial of degree (m + 1).

Independently of the choice of nodes, the degree of exactness of an n-point Newton-Cotes rule is ≤ (n − 1). In other words, this rule gives exact results for polynomials of degree ≤ (n − 1).

The quadrature rule of the form (1.1) has 2n parameters: n nodes and n weights. Hence we can hope to make it exact for all polynomials of degree (2n − 1) that have 2n coefficients. We are going to proof that the degree of exactness of an n-point Gaussian quadrature formula, where the nodes are the zeros of an orthogonal polynomial p_n(q) of degree n, is (2n − 1).

Theorem. Let {p_k(q)} ∀ k ∈ [0, n] be a set of orthogonal polynomials on [a, b] with respect to the inner product (1.12). Let q_j, j = 1, 2, …, n, be zeros of polynomial p_n(q). Then the quadrature rule given by eq. (1.1) has degree of exactness (2n − 1).

Proof. Given two univariate polynomials f(q) and p_n(q), the theorem of Euclidean division of polynomials (not proved here) asserts that there exist two unique polynomials t(q) and r(q) such that

Let f(q) be a polynomial of degree (2n − 1). Then the quotient t(q) is a polynomial of degree (n − 1) = degree(f) − degree(p_n). Multiplying both sides of eq. (1.22) by the weighting function ω(q) and integrating, yields:

The green integral in eq. (1.23) is zero because the polynomial p_n(q) is orthogonal to all polynomials of degree ≤ (n − 1), according to eq. (1.12). Since the rule (1.18) is interpolatory, it is exact for r(q). Then

Here we have used the fact that q_i's are the zeros of p_n(q), i.e. p_n(q_i) = 0. Hence rule (1.18), which is exact for r(q), is exact also for f(q), which has degree (2n −1), i.e. I(f) = Q_n(f).

It remains to prove that the rule (1.18) is not exact for all polynomials of degree 2n. Let f(q) be a polynomial of degree 2n. Then t(q) is of degree n and   <t|ω|p_n> ≠ 0. So, on one hand

i.e. I(f) ≠ I(r) = Q_n(r). On the other hand

as q_i's are the zeros of p_n(q), i.e. p_n(q_i) = 0. Hence, Q_n(f) = Q_n(r) ≠ I(f) for all polynomials f of degree 2n.

1.4 - Multi-dimensional integrals. The one-dimensional quadrature rules discussed above can be directly generalized to higher dimensions. Consider the two-dimensional integral ∫∫ f(x,y) dx dy. We can integrate out one variable as follows:

Next we approximate the integral by the product quadrature formula

Similarly, the triple integral ∭f(x,y,z) dx dy dz can be approximated by the product quadrature rule

1.5 - Transformation of coordinates. Consider a mapping function f : ℝ³ → ℝ³ which takes as input the vector u ∈ ℝ³ and produces as output the vector f(u) ∈ ℝ³. (The set ℝⁿ consists of all n-tuples of real numbers and is called the “n-dimensional real space”). It is often advantageous to use other coordinates than Cartesian x, y, z system. In general, a suitable mapping function f provides transformation equations, which specify one coordinate system in terms of the other, such as (in vector notation)

where x, y, z are Cartesian coordinates, and u₁, u₂, u₃ are some other coordinates. The functions f₁, f₂ and f₃ establish a one-to-one correspondence between the coordinate systems: they have to be continuous, have continuous partial derivatives and single valued inverse. So, a point r = (x,y,z) ∈ ℝ³ can also be expressed as r = (u₁,u₂,u₃), according to eq (1.28) which defines the so-called u-substitution.

Differentiating the vector function r = (x,y,z) yields the vector dr = (dx,dy,dz) = (df₁,df₂,df₃) = df, which can be written more explicitly and with matrix notation as

where ∂ f / ∂ u_j is a vector, whose components are the partial derivatives of x = f₁, y = f₂, and z = f₃ with respect to u_j,   ∀  j = 1, 2, 3.   Left-multiplying the the vector dr by its transpose yields the following equalities

where the Jacobian matrix J of f, also denoted by J_f and ∂(f₁,f₂,f₃) / ∂(u₁,u₂,u₃), is usually defined and arranged as follows:[N1Jacobian matrix]

or, component-wise:   J_ij = ∂ f_i / ∂ u_j .   It is worth noting that across a row of the Jacobian matrix the numerators are the same and down a column the denominators are the same. In eq (1.30) the symmetric matrix   g = J^† J   is called the fundamental (or metric) tensor of the Euclidean space in curvilinear coordinates. It has elements   g_ii = | ∂ f / ∂ u_i |2   and   g_ij = ∂ f / ∂ u_i ⋅ ∂ f / ∂ u_j   for   i,j = 1, 2, 3.   Eq (1.30) gives the so-called first fundamental form or line element

What is important for calculating integrals in 3D space, is the volume element dV = dx dy dz, which can be written as the absolute value of the mixed product of the three tangent vectors in eq (1.29):

In § 2Volume Integral eq (1.33) will be used to express the Cartesian coordinates x, y, z in terms of the spherical polar coordinates r, θ, φ.

Consider now a mapping function f : ℝ → ℝ which takes as input the scalar u ∈ ℝ and produces as output the scalar f(u) ∈ ℝ. This is the case of a single variable, for which the following substitution rule holds:

In the right hand side of eq (1.34) the variable r has been replaced by f(u), i.e. the substitution r = f(u) has been applied.  f'(u) is the first derivative of f(u) with respect to u, i.e.   f'(u) = df(u) / du.   f⁻¹(u) is the inverse of u (not to be confused with f(u)⁻¹ that is the reciprocal of f(u):   f(u)⁻¹ = 1 / f(u)).

In the following, we consider some examples of change of coordinates that allow to change the limits of integration in such a way the most appropriate rules of Gauss quadrature can be applied. Let's start with the variable transformation:

From this equation:

and

so that eq (1.34) can be rewritten as:

An equivalent expression of eq (1.38) can be obtained by inserting the factor ω(u) ω(u)⁻¹ into the integrand of its right hand side:

where it has been set:

The last integral in eq (1.39) has the proper form required by the Gauss quadrature formula:

where the u_i are either the roots of the Legendre polynomials of order n if ω(u) = 1, or the roots of the Chebyshev polynomial of the second kind if ω(u) = √1 − u², and the w_i are the associated weights. The final result is obtained by combining eqs (1.37) and (1.40) with eq (1.41):

where the transformed weights, v_i, are given by

in the case of Gauss-Legendre quadrature, and

in the case of Gauss-Chebyshev quadrature of the second kind.

Consider now a slightly different linear transformation, eq (1.45), and follows the same steps as before.

Eqs (1.45) transform an integral over the variable r ∈ [a, b] into an integral over the new variable u ∈ [0, 1]:

where it has been set:

Table 1.2 suggests applying the Gauss-Gill quadrature rule to the calculation of the integral (1.46). In this case, the weight function ω(u) = ln²u is used, and yields the transformed weights:

Integrals over a semi-infinite interval [r₀, ∞] are quite common in PAMoC. Among the Gauss quadratures reported in Table 1.2, only the Gauss-Laguerre rule is defined to cover the semi-infinite interval [0, ∞] and therefore it could be used without any transformation of the variable of integration. However, its useful to introduce the linear transformation (1.49),

which leads to integration rule:

where S is a scaling factor, ω(u) = e^−u is the weight function, Φ(u) = S e^+u φ(f(u)), and the transformed weights are

2. - Volume integral

Evaluation of atomic and molecular properties requires calculation of a volume integral, which is an integral over a three dimensional domain (i.e. a region V in ℝ³) of a function f(r) = f(x,y,z) and is usually written as a triple integral:

In the Euclidean space, the components of the position vector r are the cartesian coordinates x, y and z, and the volume element dV is given by the product of the differentials of the cartesian coordinates, dV = dx dy dz.

If function f(x,y,z) is separable along the three Cartesian coordinates, the triple integral (2.1) reduces to the product of three simple integrals, each one along a single Cartesian coordinate:

Gaussian type orbitals (GTO's), as well as their products, are an example of functions which are separable along the three Cartesian coordinates, and the corresponding one-dimensional integrals over unbounded intervals [−∞, +∞] can be estimated by the Gauss-Hermite quadrature rule or calculated analytically.

It is often advantageous to use other coordinates than Cartesian x, y, z system. The latter are usually expressed in terms of the spherical polar coordinates r, θ, φ:

with  r ∈ [0,∞],  θ ∈ [0,π],  φ ∈ [0,2π], so that the Jacobian matrix has the expression:

and the volume element is:

where we have introduced the surface element

spanning from θ to θ + dθ and φ to φ + dφ on a spherical surface at (constant) radius r, and the differential of the solid angle Ω = Ω(θ, φ):

Then, the volume integral (2.1) can be written as a three-dimensional integral in terms of the spherical polar coordinates r, θ, and φ, or as a two-dimensional integral in terms of the radial coordinate r and the solid angle Ω:

It can be calculated by the one-dimensional radial integration:

of a function g(r), which provides the spherical average of f(x, y, z) through the evaluation of a two-dimensional surface integral over the unit sphere at (constant) radius r:

The integral is numerically approximated as a sum of weighted function values, sampled over a given set of points c_i = r_i, θ_i, φ_i, Ω_i:

2.1. - Partitioning into atomic subintegrals.

In a molecule, the density ρ(r) has cusps at the positions of the nuclei, so that a simple cartesian grid for the entire molecule cannot account properly for the accumulation of density at that positions. Then, following design principles put forward by Becke in 1988, the molecular integration is broken up into separate, but overlapping atomic contributions.[2Becke, A. D.
J. Chem. Phys. 1988, 88, 2547-2553.] These single-center atomic subintegrals are then computed on atomic grids. Finally, the results are summed together to give the molecular integral (2.1).
Different partitioning schemes are available in PAMoC for experimental and/or theoretical densities:

• Stewart
• Becke
• Hirshfeld
• QTAIM
• Mulliken fifty-fifty allocation algorithm
• Stone's nearest-site allocation algorithm

2.2. - Angular integral.

The great simplicity of the three-terms product formula (N_r, N_θ, N_φ), eq. (2.13), which separates the integral into three one-dimensional quadratures over the coordinates r, θ, and φ, has led Murray, Handy and Laming in 1993 to recommend it.[3Murray, C. W.; Handy, N. C.; Laming, G. J.
Mol. Phys. 1993, 78, 997-1014.] In this case, a standard choice for integration of the angular part is a spherical product formula which uses the trapezoidal rule along φ direction and Gauss-Legendre quadrature along θ direction. Usually, the angular grid is characterized by a single integer L, which describes the maximum degree of spherical harmonic which may be integrated exactly; this implies N_φ = L + 1 equally spaced points in φ space and N_θ = (L + 1)/2 Gauss-Legendre points in θ space, for a total of N_θ⋅N_φ = (L + 1)²/2 angular quadrature points. This quadrature rule is still available in PAMoC for historical reasons, since it was used in some old codes from which PAMoC evolved.

There are highly efficient algorithms for numerically integrating on the surface of a sphere that clearly outperform the alternative integration over the individual angular coordinates. Indeed, Becke in 1988[2Becke, A. D.
J. Chem. Phys. 1988, 88, 2547-2553.] and Gill, Johnson and Pople (GJP) in 1993[4Gill, P. M. W.; Johnson, B. G.; Pople, J. A.
Chem. Phys. Lett. 1993, 209, 506-512.] found that a two-terms product formula (N_r, N_Ω), which separates the integral (2.13) into two one-dimensional quadratures on the radial coordinate r and the angular coordinate Ω is more efficient than separating also the two angular integrations on θ and φ. Since then, the general consensus seems to have focused on the so-called two-dimensional Lebedev grids,[5(a) Lebedev, V.I.
USSR Comp. Math. and Math. Phys. 1975, 15(1), 44-51. Zh. vychisl. Mat. mat. Fiz. 1975, 15(1), 48-54.
(b) Lebedev, V.I.
USSR Comp. Math. and Math. Phys. 1976, 16(2), 10-24. Zh. vychisl. Mat. mat. Fiz. 1976, 16(2), 293-306.
(c) Lebedev, V.I.
Siberian. Math. J. 1977, 18(1), 99-107. Sibirskii Matematicheskii Zhurnal 1977, 18(1), 132-142.
(d) Lebedev, V. I.; Skorokhodov, A. L.
Russian Acad. Sci. Dokl. Math. 1992, 45, 587-592.
(e) Lebedev, V. I.
Russian Acad. Sci. Dokl. Math. 1995, 50, 283-286.
(f) Lebedev, V.I.; Laikov, D.N.
Dokl. Math. 1999, 59, 477-481.] which are the PAMoC method of choice to evaluate the spherical integral of Eq. (2.12). Each of Lebedev grids are characterized by the fixed number of angular quadrature points N_Ω ≈ (L + 1)²/3, which exactly integrate the spherical harmonics having up to a maximum degree L.

2.3. - Radial integral and transformations of the radial coordinate.

The general methodology to calculate the integral (2.11) is to map the interval [0, ∞] to a finite interval [a, b] of the variable q by a suitable variable transformation r = r(q)  ⇒  dr = r'(q) dq, and then to use a quadrature formula defined on this interval. The integral (2.11) takes then the form

where we have set G(q) = ω⁻¹(q) r²(q) r'(q) g[r(q)] and v_i = w_i ω⁻¹(q_i) r²(q_i) r'(q_i) . The choice of the mapping is crucial. The mapping determines how the radial points are distributed in the molecular space and if the core and the chemical bonding regions are appropriately represented in the integration. The function r(q) has to be chosen in a way that a high resolution is provided in the core and, to a lesser extent, in the valence region of the atom to ensure an accurate integration but sparse outside the atomic region to allow an economical calculation. The resolution is controlled by the value of the derivative r'(q): a small value of r'(q) indicates that a given interval for r is mapped onto a large interval of q, which implies a high resolution; consequently, a large value of r'(q) implies a low resolution. In summary, r'(q) should be small where r(q) is small and large where r(q) is large.

Since it is hard to define the detailed prescriptions a mapping function should obey, in order for the radial quadrature to be efficient, many different mapping schemes have been introduced in literature so far.[2Becke, A. D.
J. Chem. Phys. 1988, 88, 2547-2553.; 3Murray, C. W.; Handy, N. C.; Laming, G. J.
Mol. Phys. 1993, 78, 997-1014.; 4Gill, P. M. W.; Johnson, B. G.; Pople, J. A.
Chem. Phys. Lett. 1993, 209, 506-512.; 6Gill, P.M.W.; Johnson, B.J.; Pople, J.A.; Frisch, M.J.
Intern. J. Quantum Chem. Symp. 1992, 26, 319-331.; 7Treutler, O.; Ahlrichs, R.
J. Chem. Phys. 1995, 102, 346-354.; 8Mura, M. E.; Knowles, P. J.
J. Chem. Phys. 1996, 104, 9848-9858.; 9Lindh, R.; Malmqvist, P.-A.; Gagliardi, L.
Theor. Chem. Acc. 2001, 106, 178.; 10Gill, P. M. W.; Chien, S.-H.
J. Comput. Chem. 2003,24, 732-740.] Finding a mapping function consists in transforming the integral by transition to another variable of integration. Such a procedure is one of the methods for calculating an integral in one real variable, which is known as integration by substitution.

The transformed coordinate q_i can be sampled over a set of equally spaced points and used in combinations with extended forms of Newton-Cotes quadrature formulas (cf Table 1.1). On the other hand, the evaluation points q_i can be chosen as the roots of a polynomial which is orthogonal in the same interval of the variable q with respect to a weighting function ω(q). A representative set of these orthogonal polynomials is given in Table 1.2.

A number of different mappings of the radial coordinate r, defined over the semi-infinite interval [r₀, ∞] or the finite interval [r₀, r_max], onto a new coordinate q, defined over the interval [a, b] with values chosen as either equispaced points or roots of polynomials orthogonal in the same interval [a, b], is listed in Table 2.1. Additional information can be found in Sections 3−8 and 10.

2.3.4. - Standardization of quadrature grids. In order to compare the roots and weights from the various quadrature schemes, Gill and Chien[10Gill, P. M. W.; Chien, S.-H.
J. Comput. Chem. 2003, 24, 732-740.] chose the R values so that the middle point of each quadrature is unity, i.e. r(q_k) − r₀ = 1 for k = (n+1)/2. They called the resulting values the standardized roots and weights, and the standardizing factor was defined as:

However, such a definition makes R dependent on n and is ambigous when n is even. For this reason, and in a more general way, PAMoC uses a standardization factor σ, which is defined in terms of the value of the coordinate q at the center of the interval spanned by the variable q itself:

As a consequence, the standardizing factor, σ, depends only on the type of transformation, and not on how the points q_i are chosen nor on their number n. In addition, PAMoC defines R as the product of σ and the parameter R_cov, which eventually can be used to adjust the grid to fit different atomic shapes:

Typically, R_cov is chosen as the covalent radius of the atom (see next section).

As shown in Table 2.1, MultiExp, Knowles, and Handy transformations map the radial coordinate r ∈ [r₀, ∞] into the finite interval [0, 1], whose midpoint is q_center = ¹⁄₂. Among the qudrature rules available in PAMoC, the Gauss-Gill and the extended trapezoidal rules have roots in the range [0, 1], and therefore can be used in combination with the MultiExp, Knowles and Handy grids, which are standardized by the following factors:

Becke and Ahlrichs transformations map the radial coordinate into the finite interval [−1, +1], whose midpoint is q_center = 0 and the standardization factor is σ_Becke = σ_Ahlrichs = 1. The roots of Legendre and Chebyshev polynomials fall in the same interval, occur symmetrically about the origin and include the origin itself, so that combining the Gauss-Legendre and Gauss-Chebyshev quadratures with Becke and Ahlrichs transformations yield grids a which are standardized implicitly.

The Laguerre grid spans a semi-infinite interval for which q_center cannot be defined. Therefore its standardizing factor is calculated by eq. (2.15) and depends on the value of the middle root and n:

2.4. - Scale factor of radial quadrature points: the atomic radius. The nodes of the quadrature rules are proportional to arbitrary scale factors R (not necessarily the same) and the corresponding weights are proportional to R³. The value of this parameter corresponds to the midpoint of the integration interval at q = 0. The parameter R is an atomic radius which measures the radial extent of the atom in question. In Becke's grid R is chosen as half of the Bragg-Slater radius[16Bragg, W. L.
Philos. Mag. 1920, 40, 169-189., 17Slater, J. C.
J. Chem. Phys. 1964, 41, 3199-3204.] of the respective atom, except for hydrogen in which case the factor 1/2 is not applied. Original values of Bragg-Slater radii[17Slater, J. C.
J. Chem. Phys. 1964, 41, 3199-3204.] are reported in Table 2.2. Revisited values[18Cordero, B.; Gómez, V.; Platero-Prats, A. E.; Revés, M.; Echeverría, J.; Cremades, E.; Barragán, F.; Alvarez, S.
Dalton Trans. 2008, 2832-2838.] are given in Table 2.3.

In the so-called SG-1 grid by GJP,[4Gill, P. M. W.; Johnson, B. G.; Pople, J. A.
Chem. Phys. Lett. 1993, 209, 506-512.] the parameter R was chosen as the abscissa of the maximum of the radial probability function 4 πr²χ²(r) of the valence atomic orbital χ(r) given by Slater's well known rules.[19Slater, J. C.
Phys. Rev. 1930, 36, 57-64.] The atomic radii which follow from this definition are given in Table 2.4 for the atoms H to Ar. The atomic radii for the atoms K to Rn are those determined by Clementi[20Clementi, E.; Raimondi, D. L.; Reinhardt, W. P.
J. Chem. Phys. 1967, 47, 1300-1307.] from minimal-basis-set SCF functions and reported in Table 2.5. The atomic radii reported in Tables 2.2 - 2.5 are compared in Figure 2.1.

Figure 2.1: Covalent radius (bohr) as a function of the atomic number: (a) on mouse out, for atoms H to Ar, (b) on mouse over, for atoms H to Cm.

2.5. - Grid pruning. In order to cut down the number of grid points and hence to increase the efficiency of the numerical integration, a technique called grid pruning is frequently used. The underlying idea is that as one approaches the nucleus, the electron density loses more and more its angular structure and becomes increasingly spherically symmetric. Hence, for spheres at small distances from the nucleus a progressively smaller amount of angular grid points should suffice. Similar arguments apply if we analyze the situation at large values of r. Here, the actual magnitude of ρ(r) becomes so small that again we can get away with much less sophisticated angular grids without loosing any significant accuracy. In a pruned grid one exploits these observations and the space around each atom is partitioned into various regions. Within these regions, whose sizes of course depend on the actual atom, Lebedev grids of varying density are used. Close to the nucleus fairly coarse grids are sufficient, while dense ones are employed at intermediate distances and again coarser grids as we move farther away from the nucleus. This technique is used in most modern density functional programs.

For example, in 1993 Gill, Johnson and Pople proposed that the quadrature grids obtained by the combination of the N_r-point Euler-Maclaurin radial grid and the N_Ω-point Lebedev angular grid, viz. EML(N_r,N_Ω), or the N_θ×N_φ-point spherical product grid, viz. EMSP(N_r,N_θ,N_φ) could be reduced in size, while maintaining their effectiveness, by allowing the number of points N_Ω or N_θ×N_φ to become dependent on r, so that different angular grids can be used on different concentric spheres.[4Gill, P. M. W.; Johnson, B. G.; Pople, J. A.
Chem. Phys. Lett. 1993, 209, 506-512.] Given four numbers {α₁, α₂, α₃, α₄} and an atomic radius R, the four spheres of radius α₁R, α₂R, α₃R, and α₄R obviously partition an atom into five regions. The size of each region depends on the central atom. The angular grids used in these regions, together with the appropriate α-parameters determined for the atoms H to Kr, are given in Table 2.6. The so-called standard grid SG-1 was designed to give numerical integration errors of about 0.2 kcal mol⁻¹ for medium sized molecules, while using as few grid points as possible. The grid is derived from the EML(50,194) grid, which has 50 radial points, given by the Euler-Maclaurin rules, and 194 angular points positioned by the Lebedev rules. The SG-1 grid partitions space into five spherical regions around the atom and then uses Lebedev grids with 6, 38, 86, 194 and 86 points respectively. This produces about 2500 points per atom, roughly a quarter the size of the EML(50,194) grid, yet yielding similar (within a few μ-Hartree) accuracy. Similarly, the fairly dense integration mesh EML(75,302) contains 75 radial shells and 302 angular points per shell. Due to the pruning the actual number of integration points per atom is reduced to only around 7000, just a third of the regular size of 75 × 302 = 22650 points.

More recently, in 2006 Chien and Gill reported "the development of a new standard quadrature grid, SG-0, designed to be approximately half as large as, and to provide approximately half the accuracy of, the established SG-1 grid".[26Chien, S.-H.; Gill, P. M. W.
J. Comput. Chem. 2006, 27, 730-739.] It is based on MultiExp[10Gill, P. M. W.; Chien, S.-H.
J. Comput. Chem. 2003, 24, 732-740.] and Lebedev quadrature for radial and angular coordinates, respectively, and can be referred to as MEL(23,170) for first-row atoms and MEL(26,170) for the second row elements. The pruning process is more highly resolved in SG-0 than in SG-1. "The SG-0 grid for each of the first- and second-row elements studied is given in Table 2.7. Each grid is defined by the number N_r of radial grid points, the radial scale factor R, and the number of angular grid points at each of the radial points. The total number N_tot of grid points on each atom is also listed. The SG-0 grid for any element not listed in Table 2.7 is defined to be the same as the SG-1 grid for that element. In the second column of Table 2.7, the notation x^y indicates that the x-point Lebedev grid is used at y successive radial points. For example, the SG-0 grid for the hydrogen atom has 23 radial points and the configuration 6⁶, 18³, 26¹, 38¹, 74¹, 110¹, 146⁶, 86¹, 50¹, 38¹, 18¹ indicates that a six-point Lebedev angular grid is used at each of the six innermost radial points, an 18-point angular grid at each of the next three radial points, followed by a single 26-point grid, a single 38-point grid, and so forth."

3. Lebedev angular quadrature

This section briefly discusses the evaluation of the angular integral (2.12) using Lebedev quadrature:

where the particular grid points (θ_i, φ_i, or Ω_i) and grid weights, w_i, are to be determined. It's worth recalling that the angular coordinate Ω is related to cartesian coordinates r ≡ {x, y, z}, via Ω = r / |r|. In other words, Ω represents the Cartesian coordinates of a vector pointing to the surface of a unit sphere. Angular integration, therefore, is the integration of a function living on the surface of a unit sphere (centered at a constant point r). The orthogonal (with respect to the weight function ω(Ω) = 1) polynomials living on the surface of a unit sphere are special functions called spherical harmonics. A set of quadrature weights and roots for spherical harmonics could be obtained by the usual procedure to generate Gauss quadrature formulas.[24 Anderson, D.G.
Math. Comp. 1965, 19, 477-481.; 25Golub, G.H.; Welsch, J.H.
Math. Comp. 1969, 23, 221-230.] For an ideal formula, the integration points should obviously be distributed over the surface of the sphere as uniformly as possible. A perfectly uniform distribution is provided by the vertices or the face centroids of a regular polyhedron. Lebedev and Laikov developed a highly efficient set of quadrature points which are invariant under the octahedral rotation group with inversion. The resulting sets of weights and points are referred to as Lebedev angular quadrature grids and are generally regarded as being superior (number of points vs. accuracy) to other angular quadrature schemes.[5(a) Lebedev, V.I.
USSR Comp. Math. and Math. Phys. 1975, 15(1), 44-51. Zh. vychisl. Mat. mat. Fiz. 1975, 15(1), 48-54.
(b) Lebedev, V.I.
USSR Comp. Math. and Math. Phys. 1976, 16(2), 10-24. Zh. vychisl. Mat. mat. Fiz. 1976, 16(2), 293-306.
(c) Lebedev, V.I.
Siberian. Math. J. 1977, 18(1), 99-107. Sibirskii Matematicheskii Zhurnal 1977, 18(1), 132-142.
(d) Lebedev, V. I.; Skorokhodov, A. L.
Russian Acad. Sci. Dokl. Math. 1992, 45, 587-592.
(e) Lebedev, V. I.
Russian Acad. Sci. Dokl. Math. 1995, 50, 283-286.
(f) Lebedev, V.I.; Laikov, D.N.
Dokl. Math. 1999, 59, 477-481.] Each of Lebedev grids are characterized by the fixed number of angular quadrature points N_Ω ≈ (L + 1)²/3, which exactly integrate the spherical harmonics having up to a maximum degree L. The sets of fixed sizes, called rules, which are available in PAMoC, are listed in Table 9.1. The corresponding code is distributed through CCL (fortran source).

3. Gauss-Legendre quadrature

Gauss-Legendre quadrature provides an approximate solution to the definite integral of a function f(q) over the interval [−1, 1] as a weighted sum of function values at specified points within the domain of integration:

A n-point quadrature rule as above will only produce accurate results if the function f(q) is well approximated by a polynomial function within the range [−1, 1]. It is constructed to yield an exact result for polynomials of degree (2n − 1) or less, by a suitable choice of the points q_i and weights w_i for i = 1, ..., n. The evaluation points q_i for quadrature order n are just the roots of the Legendre polynomials P_n(q), which occur symmetrically about 0. In other words, with the n-th polynomial normalized to give P_n(1) = 1, the i-th node, q_i, is the i-th root of P_n, with associated weight w_i.

The integral of a function f(r) over the interval [a, b] can still be evaluated by the Gauss-Legendre formula, after a suitable variable substitution like the following:

which changes the integral over [a, b] into an integral over [-1, 1] in the following way:

where n is the number of quadrature points used, the q_i are the roots of the Legendre polynomial P_n(q) with associated weights w_i, and the   ν_i = b − a / 2 w_i   are the transformed weights.

The following paragraphs describe two main applications of the Gauss-Legendre quadrature rule in PAMoC.

theta quadrature over the interval [0, π]. In the case of explicit user request to estimate the angular integral of eq. (2.12) by the two-term (N_θ, N_φ) product formula ("intgrid 8 | 9 | 10" and "intgrid -mmmlll"), PAMoC employes the Gauss-Legendre formula for the theta quadrature over the interval [0, π]. To this aim, the variable substitution

changes the integral over θ ∈ [0, π] into an integral over q ∈ [−1, 1] in the following way

Radial integral. The Gauss-Legendre quadrature rule can also be used by PAMoC to estimate the radial integral (2.11), in conjunction with a suitable variable transformation to map the semi-infinite interval spanned by the radial coordinate r into the finite interval [−1, +1] of the variable q. This is accomplished by two transformations of the radial coordinate, due to Becke[2Becke, A. D.
J. Chem. Phys. 1988, 88, 2547-2553.] and Ahlrichs,[7Treutler, O.; Ahlrichs, R.
J. Chem. Phys. 1995, 102, 346-354.] which are described in section 7 with more details. In addition, the linear transformation (3.2) is also available. The appropriate weights for the Gauss-Legendre quadrature are given by eqs. (3.6), (3.7), and (3.8), respectively,

where q_i and w_i are roots and weights of the Gauss-Legendre quadrature rule, and r_i is defined by eqs. (8.4), (8.8), and (3.2), respectively.

The n-point Gauss-Legendre quadrature rule integrates exactly all spherical harmonics of degree L = 2n − 1 or less, so that, for a given L, a number of n = (L + 1)/2 integration points is needed.
In PAMoC, radial integration by the Gauss-Legendre quadrature rule combined with Becke and Ahlrichs transformations of the radial coordinate, can be selected by the keywords "radrule 15", "radrule 16", and "radrule 17", respectively.

Becke, Ahlrichs, and linear grids, evaluated using the roots (and the associated weights) of the Legendre polynomials, are reported in Tables 3.2 and 3.3 for n = 11, and are compared in Figure 8.1 with the analogous grids obtained using the roots (and the associated weights) of the Chebyshev polynomials of the second kind.