A survey of the numerical integration methods used in PAMoC is presented in this section of the manual.
The basic problem in numerical integration is to compute an approximate solution to a definite integral of a function f(q) as a weighted sum of function values at specified points within the domain of integration:
b ∫ a f(q) dq ≈ n ∑ i=1 w_{i} f(q_{i}) | (1.1) |
A natural way to approximate the integral (1.1) is to find a polynomial L_{n}(f; q_{1}, q_{2}, ..., q_{n}; q) of degree at most (n − 1) that interpolates f(q) in n distinct points q_{1}, q_{2}, ..., q_{n} in [a, b], i.e. such that
L_{n}(q_{i}) = f(q_{i}) ∀ i = 1, 2, ..., n | (1.2) |
and then write
b ∫ a f(q) dq ≈ b ∫ a L_{n}(f; q_{1}, q_{2}, ..., q_{n}; q) dq | (1.3) |
Given a basis l_{1}, l_{2}, ..., l_{n} of the space of polynomials of degree less or equal to (n − 1), we write
L_{n}(q) = a_{1} l_{1}(q) + a_{2} l_{2}(q) + ... + a_{n} l_{n}(q) | (1.4) |
and find the coefficients a_{1}, a_{2}, ..., a_{n} by solving the linear system
| (1.5) |
We choose, as a basis, the Lagrange polynomials of degree (n − 1)
l_{i}(q) =
n
∏
k=1 k≠i q − q_{k} q_{i} − q_{k} | (1.6) |
which satisfy the condition
l_{i}(q_{j}) = δ_{ij} | (1.7) |
where δ_{ij} is the Kronecker delta. Because of this condition, the matrix in eq. (1.5) is equal to the identity matrix and a_{i} = f(q_{i}) ∀ i. Therefore, the interpolating polynomial (1.4) is given by
L_{n}(q) = n ∑ i=1 l_{i}(q) f(q_{i}) | (1.8) |
This makes L_{n}(q) a polynomial of degree (n − 1) with L_{n}(q_{j}) = n ∑ i=1 δ_{ij} f(q_{i}) = f(q_{j}). Inserting the expression (1.8) for L_{n}(q) into eq. (1.3), gives the quadrature rule (1.1), with
w_{i} = b ∫ a l_{i}(q) dq | (1.9) |
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_{1}, b = q_{n}] are used, "open" if only the points [q_{2}, 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.
n | L_{n}(q) = n ∑ i=1 l_{i}(q) f_{i} | b ∫ a L_{n}(q) dq = n ∑ i=1 f_{i} b ∫ a l_{i}(q) dq | formula type |
---|---|---|---|
2 | q − q_{2} q_{1} − q_{2} f_{1} + q − q_{1} q_{2} − q_{1} f_{2} = f_{2} − f_{1} q_{2} − q_{1} q + q_{2} f_{1} − q_{1} f_{2} q_{2} − q_{1} | 1 2 f_{2} − f_{1} q_{2} − q_{1} (q_{2}^{2} − q_{1}^{2}) + q_{2} f_{1} − q_{1} f_{2} q_{2} − q_{1} (q_{2} − q_{1}) = 1 2 (q_{2} − q_{1}) (f_{1} + f_{2}) = 1 2 h (f_{1} + f_{2}) | trapezoidal rule |
3 | (q − q_{2}) (q − q_{3}) (q_{1} − q_{2}) (q_{1} − q_{3}) f_{1} + (q − q_{1}) (q − q_{3}) (q_{2} − q_{1}) (q_{2} − q_{3}) f_{2} + (q − q_{1}) (q − q_{2}) (q_{3} − q_{1}) (q_{3} − q_{2}) f_{3} | 1 3 h (f_{1} + 4 f_{2} + f_{3}) | Simpson's rule |
4 | (q − q_{2}) (q − q_{3}) (q − q_{4}) (q_{1} − q_{2}) (q_{1} − q_{3}) (q_{1} − q_{4}) f_{1} + (q − q_{1}) (q − q_{3}) (q − q_{4}) (q_{2} − q_{1}) (q_{2} − q_{3}) (q_{2} − q_{4}) f_{2} + (q − q_{1}) (q − q_{2}) (q − q_{4}) (q_{3} − q_{1}) (q_{3} − q_{2}) (q_{3} − q_{4}) f_{3} + (q − q_{1}) (q − q_{2}) (q − q_{3}) (q_{4} − q_{1}) (q_{4} − q_{2}) (q_{4} − q_{3}) f_{4} | 3 8 h (f_{1} + 3 f_{2} + 3 f_{3} + f_{4}) | Simpson's 3/8 rule |
5 | (q − q_{2}) (q − q_{3}) (q − q_{4}) (q − q_{5}) (q_{1} − q_{2}) (q_{1} − q_{3}) (q_{1} − q_{4}) (q_{1} − q_{5}) f_{1} + (q − q_{1}) (q − q_{3}) (q − q_{4}) (q − q_{5}) (q_{2} − q_{1}) (q_{2} − q_{3}) (q_{2} − q_{4}) (q_{2} − q_{5}) f_{2} + (q − q_{1}) (q − q_{2}) (q − q_{4}) (q − q_{5}) (q_{3} − q_{1}) (q_{3} − q_{2}) (q_{3} − q_{4}) (q_{3} − q_{5}) f_{3} + (q − q_{1}) (q − q_{2}) (q − q_{3}) (q − q_{5}) (q_{4} − q_{1}) (q_{4} − q_{2}) (q_{4} − q_{3}) (q_{4} − q_{5}) f_{4} + (q − q_{1}) (q − q_{2}) (q − q_{3}) (q − q_{4}) (q_{5} − q_{1}) (q_{5} − q_{2}) (q_{5} − q_{3}) (q_{5} − q_{4}) f_{5} | 2 45 h (7 f_{1} + 32 f_{2} + 12 f_{3} + 32 f_{4} + 7 f_{5}) | Boole's rule |
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
b = q_{n} ∫ a = q_{1} f(q) dq = n−1 ∑ i=1 q_{i+1} ∫ q_{i} f(q) dq ≈ 1 2 n−1 ∑ i=1 (q_{i+1} − q_{i}) (f_{i} + f_{i+1}) = h ( 1 2 f_{1} + n−1 ∑ i=2 f_{i} + 1 2 f_{n} ) | (1.10) |
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
b = q_{n} ∫ a = q_{1} f(q) dq − h ( 1 2 f_{1} + n−1 ∑ i=2 f_{i} + 1 2 f_{n} ) = m ∑ k=1 B_{2k} (2k)! h^{2k} [ f^{ (2k−1)}(a) − f^{ (2k−1)}(b) ] + E_{m}(f) | (1.11) |
The numbers B_{j} that appear in eq. (1.11) are Bernoulli numbers (B_{0} = 1, B_{1} = −^{1}⁄_{2}, B_{2} = ^{1}⁄_{6}, B_{4} = −^{1}⁄_{30}, B_{6} = ^{1}⁄_{42}, B_{8} = −^{1}⁄_{30}, ... ; B_{3} = B_{5} = B_{7} = B_{9} = ... = 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_{0}(q), p_{1}(q), …, p_{n−1}(q), defined in some interval [a, b] and such that
<p_{i}|ω|p_{j}> ≔ b ∫ a ω(q) p_{i}(q) p_{j}(q) dq = δ_{ij} b ∫ a ω(q) p_{i}^{2}(q) dq ≕ δ_{ij} <p_{i}|ω|p_{i}> | (1.12) |
In addition, if <p_{i}|ω|p_{i}> = 1, the polynomials are orthonormal. In particular, eq. (1.12) states that
b ∫ a ω(q) p_{i}(q) p_{0}(q) dq = 0 |
Since p_{0}(q) is a polynomial of degree zero, it must be a constant and does not affect the integral, so that we may write
b ∫ a ω(q) p_{i}(q) dq = 0 ∀ i > 0. | (1.13) |
Consider the set of roots q_{1}, q_{2}, ..., 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
b ∫ a ω(q) p_{k}(q) dq = n ∑ i=1 v_{i} p_{k}(q_{i}) = 0. | (1.14) |
This is done by solving the simultaneous set of equations
| (1.15) |
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_{0}(q), p_{1}(q), …, p_{n−1}(q)
f(q) = ω(q) n ∑ k=0 c_{k} p_{k}(q) | (1.16) |
where
c_{k} = b ∫ a f(q) p_{k}(q) dq b ∫ a ω(q) p_{k}^{2}(q) dq | (1.17) |
is an expansion coefficient. Integrating both sides of eq. (1.16) over the interval [a, b] yields:
b ∫ a f(q) dq = n ∑ k=0 c_{k} b ∫ a ω(q) p_{k}(q) dq eq. (1.14) = n ∑ k=0 c_{k} n ∑ i=1 v_{i} p_{k}(q_{i}) = n ∑ i=1 v_{i} n ∑ k=0 c_{k} p_{k}(q_{i}) eq. (1.16) = n ∑ i=1 w_{i} f(q_{i}) | (1.18) |
where
w_{i} = v_{i} ω(q_{i}) | (1.19) |
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
f(q) = ω(q) g(q) | (1.20) |
eq. (1.18) becomes
b ∫ a ω(q) g(q) dq ≈ n ∑ i=1 v_{i} g(q_{i}) | (1.21) |
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.
ω(q) | Interval [a,b] |
Orthogonal polynomials | weights | radrule |
---|---|---|---|---|
1 | [-1, 1] | n-th Legendre polynomial P_{n}(q) = 1 2^{n} ⌊n/2⌋ ∑ k=0 (−1)^{k} (2n − 2k)! k! (n − k)! (n − 2k)! q^{n−2k} |
w_{i} =
2
(1 − q_{i}^{2})
[P'_{n}(q_{i})]^{2}
(i = 1,2,…,n) |
10 |
(1 − q^{2})^{½} | [-1, 1] | n-th Chebyshev polynomial of the second kind
U_{n}(q) = ⌊n/2⌋ ∑ k=0 (-1)^{k} ( n − k k ) (2q)^{n−2k} with analytical roots: q_{i} = cos ( i π n + 1 ) (i = 1,2,…, n) |
w_{i} =
π
n + 1
(1 − q_{i}^{2})
(i = 1,2,…, n) |
20 |
ln^{2}q | [0, 1] | Gill polynomials Q_{0}(q) = 1, Q_{1}(q) = 8q − 1, Q_{2}(q) = 7992q^{2} − 4104q + 217, … |
w_{i} (i = 1,2,…,n) |
30 |
e^{−q} | [0, ∞] | n-th Laguerre polynomial L_{n}(q) = n ∑ k=0 (−1)^{k} k! ( n k ) q^{k} |
w_{i} =
q_{i}
(n + 1)^{2}
[L_{n+1}(q_{i})
]^{2}
(i = 1,2,…, n) |
40 |
q^{α} e^{−q}, α > −1 | [0, ∞] | n-th Generalized Laguerre polynomial L_{n}^{(α)}(q) = n ∑ k=0 ( n + α n − k) (−1)^{k} k! q^{k} |
w_{i} (i = 1,2,…,n) |
50 |
e^{−q2} | [−∞, ∞] | n-th Hermite polynomial H_{n}(q) = (−1)^{n} e^{q²} d^{n} dq^{n} e^{−q²} |
w_{i} =
2^{n−1} n!
√π
n^{2}
[H_{n−1}(q_{i})
]^{2}
(i = 1,2,…, n) |
60 |
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
f(q) = p_{n}(q) t(q) + r(q) and degree(r) ≤ (n − 1) | (1.22) |
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:
b ∫ a f(q) ω(q) dq = b ∫ a p_{n}(q) t(q) ω(q) dq + b ∫ a r(q) ω(q) dq = b ∫ a r(q) ω(q) dq | (1.23) |
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
b ∫ a r(q) ω(q) dq = n ∑ i=1 w_{i} r(q_{i}) = n ∑ i=1 w_{i} p_{n}(q_{i}) t(q_{i}) + n ∑ i=1 w_{i} r(q_{i}) = n ∑ i=1 w_{i} [p_{n}(q_{i}) t(q_{i}) + r(q_{i})] = n ∑ i=1 w_{i} f(q_{i}) | (1.24) |
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
b ∫ a f(q) ω(q) dq = b ∫ a p_{n}(q) t(q) ω(q) dq + b ∫ a r(q) ω(q) dq ≠ b ∫ a r(q) ω(q) dq |
i.e. I(f) ≠ I(r) = Q_{n}(r). On the other hand
n ∑ i=1 w_{i} f(q_{i}) = n ∑ i=1 w_{i} [p_{n}(q_{i}) t(q_{i}) + r(q_{i})] = n ∑ i=1 w_{i} r(q_{i}) |
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:
∫ f(x,y) dy ≈ n_{y} ∑ j=1 w_{y,j} f(x,y_{j}) | (1.25) |
Next we approximate the integral by the product quadrature formula
∫∫ f(x,y) dx dy ≈ n_{x} ∑ i=1 n_{y} ∑ j=1 w_{x,i} w_{y,j} f(x_{i},y_{j}) | (1.26) |
Similarly, the triple integral ∭f(x,y,z) dx dy dz can be approximated by the product quadrature rule
∫∫∫ f(x,y,z) dx dy dz ≈ n_{x} ∑ i=1 n_{y} ∑ j=1 n_{z} ∑ k=1 w_{x,i} w_{y,j} w_{z,k} f(x_{i},y_{j},z_{k}) | (1.27) |
1.5 - Transformation of coordinates. Consider a mapping function f : ℝ^{3} → ℝ^{3} which takes as input the vector u ∈ ℝ^{3} and produces as output the vector f(u) ∈ ℝ^{3}. (The set ℝ^{n} 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)
r(x,y,z) =
| (1.28) |
where x, y, z are Cartesian coordinates, and u_{1}, u_{2}, u_{3} are some other coordinates. The functions f_{1}, f_{2} and f_{3} 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) ∈ ℝ^{3} can also be expressed as r = (u_{1},u_{2},u_{3}), 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_{1},df_{2},df_{3}) = df, which can be written more explicitly and with matrix notation as
dr =
| (1.29) |
where ∂ f ∂ u_{j} is a vector, whose components are the partial derivatives of x = f_{1}, y = f_{2}, and z = f_{3} with respect to u_{j}, ∀ j = 1, 2, 3. Left-multiplying the the vector dr by its transpose yields the following equalities
dr^{†} dr = dr ⋅ dr = dr^{2} = du^{†}( ∂ f ∂ u )† ( ∂ f ∂ u ) du = du^{†} J^{†} J du = du^{†} g du | (1.30) |
where the Jacobian matrix J of f, also denoted by J_{f} and ∂(f_{1},f_{2},f_{3}) ∂(u_{1},u_{2},u_{3}), is usually defined and arranged as follows:[N1Jacobian matrix]
J =
∂ f
∂ u
=
(
∂ f
∂ u_{1}
∂ f
∂ u_{2}
∂ f
∂ u_{3}
) =
| (1.31) |
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
dr^{2} = dx^{2} + dy^{2} +
dz^{2}
= g_{11} du_{1}^{2} + g_{22} du_{2}^{2} + g_{33} du_{3}^{2} + 2 g_{12} du_{1} du_{2} + 2 g_{13} du_{1} du_{3} + 2 g_{23} du_{2} du_{3} | (1.32) |
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):
dV = | ∂ f ∂ u_{1} ⋅ ( ∂ f ∂ u_{2} × ∂ f ∂ u_{3} ) | du_{1} du_{2} du_{3} = det( J ) du_{1} du_{2} du_{3} = √det( g ) du_{1} du_{2} du_{3} | (1.33) |
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:
b ∫ a φ(r) dr = f^{−1}(b) ∫ f^{−1}(a) φ(f(u)) f'(u) du | (1.34) |
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^{−1}(u) is the inverse of u (not to be confused with f(u)^{−1} that is the reciprocal of f(u): f(u)^{−1} = 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:
r = b − a 2 u + a + b 2 ≡ f(u) | (1.35) |
From this equation:
f'(u) = df(u) du = b − a 2 | (1.36) |
and
f^{−1}(r) = 2r − a − b b − a ≡ u ⇒ f^{−1}(a) = − 1, f^{−1}(b) = + 1 | (1.37) |
so that eq (1.34) can be rewritten as:
b ∫ a φ(r) dr = +1 ∫ −1 b − a 2 φ(f(u)) du | (1.38) |
An equivalent expression of eq (1.38) can be obtained by inserting the factor ω(u) ω(u)^{−1} into the integrand of its right hand side:
b ∫ a φ(r) dr = +1 ∫ −1 b − a 2 ω(u) ω(u)^{−1} φ(f(u)) du = +1 ∫ −1 ω(u) Φ(u) du | (1.39) |
where it has been set:
Φ(u) = b − a 2 ω(u)^{−1} φ(f(u)) | (1.40) |
The last integral in eq (1.39) has the proper form required by the Gauss quadrature formula:
b ∫ a φ(r) dr ≈ n ∑ i=1 w_{i} Φ(u_{i}) | (1.41) |
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^{2}, 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):
b ∫ a φ(r) dr ≈ n ∑ i=1 v_{i} φ(r_{i}) | (1.42) |
where the transformed weights, v_{i}, are given by
v_{i} = b − a 2 w_{i} | (1.43) |
in the case of Gauss-Legendre quadrature, and
v_{i} = (b − a)^{2} 4[(r − a)(b − r)]^{1/2} w_{i} | (1.44) |
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.
r = a + (b − a) u ≡
f(u) ⇓ f'(u) = b − a, f^{−1}(r) = r − a b − a ≡ u, f^{−1}(a) = 0, f^{−1}(b) = 1 | (1.45) |
Eqs (1.45) transform an integral over the variable r ∈ [a, b] into an integral over the new variable u ∈ [0, 1]:
b ∫ a φ(r) dr = 1 ∫ 0 (b − a) φ(f(u)) du = 1 ∫ 0 (b − a) ω(u) ω(u)^{−1} φ(f(u)) du = 1 ∫ 0 ω(u) Φ(u) du ≈ n ∑ i=1 w_{i} Φ(u_{i}) = n ∑ i=1 v_{i} φ(r_{i}) | (1.46) |
where it has been set:
Φ(u) = (b − a) ω(u)^{−1} φ(f(u)) | (1.47) |
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^{2}u is used, and yields the transformed weights:
v_{i} = (b − a) (ln r_{i} − a b − a )^{−2} w_{i} | (1.48) |
Integrals over a semi-infinite interval [r_{0}, ∞] 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),
r = r_{0} + S u ≡
f(u) ⇓ f'(u) = S, f^{−1}(r) = r − r_{0} S ≡ u, f^{−1}(r_{0}) = 0, f^{−1}(∞) = ∞ | (1.49) |
which leads to integration rule:
∞ ∫ r_{0} φ(r) dr = ∞ ∫ 0 S φ(f(u)) du = ∞ ∫ 0 e^{−u} e^{+u} S φ(f(u)) du = ∞ ∫ 0 e^{−u} Φ(u) du ≈ n ∑ i=1 w_{i} Φ(u_{i}) = n ∑ i=1 v_{i} f(r_{i}) | (1.50) |
where S is a scaling factor, ω(u) = e^{−u} is the weight function, Φ(u) = S e^{+u} φ(f(u)), and the transformed weights are
v_{i} = S exp( r_{i} − r_{0} S ) w_{i} | (1.51) |
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 ℝ^{3}) of a function f(r) = f(x,y,z) and is usually written as a triple integral:
I = ∫ V f(r) d^{3}r = ∫ V f(r) dV = ∭ V f(x,y,z) dx dy dz | (2.1) |
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:
∀ f(x,y,z) = g_{X}(x) ⋅
g_{Y}(y) ⋅ g_{Z}(z)
⇓ ∭ V f(x,y,z) dx dy dz = x_{2} ∫ x_{1} g_{X}(x) dx ⋅ y_{2} ∫ y_{1} g_{Y}(y) dy ⋅ z_{2} ∫ z_{1} g_{Z}(z) dz | (2.2) |
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, θ, φ:
x = r cosθ sinφ, y = r sinθ sinφ, z = r cosθ | (2.5) |
with r ∈ [0,∞], θ ∈ [0,π], φ ∈ [0,2π], so that the Jacobian matrix has the expression:
J(r,θ,φ) =
| (2.6) |
and the volume element is:
dV = r^{2} sinθ dr dθ dφ = dr dS = r^{2} dr dΩ | (2.7) |
where we have introduced the surface element
dS = r^{2} sinθ dθ dφ | (2.8) |
spanning from θ to θ + dθ and φ to φ + dφ on a spherical surface at (constant) radius r, and the differential of the solid angle Ω = Ω(θ, φ):
dΩ = dS r^{2} = sinθ dθ dφ | (2.9) |
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 Ω:
I = +∞ ∫ 0 r^{2} dr π ∫ 0 sinθ dθ 2π ∫ 0 f(r,θ,φ) dφ = +∞ ∫ 0 r^{2} dr ∫ Ω f(r,Ω) dΩ | (2.10) |
It can be calculated by the one-dimensional radial integration:
I = +∞ ∫ 0 r^{2} g(r) dr | (2.11) |
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:
g(r) = π ∫ 0 sinθ dθ 2π ∫ 0 f(r,θ,φ) dφ = ∫ Ω f(r,Ω) dΩ | (2.12) |
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}:
∫ V f(r) dV ≈ N_{r} ∑ i=1 w_{r,i} N_{θ} ∑ j=1 w_{θ,j} N_{φ} ∑ k=1 w_{φ,k} f(r_{i},θ_{j},φ_{k}) ≈ N_{r} ∑ i=1 w_{r,i} N_{Ω} ∑ j=1 w_{Ω,j} f(r_{i},Ω_{j}) | (2.13) |
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:
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}/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)^{2}/3, which exactly integrate the
spherical harmonics having up to a maximum degree L.
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
+∞
∫
0
r^{2} g(r) dr =
b
∫
a
ω(q) ω^{−1}(q)
r^{2}(q) r'(q)
g[r(q)] dq =
b
∫
a
ω(q) G(q) dq ≈ n ∑ i=1 w_{i} G(q_{i}) = n ∑ i=1 w_{i} ω^{−1}(q_{i}) r^{2}(q_{i}) r'(q_{i}) g[r(q_{i})] = n ∑ i=1 v_{i} g(q_{i}) | (2.14) |
where we have set G(q) = ω^{−1}(q) r^{2}(q) r'(q) g[r(q)] and v_{i} = w_{i} ω^{−1}(q_{i}) r^{2}(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_{0}, ∞] or the finite interval [r_{0}, 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.
{r_{i}} | grid | r_{i} − r_{0} | q_{i} | {q_{i}} | radrule |
---|---|---|---|---|---|
[r_{0}, ∞] | MultiExp | − R ln q_{i} | exp [−(r_{i} − r_{0})/R] | [0, 1] | 31, 71 |
[r_{0}, ∞] | Knowles | − R ln (1 − q_{i}^{k}) | ^{k}√1 − exp [−(r_{i} − r_{0})/R] | [0, 1] | 32, 72 |
[r_{0}, ∞] | Handy | R q_{i}^{m} (1 − q_{i})^{m} | ^{m}√r_{i} − r_{0} ^{m}√R − ^{m}√r_{i} − r_{0} | [0, 1] | 33, 73 |
[r_{0}, r_{max}] | Handy | (r_{max} − r_{0}) q_{i}^{m} 1 + (r_{max} − r_{0} − 2^{m})(1 − q_{i})^{m} | m√ (r_{i} − r_{0}) (r_{max} − r_{0} − 2^{m} + 1) (r_{i} − r_{0})(r_{max} − r_{0} − 2^{m}) + r_{max} − r_{0} | [0, 1] | 34, 74 |
[r_{0}, ∞] | Becke | R 1 + q_{i} 1 − q_{i} | r_{i} − r_{0} − R r_{i} − r_{0} + R | [−1, +1] | 15, 25, 75 |
[r_{0}, ∞] | Ahlrichs | − R (1 + q_{i})^{α} ln 2 ln 1 − q_{i} 2 | [−1, +1] | 16, 26, 76 | |
[r_{0}, r_{max}] | linear | r_{max} − r_{0} 2(1 + q_{i}) | 2 r_{i} − (r_{max} + r_{0}) r_{max} − r_{0} | [−1, +1] | 17, 27, 77 |
[r_{0}, r_{max}] | linear | (r_{max} − r_{0}) q_{i} | r_{i} − r_{0} r_{max} − r_{0} | [0, 1] | 37, 77 |
[r_{0}, ∞] | linear | R q_{i} | r_{i} − r_{0} R | [0, ∞] | 47, 57, 77 |
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_{0} = 1 for k = (n+1)/2.
They called the resulting values the standardized roots and weights, and
the standardizing factor was defined as:
R = 1 r[q_{(n+1)/2}] − r_{0} | (2.15) |
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:
σ = 1 r(q_{center}) − r_{0} | (2.16) |
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:
R = σ R_{cov} | (2.17) |
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_{0}, ∞] into the finite interval [0, 1], whose midpoint is q_{center} = ^{1}⁄_{2}. 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:
σ_{MutiExp} = − 1 ln q_{center} = 1 ln 2 = 1.442695 | |
σ_{Knowles} = − 1 ln (1 − q^{k}_{center}) (e.g.: σ_{Log1} = 1.442695 and σ_{Log3} = 7.488876) | |
σ_{Handy} = ( 1 − q_{center} q_{center})^{m} = 1 |
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:
σ_{Laguerre} = 1 q_{(n+1)/2} |
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^{3}. 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.
H | He | ||||||||||||||||
0.25 | |||||||||||||||||
Li | Be | B | C | N | O | F | Ne | ||||||||||
1.45 | 1.05 | 0.85 | 0.70 | 0.65 | 0.60 | 0.50 | |||||||||||
Na | Mg | Al | Si | P | S | Cl | Ar | ||||||||||
1.80 | 1.50 | 1.25 | 1.10 | 1.00 | 1.00 | 1.00 | |||||||||||
K | Ca | Sc | Ti | V | Cr | Mn | Fe | Co | Ni | Cu | Zn | Ga | Ge | As | Se | Br | Kr |
2.20 | 1.80 | 1.60 | 1.40 | 1.35 | 1.40 | 1.40 | 1.40 | 1.35 | 1.35 | 1.35 | 1.35 | 1.30 | 1.25 | 1.15 | 1.15 | 1.15 | |
Rb | Sr | Y | Zr | Nb | Mo | Tc | Ru | Rh | Pd | Ag | Cd | In | Sn | Sb | Te | I | Xe |
2.35 | 2.00 | 1.80 | 1.55 | 1.45 | 1.45 | 1.35 | 1.30 | 1.35 | 1.40 | 1.60 | 1.55 | 1.55 | 1.45 | 1.45 | 1.40 | 1.40 | |
Cs | Ba | La | Hf | Ta | W | Re | Os | Ir | Pt | Au | Hg | Tl | Pb | Bi | Po | At | Rn |
2.60 | 2.15 | 1.95 | 1.55 | 1.45 | 1.35 | 1.35 | 1.30 | 1.35 | 1.35 | 1.35 | 1.50 | 1.90 | 1.80 | 1.60 | 1.90 | ||
↓ | |||||||||||||||||
Ce | Pr | Nd | Pm | Sm | Eu | Gd | Tb | Dy | Ho | Er | Tm | Yb | Lu | ||||
1.85 | 1.85 | 1.85 | 1.85 | 1.85 | 1.85 | 1.80 | 1.75 | 1.75 | 1.75 | 1.75 | 1.75 | 1.75 | 1.75 |
H | He | ||||||||||||||||
0.31 | 0.28 | ||||||||||||||||
Li | Be | B | C | N | O | F | Ne | ||||||||||
1.28 | 0.96 | 0.84 | 0.76 | 0.71 | 0.66 | 0.57 | 0.58 | ||||||||||
Na | Mg | Al | Si | P | S | Cl | Ar | ||||||||||
1.66 | 1.41 | 1.21 | 1.11 | 1.07 | 1.05 | 1.02 | 1.06 | ||||||||||
K | Ca | Sc | Ti | V | Cr | Mn | Fe | Co | Ni | Cu | Zn | Ga | Ge | As | Se | Br | Kr |
2.03 | 1.76 | 1.70 | 1.60 | 1.53 | 1.39 | 1.39 | 1.32 | 1.26 | 1.24 | 1.32 | 1.22 | 1.22 | 1.20 | 1.19 | 1.20 | 1.20 | 1.16 |
Rb | Sr | Y | Zr | Nb | Mo | Tc | Ru | Rh | Pd | Ag | Cd | In | Sn | Sb | Te | I | Xe |
2.20 | 1.95 | 1.90 | 1.75 | 1.64 | 1.54 | 1.47 | 1.46 | 1.42 | 1.39 | 1.45 | 1.44 | 1.42 | 1.39 | 1.39 | 1.38 | 1.39 | 1.40 |
Cs | Ba | La | Hf | Ta | W | Re | Os | Ir | Pt | Au | Hg | Tl | Pb | Bi | Po | At | Rn |
2.44 | 2.15 | 2.07 | 1.75 | 1.70 | 1.62 | 1.51 | 1.44 | 1.41 | 1.36 | 1.36 | 1.32 | 1.45 | 1.46 | 1.48 | 1.40 | 1.50 | 1.50 |
↓ | |||||||||||||||||
Ce | Pr | Nd | Pm | Sm | Eu | Gd | Tb | Dy | Ho | Er | Tm | Yb | Lu | ||||
2.04 | 2.03 | 2.01 | 1.99 | 1.98 | 1.98 | 1.96 | 1.94 | 1.92 | 1.92 | 1.89 | 1.90 | 1.87 | 1.87 |
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^{2}χ^{2}(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.
H | He | ||||||
1.0000 | 0.5882 | ||||||
Li | Be | B | C | N | O | F | Ne |
3.0769 | 2.0513 | 1.5385 | 1.2308 | 1.0256 | 0.8791 | 0.7692 | 0.6838 |
Na | Mg | Al | Si | P | S | Cl | Ar |
4.0909 | 3.1579 | 2.5714 | 2.1687 | 1.8750 | 1.6514 | 1.4754 | 1.3333 |
H | He | ||||||||||||||||
a_{0} | 0.31 | ||||||||||||||||
Li | Be | B | C | N | O | F | Ne | ||||||||||
1.67 | 1.12 | 0.87 | 0.67 | 0.56 | 0.48 | 0.42 | 0.38 | ||||||||||
Na | Mg | Al | Si | P | S | Cl | Ar | ||||||||||
1.90 | 1.45 | 1.18 | 1.11 | 0.98 | 0.88 | 0.79 | 0.71 | ||||||||||
K | Ca | Sc | Ti | V | Cr | Mn | Fe | Co | Ni | Cu | Zn | Ga | Ge | As | Se | Br | Kr |
2.43 | 1.94 | 1.84 | 1.76 | 1.71 | 1.66 | 1.61 | 1.56 | 1.52 | 1.49 | 1.45 | 1.42 | 1.36 | 1.25 | 1.14 | 1.63 | 0.94 | 0.88 |
Rb | Sr | Y | Zr | Nb | Mo | Tc | Ru | Rh | Pd | Ag | Cd | In | Sn | Sb | Te | I | Xe |
2.65 | 2.19 | 2.12 | 2.06 | 1.98 | 1.90 | 1.83 | 1.78 | 1.73 | 1.69 | 1.65 | 1.61 | 1.56 | 1.45 | 1.33 | 1.23 | 1.15 | 1.08 |
Cs | Ba | La | Hf | Ta | W | Re | Os | Ir | Pt | Au | Hg | Tl | Pb | Bi | Po | At | Rn |
2.98 | 2.53 | 2.50 | 2.08 | 2.00 | 1.93 | 1.88 | 1.85 | 1.80 | 1.77 | 1.74 | 1.71 | 1.56 | 1.54 | 1.43 | 1.35 | 1.27 | 1.20 |
↓ | |||||||||||||||||
Ce | Pr | Nd | Pm | Sm | Eu | Gd | Tb | Dy | Ho | Er | Tm | Yb | Lu | ||||
2.50 | 2.08 | 2.00 | 1.93 | 1.88 | 1.85 | 1.80 | 1.77 | 1.74 | 1.71 | 1.56 | 1.54 | 1.43 | 1.35 |
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. |
Size of atomic regions | ||||||
---|---|---|---|---|---|---|
Atom | α_{1} | α_{2} | α_{3} | α_{4} | outer region | |
H-He | 0.2500 | 0.5000 | 1.0000 | 4.5000 | ||
Li-Ne | 0.1667 | 0.5000 | 0.9000 | 3.5000 | ||
Na-Ar | 0.1000 | 0.4000 | 0.8000 | 2.5000 | ||
K-Kr | 0.0200 | 0.1000 | 0.2000 | 3.5000 | ||
Lebedev grid | ||||||
(N_{r},N_{Ω}) | L | α_{1} | α_{2} | α_{3} | α_{4} | outer region |
(50,194) | 23 | 6 | 38 | 86 | 194 | 86 |
(75,302) | 29 | 26 | 86 | 146 | 302 | 146 |
(99,590) | 41 | 86 | 194 | 302 | 590 | 302 |
(110,770) | 47 | 110 | 230 | 350 | 770 | 350 |
(150,974) | 53 | 146 | 266 | 434 | 974 | 434 |
(185,1202) | 59 | 170 | 350 | 590 | 1202 | 590 |
(225,5720) | 131 | 434 | 590 | 1202 | 5720 | 1202 |
Spherical product grid | ||||||
(N_{r},N_{θ},N_{φ}) | L | α_{1} | α_{2} | α_{3} | α_{4} | outer region |
(50,15,30) | 29 | 50 | 128 | 200 | 450 | 200 |
(75,25,50) | 49 | 128 | 200 | 450 | 1250 | 450 |
(99,48,96) | 95 | 200 | 450 | 1250 | 4608 | 1250 |
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 {α_{1}, α_{2},
α_{3}, α_{4}} and an atomic radius
R, the four spheres of radius α_{1}R,
α_{2}R, α_{3}R, and
α_{4}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^{−1} 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.
Element | Lebedev partition | N_{r} | R | N_{tot} |
---|---|---|---|---|
H | 6^{6} 18^{3} 26^{1} 38^{1} 74^{1} 110^{1} 146^{6} 86^{1} 50^{1} 38^{1} 18^{1} | 23 | 1.30 | 1406 |
Li | 6^{6} 18^{3} 26^{1} 38^{1} 74^{1} 110^{1} 146^{6} 86^{1} 50^{1} 38^{1} 18^{1} | 23 | 1.95 | 1406 |
Be | 6^{4} 18^{2} 26^{1} 38^{2} 74^{1} 86^{1} 110^{2} 146^{5} 50^{1} 38^{1} 18^{1} 6^{2} | 23 | 2.20 | 1390 |
B | 6^{4} 26^{4} 38^{3} 86^{3} 146^{6} 38^{1} 6^{2} | 23 | 1.45 | 1426 |
C | 6^{6} 18^{2} 26^{1} 38^{2} 50^{2} 86^{1} 110^{1} 146^{1} 170^{2} 146^{2} 86^{1} 38^{1} 18^{1} | 23 | 1.20 | 1390 |
N | 6^{6} 18^{3} 26^{1} 38^{2} 74^{2} 110^{1} 170^{2} 146^{3} 86^{1} 50^{2} | 23 | 1.10 | 1414 |
O | 6^{5} 18^{1} 26^{2} 38^{1} 50^{4} 86^{1} 110^{5} 86^{1} 50^{1} 38^{1} 6^{1} | 23 | 1.10 | 1154 |
F | 6^{4} 38^{2} 50^{4} 74^{2} 110^{2} 146^{2} 110^{2} 86^{3} 50^{1} 6^{1} | 23 | 1.20 | 1494 |
Na | 6^{6} 18^{2} 26^{3} 38^{1} 50^{2} 110^{8} 74^{2} 6^{2} | 26 | 2.30 | 1328 |
Mg | 6^{5} 18^{2} 26^{2} 38^{2} 50^{2} 74^{1} 110^{2} 146^{4} 110^{1} 86^{1} 38^{2} 18^{1} 6^{1} | 26 | 2.20 | 1492 |
Al | 6^{6} 18^{2} 26^{1} 38^{2} 50^{2} 74^{1} 86^{1} 146^{2} 170^{2} 110^{2} 86^{1} 74^{1} 26^{1} 18^{1} 6^{1} | 26 | 2.10 | 1496 |
Si | 6^{5} 18^{4} 38^{4} 50^{3} 74^{1} 110^{2} 146^{1} 170^{3} 86^{1} 50^{1} 6^{1} | 26 | 1.30 | 1496 |
P | 6^{5} 18^{4} 38^{4} 50^{3} 74^{1} 110^{2} 146^{1} 170^{3} 86^{1} 50^{1} 6^{1} | 26 | 1.30 | 1496 |
S | 6^{4} 18^{1} 26^{8} 38^{2} 50^{1} 74^{2} 110^{1} 170^{3} 146^{1} 110^{1} 50^{1} 6^{1} | 26 | 1.10 | 1456 |
Cl | 6^{4} 18^{7} 26^{2} 38^{2} 50^{1} 74^{1} 110^{2} 170^{3} 146^{1} 110^{1} 86^{1} 6^{1} | 26 | 1.45 | 1480 |
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^{6}, 18^{3}, 26^{1},
38^{1}, 74^{1}, 110^{1}, 146^{6},
86^{1}, 50^{1}, 38^{1}, 18^{1}
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."
This section briefly discusses the evaluation of the angular integral (2.12) using Lebedev quadrature:
g(r) = π ∫ 0 sinθ dθ 2π ∫ 0 f(r,θ,φ) dφ = ∫ Ω f(r,Ω) dΩ ≈ n ∑ i=1 w_{i} f(r, Ω_{i}) = n ∑ i=1 w_{i} f(r, θ_{i}, φ_{i}) | (9.1) |
L | N_{Ω} | L | N_{Ω} | L | N_{Ω} | L | N_{Ω} |
---|---|---|---|---|---|---|---|
3 | 6 | 19 | 146 | 41 | 590 | 89 | 2702 |
5 | 14 | 21 | 170 | 47 | 770 | 95 | 3074 |
7 | 26 | 23 | 194 | 53 | 974 | 101 | 3470 |
9 | 38 | 25 | 230 | 59 | 1202 | 107 | 3890 |
11 | 50 | 27 | 266 | 65 | 1454 | 113 | 4334 |
13 | 74 | 29 | 302 | 71 | 1730 | 119 | 4802 |
15 | 86 | 31 | 350 | 77 | 2030 | 125 | 5294 |
17 | 110 | 35 | 434 | 83 | 2354 | 131 | 5810 |
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)^{2}/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).
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:
+1 ∫ −1 f(q) dq ≈ n ∑ i=1 w_{i} f(q_{i}) | (3.1) |
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}.
P_{n}(q) = 1 2^{n} ⌊n/2⌋ ∑ k=0 (−1)^{k} (2n − 2k)! k! (n − k)! (n − 2k)! q^{n−2k} | q_{i} | w_{i} = 2 (1 − q_{i}^{2}) [P'_{n}(q_{i})]^{2} |
---|---|---|
P_{0}(q) = 1 | ||
P_{1}(q) = q | q_{1} = 0 | w_{1} = 2 |
P_{2}(q) = (3q^{2} − 1)/2 | q_{1,2} = ∓√⅓ = ∓0.577350 | w_{1,2} = 1 |
P_{3}(q) = (5q^{3} − 3q)/2 | q_{1,3} = ∓√⅗ = ∓0.774597 | w_{1,3} = ^{5}⁄_{9} = 0.555556 |
q_{2} = 0 | w_{2} = ^{8}⁄_{9} = 0.888889 | |
P_{4}(q) = (35q^{4} − 30q^{2} + 3)/8 | q_{1,4} = ∓√^{3}⁄_{7} + ^{2}⁄_{7}√^{6}⁄_{5} = ∓0.861136 | w_{1,4} = (18 − √30)/36 = 0.347855 |
q_{2,3} = ∓√^{3}⁄_{7} − ^{2}⁄_{7}√^{6}⁄_{5} = ∓0.339981 | w_{2,3} = (18 + √30)/36 = 0.652145 | |
P_{5}(q) = (63q^{5} − 70q^{3} + 15q)/8 | q_{1,5} = ∓⅓√5 + 2√^{10}⁄_{7} = ∓0.906180 | w_{1,5} = (322 − 13√70)/900 = 0.236927 |
q_{2,4} = ∓⅓√5 − 2√^{10}⁄_{7} = ∓0.538469 | w_{2,4} = (322 + 13√70)/900 = 0.478629 | |
q_{3} = 0 | w_{3} = 128/225 = 0.568889 | |
P_{6}(q) = (231q^{6} − 315q^{4} + 105q^{2} − 5)/16 | q_{1,6} = ∓0.932470 | w_{1,6} = 0.171324 |
q_{2,5} = ∓0.661209 | w_{2,5} = 0.360762 | |
q_{3,4} = ∓0.238619 | w_{3,4} = 0.467914 | |
The floor function ⌊n/2⌋ specifies the largest integer less than or equal to n/2. |
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:
r = b − a 2 q + a + b 2 ≡ r(q) ⇒ dr = b − a 2 dq | (3.2) |
which changes the integral over [a, b] into an integral over [-1, 1] in the following way:
b ∫ a f(r) dr = +1 ∫ −1 b − a 2 f(r) dq ≈ n ∑ i=1 b − a 2 w_{i} f(r_{i}) = n ∑ i=1 ν_{i} f(r_{i}) | (3.3) |
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
q = cosθ ⇒ dq = − sinθ dθ | (3.4) |
changes the integral over θ ∈ [0, π] into an integral over q ∈ [−1, 1] in the following way
π ∫ 0 f(θ) sinθ dθ = − −1 ∫ 1 f(q) dq = 1 ∫ −1 f(q) dq ≈ n ∑ j=1 w_{j} f(q_{j}) | (3.5) |
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,
Becke: ν_{i} = 2 1 − q_{i}^{2} (r_{i} − r_{0}) r_{i}^{2} w_{i} | (3.6) |
Ahlrichs: ν_{i} = (r_{i} − r_{0}) r_{i}^{2} ( α 1 + q_{i} − 1 1 − q_{i} ln^{−1} 1 − q_{i} 2) w_{i} | (3.7) |
Linear: ν_{i} = b − a 2 r_{i}^{2} w_{i} | (3.8) |
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.
Table 3.2 − Roots of the n-point Gauss-Legendre quadrature rule, q_{i}, combined with Becke, Ahlrichs, and linear grids, r_{i} ≡ r(q_{i}). | Table 3.3 − Weights^{(a)} of the n-point Gauss-Legendre quadrature rule, w_{i}, combined with Becke, Ahlrichs, and linear grids, ν_{i}. | |||||||||
i | Gauss-Legendre | Becke^{(a)} | Ahlrichs^{(a)} | Linear^{(b)} | i | Gauss-Legendre | Becke^{(b)} | Ahlrichs^{(b)} | Linear^{(c)} | |
---|---|---|---|---|---|---|---|---|---|---|
1 | âˆ’0.9782 | 0.0110 | 0.0016 | 0.1089 | 1 | (-2)5.567 | (-6)3.450 | (-8)2.000 | (-3)3.298 | |
2 | âˆ’0.8871 | 0.0598 | 0.0227 | 0.5647 | 2 | (-1)1.256 | (-4)2.526 | (-5)2.108 | (-1)2.002 | |
3 | âˆ’0.7302 | 0.1560 | 0.0953 | 1.3492 | 3 | (-1)1.863 | (-3)3.028 | (-3)1.001 | (+0)1.696 | |
4 | âˆ’0.5191 | 0.3166 | 0.2557 | 2.4045 | 4 | (-1)2.332 | (-2)2.025 | (-2)1.420 | (+0)6.741 | |
5 | âˆ’0.2695 | 0.5754 | 0.5431 | 3.6523 | 5 | (-1)2.628 | (-1)1.080 | (-1)1.075 | (+1)1.753 | |
6 | 0.0000 | 1.0000 | 1.0000 | 5.0000 | 6 | (-1)2.729 | (-1)5.459 | (-1)5.575 | (+1)3.412 | |
7 | 0.2695 | 1.7380 | 1.6768 | 6.3477 | 7 | (-1)2.628 | (+0)2.976 | (+0)2.270 | (+1)5.295 | |
8 | 0.5191 | 3.1588 | 2.6425 | 7.5955 | 8 | (-1)2.332 | (+1)2.012 | (+0)7.977 | (+1)6.727 | |
9 | 0.7302 | 6.4116 | 4.0153 | 8.6508 | 9 | (-1)1.863 | (+2)2.103 | (+1)2.649 | (+1)6.971 | |
10 | 0.8871 | 16.7089 | 6.0694 | 9.4353 | 10 | (-1)1.256 | (+3)5.497 | (+1)9.543 | (+1)5.590 | |
11 | 0.9782 | 90.8639 | 8.8199 | 9.8911 | 11 | (-2)5.567 | (+6)1.939 | (+2)5.516 | (+1)2.723 | |
(a) R = 1 and r_{0} = 0. (b) a = 0 and b = 10. | (a) The weights are given in scientific notation, the power of ten being shown in parentheses. (b) R = 1 and r_{0} = 0. (c) a = 0 and b = 10. |
Gauss-Laguerre quadrature provides an approximate solution to the definite integral of a function e^{−q} f(q) over the interval [0, +∞] as a weighted sum of function values at specified points within the domain of integration:
+∞ ∫ 0 e^{−q} f(q) dq ≈ n ∑ i=1 w_{i} f(q_{i}) | (4.1) |
It fits all polynomials of degree (2n − 1). The evaluation points q_{i} for quadrature order n are just the roots of the Laguerre polynomials L_{n}(q), with associated weight w_{i}.
L_{n}(q) = n ∑ k=0 ( n k ) (−1)^{k} k! q^{k} | q_{i} | w_{i} = q_{i} (n + 1)^{2} [ L_{n+1}(q_{i}) ]^{2} |
---|---|---|
L_{0}(q) = 1 | ||
L_{1}(q) = −q + 1 | q_{1} = 1 | w_{1} = 1 |
L_{2}(q) = ½(q^{2} − 4q + 2) | q_{1} = 2 − √2 = 0.585786 | w_{1} = ¼(2 + √2) = 0.853553 |
q_{2} = 2 + √2 = 3.41421 | w_{2} = ¼(2 − √2) = 0.146447 | |
L_{3}(q) = ⅙(−q^{3} + 9q^{2} − 18q + 6) | q_{1} = 0.415775 | w_{1} = 0.711093 |
q_{2} = 2.29428 | w_{2} = 0.278518 | |
q_{3} = 6.28995 | w_{3} = 0.010389 | |
L_{4}(q) = 1 24 (q^{4} − 16q^{3} + 72q^{2} − 96q + 24) | q_{1} = 0.322548 | w_{1} = 0.603154 |
q_{2} = 1.74576 | w_{2} = 0.357419 | |
q_{3} = 4.53662 | w_{3} = 0.038888 | |
q_{4} = 9.39507 | w_{4} = 0.000539 | |
L_{5}(q) = 1 120 (−q^{5} + 25q^{4} − 200q^{3} + 600q^{2} − 600q + 120) | q_{1} = 0.26356 | w_{1} = 0.521756 |
q_{2} = 1.4134 | w_{2} = 0.398667 | |
q_{3} = 3.59643 | w_{3} = 0.075942 | |
q_{4} = 7.08581 | w_{4} = 0.003612 | |
q_{5} = 12.6408 | w_{5} = 0.000023 | |
L_{6}(q) = 1 720 (q^{6} − 36q^{5} + 450q^{4} − 2400q^{3} + 5400q^{2} − 4320q + 720) | q_{1} = 0.222847 | w_{1} = 0.458965 |
q_{2} = 1.188932 | w_{2} = 0.417001 | |
q_{3} = 2.992736 | w_{3} = 0.113373 | |
q_{4} = 5.775144 | w_{4} = 0.010399 | |
q_{5} = 9.837467 | w_{5} = 0.000261 | |
q_{6} = 15.98287 | w_{6} = 0.000001 |
Radial integral. PAMoC can use the Gauss-Laguerre rule to estimate the radial integral (2.11), provided that the keyword "radrule" is set to 47. Expressing the radial coordinate as a linear function of an adjustable scale factor R,
r = r_{0} + R q ≡ r(q) ⇒ dr = R dq, | (4.2) |
yields:
+∞
∫
0
r^{2} g(r) dr =
+∞
∫
0
R r^{2} e^{−q} e^{+q} g[r(q)] dq =
+∞
∫
0
e^{−q} G(q) dq ≈ ≈ n ∑ i=1 w_{i} G(q_{i}) = n ∑ i=1 w_{i} R r_{i}^{2} e^{qi} g(r_{i}) = n ∑ i=1 ν_{i} g(r_{i}) | (4.3) |
where G(q) = R r^{2} e^{q} g(r), and r_{i} and v_{i} define the Laguerre grid:
r_{i} = r_{0} + R q_{i} ν_{i} = R r_{i}^{2} e^{qi} w_{i} = R r_{i}^{2} q_{i} e^{qi} (n + 1)^{2} [L_{n+1}(q_{i})]^{2} | (4.4) |
As the number of quadrature points n increases, the value of
r_{n} increases accordingly (i.e.
r_{n} ≈ −12.58 + 3.87 n for
n > 20, cf Figure 4.1). As discussed in § 10 about the
Handy grid, also the
roots of the Laguerre polynomial comprise a number of points which are very
far from the nucleus, and which probably contribute very little to the radial
integral. Usually a distance of 10 bohr is large enough to be a good
approximation to infinity, so that points beyond this value can be neglected.
The percentage of roots lower than 10 bohr is about
[539 ln(n) − 616]/n
(that is less than 50% for n = 21) and decreases approximately as
2633/(37 + n).
This waste of evaluation points can be avoided by using a scaling factor in the
manner described in § 2.3.4.
PAMoC expresses the scaling parameter R as the product of
two factors, namely, the covalent radius of the atom, R_{cov},
times a parameter σ, the value of which is chosen so that the
the middle root of the quadrature is equal to R_{cov}. The
resulting values, for R_{cov} =
1, are the standardized roots and weights defined by Gill and
Chien.[10Gill, P. M. W.; Chien, S.-H.
J. Comput. Chem. 2003,24, 732-740.]
Figure 4.1 shows that also the value of the highest standardized root
r_{n} increases with n, but it seems to reach an
asymptotic value r_{asymp} of about 5.7 units for
n > 80. This fact suggests to scale the radial coordinate furtherly
by the factor r_{∞}/r_{asymp},
that moves the asymptote closer to r_{∞}, i.e. a
value which is large enough to be a good approximation to infinity.
In PAMoC the type of standardization of the radial coordinate is
controlled by the keyword "stdtyp", which is 0 for no standardization, 1 for
middle-root standardization, and 3 for middle-root standardization with
additional correction by the nearest integer of the multiplying factor
r_{∞}/r_{asymp}.
Gill and Chien pointed out that the Laguerre grid (4.4) is exact if G(q) = L_{2n−1}(q), which implies that
g(r) = R^{−1} r^{−2} exp(− r − r_{0} R ) L_{2n−1}( r − r_{0} R ) | (4.5) |
"The grid also works well when eq. (4.4) is approximately true, but
it becomes less satisfactory when g(r) either decays more slowly
than an exponential or contains multiple exponential components
with significantly different rates of decay."[10Gill, P. M. W.; Chien, S.-H.
J. Comput. Chem. 2003, 24, 732-740.]
GeneralizedGauss-Laguerre quadrature provides an approximate solution to the definite integral of a function q^{α} e^{−q} f(q) over the interval [0, +∞] and for some real number α > −1 as a weighted sum of function values at specified points within the domain of integration:
+∞ ∫ 0 q^{α} e^{−q} f(q) dq ≈ n ∑ i=1 w_{i} f(q_{i}) | (5.1) |
It fits all polynomials of degree (2n − 1). The evaluation points q_{i} for quadrature order n are just the roots of the generalized Laguerre polynomial L_{n}^{(α)}(q), with associated weight w_{i}.
L_{n}^{(α)}(q) = n ∑ k=0 ( n + α n − k ) (−1)^{k} k! q^{k} | q_{i} |
---|---|
L_{0}^{(α)}(q) = 1 | |
L_{1}^{(α)}(q) = −q + α + 1 | q_{1} = α + 1 |
L_{2}^{(α)}(q) = ½[q^{2} − (α + 2)(2q − α − 1)] | q_{1} = (α + 2) − √α + 2 |
q_{2} = (α + 2) + √α + 2 | |
L_{3}^{(α)}(q) = ⅙{−q^{3} + (α + 3) [3q^{2} − 3(α + 2)q + (α + 1)(α + 2)]} | … |
Radial integral. PAMoC can use the Generalized Gauss-Laguerre rule to estimate the radial integral (2.11), provided that the keyword "radrule" is set to 57. Expressing the radial coordinate as a function of an adjustable scale factor R, as in eq. (4.2), yields:
+∞
∫
0
r^{2} g(r) dr =
+∞
∫
0
R r^{2} q^{−α}
q^{α}
e^{−q} e^{+q}
g(r) dq =
+∞
∫
0
q^{α} e^{−q} G(q) dq ≈ ≈ n ∑ i=1 w_{i} G(q_{i}) = n ∑ i=1 w_{i} R r_{i}^{2} q_{i}^{−α} e^{qi} g(r_{i}) = n ∑ i=1 ν_{i} g(r_{i}) | (5.2) |
where G(q) = R r^{2} q^{−α} e^{q} g(r_{i}), and r_{i} and ν_{i} define the generalized Laguerre grid:
r_{i} = r_{0} + R q_{i} ν_{i} = R r_{i}^{2} q_{i}^{−α} e^{qi} w_{i} | (5.3) |
It is worth noting that, for r_{0} = 0 and α = 2, the expression of the quadrature weights simplifies to ν_{i} = R^{3} e^{qi} w_{i}. In this case, the generalized Laguerre grid (5.3) is exact if g(r) = R^{−3} exp(−r/R) L_{2n−1}^{(2)}(r/R).
Table 5.2 − Roots^{(a)} of the 11-point Gauss-Laguerre and Generalized Gauss-Laguerre (α = 2) quadrature rules, q_{i} and standardized r_{i} ≡ r(q_{i}). | Table 5.3 − Weights^{(a,b)} of the 11-point Gauss-Laguerre and Generalized Gauss-Laguerre (α = 2) quadrature rules, w_{i} and standardized ν_{i}. | |||||||||||||
Gauss-Laguerre | Generalized Gauss-Laguerre | Gauss-Laguerre | Generalized Gauss-Laguerre | |||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
q_{i} and r_{i} (σ = 1) |
σ = ^{1}⁄_{q6} | σ = ^{2}⁄_{q6} | q_{i} and r_{i} (σ = 1) |
σ = ^{1}⁄_{q6} | σ = ^{2}⁄_{q6} | w_{i} | σ = 1 | σ = ^{1}⁄_{q6} | σ = ^{2}⁄_{q6} | w_{i} | σ = 1 | σ = ^{1}⁄_{q6} | σ = ^{2}⁄_{q6} | |
0.1258 | 0.0168 | 0.0335 | 0.5298 | 0.0555 | 0.1666 | (-1)2.85 | (-3)5.11 | (-5)1.21 | (-5)9.66 | (-1)2.85 | (-3)5.11 | (-5)1.21 | (-5)9.66 | |
0.6654 | 0.0886 | 0.1772 | 1.4318 | 0.1501 | 0.4502 | (-1)3.90 | (-1)3.36 | (-4)7.93 | (-3)6.34 | (-1)3.90 | (-1)3.36 | (-4)7.93 | (-3)6.34 | |
1.6472 | 0.2193 | 0.4387 | 2.7533 | 0.2886 | 0.8657 | (-1)2.33 | (+0)3.28 | (-3)7.74 | (-2)6.19 | (-1)2.33 | (+0)3.28 | (-3)7.74 | (-2)6.19 | |
3.0911 | 0.4116 | 0.8232 | 4.5189 | 0.4736 | 1.4208 | (-2)7.66 | (+1)1.61 | (-2)3.80 | (-1)3.04 | (-2)7.66 | (+1)1.60 | (-2)3.80 | (-1)3.04 | |
5.0293 | 0.6697 | 1.3394 | 6.7643 | 0.7010 | 2.1268 | (-2)1.44 | (+1)5.56 | (-1)1.31 | (+0)1.05 | (-2)1.44 | (+1)5.56 | (-1)1.31 | (+0)1.05 | |
7.5099 | 1.0000 | 2.0000 | 9.5412 | 1.0000 | 3.0000 | (-3)1.52 | (+2)1.56 | (-1)3.69 | (+0)2.95 | (-3)1.52 | (+2)1.56 | (-1)3.69 | (+0)2.95 | |
10.6060 | 1.4123 | 2.8245 | 12.9259 | 1.3548 | 4.0649 | (-5)8.51 | (+2)3.87 | (-1)9.13 | (+0)7.30 | (-5)8.51 | (+2)3.87 | (-1)9.13 | (+0)7.30 | |
14.4316 | 1.9217 | 3.8434 | 17.0367 | 1.7856 | 5.3568 | (-6)2.29 | (+2)8.84 | (+0)2.09 | (+1)1.67 | (-6)2.29 | (+2)8.84 | (+0)2.09 | (+1)1.67 | |
19.1789 | 2.5538 | 5.1076 | 22.0710 | 2.3132 | 6.9397 | (-8)2.00 | (+3)1.95 | (+0)4.61 | (+1)3.69 | (-8)2.00 | (+3)1.95 | (+0)4.61 | (+1)3.69 | |
25.2177 | 3.3579 | 6.7159 | 28.4079 | 2.9774 | 8.9322 | 0 | (+3)4.39 | (+1)1.04 | (+1)8.29 | 0 | (+3)4.39 | (+1)1.04 | (+1)8.29 | |
33.4972 | 4.4604 | 8.9208 | 37.0190 | 3.8799 | 11.6398 | 0 | (+4)1.14 | (+1)2.70 | (+2)2.16 | 0 | (+4)1.14 | (+1)2.70 | (+2)2.16 | |
(a) R = 1 and r_{0} = 0. | (a) R = 1 and r_{0} = 0. (b) The weights are given in scientific notation, the power of ten being shown in parentheses. |
Gauss-Hermite quadrature provides an approximate solution to the integral of a function e^{-q²} f(q) over the interval [−∞, +∞] as a weighted sum of function values at specified points
+∞ ∫ −∞ e^{-q²} f(q) dq ≈ n ∑ i=1 w_{i} f(q_{i}) | (6.1) |
where n is the number of sample points used. The q_{i} are the roots of the Hermite polynomial H_{n}(q) with associated weights w_{i}.
H_{n}(q) = (−1)^{n}
e^{q²}
d^{n}
dq^{n}
e^{−q²} = = n! ⌊n/2⌋ ∑ k=0 (−1)^{k} (2 q)^{n−2k} k! (n − 2k)! |
q_{i} | w_{i} = 2^{n−1} n! √π n^{2} [ H_{n−1}(q_{i}) ]^{2} |
---|---|---|
H_{0}(q) = 1 | ||
H_{1}(q) = 2q | q_{1} = 0 | w_{1} = √π = 1.772454 |
H_{2}(q) = 4q^{2} − 2 | q_{1,2} = ∓√½ = ∓0.707107 | w_{1,2} = ½√π = 0.886002 |
H_{3}(q) = 8q^{3} − 12q | q_{1,3} = ∓√(^{3}⁄_{2}) = ∓1.224745 | w_{1,3} = 0.295409 |
q_{2} = 0 | w_{2} = 1.181636 | |
H_{4}(q) = 16q^{4} − 48q^{2} + 12 | q_{1,3} = ∓√[(3+√6)/2] = ∓1.650680 | w_{1,3} = 0.081313 |
q_{2,4} = ∓√[(3-√6)/2] = ∓0.524648 | w_{2,4} = 0.804914 |
Multipole moments of the overlap density between two GTO's.
Stone pointed out that the unabridged Cartesian multipole moments
of the overlap density between two cartesian Gaussian type orbitals (GTO's)
can be efficiently and accurately estimated using Gauss-Hermite
quadrature.[22Stone, A. J.
J. Chem. Theory Comput. 2005, 1, 1128-1132.]
The following integrals are needed:
S_{μ∊A,ν∊B}^{(l,m,n)}(A) = <χ_{μ}^{A}|x_{A}^{l} y_{A}^{m} z_{A}^{n}|χ_{ν}^{B}> = +∞ ∫ −∞ x_{A}^{l} y_{A}^{m} z_{A}^{n} χ^{*}_{μ}^{A}(r) χ_{ν}^{B}(r) d^{3}r, | (6.2)) |
where
χ_{μ}^{A}(r) = N_{1} x_{A}^{l1} y_{A}^{m1} z_{A}^{n1} exp(-ζ_{1} r_{A}^{2}) χ_{ν}^{B}(r) = N_{2} x_{B}^{l2} y_{B}^{m2} z_{B}^{n2} exp(-ζ_{2} r_{B}^{2}) | (6.3) |
r_{A} = r − A and r_{B} = r − B are the electron positions relative to the positions of atoms A and B on which the GTO's are centered, ζ_{1} and ζ_{2} are the exponents, and N_{1} and N_{2} are normalization constants.
Boys[23Boys, S. F.
Proc. Roy. Soc. A 1950, 200, 542-554.]
showed that the product of two 1s GTO's, one centered at A and
the other centered at B, is another Gaussian centered at
P =
ζ_{1}A +
ζ_{2}B
ζ_{1} + ζ_{2}
with exponent ζ_{3} =
ζ_{1} + ζ_{2}:
exp(−ζ_{1} r_{A}^{2}) ⋅ exp(−ζ_{2} r_{B}^{2}) = exp[− ζ_{1} ζ_{2} ζ_{3} (A − B)^{2}] ⋅ exp(− ζ_{3} r_{P}^{2}) ≡ K_{AB} exp(− ζ_{3} r_{P}^{2}), | (6.4) |
where K_{AB} = exp[− ζ_{1} ζ_{2} ζ_{3} (A − B)^{2}] and r_{P} = r − P. Since 1s GTO's are separable along the three Cartesian coordinates:
exp(− ζ_{3} r_{P}^{2}) = exp(− ζ_{3} x_{P}^{2}) ⋅ exp(− ζ_{3} y_{P}^{2}) ⋅ exp(− ζ_{3} z_{P}^{2}), | (6.5) |
the integral S_{μ∊A,ν∊B}^{(l,m,n)}(A) can be written as:
S_{μ∊A,ν∊B}^{(l,m,n)}(A) = N_{1} N_{2} K_{AB} ∏ (t; k, k_{1}, k_{2}) +∞ ∫ −∞ G_{A}(t; k, k_{1}, k_{2}) dt, | (6.6) |
where
G_{A}(t; k, k_{1}, k_{2}) = (t − A_{t})^{k+k1} (t − B_{t})^{k2} exp[−ζ_{3}(t − P_{t})^{2} = t_{A}^{k+k1} t_{B}^{k2} exp[−ζ_{3}t_{P}^{2}], | (6.7) |
for (t; k, k_{1}, k_{2}) = (x; l, l_{1}, l_{2}), (y; m, m_{1}, m_{2}), and (z; n, n_{1}, n_{2}). As the integrand G_{A}(t) doesn't exactly correspond to the integrand of the quadrature formula, we need a change of variable
q = √ζ_{3} (t − P_{t}) ⇔ t = ^{q}⁄_{√ζ3} + P_{t} ⇒ dt = (^{1}⁄_{√ζ3}) dq, | (6.8) |
which, coupled with the integration by substitution, yields for each integral in eq. (6.6):
+∞
∫
−∞
G_{A}(t) dt =
+∞
∫
−∞
(^{1}⁄_{√ζ3})
(^{q}⁄_{√ζ3}
+ P_{t} - A_{t})^{k+k1}
(^{q}⁄_{√ζ3}
+ P_{t} - B_{t})^{k2}
exp(−q^{2}) dq = = +∞ ∫ −∞ F_{A}(q) exp(−q^{2}) dq ≈ n ∑ i=1 w_{i} F_{A}(q_{i}) = = n ∑ i=1 (^{wi}⁄_{√ζ3}) (^{qi}⁄_{√ζ3} + P_{t} - A_{t})^{k+k1} (^{qi}⁄_{√ζ3} + P_{t} - B_{t})^{k2} = n ∑ i=1 (^{wi}⁄_{√ζ3}) (t_{i} − A_{t})^{k+k1} (t_{i} − B_{t})^{k2} | (6.9) |
where q_{i} are the roots of the Hermite polynomial of order n and w_{i} are the associated weights. The coordinates t_{i} are defined in terms of the roots q_{i}, according to transformation (6.8).
The n-point Gauss-Hermite rule integrates exactly all polynomials of degree (2 n − 1) or less. Since the integrand F_{A}(q) in eq. (6.9) is a polynomial of degree (k + k_{1} + k_{2}), a number n = (k + k_{1} + k_{2})/2 is needed.
It is worth noting that in the case of monocentric integrals we have
S_{μ∊A,ν∊A}^{(l,m,n)}(A) = <χ_{μ}^{A}|x_{A}^{l} y_{A}^{m} z_{A}^{n}|χ_{ν}^{A}> = N_{1} N_{2} ∏ t=x,y,z +∞ ∫ −∞ t_{A}^{k+k1+k2} exp[-ζ_{3}t_{A}^{2}] dt, | (6.10) |
where (k+k_{1}+k_{2}) = (l+l_{1}+l_{2}), (m+m_{1}+m_{2}), (n+n_{1}+n_{2}) for t = x, y, z, respectively. For even (k+k_{1}+k_{2}) an analytical solution is available for the integral
+∞ ∫ −∞ t_{A}^{k+k1+k2} exp[-ζ_{3}t_{A}^{2}] dt = ( π ζ_{3} )^{½} (k+k_{1}+k_{2}−1)!! (2 ζ_{3})^{(k+k1+k2)/2} ^{,} | (6.11) |
whereas for odd (k+k_{1}+k_{2}) the integral is equal to zero. The final expression for the Mulliken unabridged cartesian multipole moments of atom A is
m_{A}^{(l,m,n)} = ∑ μ∊A ∑ B ∑ ν∊B D_{μν}^{AB} S_{μ∊A,ν∊B}^{(l,m,n)}(A). | (6.12) |
The procedure outlined above is implemented in PAMoC routines muldma.f and stodma.f.
Radial integral. PAMoC can use the Gauss-Hermite rule to estimate the radial integral (2.11), provided that the keyword "radrule" is set to 67. Expressing the radial coordinate as a function of an adjustable scale factor R, as in eq. (4.2), yields:
^{1}⁄_{2}
+∞
∫
−∞
r^{2} g(r) dr =
+∞
∫
−∞
^{1}⁄_{2} R r^{2}
e^{−q2} e^{+q2}
g(r) dq =
+∞
∫
−∞
e^{−q2} G(q) dq
≈ ≈ n ∑ i=1 w_{i} G(q_{i}) = n ∑ i=1 ^{1}⁄_{2} w_{i} R r_{i}^{2} e^{qi2} g(r_{i}) = n ∑ i=1 ν_{i} g(r_{i}) | (6.13) |
where G(q) = ^{1}⁄_{2} R r^{2} e^{q2} g(r), and r_{i} and ν_{i} define the Hermite grid:
r_{i} = r_{0} + R q_{i} ν_{i} = ^{1}⁄_{2} R r_{i}^{2} e^{qi2} w_{i} | (6.14) |
It is worth noting that, for r_{0} = 0, the expression of the quadrature weights simplifies to ν_{i} = ^{1}⁄_{2} R^{3} e^{qi2} w_{i}. In this case, the Hermite grid (6.14) is exact if g(r) = 2 R^{−3} exp(−r^{2}/R) H_{2n−1}(r/R).
Gauss-Chebyshev quadrature of the second kind provides an approximate solution to the definite integral of a function √1 − q^{2} f(q) over the interval [−1, +1] as a weighted sum of function values at specified points within the domain of integration:
+1 ∫ −1 √1 − q^{2} f(q) dq ≈ n ∑ i=1 w_{i} f(q_{i}) | (8.1) |
It fits all polynomials of degree (2n − 1). The evaluation points q_{i} for quadrature order n are just the roots of the Chebyshev polynomial of the second kind U_{n}(q), with associated weight w_{i}.
U_{n}(q) = ⌊n/2⌋ ∑ k=0(-1)^{k} ( n − k k ) (2q)^{n−2k} | q_{i} = cos ( i n + 1π ) | w_{i} =
π
n + 1
sin^{2} (
i
n + 1π
) = π n + 1 (1 − q_{i}^{2}) |
---|---|---|
U_{0}(q) = 1 | ||
U_{1}(q) = 2 q | q_{1} = 0 | w_{1} = π/2 |
U_{2}(q) = 4 q^{2} − 1 | q_{1,2} = ±^{1}⁄_{2} | w_{1,2} = π/4 |
U_{3}(q) = 8 q^{3} − 4 q | q_{2} = 0 q_{1,3} = ±1/√2 | w_{2} = π/4 w_{1,3} = π/8 |
U_{4}(q) = 16 q^{4} − 12 q^{2} + 1 | q_{1,2,3,4} = ±√(3 ± √5)/8 | w_{1,2} = π(5 − √5)/40 w_{3,4} = π(5 + √5)/40 |
U_{5}(q) = 32 q^{5} − 32 q^{3} + 6 q | q_{3} = 0 q_{1,5} = ±0.866025 q_{2,4} = ±0.5 |
w_{3} = 0.523599 w_{1,5} = 0.130900 w_{2,4} = 0.392699 |
U_{6}(q) = 64 q^{6} − 80 q^{4} + 24 q^{2} − 1 | q_{1,6} = ±0.900969 q_{2,5} = ±0.623490 q_{3,4} = ±0.222521 |
w_{1,6} = 0.084489 w_{2,5} = 0.274333 w_{3,4} = 0.426576 |
U_{7}(q) = 128 q^{7} − 192 q^{5} + 80 q^{3} − 8 q | q_{4} = 0 q_{1,7} = ±0.923880 q_{2,6} = ±0.707107 q_{3,5} = ±0.382683 |
w_{4} = 0.392699 w_{1,7} = 0.057509 w_{2,6} = 0.196350 w_{3,5} = 0.335190 |
U_{8}(q) = 256 q^{8} − 448 q^{6} + 240 q^{4} − 40 q^{2} + 1 | q_{1,8} = ±0.939693 q_{2,7} = ±0.766044 q_{3,6} = ±0.5 q_{4,5} = ±0.173648 |
w_{1,8} = 0.040833 w_{2,7} = 0.144226 w_{3,6} = 0.261799 w_{4,5} = 0.338540 |
U_{9}(q) = 512 q^{9} − 1024 q^{7} + 672 q^{5} − 160 q^{3} + 10 q | q_{5} = 0 q_{1,9} = ±0.951057 q_{2,8} = ±0.809017 q_{3,7} = ±0.587785 q_{4,6} = ±0.309017 |
w_{5} = 0.314159 w_{1,9} = 0.029999 w_{2,8} = 0.108539 w_{3,7} = 0.205620 w_{4,6} = 0.284160 |
The floor function ⌊n/2⌋ specifies the largest integer less than or equal to n/2. |
The Chebyshev polynomials of the second kind are orthogonal with respect to the weight function √1 − q^{2} over the interval [−1, +1], i.e.:
+1 ∫ −1 U_{n}(q) U_{m}(q) √1 − q^{2} dq = π 2 δ_{nm} | (8.2) |
As a consequence:
n ∑ i=1 w_{i} = π 2 | (8.3) |
Radial integral. The Gauss-Chebyshev second kind quadrature can be used to estimate the radial integral (2.11), provided that a suitable variable transformation is given to map the improper interval [0, ∞] to a finite interval [−1, +1].
Becke's grid. In 1988, Becke introduced the
transformation:[2Becke, A. D.
J. Chem. Phys. 1988, 88, 2547-2553.]
r = r_{0} + R 1 + q 1 − q ≡ r(q) ⇒ dr = 2 R (1 − q)^{2} dq = 2 r − r_{0} 1 − q^{2} dq | (8.4) |
The displacement of the origin from 0 to r_{0} was not considered in the original paper, but has been added here for completeness. At q = − 1, r(q) = r_{0}, and as q → 1, the sample points tend toward infinity. Under this map, the radial integral becomes:
+∞
∫
0
r^{2} g(r) dr =
+1
∫
−1
2
r^{2}
(r − r_{0})
1
− q^{2}
g(r) dq =
+1
∫
−1
√1 − q^{2}
G(q) dq ≈ ≈ n ∑ i=1 w_{i} G(q_{i}) = n ∑ i=1 2 w_{i} r_{i}^{2} (r_{i} − r_{0}) (1 − q_{i}^{2})^{3/2} g(r) = n ∑ i=1 ν_{i} g(r_{i}) | (8.5) |
where G(q) = 2 r^{2} (r − r_{0}) (1 − q^{2})^{3/2} g(r), and r_{i} and ν_{i} define the Becke grid:
r_{i} = r_{0} + R 1 + q_{i} 1 − q_{i} ν_{i} = 2 r_{i}^{2} (r_{i} − r_{0}) (1 − q_{i}^{2})^{3/2} w_{i} = 2π n + 1 r_{i}^{2} (r_{i} − r_{0}) √1 − q_{i}^{2} q_{i} = cos ( i n + 1 π) ∀ i = 1, 2, ..., n | (8.6) |
At q = 0, r(0) − r_{0} = R
and, in particular,
q_{(n+1)/2} = 0 ∀ n odd,
i.e. the middle root of the Becke grid is equal to R. Thus, for
R = 1, eq. (8.6) gives
standardized roots and weights, as defined by Gill and
Chien.[10Gill, P. M. W.; Chien, S.-H.
J. Comput. Chem. 2003, 24, 732-740.].
The value of the scale factor R controls the extent of the
integrand.
Becke chose R as half of the Bragg-Slater radius
of the respective atom, except for hydrogen in which case the factor 1/2 is
not applied.[2Becke, A. D.
J. Chem. Phys. 1988, 88, 2547-2553.]
The Becke grid (8.6) is exact if G(q)
= U_{2n−1}(q), which implies
g(q) = R^{3} 4 r^{3/2} (r + R)^{3} U_{2n−1}( r − R r + R ) | (8.7) |
In general, such g(r) possess a 1/r^{3/2} singularity at the origin and decay as 1/r^{9/2} at large r.
Ahlrichs's grid. In 1995, Treutler and
Ahlrichs introduced the transformation:[7Treutler, O.; Ahlrichs, R.
J. Chem. Phys. 1995, 102, 346-354.]
r = r_{0} − R ln 2 (1 + q)^{α} ln 1 − q 2 ≡ r(q) ⇒ dr = (r − r_{0}) ( α 1 + q − 1 1 − q ln^{−1} 1 − q 2 ) dq | (8.8) |
which is normalized to map the center of the q interval onto r(0) - r = R. At q = −1, r(q) − r_{0} = 0, and as q → 1, the sample points tend toward infinity. Under this map, the radial integral becomes:
+∞
∫
0
r^{2} g(r) dr =
+1
∫
−1
r^{2} (r − r_{0})
(
α
1 + q −
1
1 − q ln^{−1}
1 − q
2
)
√1 − q^{2}
√1 −
q^{2}
g(r) dq =
+1
∫
−1
√1 − q^{2}
G(q) dq ≈ n ∑ i=1 w_{i} G(q_{i}) = n ∑ i=1 π n + 1 (1 − q_{i}^{2}) G(q_{i}) = n ∑ i=1 π n + 1 r_{i}^{2} (r_{i} − r_{0}) ( α√ 1 − q_{i} 1 + q_{i} − √ 1 + q_{i} 1 − q_{i} ln^{−1} 1 − q_{i} 2 ) g(r) = n ∑ i=1 ν_{i} g(r_{i}) | (8.9) |
where r_{i} and ν_{i}, with α = 0.6, define the Ahlrichs grid:
r_{i} = r_{0} −
R
ln 2
(1 + q_{i})^{α} ln
1 − q_{i}
2
ν_{i} = r_{i}^{2}
(r_{i} − r_{0})
(
α
1 + q_{i}
−
1
1 − q_{i}
ln^{−1}
1 − q_{i}
2
)
w_{i}
√1 − q_{i}^{2}
= π n + 1 r_{i}^{2} (r_{i} − r_{0}) [ α √ 1 − q_{i} 1 + q_{i} − √ 1 + q_{i} 1 − q_{i} ln^{−1} 1 − q_{i} 2 ] q_{i} = cos ( i n + 1 π) ∀ i = 1, 2, ..., n | (8.10) |
Linear grid. Although it is rather inappropriate for radial integration, the linear transformation (3.2) can also be used in combination with the 2nd kind Gauss-Chebyshev quadrature rule, yielding the following quadrature points and weights:
r_{i} = b − a 2 q_{i} + a + b 2 | |
ν_{i} = b − a 2 √1 − q^{2} r_{i}^{2} w_{i} = π (b − a) 2 (n + 1) r_{i}^{2} √1 − q^{2} ∀ i = 1, 2, ..., n | (8.11) |
with q_{i} = cos ( i n + 1π) |
Becke, Ahlrichs, and linear grids, evaluated using the roots (and the associated weights) of the Chebyshev polynomials of the second kind, are reported in Tables 8.2 and 8.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 Legendre polynomials.
Table 8.2 − Roots of the n-point Gauss-Chebyshev quadrature rule of the second kind, q_{i}, combined with Becke, Ahlrichs, and linear grids, r_{i} ≡ r(q_{i}). | Table 8.3 − Weights^{(a)} of the n-point Gauss-Chebyshev quadrature rule of the second kind, w_{i}, combined with Becke, Ahlrichs, and linear grids, ν_{i}. | |||||||||
i | 2nd kind Gauss-Chebyshev |
Becke^{(a)} | Ahlrichs^{(a)} | Linear^{(b)} | i | 2nd kind Gauss-Chebyshev |
Becke^{(b)} | Ahlrichs^{(b)} | Linear^{(c)} | |
---|---|---|---|---|---|---|---|---|---|---|
1 | âˆ’0.9659 | 0.0173 | 0.0033 | 0.1704 | 1 | (-2)1.754 | (-5)1.053 | (-7)1.100 | (-3)9.834 | |
2 | âˆ’0.8660 | 0.0718 | 0.0299 | 0.6699 | 2 | (-2)6.545 | (-4)3.876 | (-5)4.292 | (-1)2.937 | |
3 | âˆ’0.7071 | 0.1716 | 0.1093 | 1.4645 | 3 | (-1)1.309 | (-3)3.740 | (-3)1.391 | (+0)1.985 | |
4 | âˆ’0.5000 | 0.3333 | 0.2738 | 2.5000 | 4 | (-1)1.963 | (-2)2.239 | (-2)1.637 | (+0)7.085 | |
5 | âˆ’0.2588 | 0.5888 | 0.5581 | 3.7059 | 5 | (-1)2.443 | (-1)1.106 | (-1)1.110 | (+1)1.736 | |
6 | 0.0000 | 1.0000 | 1.0000 | 5.0000 | 6 | (-1)2.618 | (-1)5.236 | (-1)5.348 | (+1)3.272 | |
7 | 0.2588 | 1.6984 | 1.6442 | 6.2941 | 7 | (-1)2.443 | (+0)2.656 | (+0)2.063 | (+1)5.009 | |
8 | 0.5000 | 3.0000 | 2.5508 | 7.5000 | 8 | (-1)1.963 | (+1)1.632 | (+0)6.934 | (+1)6.377 | |
9 | 0.7071 | 5.8284 | 3.8201 | 8.5355 | 9 | (-1)1.309 | (+2)1.466 | (+1)2.197 | (+1)6.743 | |
10 | 0.8660 | 13.9282 | 5.6704 | 9.3301 | 10 | (-2)6.545 | (+3)2.830 | (+1)7.357 | (+1)5.697 | |
11 | 0.9659 | 57.6955 | 8.8138 | 9.8296 | 11 | (-2)1.754 | (+5)3.885 | (+2)3.485 | (+1)3.273 | |
(a) R = 1 and r_{0} = 0. (b) a = 0 and b = 10. | (a) The weights are given in scientific notation, the power of ten being shown in parentheses. (b) R = 1 and r_{0} = 0. (c) a = 0 and b = 10. |
In PAMoC, radial integration by the 2nd kind Gauss-Chebyshev quadrature rule combined with Becke, Ahlrichs, and linear transformations of the radial coordinate, can be selected by the keywords "radrule 25", "radrule 26", and "radrule 27", respectively.
The Gauss-Gill quadrature provides an approximate solution to the definite
integral of a function ln^{2} q
f(q) over the interval [0, 1] as a weighted sum of
function values at specified points within the domain of integration:[10Gill, P. M. W.; Chien, S.-H.
J. Comput. Chem. 2003, 24, 732-740.]
1 ∫ 0 ln^{2} q f(q) dq ≈ n ∑ i=1 w_{i} f(q_{i}) | (9.1) |
It fits all polynomials of degree (2n − 1). The evaluation points q_{i} for quadrature order n are the roots of the polynomial Q_{n}(q) orthogonal on the interval [0, 1] with respect to the weight function ln^{2} q, and w_{i} are the associated weights.
n | Q_{n}(q) | q_{i} | w_{i} |
---|---|---|---|
1 | 8q − 1 | q_{1} = 1/8 | w_{1} = 2.000000000000000 |
2 | 7992q^{2} − 4104q + 217 | q_{1} = 0.059850992523974 q_{2} = 0.453662520989539 |
w_{1} = 1.669136108179106 w_{2} = 0.330863891820894 |
3 | … | q_{1} = 0.036263311146964 q_{2} = 0.273148602374171 q_{3} = 0.653711089636059 |
w_{1} = 1.363830383647107 w_{2} = 0.565815459643824 w_{3} = 0.070354156709070 |
Roots and weights for n = 1-26, 30, 35, 40, 45, 50 can be obtained from Gill's website at the Australian National University (ANU). Their values for n = 1-3 are reported in Table 9.1
MultiExp grid for radial quadrature.
The Gauss-Gill "log-squared" quadrature rule (9.1) was devised to
estimate the radial integral (2.11), provided that a suitable variable
transformation is given to map the improper interval [0, ∞] to a finite
interval [0,1]. Gill and Chien proposed the transformation:[10Gill, P. M. W.; Chien, S.-H.
J. Comput. Chem. 2003, 24, 732-740.]
r = r_{0} − R ln q ⇒ dr = − R q dq | (9.2) |
At q = 1, r(q) = r_{0}, and as q → 0, the sample points tend toward infinity. Under this map, the radial integral becomes:
+∞
∫
0
r^{2} g(r) dr = R
1
∫
0
r^{2} q^{−1}
ln^{−2}q ln^{2}q
g(r) dq =
1
∫
0
ln^{2}q G(q) dq
≈ ≈ n ∑ i=1 w_{i} G(q_{i}) = R n ∑ i=1 r^{2} q_{i}^{−1} ln^{−2}q_{i} w_{i} g(r_{i}) = n ∑ i=1 ν_{i} g(r_{i}) | (9.3) |
where
G(q) = R (r_{0} − R ln q)^{2} q^{−1} ln^{−2}q g(r) = R^{3} r^{2} (r − r_{0})^{−2} e^{(r−r0)/R} g(r) |
and r_{i} and ν_{i} define the MultiExp grid:
r_{i} = r_{0} − R ln q_{i} ν_{i} = R (r_{0} − R ln q_{i})^{2} q_{i}^{−1} ln^{−2} q_{i} w_{i} = R^{3} r^{2} (r − r_{0})^{−2} e^{(r−r0)/R} w_{i} | (9.4) |
The standardized roots and weights introduced by Gill and
Chien[10Gill, P. M. W.; Chien, S.-H.
J. Comput. Chem. 2003, 24, 732-740.]
are obtained by using R =
R_{cov} σ_{Gill-Chien,}
with R_{cov} = 1
and σ_{Gill-Chien,} =
− 1/ln q_{(n+1)/2}. For
example, the standardizing factor for the 3-point quadrature formula is
σ_{Gill-Chien} = −
1/ln q_{2} = 0.770570791. The values of
σ_{Gill-Chien} depend
on the value of the middle root of the quadrature,
q_{(n+1)/2}, as well
as on the number n of evaluation points. In order to avoid such a
dependence, the value of q at the center of the interval where q
itself is defined, i.e. q_{center} =
^{1}⁄_{2}, can be used to compute the standardizing
factor as σ_{center} =
− ln^{−1}(q_{center}) =
ln^{−1}(2) = 1.442695.
The MultiExp grid (9.4) is exact if G(q) = Q_{2n−1}(q), which implies
g(r) = R^{−3} r^{−2} (r − r_{0})^{2} e^{−(r−r0)/R} Q_{2n−1}(e^{−(r−r0)/R}) | (9.5) |
For r_{0} = 0, eq. (9.5) can be written as
g(r) = R^{−3} 2n−1 ∑ k=0 c_{k} e^{−(k+1)r/R} | (9.6) |
Obviously, such g(r) are linear combinations of exponential
functions, and for this reason Gill and Chien named "MultiExp" the grid
defined by eq. (9.4). These authors concluded that "The accuracy with which
a quadrature scheme estimates the
value of an integral depends on how close the integrand is to the
set of functions for which the quadrature is exact. Our principal
interest lies in the development of efficient quadrature schemes for
the types of integral that arise in atomic and molecular calculations
using density functional theory. Because of shell structure, the
radial electron density around an atom can usually be modeled
well by a sum of exponential functions with a variety of decay
constants. It has previously been noted[3Murray, C. W.; Handy, N. C.; Laming, G.
J.
Mol. Phys. 1993, 78, 997-1014.; 7Treutler, O.; Ahlrichs, R.
J. Chem. Phys. 1995, 102, 346-354.]
that such multiexponential character may limit the usefulness of the
Laguerre grid, but it is clear that the MultiExp grid is well suited to such
problems."[10Gill, P. M. W.;
Chien, S.-H.
J. Comput. Chem. 2003, 24, 732-740.]
Gauss-Gill quadrature with Knowles, Handy, and
linear grids.
The Gauss-Gill quadrature can also be used in conjunction with any variable
transformation which maps the semi-infinite interval spanned by the radial
coordinate r into the finite interval [0, 1] of the variable q.
PAMoC can use the transformations due to Knowles[8Mura, M. E.; Knowles, P. J.
J. Chem. Phys.
1996, 104, 9848-9858.] and Handy,[3Murray, C. W.; Handy, N. C.; Laming, G.
J.
Mol. Phys. 1993, 78, 997-1014.]
which are described in section 10
(eqs. 10.14 and 10.10, respectively) with more details, and of course some
linear transformation like eq. (10.23).
The appropriate weights for the Gauss-Gill quadrature are given by equations
from (9.7) through (9.10), respectively,
Knowles: ν_{i} = k R r_{i}^{2} q_{i}^{k−1} (1 − q_{i}^{k}) ln^{2} q_{i} w_{i} | (9.7) |
Handy: ν_{i} = m (r_{i} − r_{0}) r_{i}^{2} q_{i} (1 − q_{i}) ln^{2} q_{i} w_{i} | (9.8) |
modified Handy: ν_{i} = m (r_{max} − r_{0}) r_{i}^{2} [1 + (r_{max} − r_{0} − 2^{m})(1 − q_{i})^{m−1}] q_{i}^{m−1} [1 + (r_{max} − r_{0} − 2^{m}) (1 − q_{i})^{m}]^{2} ln^{2} q_{i} w_{i} | (9.9) |
Linear: ν_{i} = (r_{max} − r_{0}) r_{i}^{2} ln^{2} q_{i} w_{i} | (9.10) |
where q_{i} and w_{i} are roots and weights of the Gauss-Gill quadrature rule, and r_{i} is defined by eqs. (10.14), (10.10), (10.20), and (10.23), respectively.
In PAMoC, radial integration by the MultiExp
quadrature is available through the keyword
"radrule 31",
which actually is the default choice.
Radial integration using Knowles, Handy, and linear grids with the roots and
weights of the Gauss-Gill quadrature rule, is afforded through the keywords
"radrule 32",
"radrule 33",
"radrule 34", and
"radrule 37",
respectively.
Roots and weights of Gauss-Gill quadrature rule combined with the 11-point MultiExp, Knowles' Log3 Handy and linear transformations of the radial coordinate are compared in Tables 9.2 and 9.3, respectively.
MultiExp^{(b)} | Knowles' Log3^{(b)} | Handy^{(c)} | ||||||||
---|---|---|---|---|---|---|---|---|---|---|
i | σ = 1 | σ = 1.1374 | σ = 1.4427 | σ = 1 | σ = 13.4747 | σ = 7.4889 | σ = 1 | σ = 1.9854 | modified Handy^{(d)} |
Linear^{(d)} |
1 | 0.0455 | 0.0517 | 0.0656 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0001 | 0.0000 | 0.0530 |
2 | 0.1237 | 0.1407 | 0.1785 | 0.0000 | 0.0007 | 0.0004 | 0.0015 | 0.0029 | 0.0021 | 0.3672 |
3 | 0.2402 | 0.2732 | 0.3465 | 0.0009 | 0.0124 | 0.0069 | 0.0116 | 0.0231 | 0.0161 | 0.9732 |
4 | 0.3995 | 0.4544 | 0.5763 | 0.0063 | 0.0846 | 0.0470 | 0.0511 | 0.1014 | 0.0680 | 1.8430 |
5 | 0.6088 | 0.6924 | 0.8783 | 0.0254 | 0.3416 | 0.1898 | 0.1710 | 0.3394 | 0.2137 | 2.9252 |
6 | 0.8792 | 1.0000 | 1.2685 | 0.0742 | 1.0000 | 0.5558 | 0.5037 | 1.0000 | 0.5645 | 4.1510 |
7 | 1.2292 | 1.3981 | 1.7734 | 0.1756 | 2.3655 | 1.3147 | 1.4234 | 2.8260 | 1.3168 | 5.4401 |
8 | 1.6912 | 1.9235 | 2.4399 | 0.3590 | 4.8379 | 2.6888 | 4.1468 | 8.2331 | 2.7247 | 6.7066 |
9 | 2.3297 | 2.6497 | 3.3611 | 0.6665 | 8.9806 | 4.9912 | 13.5684 | 26.9387 | 4.8570 | 7.8649 |
10 | 3.3044 | 3.7582 | 4.7672 | 1.1710 | 15.7793 | 8.7697 | 57.6650 | 114.4882 | 7.2214 | 8.8364 |
11 | 5.2406 | 5.9604 | 7.5606 | 2.0593 | 27.7479 | 15.4216 | 461.8325 | 916.9235 | 9.0235 | 9.5554 |
(a) R = 1, r_{0} = 0, and r_{max} = 10. (b) The first column (σ = 1) refers to the non-standardized grid, the second column refers to the grid standardized by means of the middle root q_{(n+1)/2}, and the third column refers to the grid standardized by means of the central point of the interval q_{center} = ^{1}⁄_{2}. (c) The first column refers to the grid standardized by means of the central point of the interval q_{center} = ^{1}⁄_{2}, which yields σ = 1 by definition. The second column refers to the grid standardized by means of the middle root q_{(n+1)/2}. (d) Non-standardized grid. |
The 11-point grids of Table 9.2 are illustrated in Figure 9.1. The effect of standardization is to increase the resolution. This effect is relatively moderate in the case of the MultiExp grid, which has the narrowest range of the grids reported in Table 9.2. The Handy grid is the most extended, with three evaluation points (27%) which are too far to contribute significantly to the integral.
MultiExp^{(b)} | Knowles' Log3^{(b)} | Handy^{(c)} | ||||||||
---|---|---|---|---|---|---|---|---|---|---|
i | σ = 1 | σ = 1.1374 | σ = 1.4427 | σ = 1 | σ = 13.4747 | σ = 7.4889 | σ = 1 | σ = 1.9854 | modified Handy^{(d)} |
Linear^{(d)} |
1 | (-4)1.25 | (-4)1.83 | (-1)3.74 | 0 | 0 | 0 | 0 | 0 | 0 | (-4)4.69 |
2 | (-3)1.48 | (-3)2.18 | (-1)4.45 | 0 | 0 | 0 | (-8)1 | (-8)6 | (-8)2 | (-2)6.24 |
3 | (-3)7.90 | (-2)1.16 | (-1)2.37 | 0 | (-6)4.41 | (-7)7.6 | (-6)2.66 | (-5)2.08 | (-6)6.94 | (-1)7.05 |
4 | (-2)2.92 | (-2)4.29 | (-1)8.76 | (-7)4.0 | (-4)9.76 | (-4)1.68 | (-4)1.75 | (-3)1.37 | (-4)3.98 | (+0)3.35 |
5 | (-1)8.80 | (-1)1.30 | (-1)2.64 | (-5)1.97 | (-2)4.83 | (-3)8.29 | (-3)5.63 | (-2)4.41 | (-2)1.02 | (+0)9.98 |
6 | (-1)2.37 | (-1)3.48 | (-1)7.11 | (-4)3.90 | (-1)9.54 | (-1)1.64 | (-1)1.34 | (+0)1.05 | (-1)1.63 | (+1)2.19 |
7 | (-1)6.04 | (-1)8.87 | (+0)1.81 | (-3)4.21 | (+1)1.03 | (+0)1.77 | (+0)3.00 | (+1)2.35 | (+0)1.80 | (+1)3.82 |
8 | (+0)1.53 | (+0)2.25 | (+0)4.60 | (-2)3.05 | (+1)7.47 | (+1)1.28 | (+1)7.92 | (+2)6.20 | (+1)1.33 | (+1)5.52 |
9 | (+0)4.15 | (+0)6.11 | (+1)1.25 | (-1)1.73 | (+2)4.23 | (+1)7.26 | (+3)3.20 | (+4)2.51 | (+1)5.62 | (+1)6.66 |
10 | (+1)1.37 | (+1)2.02 | (+1)4.13 | (-1)8.86 | (+3)2.17 | (+2)3.72 | (+5)3.19 | (+6)2.50 | (+2)1.14 | (+1)6.67 |
11 | (+1)8.67 | (+2)1.28 | (+2)2.60 | (+4)1.28 | (+4)1.28 | (+3)2.20 | (+8)2.67 | (+9)2.09 | (+2)1.11 | (+1)5.25 |
(a) The weights are given in scientific notation, the power of ten being shown in parentheses. (b) The first column (σ = 1) refers to the non-standardized grid, the second column refers to the grid standardized by means of the middle root q_{(n+1)/2}, and the third column refers to the grid standardized by means of the central point of the interval q_{center} = ^{1}⁄_{2}. (c) The first column refers to the grid standardized by means of the central point of the interval q_{center} = ^{1}⁄_{2}, which yields σ = 1 by definition. The second column refers to the grid standardized by means of the middle root q_{(n+1)/2}. (d) Non-standardized grid. |
The definite integral of a function f(q) over some interval [a, b] can be approximated by the trapezoidal rule (also known as the trapezoid rule or trapezium rule, see left side of Figure 10.1), which estimates the area under the curve in the interval [a, b] by a right trapezoid with vertical bases f(a) and f(b) and horizontal height h = b − a:
b ∫ a f(q)dq ≈ b − a 2 [ f(a) + f(b) ] | (10.1) |
Figure 10.1: Trapezoidal rule (on the left) and 4-point closed extended trapezoidal rule (on the right). |
As shown in the right side of Figure 10.1, a better approxination is obtained by dividing the interval into a number of subintervals of equal length. Then the trapezoidal rule is used to estimate the area under the curve in each subinterval, and the results for all subintervals are summed together. For n equispaced points, using the trapezoidal rule (n - 1) times and adding the results gives the "extended" trapezoidal rule:
b = q_{n} ∫ a = q_{1} f(q) dq = n−1 ∑ i=1 q_{i+1} ∫ q_{i} f(q) dq ≈ 1 2 n−1 ∑ i=1 (q_{i+1} − q_{i}) (f_{i} + f_{i+1}) = h ( 1 2 f_{1} + n−1 ∑ i=2 f_{i} + 1 2 f_{n} ) = n ∑ i=1 w_{i} f(q_{i}) | (10.2) |
where the evaluation points, q_{i}, and the associated weights, w_{i}, are defined by
q_{i} = a + (i − 1) h (∀ i = 1, 2, ..., n) | |
w_{1} = w_{n} = h 2 and w_{2} = w_{3} = … = w_{n − 1} = h | (10.3) |
h = b − a n − 1 |
and f(q_{i}) ≡ f_{i}, f(a) ≡ f(q_{1}) ≡ f_{1}, f(b) ≡ f(q_{n}) ≡ f_{n}. This n-point closed extended trapezoidal rule belongs to a large family of numerical integration techniques for equispaced points, known as Newton-Cotes formulas. The term "closed" means that points a and b are included in the formula, whereas "open" formulas consider only points belonging to the interval [q_{2}, q_{n−1}] and "semi-open" formulas refer to points in the intervals [q_{1}, q_{n−1}] and [q_{2}, q_{n}]. For instance, if f(q) has a singularity at the endpoint q_{n}, the last integral in the summation of Eq. (10.2) can be approximated by:
q_{n} ∫ q_{n−1} f(q) dq ≈ h f(q_{n−1}) ≡ h f_{n−1} | (10.4) |
and Eq. (10.2) yields the n-point semi-open trapezoidal formula:
b = q_{n} ∫ a = q_{1} f(q) dq ≈ h ( 1 2 f_{1} + n−2 ∑ i=2 f_{i} + 3 2 f_{n−1} ) | (10.5) |
which is de facto reduced to a (n − 1)-point formula. In order that the semi-open trapezoidal formula has n evaluation points effectively, a (n + 1)-point formula is needed with n equal intervals of length h = (b − a)/n, and weights w_{1} = ^{1}⁄_{2}h, w_{2} = w_{3} = … = w_{n−1} = h, and w_{n} = ^{3}⁄_{2}h.
The n-point trapezoid rule integrates exactly all spherical harmonics of degree L = n − 1 or less, so that, for a given L, a number of n = (L + 1) integration points is needed.
The following paragraphs describe two main applications of the extended trapezoidal rule in PAMoC.
φ quadrature over the interval [0, 2π]. In the case of explicit user request to estimate the angular integral by the two-term (N_{θ}, N_{φ}) product formula ("intgrid 8 | 9 | 10" and "intgrid -mmmlll"), PAMoC uses the closed extended trapezoidal rule for the φ quadrature over the interval [0, 2π]. In this case, using (n + 1) equispaced points, from φ_{1} = 0 to φ_{n+1} = 2π with step-size h = 2π/n, and observing that f(0) = f(2π), we have:
2π ∫ 0 f(φ)dφ ≈ 2π n n ∑ i=1 f(φ_{i}), | (10.6) |
where φ_{i} = (i − 1) (2π/n), ∀ i = 1, 2, …, n.
Radial integral. The
trapezoidal rule is also used by PAMoC to estimate the radial integral
(2.11), in conjunction with a suitable variable transformation to map the
improper interval [r_{0}, ∞] of the radial coordinate
r to the finite interval [a, b] of the variable q. Two
transformations were just designed by Handy[3Murray, C. W.; Handy, N. C.; Laming, G. J.
Mol. Phys. 1993, 78, 997-1014.] and
Knowles.[8Mura, M. E.; Knowles, P. J.
J. Chem. Phys. 1996, 104, 9848-9858.] for use
in combination with the composite trapezoid rule.
However, also Becke, Ahlrichs, and Gill-Chien transformations can be used
as well.
If the radial coordinate r is expressed as a function of the variable q ∈ [a, b], i.e. r ≡ r(q), so that dr = r'(q) dq, then the radial integral (2.11) can be solved by substitution, yielding
+∞ ∫ 0 r^{2} g(r) dr = b ∫ a r^{2}(q) r'(q) g(r) dq ≡ b ∫ a G(q) dq | (10.7) |
where the integrand is now given by G(q) = r^{2}(q) r'(q) g(r). Applying the quadrature rule (10.2), yields
b ∫ a G(q) dq ≈ n ∑ i=1 w_{i} G(q_{i}) = n ∑ i=1 ν_{i} g(r_{i}) | (10.8) |
where the ν_{i} are the transformed weights
ν_{i} = r^{2}(q_{i}) r'(q_{i}) w_{i} | (10.9) |
and the q_{i} and w_{i} have expressions like those of eq. (10.3).
Handy's grid. In 1993, Murray, Handy and Laming
introduced the transformation:[3Murray, C. W.;
Handy, N. C.; Laming, G. J.
Mol. Phys. 1993, 78,
997-1014.]
r = r_{0} + R q^{m} (1 − q)^{m} ⇒ dr = m R q^{m−1} (1 − q)^{m+1} dq | (10.10) |
The displacement of the origin from 0 to r_{0} was not considered in the original paper, but has been added here for completeness. At the origin is q = 0, r(0) = r_{0}, r'(0) = 0, and G(0) = r^{2}(0) r'(0) g[r(0)] = 0. On the other hand, as q → 1^{−}, the sample points eventually tend toward infinity and G(q) becomes singular. For this reason, a semi-open extended trapezoidal rule is needed, and the variable q is conveniently sampled over (n + 2) equispaced points, according to the relationship
q_{i} = i h = i n + 1 ∀ i = 0, 1, 2, …, n, n+1 | (10.11) |
which satisfies the conditions q_{0} = 0 and q_{n+1} = 1. However the quadrature rule comprises only n points q_{i} (i = 1, 2, …, n), because the first point does not contribute to the summation as r'(q_{0}) and G(q_{0}) are always zero (for m > 1), and the endpoint is discarded, since G(q_{n+1}) is singular. So, the n evaluation points r_{i} and the associated transformed weights ν_{i}, which define the Handy grid, are
r_{i} = r_{0} + R q_{i}^{m} (1 − q_{i})^{m} = r_{0} + R i^{m} (n − i + 1)^{m} | |
ν_{i} = w_{i} m q_{i} (1 − q_{i}) (r_{i} − r_{0}) r_{i}^{2} = w_{i} m (n + 1)^{2} i (n − i + 1) (r_{i} − r_{0}) r_{i}^{2} ∀ i = 1, 2, …, n | (10.12) |
q_{i} = i h, w_{1} = w_{2} = … = w_{n − 1} = h, w_{n} = 3 2 h, h = 1 n + 1, m ≥ 2 |
At the grid mid-point, if n is odd, we have i =
(n + 1)/2,
q_{(n + 1)/2} = ^{1}⁄_{2}, and
r_{(n + 1)/2} - r_{0}
R = 1, which states that the Handy grid is
standardized, according to the definition by Gill and Chien.[10Gill, P. M. W.; Chien, S.-H.
J. Comput. Chem. 2003, 24, 732-740.]
If m is equal to one, then r'(q_{0}) is not null. Indeed in this case r'(q) = R/(1 − q), r'(q_{0}) = R, and ν_{0} is proportional to R r_{0}^{2}. Then, if r_{0} is not null, an (n+1)-point semi-open extended trapezoid quadrature provides the n evaluation points q_{i} = (i − 1)h ∀ i = 1, 2, …, n, and the associated weights w_{1} = ^{1}⁄_{2} h, w_{2} = w_{3} = … = w_{n − 1} = h, w_{n} = ^{3}⁄_{2} h, with h = ^{1}⁄_{n}.
Gill and Chien showed[10Gill,
P. M. W.; Chien, S.-H.
J. Comput. Chem. 2003, 24, 732-740.]
that the Handy grid is exact if G(q) =
m R^{3}
q^{3m−1}
(1 − q)^{3m+1}
g[r(q)] equals
L(q), a continuous function that is linear
between the roots q_{i} = i/(n + 1) and which
vanishes at q = 0 and q = 1. This implies, for
r_{0} = 0,
g(r) = ^{m}√R m r^{(3m−1)/m} ( ^{m}√R + ^{m}√r )^{2} L ( ^{m}√r ^{m}√R + ^{m}√r) ) | (10.13) |
In general, such g(r) possess a
1/r^{(3m−2)/m} singularity at the origin
and decay as 1/r^{(3m+1)/m} for large
r.[10Gill, P. M. W.;
Chien, S.-H.
J. Comput. Chem. 2003, 24, 732-740.]
In PAMoC, radial integration by the semi-open extended trapezoidal rule, combined with the Handy grid with m = 2, can be selected by the keyword "radrule 73". This is the only option available in all versions of PAMoC before 2013.
Knowles' grid. In 1996, Mura and Knowles
introduced the transformation:[8Mura, M. E.;
Knowles, P. J.
J. Chem. Phys. 1996, 104,
9848-9858.]
r = r_{0} − R ln (1 − q^{k}) ⇒ dr = k R q^{k−1} 1 − q^{k} dq | (10.14) |
which, at both the origin and the endpoint, behaves as the Handy transformation. Therefore, sampling q over (n + 2) equispaced points according to eq. (10.11), leads to a (n + 2)-point semi-open quadrature formula, which de facto reduces to a n-point formula like (10.10), whose evaluation points and weights define the Knowles grid:
r_{i} = r_{0} − R ln (1 − q_{i}^{k}) | |
ν_{i} = k R q_{i}^{k−1} 1 − q_{i}^{k} w_{i} r_{i}^{2} ∀ i = 1, 2, …, n | (10.15) |
q_{i} = i h, w_{1} = w_{2} = … = w_{n − 1} = h, w_{n} = 3 2 h, h = 1 n + 1, k ≥ 2 |
For i = n + 1 (i.e. q_{n+1} = 1) both r_{i} and ν_{i} are singular, so that the last evaluation point is discarded and the weight of the n-th point is multiplied by ^{3}⁄_{2}. For i = 0 (i.e. q_{0} = 0) it is r_{i} = r_{0} and ν_{0} = 0, so that the first evaluation point, q_{0}, doesn't contribute to the integral and therefore has been excluded from the grid. However, if k = 1, q_{0} can be excluded from the grid only if r_{0} = 0, as ν_{0} is proportional to R r_{0}^{2}. So, as in the case of the Handy grid for m = 1, if k = 1 and r_{0} is not null, an (n+1)-point semi-open extended trapezoid quadrature provides the n evaluation points q_{i} = (i − 1)h ∀ i = 1, 2, …, n, and the associated weights w_{1} = ^{1}⁄_{2} h, w_{2} = w_{3} = … = w_{n − 1} = h, w_{n} = ^{3}⁄_{2} h, with h = ^{1}⁄_{n}.
At the grid mid-point, if n is odd, we have i = (n + 1)/2, q_{(n + 1)/2} = ^{1}⁄_{2}, and r_{(n + 1)/2} - r_{0} R = ln 2^{k} 2^{k} − 1 = 1 σ_{k}. Here, σ_{k} = ln^{−1} 2^{k} 2^{k} − 1 is the standardization factor for the Knowles grid (cf § 2.3.4). In particular, we have σ_{1} = ln^{−1}(2) = 1.442695, σ_{2} = ln^{−1}(^{4}⁄_{3}) = 3.476059, σ_{3} = ln^{−1}(^{8}⁄_{7}) = 7.488876, and σ_{4} = ln^{−1}(^{16}⁄_{15}) = 15.494622, for k = 1, 2, 3, and 4, respectively. Mura and Knowles called the transformation (10.14) the “Log-k” transformation. They found that the Log3 grid was the most efficient in integrating the atomic density and common exchange correlation functionals, with scaling parameter values of R = 5.0 (for H-He, B-Ne, Al-Zr, and Sc-Zn) and R = 7.0 (for Li-Be, Na-Mg, and K-Ca), which are in the range of the value σ_{3} = 7.488876 adopted in PAMoC.
For r_{0} = 0 and R = 1, the Knowles grid is exact if G(q) = k q^{k−1} ln^{2}(1 − q^{k}) 1 − q^{k} g[r(q)] equals L(q), which implies
g(r) = exp(−r) k r^{2} [1 − exp(−r)]^{(k−1)/k} L( ^{k}√1 − exp(−r) ) | (10.16) |
In general, such g(r) possess a
1/r^{(3k−2)/k} singularity at the origin
and decay exponentially.[10Gill, P. M. W.; Chien, S.-H.
J. Comput. Chem. 2003, 24, 732-740.]
In PAMoC, radial integration by the semi-open extended trapezoidal rule, combined with Knowles' Log3 grid, is selected by the keyword "radrule 72".
Gill and Chien transformation of the
radial coordinate.
The (n+1)-point semi-open extended trapezoid rule can be used in
combination with the transformation (9.2), introduced by Gill and Chien,[10Gill, P. M. W.; Chien, S.-H.
J. Comput. Chem. 2003,24, 732-740.] to define
the grid:
r_{i} = r_{0} − R ln q_{i} | |
ν_{i} = − R r_{i}^{2} w_{i} q_{i} ∀ i = 1, 2, …, n | (10.17) |
q_{i} = (n − i + 1) h,
h =
1
n and w_{1} = 1 2 h, w_{2} = w_{3} = … = w_{n − 1} = h, w_{n} = 3 2 h, r_{0} ≠ 0. |
Indeed, this transformation of the radial coordinate has been thought to
be used in conjunction with the Gauss-Gill “log-squared”
quadrature rule (9.1).[10Gill, P. M. W.;
Chien, S.-H.
J. Comput. Chem. 2003,24,
732-740.] However, nothing prevents to apply it to a set of
evenly spaced points.
For i = n + 1 (i.e. q_{n+1} = 0) both r_{i} and ν_{i} are singular, so that the (n+1)-th evaluation point is discarded and the weight of the n-th point is multiplied by ^{3}⁄_{2}. For i = 1 (i.e. q_{1} = 1) it is r_{1} = r_{0} and ν_{1} is proportional to R r_{0}^{2}, so that the first evaluation point contributes to the integral only if r_{0} ≠ 0. Otherwise, if r_{0} = 0, the n evaluation points may be calculates as q_{i} = (n − i + 1)/(n + 1) ∀ i = 1, 2, …, n, and the associated weights as w_{1} = w_{2} = … = w_{n−1} = h, w_{n} = ^{3}⁄_{2} h, with h = 1/(n + 1).
At the center of the interval, it is q = ^{1}⁄_{2}, i = ^{(n+2)}⁄_{2}, and r = r_{0} + R ln 2, so that the standardization factor can be defined as either σ_{PAMoC} = 1/ln 2 = 1.4427 or σ_{Gill-Chien} = 1/q_{(n+2)/2}. The two values are rather close, depending on n, and become closer and closer as n increases. In PAMoC, grid (10.17) can be used for radial integration through the keyword "radrule 71".
Becke's grid.
In Becke's transformation of the radial coordinate,[2Becke, A. D.
J. Chem. Phys. 1988,
88, 2547-2553.] defined by eq. (8.4), the variable q
is defined in the interval [−1, +1], so that r = r_{0} at
q = −1, and is singular at q = 1. So, an (n+1)-point
semi-open extended trapezoid rule can be used to obtain the grid:
r_{i} = r_{0} + R 1 + q_{i} 1 − q_{i} = r_{0} + R i − 1 n − i + 1 | |
ν_{i} = 2 R r_{i}^{2} (1 − q_{i})^{2} w_{i} ∀ i = 1, 2, …, n | (10.18) |
q_{i} =
2 i − n
− 2
n, h =
2
n and w_{1} = 1 2 h, w_{2} = w_{3} = … = w_{n − 1} = h, w_{n} = 3 2 h, r_{0} ≠ 0 |
For i = n + 1 (i.e. q_{n+1} = 1) both r_{i} and ν_{i} are singular, so that the (n+1)-th evaluation point is discarded and the weight of the n-th point is multiplied by ^{3}⁄_{2}. For i = 1 (i.e. q_{1} = − 1) it is r_{1} = r_{0} and ν_{1} is proportional to R r_{0}^{2}, so that the first evaluation point contributes to the integral only if r_{0} ≠ 0. Otherwise, if r_{0} = 0, the n evaluation points may be calculates as q_{i} = (2 i − n − 1)/(n + 1) ∀ i = 1, 2, …, n, and the associated weights as w_{1} = w_{2} = … = w_{n−1} = h, w_{n} = ^{3}⁄_{2} h, with h = 2/(n + 1).
At the center of the interval, it is q = 0, i = (n+2)/2, and r = r_{0} + R, so that the standardization factor is equal one and the grid doesn't need to be standardized. In PAMoC, grid (10.18) can be used for radial integration through the keyword "radrule 75".
Ahlrichs' grid.
In Ahlrichs' transformation,[7Treutler, O.; Ahlrichs, R.
J. Chem. Phys. 1995, 102, 346-354.]
defined by eq. (8.8), the variable q
is defined in the interval [−1, +1], so that r = r_{0} at
q = −1, and is singular at q = 1. So, an (n+2)-point
semi-open extended trapezoid rule can be used to obtain the grid:
r_{i} = r_{0} − R ln 2 (1 + q_{i})^{α} ln 1 − q_{i} 2 | |
ν_{i} = r_{i}^{2} R ln 2 (1 + q_{i})^{α−1} ( 1 + q_{i} 1 − q_{i} − α ln 1 − q_{i} 2) w_{i} ∀ i = 1, 2, …, n | (10.19) |
q_{i} =
2 i − n
− 1
n + 1, h =
2
n + 1 and w_{1} = w_{2} = … = w_{n − 1} = h, w_{n} = 3 2 h |
For i = n + 1 (i.e. q_{n+1} = 1) both r_{i} and ν_{i} are singular, so that the (n+2)-th evaluation point is discarded and the weight of the n-th point is multiplied by ^{3}⁄_{2}. For i = 0 (i.e. q_{0} = − 1) it is r_{i} = r_{0} and ν_{0} = 0, so that the first evaluation point, q_{0}, doesn't contribute to the integral ond therefore has been excluded from the grid. At the center of the interval, it is q = 0, i = (n+1)/2, and r = r_{0} + R, so that the standardization factor is equal one and the grid doesn't need to be standardized. In PAMoC, grid (10.19) can be used for radial integration through the keyword "radrule 76".
Modified Handy's grid for finite
intervals of the radial coordinate r.
Figure 10.4 shows that Handy's grid comprises a number of points which are
very far from the nucleus, and which probably contribute very little to the
integral of typically exponentially (or Gaussian) decaying integrands.
Such points can, of course, be systematically neglected within computer programs
by dynamic pruning of the grid, which is reduced by about 20% of its original
size.[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., 10Gill,
P. M. W.; Chien, S.-H.
J. Comput. Chem. 2003, 24,
732-740.] On this line, PAMoC truncates
the summation in the right hand side of eq. (10.10) at a value
N_{r} < n such that
q_{Nr} ≤
r_{max}, where r_{max} is large enough
(usually 10 bohr) to be a good approximation to infinity.
However. Mura and Knowles pointed out that it is preferable to use more
appropriate grids from the outset[8Mura, M. E.;
Knowles, P. J.
J. Chem. Phys. 1996, 104, 9848-9858.]
and indeed the Knowles grid is relatively compact, as shown in Figure 10.4.
Therefore, a modified form of the Murray, Handy and Laming
transformation[3Murray, C. W.; Handy, N. C.;
Laming, G. J.
Mol.
Phys. 1993, 78, 997-1014.] is made available in
PAMoC to map the interval [r_{0},
r_{max}] to the interval [0, 1]:
r(q) = r_{0} +
(r_{max}
− r_{0}) q^{m}
1 +
(r_{max} − r_{0}
− 2^{m})(1 −
q)^{m}
⇒ dr = m (r_{max} − r_{0}) [1 + (r_{max} − r_{0} − 2^{m})(1 − q)^{m−1}] q^{m−1} [1 + (r_{max} − r_{0} − 2^{m}) (1 − q)^{m}]^{2} dq | (10.20) |
which provides a distribution of the sample points concentrated at the origin r_{0} and with the furthest points at a distance from the nucleus which does not exceed r_{max}. At the origin and for m ≥ 2, it is q = 0, r(0) = r_{0}, r'(0) = 0, and G(0) = r^{2}(0) r'(0) g[r(0)] = 0. At the endpoint it is q = 1, r(1) = r_{max}, r'(1) = m (r_{max} − r_{0}), and G(1) = r^{2}(1) r'(1) g[r(1)] = m r_{max}^{2} (r_{max} − r_{0}) g(r_{max}). At the grid midpoint it is q = ^{1}⁄_{2} and r(^{1}⁄_{2}) = 1, so that the grid is standardized by construction. Sampling q over n equispaced points, i.e. q = ^{(i−1)}⁄_{(n−1)} ∀ i = 1, 2, …, n, leads to a n-point closed quadrature formula, whose leading term is null because r'(0) and then ν_{1} are zero by definition. Thus, the change of variable of eq. (10.20) coupled with the integration by substitution and with the trapezoidal rule, yields:
r_{max}
∫
r_{0}
r^{2} g(r) dr =
1
∫
0
r^{2}(q) r'(q) g[r(q)] dq ≡
1
∫
0
G(q) dq ≈ 1 n − 1 ( 1 2 G_{1} + n−1 ∑ i=2 G_{i} + 1 2 G_{n} ) = n ∑ i=1 w_{i} G_{i} = n ∑ i=1 ν_{i} g(r_{i}) | (10.21) |
The evalution points r_{i} and their associated weights ν_{i} and w_{i} define the modified Handy grid for the finite interval [r_{0}, r_{max}]:
r_{i} = r_{0} + (r_{max} − r_{0}) q_{i}^{m} 1 + (r_{max} − r_{0} − 2^{m})(1 − q_{i})^{m} | |
ν_{i} = m (r_{max} − r_{0}) w_{i} r_{i}^{2} [1 + (r_{max} − r_{0} − 2^{m})(1 − q_{i})^{m−1}] q_{i}^{m−1} [1 + (r_{max} − r_{0} − 2^{m}) (1 − q_{i})^{m}]^{2} ∀ i = 1, 2, ..., n | (10.22) |
q_{i} = (i − 1) h, w_{1} = w_{n} = 1 2 h, w_{2} = w_{3} = … = w_{n − 1} = h, h = 1 n − 1, m ≥ 2 |
If m is equal to one, then r'(q_{0}) is not null and ν_{1} is proportional to (r_{max} − r_{0}) r_{0}^{2} / (r_{max} − r_{0} − 1). Then, if r_{0} is not null, an n-point closed extended trapezoid quadrature provides the n evaluation points q_{i} = (i − 1)h ∀ i = 1, 2, …, n, and the associated weights w_{1} = w_{n} = ^{1}⁄_{2} h, w_{2} = w_{3} = … = w_{n−1} = h, h = 1/(n − 1), m = 1, r_{0} ≠ 0.
This grid is a good approximation for the integral of a radial function which changes rapidly in proximity of the origin and smootly vanishes at large distances. It is available in PAMoC through the keyword "radrule 74".
Linear grid for finite intervals of the radial coordinate r. For sake of completeness, is worth noting that if the radial coordinate r ∈ [r_{0}, r_{max}] is given by evenly spaced points, then the corresponding grid is:
r_{i} = r_{0} + q_{i} | |
ν_{i} = w_{i} r_{i}^{2} ∀ i = 1, 2, …, n | (10.23) |
q_{i} = (i − 1) h, w_{1} = w_{n} = 1 2 h, w_{2} = w_{3} = … = w_{n − 1} = h, and h = r_{max} − r_{0} n − 1 |
Such a grid provides sample points equally spaced from the origin r_{0} to the furthest point at a distance r_{max} from the origin. At the origin it is q_{1} = 0, r_{1} = r_{0}, r'_{1} = 1, and G_{1} = r_{0}^{2} g(r_{0}). At the endpoint it is q_{n} = r_{max} − r_{0}, r_{n} = r_{max}, r_{n}' = 1, and G_{n} = r_{max}^{2} g(r_{max}). At the grid midpoint it is q_{(n+1)/2} = r_{max} − r_{0}, and r_{(n+1)/2} = (r_{0} + r_{max})/2, so that the grid is not standardized.
It goes without saying that this is the worst grid that can be used for the radial quadrature. However PAMoC makes it available to masochists through keyword "radrule 77".
In order to compare the roots and weights from various transformations of the radial coordinate, coupled with the extended trapezoid quadrature scheme, the roots and weights for the 11-point Gill-Chien, Knowles, Handy, Becke, Ahlrichs, and linear grids are listed in Tables 10.2 and 10.3, respectively.
Gill-Chien and Knowles' Log1^{(b)} |
Knowles' Log3^{(b)} | ||||||||
---|---|---|---|---|---|---|---|---|---|
i | σ = 1 | σ = 1.4427 | σ = 1 | σ = 7.4889 | Handy^{(c)} | modified Handy^{(c)} |
Becke^{(c)} | Ahlrichs^{(c)} | Linear^{(d)} |
1 | 0.0870 | 0.1255 | 0.0006 | 0.0043 | 0.0083 | 0.0139 | 0.0909 | 0.0428 | 0.9091 |
2 | 0.1823 | 0.2630 | 0.0046 | 0.0348 | 0.0400 | 0.0659 | 0.2000 | 0.1361 | 1.8182 |
3 | 0.2877 | 0.4150 | 0.0157 | 0.1179 | 0.1111 | 0.1782 | 0.3333 | 0.2738 | 2.7273 |
4 | 0.4055 | 0.5850 | 0.0377 | 0.2826 | 0.2500 | 0.3855 | 0.5000 | 0.4586 | 3.6364 |
5 | 0.5390 | 0.7776 | 0.0751 | 0.5623 | 0.5102 | 0.7418 | 0.7143 | 0.6970 | 4.5455 |
6 | 0.6931 | 1.0000 | 0.1335 | 1.0000 | 1.0000 | 1.3284 | 1.0000 | 1.0000 | 5.4545 |
7 | 0.8755 | 1.2630 | 0.2213 | 1.6570 | 1.9600 | 2.2581 | 1.4000 | 1.3854 | 6.3636 |
8 | 1.0986 | 1.5850 | 0.3514 | 2.6316 | 4.0000 | 3.6571 | 2.0000 | 1.8836 | 7.2727 |
9 | 1.3863 | 2.0000 | 0.5480 | 4.1036 | 9.0000 | 5.5862 | 3.0000 | 2.5508 | 8.1818 |
10 | 1.7918 | 2.5850 | 0.8644 | 6.4735 | 25.0000 | 7.8740 | 5.0000 | 3.5121 | 9.0909 |
11 | 2.4849 | 3.5850 | 1.4708 | 11.0145 | 121.0000 | 10.0000 | 11.0000 | 5.1574 | 10.0000 |
(a) R = 1, r_{0} = 0, and r_{max} = 10. (b) Both standardized (σ = 1.4427) and non-standardized (σ = 1) grids are reported. (c) This grid is standardized by construction. (d) Non-standardized grid. |
Gill-Chien and Knowles' Log1^{(c)} |
Knowles' Log3^{(c)} | ||||||||
---|---|---|---|---|---|---|---|---|---|
i | σ = 1 | σ = 1.4427 | σ = 1 | σ = 7.4889 | Handy^{(d)} | modified Handy^{(d)} |
Becke^{(d)} | Ahlrichs^{(d)} | Linear^{(e)} |
1 | (-4)6.88 | (-3)2.07 | 0 | (-7)2.4 | (-6)1.23 | (-6)5.78 | (-4)8.20 | (-4)1.29 | (-1)7.51 |
2 | (-3)3.32 | (-3)9.98 | (-7)1.5 | (-5)6.31 | (-5)7.68 | (-4)3.37 | (-3)4.80 | (-3)2.14 | (+0)3.01 |
3 | (-3)9.20 | (-2)2.76 | (-6)3.94 | (-3)1.65 | (-3)1.22 | (-3)4.85 | (-2)1.65 | (-2)1.20 | (+0)6.76 |
4 | (-2)2.06 | (-2)6.17 | (-5)4.11 | (-2)1.73 | (-2)1.17 | (-2)4.03 | (-2)4.69 | (-2)4.42 | (+1)1.20 |
5 | (-2)4.15 | (-1)1.25 | (-4)2.64 | (-1)1.11 | (-2)9.11 | (-1)2.51 | (-1)1.25 | (-1)1.30 | (+1)1.88 |
6 | (-2)8.01 | (-1)2.40 | (-3)1.27 | (-1)5.35 | (-1)6.67 | (+0)1.30 | (-1)3.33 | (-1)3.40 | (+1)2.70 |
7 | (-1)1.53 | (-1)4.60 | (-3)5.20 | (+0)2.18 | (+0)5.16 | (+0)5.84 | (-1)9.41 | (-1)8.35 | (+1)3.68 |
8 | (-1)3.02 | (-1)9.06 | (-2)1.95 | (+0)8.19 | (+1)4.80 | (+1)2.23 | (+0)3.00 | (+0)2.02 | (+1)4.81 |
9 | (-1)6.41 | (+0)1.92 | (-2)7.30 | (+1)3.07 | (+2)6.48 | (+1)6.76 | (+1)1.20 | (+0)5.10 | (+1)6.09 |
10 | (+0)1.61 | (+0)4.82 | (-1)3.08 | (+2)1.29 | (+4)1.88 | (+2)1.44 | (+1)7.50 | (+1)1.47 | (+1)7.51 |
11 | (+0)9.26 | (+1)2.78 | (+0)2.97 | (+3)1.25 | (+6)5.80 | (+1)9.09 | (+3)2.18 | (+1)9.40 | (+1)4.55 |
(a) The weights are given in scientific notation, the power of ten being shown in parentheses. (b) R = 1, r_{0} = 0, and r_{max} = 10. (c) Both standardized (σ = 1.4427) and non-standardized (σ = 1) grids are reported. (d) This grid is standardized by construction. (e) Non-standardized grid. |
The 11-point grids of Table 10.2 are illustrated in Figure 10.4. The 11-point Gill-Chien and Knowles' Log1 grids share the same roots (and weights). The only difference between them is that the former gives r_{i} values which increase as q_{i} values decrease, whereas the latter gives r_{i} values which increase with q_{i} values. The effect of standardization is to increase the resolution. This effect is clearly shown in Figure 10.4 for the Gill-Chien and Knowles grids, which, together with the linear, grid are not standardized by construction. Standardization ensures that the curves coincide at the middle root (r_{1} = 1), but the Handy, Knowles and Becke grids extend further out, and the Ahlrichs and Knowles grids extend further in than the others. The Gill-Chien grid has the narrowest range of the seven. In all cases, resolution increases with the number of roots n.
Euler-Maclaurin quadrature.
The extended trapezoid quadrature is closely related to the Euler-Maclaurin
summation formula, which is an important tool in numerical analysis
to estimate sums or integrals of functions. For an integral in the interval
[a, b] and a step length h =
(b − a) / (n − 1), it takes the
form:[16Smith, F. J.
Numer. Math.
1965, 7, 406-411., 17Graham, R. L.; Knuth, D. E.; Patashnik, O.
"Concrete Mathematics",
Second Edition, Addison-Wesley Publishing Company, Inc., 1994, eq. (9.67), p. 469.]
b = q_{n} ∫ a = q_{1} f(q) dq = h ( 1 2 f_{1} + n−1 ∑ i=2 f_{i} + 1 2 f_{n} ) − m ∑ k=1 h^{2k} B_{2k} (2k)! [ f^{ (2k−1)}(b) − f^{ (2k−1)}(a) ] + E_{m}(f) | (10.24) |
where f(q) and its derivatives are continuous in [a, b] and q_{i} = a + h (i − 1). The remainder term E_{m}(f) has the expression:
E_{m}(f) = (a − b) h^{2m+2} B_{2m+2} (2m+2)! f^{ (2m+2)}(ξ) | (10.25) |
where ξ is some number in [a, b]. The remainder is often small and, in particular, it is zero if f^{ (2m+2)}(q) = 0. The numbers B_{j} that appear in eqs (10.24) and (10.25) are Bernoulli numbers (B_{0} = 1, B_{1} = −^{1}⁄_{2}, B_{2} = ^{1}⁄_{6}, B_{4} = −^{1}⁄_{30}, B_{6} = ^{1}⁄_{42}, B_{8} = −^{1}⁄_{30}, ... ; B_{3} = B_{5} = B_{7} = B_{9} = ... = 0). The formula (10.24) is identical with the n-point closed extended trapezoidal rule (10.2) but with an expansion of the error term which is in general an asymptotic series.
If f(q) is the integrand of a radial integral, f(q) and its derivatives vanish at the origin a = q_{0} = 0 and tend to zero at large distances, as b = q_{n} → ∞. Under these conditions, eq. (10.24) is often written as
∞ ∫ 0 f(q) dq = h n ∑ i=1 f_{i} | (10.26) |
which is a simplified form of the trapezoidal rule.
Gill and Chien warned about the opportunity “of selecting specific
numerical examples to illustrate quadrature accuracy […], for one can
always select examples to make a given scheme look as good, or as bad, as one
wishes.”[10Gill, P. M. W.;
Chien, S.-H.
J. Comput. Chem. 2003,24,
732-740.] However, taking into account their advice, and
mimicking their work, it is worth examining a few example tests, “not
because they are conclusive, but because they are interesting.”[10Gill, P. M. W.; Chien, S.-H.
J. Comput. Chem. 2003,24, 732-740.]
The radial functions reported in the
first column of Table 11.2 are used to illustrate grid accuracy.[10Gill, P. M. W.; Chien, S.-H.
J. Comput. Chem. 2003,24, 732-740.]
If the radial function g(r) is a Gaussian with unit exponent,
Table 11.2 reveals that the grids may be ranked as
Handy < modified Handy < Knowles' Log3 < Becke ≈ Linear < Gill-Chien < Ahlrichs |
Adding a second, tighter Gaussian, worsens the linear grid accuracy, yielding the ranking
Linear ≪ Handy < modified Handy < Becke ≈ Knowles' Log3 < Gill-Chien < Ahlrichs |
Adding a third, even tighter Gaussian to g(r), leads to deterioration in all grid accuracies to yield
If the same test is done using an exponential function with unit exponent, the grids are ranked as
Ahlrichs < Gill-Chien < linear < Becke ≈ Knowles' Log3 ≈ modified Handy < Handy |
Adding one and two exponential functions with exponent increasingly higher makes the linear grid accuracy the worst of all and leaves the ranking of the remaining ones unchanged.
Linear ≪ Knowles' Log3 ≈ Handy ≈ modified Handy < Ahlrichs ≈ Becke ≈ Gill-Chien |
When the grids of Tables 10.2 and 10.3 are used to estimate the slowly converging radial integral of the function 1/(1 + r^{4}) the results shown in the final row of Table 11.2 are obtained, with the grid ranking
Gill-Chien < Ahlrichs ≈ modified Handy ≈ linear < Knowles' Log3 < Handy < Becke |
favoring the grids which sample the integrand at large r.
i | g(r) |
---|---|
1 | exp(−r^{2}) |
2 | ∑_{j=0,1} 10^{j} exp(−10^{j} r^{2}) |
3 | ∑_{j=0,1,2} 10^{j} exp(−10^{j} r^{2}) |
4 | exp(−r) |
5 | ∑_{j=0,1} 10^{2j} exp(−10^{j} r) |
6 | ∑_{j=0,1,2} 10^{2j} exp(−10^{j} r) |
7 | 1/(1 + r^{4}) |
g(r) | Gill-Chien and Knowles' Log1 σ = 1.4427 |
Knowles' Log3 σ = 7.4889 |
Handy | modified Handy |
Becke | Ahlrichs | Linear | ||
---|---|---|---|---|---|---|---|---|---|
1 | 4.1 | 3.3 | 2.0 | 2.9 | 3.5 | 5.3 | 3.5 | ||
2 | 4.0 | 3.7 | 2.3 | 2.8 | 3.6 | 5.3 | 0.6 | ||
3 | 3.1 | 2.3 | 2.4 | 2.4 | 3.1 | 3.0 | 0.5 | ||
4 | 1.4 | 2.5 | 2.8 | 2.5 | 2.5 | 1.2 | 2.3 | ||
5 | 1.4 | 2.5 | 2.8 | 2.5 | 2.6 | 1.3 | 1.0 | ||
6 | 1.3 | 2.4 | 2.8 | 2.5 | 2.2 | 1.3 | 1.0 | ||
7 | 0.8 | 1.5 | 2.1 | 1.0 | 2.2 | 1.0 | 1.1 | ||
(a) As suggested by Gill and Chien,[10Gill, P. M. W.; Chien, S.-H. J. Comput. Chem. 2003,24, 732-740.] the quality of a quadrature estimation (“Approx”) of an integral is compared with the exact value (“Exact”) using the equation Accuracy = −log_{10} | Approx Exact − 1 | The resulting values are essentially the number of correct digits in the quadrature result. |
It is clear that, if r is null outside the interval [0,10], the Log1 grid should not be coupled with the Euler-Maclaurin quadrature because it covers less than half of the interval. On the other hand, the Handy grid should be discarded because two points (18%) fall outside the interval. Therefore the modified Handy grid and Knowles' Log3 grid remain as suitable choices for the Euler-Maclaurin radial quadrature.
MultiExp^{(b)} | Knowles' Log3^{(b)} | Handy^{(c)} | ||||||||
---|---|---|---|---|---|---|---|---|---|---|
i | σ = 1 | σ = 1.1374 | σ = 1.4427 | σ = 1 | σ = 13.4747 | σ = 7.4889 | σ = 1 | σ = 1.9854 | modified Handy^{(d)} |
Linear^{(d)} |
1 | 4.7 | 4.9 | 5.3 | 2.8 | 1.6 | 1.7 | 0.9 | 1.7 | 1.6 | 3.0 |
2 | 4.8 | 4.8 | 4.9 | 2.6 | 1.8 | 2.3 | 0.9 | 2.3 | 1.9 | 1.3 |
3 | 3.3 | 3.1 | 2.6 | 2.5 | 1.7 | 2.1 | 0.9 | 1.8 | 3.0 | 1.6 |
4 | 8.1 | 2.2 | 2.5 | 0.3 | 6.6 | 6.0 | 1.5 | 1.9 | 2.4 | 2.5 |
5 | 7.9 | 2.3 | 2.5 | 0.4 | 2.6 | 3.1 | 1.5 | 2.0 | 2.4 | 2.1 |
6 | 2.5 | 3.0 | 2.5 | 0.4 | 5.0 | 2.9 | 1.5 | 2.1 | 2.4 | 2.2 |
7 | 0.9 | 1.0 | 1.1 | 0.5 | 2.1 | 1.4 | 2.3 | 1.7 | 1.0 | 1.1 |
(a) See note (a) in Table 11.2 for definition of accuracy. (b-d) See the corresponding notes in Table 9.2. |
Gauss-Legendre Quadrature^{(b)} | 2nd kind Chebyshev Quadrature^{(b)} | |||||
---|---|---|---|---|---|---|
i | Becke | Ahlrichs | Linear | Becke | Ahlrichs | Linear |
1 | 2.2 | 3.4 | 2.6 | 2.3 | 3.7 | 3.6 |
2 | 2.3 | 3.5 | 1.3 | 2.4 | 3.2 | 1.2 |
3 | 2.3 | 4.0 | 1.3 | 2.5 | 2.4 | 1.4 |
4 | 2.8 | 3.9 | 2.6 | 2.5 | 3.5 | 2.5 |
5 | 2.9 | 3.9 | 2.0 | 2.5 | 3.6 | 2.4 |
6 | 3.5 | 2.9 | 1.7 | 2.5 | 2.8 | 2.3 |
7 | 3.7 | 1.2 | 1.0 | 2.6 | 1.1 | 1.0 |
(a) See note (a) in Table 11.2 for definition of accuracy. (b-d) See the corresponding notes in Table 9.2. |
None of the quadrature schemes outlined above are perfect, and each can be
criticized on both aesthetic and practical grounds.[11El-Sherbiny, A.; Poirier, R.A.
J. Comp. Chem. 2004, 25, 1378-1384.]
One can say for sure that there is no absolute favorite radial quadrature
scheme, which dominates others, so that is worth that programs for atomic
quadratures have several for users to choose from.
The quadrature methods examined so far use grid points and weights that
are known beforehand. Other schemes have been designed, in which points and
weights are chosen dynamically as the integrand is explored.[12(a) Pérez-Jordá, J. M.;
San-Fabián, E.; Moscardó, F.
Comp. Phys. Commun. 1992, 70, 271-284.
(b) Pérez-Jordá, J. M.; San-Fabián, E.
Comp. Phys. Commun. 1993, 77, 46-56.
(c) Pérez-Jordá, J. M.; Becke A. D.; San-Fabián,
E.
J. Chem. Phys. 1994, 100, 6520-6534.;
13(a) Krack, M.; Köster,
A. M.
J. Chem. Phys. 1998, 108, 3226-3234.
(b) Köster, A. M.; Flores-Moreno, R.; Reveles, J. U.
J. Chem. Phys. 2004, 121, 681-690.;
14 Gräfenstein, J.; Cremer, D.
J. Chem. Phys. 2007 127, 164113-7.;
15Kakhiani, K.; Tsereteli, K.; Tsereteli, P.
Computer Physics Communications 2009, 180, 256-268.]
the octahedron group with inversion"
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) "Quadratures on a sphere"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) "Spherical quadrature formulas exact to orders 25-29"Lebedev, V.I. Siberian. Math. J. 1977, 18(1), 99-107. Sibirskii Matematicheskii Zhurnal 1977, 18(1), 132-142.
(d) "Quadrature formulas of orders 41, 47, and 53 for the sphere"Lebedev, V. I.; Skorokhodov, A. L. Russian Acad. Sci. Dokl. Math. 1992, 45, 587-592.
(e) "A quadrature formula for the sphere of 59th algebraic order of accuracy"Lebedev, V. I. Russian Acad. Sci. Dokl. Math. 1995, 50, 283-286.
(f) "A quadrature formula for the sphere of the 131-st algebraic order of accuracy"Lebedev, V.I.; Laikov, D.N. Dokl. Math. 1999, 59, 477-481.
Pérez-Jordá, J. M.; San-Fabián, E.; Moscardó, F. Comp. Phys. Commun. 1992, 70, 271-284.
(b) "A simple, efficient and more reliable scheme for automatic numerical integration"Pérez-Jordá, J. M.; San-Fabián, E. Comp. Phys. Commun. 1993, 77, 46-56.
(c) "Automatic numerical integration techniques for polyatomic molecules"Pérez-Jordá, J. M.; Becke A. D.; San-Fabián, E. J. Chem. Phys. 1994, 100, 6520-6534.
Krack, M.; Köster, A. M. J. Chem. Phys. 1998, 108, 3226-3234.
(b) "Efficient and reliable numerical integration of exchange-correlation energies and potentials"Köster, A. M.; Flores-Moreno, R.; Reveles, J. U. J. Chem. Phys. 2004, 121, 681-690.
J =
∂ f
∂ x
=
(
∂ f
∂ x_{1}
…
∂ f
∂ x_{n}
) =
| (N1) |
As a mnemonic, 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.
If m = n, the Jacobian matrix is a square matrix, and therefore admits a determinant, the Jacobian determinant of f.
If m = 1, f is a scalar-valued function of multiple variables, f(x_{1},…,x_{n}), and the Jacobian matrix is reduced to a row vector of partial derivatives of f, i.e. the transpose of the gradient of f, ∇^{†}f (confusingly, often called "the Jacobian" as well).
The transpose of the Jacobian matrix can be viewed as the gradient of the transpose of the vector-valued function f, i.e. J^{†} = ∇f^{†}. Then, the Jacobian generalizes the gradient of a scalar-valued function of multiple variables, which itself generalizes the derivative of a scalar-valued function of a single variable. In other words, the Jacobian for a scalar-valued multivariate function is the gradient and that of a scalar-valued function of single variable is simply its derivative.
[Main source of this note: Wikipedia contributors. Jacobian matrix and determinant. Wikipedia, The Free Encyclopedia. Available at: https://en.wikipedia.org/w/index.php?title=Jacobian_matrix_and_determinant&oldid=834303707. Accessed April 13, 2018].