In the introduction it has been shown that
a definite integral of a function f(q), which is continuous over
a closed interval [a,b] of the real axis q, can be partitioned
into a sum of definite integrals over n sub-intervals
[qi,qi+1]
b∫af(q) dq =
n∑i=1qi+1∫qif(q) dq
(1)
and that each of these integrals can be approximated by the product of
the value of the function at a given point ti∈ [qi,qi+1]
and the length of the interval wi =
qi+1 − qi
qi+1∫qif(q) dq ≈
fi wi, left-hand
rule; fi+1wi,
right-hand rule; f(qi +
qi+12)wi,
midpoint rule; fi +
fi+12wi,
trapezoidal rule.
(3)
The sum of the (signed) areas of all the rectangles or trapezoids
in the interval [a,b] gives the corresponding composite (or
extended or chained or repeated) rules
When the nodes of the partition are equally spaced, as is often the case,
i.e. if wi = w =
b −
an ∀ i,
the composite rules in (4) can be simplified for calculation efficiency
If the function f(q) is given explicitly instead of simply
being tabulated at the values qi, the best numerical method
of integration is called Gaussian quadrature. By picking the intervals
at which to sample the function, this procedure produces more accurate
approximations (but is significantly more complicated to implement).
To integrate a function f(q) over some
interval [a,b], divide it into n equal parts of width h
≡ (b-a)/n such that q1 = a,
qn+1 = b, and
fi = f(qi) for i = 1, 2,
…, n, n+1. Then find polynomials which approximate the
tabulated function, and integrate them to approximate the area under the curve.
To find the fitting polynomials, use Lagrange interpolating polynomials.
The resulting formulas are called Newton-Cotes formulas.
Newton-Cotes formulas may be "closed" if all
points in the interval [a = q1,
b = qn+1] are used, "open" if only the points
[q2, qn] are used, or a variation of these
two. Let's first consider the derivation of closed formulas.
Let Ln(q) be the interpolation polynomial for
f(q) constructed from its values at the n+1 points
qi (cf Eq. (2.9) of the Introduction):
Ln(q) =
n+1∑i=1ln,i(q) f(qi)
(1.1)
Here ln,i(q) is the lagrangian basis polynomial of
the i-th node.
The integral over [a,b] of f(q) is approximately replaced
by the integral of Ln(q):
The last term in the equation above is the well known
trapezoidal rule, 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.
Quadrature formulas for n > 2 are obtained straightforwardly
in a similar way.
(♣) n = 2: 3-point closed formula
or Simpson's rule. The Lagrange interpolating polynomial is
We use integration by substitution to estimate the integral
∫b
aL2(q) dq.
Since points q1 = a < q2 <
q3 = b are equally spaced by
h = q3 − q12, we have the following
relations q2 =
q1 + q32,
q1 = q2 −
h, and q3 =
q2 + h, which lead to the expression
L2(u) =
12[u (u − 1) f1 +
2 (1 − u2) f2 +
u (u + 1) f3]
Since u = −1 for q =
q1 = q2 − h, u
= 1 for q = q3 = q2 +
h, and h du = dq, we simply
have to calculate the integral
b∫aL2(q) dq =
h1∫−1L2(u) du =
h2(f11∫−1u (u − 1) du + 2 f21∫−1
(1 − u2) du + f31∫−1u (u + 1) du) =
h2(f1[u33 −
u22]1
−1 +
2 f2[u −
u33]1
−1 +
f3[u33 +
u22]1
−1) =
h3
(f1 + 4 f2 + f3)
The last equality in the equation above is the well known Simpson's
rule, which is exact if f is a polynomial up to quadratic degree.
Because of the 1⁄3 factor Simpson's rule is also referred to as
Simpson's 1⁄3 rule.
(♥) n = 3: 4-point closed
formula or Simpson's second rule. The Lagrange interpolating polynomial
is
which is exact if f is a polynomial up to fourth degree.
Find the
(n+1)-point closed Newton-Cotes formulas for n = 1, 2, 3, 4,
and 5, using equations (1.2) and (1.9).
The values of the various factors that determine the weights
of the Newton-Cotes integration formulas according to Eq. (1.9) are reported
in the following Table for n = 1,2, ..., 5. The fifth column shows the
integral value over the interval [0,n] of the polynomial function of
degree n reported in the fourth column. The product of the factors
reported in the third and fifth columns gives the value of the weight
wi in the approximate integration formula (1.2).
Table I. Calculation of the weights of the (n+1)-point
Newton-Cotes integration formulas by means of eq. (1.9).
n
i
(−1)n−i+1(i − 1)!
(n − i + 1)!
n+1∏j=1; j≠i
(s − j + 1)
hn∫0n+1∏j=1; j≠i
(s − j + 1) ds
wi ≡ αi h
1
1
−1
s−1
−12h
12h
2
1
s
12h
12h
2
1
12
(s−1)(s−2)
23h
13h
2
−1
s(s−2)
−43h
43h
3
12
s(s−1)
23h
13h
3
1
−16
(s−1)(s−2)(s−3)
−94h
38h
2
12
s(s−2)(s−3)
94h
98h
3
−12
s(s−1)(s−3)
−94h
98h
4
16
s(s−1)(s−2)
94h
38h
4
1
124
(s−1)(s−2)(s−3)
(s−4)
11215h
1445h
2
−16
s(s−2)(s−3)(s−4)
−12815h
6445h
3
14
s(s−1)(s−3)(s−4)
3215h
815h
4
−16
s(s−1)(s−2)(s−4)
−12815h
6445h
5
124
s(s−1)(s−2)(s−3)
11215h
1445h
5
1
−1120
(s−1)(s−2)(s−3)
(s−4)(s−5)
−47512h
95288h
2
124
s(s−2)(s−3)(s−4)
(s−5)
1254h
12596h
3
−112
s(s−1)(s−3)(s−4)
(s−5)
−12512h
125144h
4
112
s(s−1)(s−2)(s−4)
(s−5)
12512h
125144h
5
−124
s(s−1)(s−2)(s−3)
(s−5)
−1254h
12596h
6
1120
s(s−1)(s−2)(s−3)
(s−4)
47512h
95288h
In the last column of Table I the weights αi
= wi / h have been introduced. These weights are also known
as Cotes numbers. They depend on i and n, but not on nodes
qi and have two interesting properties, i.e. they are such
that αi = αn−i+2
and ∑n+1 i=1αi =
n, due to the property of the Lagrange coefficients
li(q) to form a partition of the unit.
Table 1, in the chapter text, lists some of the Newton–Cotes
formulas of the closed type (from 2- to 6-point).
Higher order rules, from 7-point to 11-point, are reported in
reference [L7Weisstein, Eric W.
"Newton-Cotes Formulas." From MathWorld − A Wolfram Web
Resource.]. All formulas can be concisely represented by the
triangle of coefficients shown in the following table
Table II. Cotes numbers which determine the coefficients of
the (n+1)-point Newton-Cotes integration formulas (see Table 1 in the
chapter text).
↓n \ i→
1
2
3
4
5
6
7
8
9
10
11
1
12
12
2
13
43
13
3
38
98
98
38
4
1445
6445
815
6445
1445
5
95288
12596
125144
125144
12596
95288
6
41140
216140
27140
272140
27140
216140
41140
7
525717280
2503917280
962117280
2092317280
2092317280
962117280
2503917280
525717280
8
395614175
2325214175
−371214175
4198414175
−1816014175
4198414175
−371214175
2325214175
395614175
9
2571389600
14166989600
972089600
17409689600
5200289600
5200289600
17409689600
972089600
14166989600
2571389600
10
80335299376
531500299376
−242625299376
1362000299376
−1302750299376
2136840299376
−1302750299376
1362000299376
−242625299376
531500299376
80335299376
Note that the sum of the elements of a row gives n.
Let's now turn to (n+1)-point open
Newton-Cotes formulas. In this case function evaluation at the endpoints
of the interval is excluded from the quadrature rule.
This implies that there are only (n−1) points that can be
interpolated and that the degree of the Lagrange interpolating polynomial is
at most (n−2).
The step size is still h = b − an because n is the
number of intervals in a sequence of (n+1) points.
Eqs (1.2) and (1.9) still hold, but i and j range between 2
and n. That is
b∫af(q) dq ≈
n∑i=2wif(qi)
(1.10)
and
wi =
(−1)n−i+1hn!(ni − 1)n∫0(n+1∏j=1; j≠i
(s − j + 1)
)
ds
∀ i = 2, 3, …, n
(1.11)
An open formula is often referred to by the number of points where the
function is known, i.e. m = n − 1. Then the degree of the Lagrange
interpolating polynomial is at most m − 1 and the step size becomes
h = b −
am + 1. In addition,
the numbering of points starts from 1, and each point is indicated by
qi = a + i h
∀ i = 1, 2, …, m.
Find
the (n + 1)-point open Newton-Cotes formulas for n = 2, 3, 4, and 5.
(♠) n = 2: 3-point open
formula or midpoint rule. Given a sequence of three equally spaced points
q1 < q2 < q3,
find an approximate value of the integral of f(q) from
q1 to q3 when the value of
f(q) is known only at point q2.
The Lagrange interpolating polynomial of order n = 2,
L2(q), contains three terms. We retain only that
corresponding to the known point [q2,
f(q2)]:
Since the points are equally spaced by
h = q3 − q12 =
b − a2, we have the following
relations q2 =
q1 + q32 =
a + b2,
q1 = q2 −
h, and
q3 =
q2 + h, which lead to the expression
l2,2(q) =
−h−2 (q − a)
(q − b). Since Eq. (1.6) allows us to replace
q with a + sh, we have
l2,2(q) =
s (1 − s) = l2,2(s).
The following integral provides the final solution
Ultimately, the Newton-Cotes formula for n = 2
(3-point open formula) is the so-called rectangle or midpoint rule.
This formula approximates the integral (and therefore the area subtended by
the function) as a rectangle of base (b − a) and height
f(a + b2),
where a and b are the extremes of integration and
a
+ b2
is the midpoint of the interval.
(♣) n = 3: 4-point open formula
or trapezoidal rule. Given a sequence of four equally spaced points
q1 < q2 < q3 <
q4, find an approximate value of the integral of
f(q) from q1 to q4 when the
value of f(q) is known only at points q2 and
q3.
There are two points to interpolate with a polynomial of first degree.
Let the step size be h =
q4 −
q13 =
b − a3 and q1 =
a, q2 = a + h, q3 = a +
2h, q4 = a + 3h = b.
Use the Lagrange method to interpolate the two points (q2,
f2) and (q3, f3), where
the notation fi is a shorthand for
f(qi). According to eqs (1.1) and (1.3), the
interpolating polynomial is
The following material will be moved to other
sections (web pages)
1.1. - Newton-Cotes quadrature formulas.
The interpolatory procedure outlined above leads to an extremely useful and
straightforward family of numerical integration techniques, which are known
as Newton-Cotes formulas and apply to n equispaced quadrature points
or nodes in the range [a, b], such that
qk = 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 = q1, b = qn] are used, "open"
if only the points [q2, qn−1]
are used, or a variation of these two.
A few closed Newton-Cotes quadrature formulas for small values of n are
reported in Table 1.1.
Using polynomial interpolation with polynomials of high degree over a set of
equispaced points may result in an oscillating pattern that magnifies near the
ends of the interpolation points and is known as Runge's
phenomenon.[1Runge, C. Zeitschrift für Mathematik und Physik1901, 46,
224-243.]
This implies that going to higher degrees does not always improve accuracy.
Runge's phenomenon can be avoided by using a composite or
extended rule.
Closed "extended" rules use multiple copies of lower
order closed rules to build up higher order rules. For n tabulated
points, using the 2-point trapezoidal rule (n − 1) times and
adding the results gives
where h ≡ (b − a) / (n − 1) and
qi = a + h (i − 1).
It can be shown that the error associated with the n-point
closed extended trapezoidal rule can be expressed in terms of derivatives
at the endpoints, according to the Euler-Maclaurin formula
The numbers Bj that appear in eq. (1.11) are Bernoulli
numbers
(B0 = 1,
B1 = −1⁄2,
B2 = 1⁄6,
B4 = −1⁄30,
B6 = 1⁄42,
B8 = −1⁄30, ... ;
B3 = B5 = B7 = B9 = ... = 0),
and the remainder Em(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 p0(q),
p1(q), …,
pn−1(q), defined in some
interval [a, b] and such that
pi(q) is a polynomial of degree i;
the inner product of any two polynomials
pi(q) and pj(q)
with respect to a positive function ω(q), called a
weighting function, is given by
In addition, if
<pi|ω|pi> = 1,
the polynomials are orthonormal.
In particular, eq. (1.12) states that
b∫aω(q) pi(q) p0(q) dq = 0
Since p0(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) pi(q) dq = 0
∀ i > 0.
(1.13)
Consider the set of roots q1, q2, ...,
qn of the polynomial pn(q).
It is trivial to show that
n∑i=1vi pn(qi) = 0
for any set of vi's, because
pn(qi) =
0 for all qi.
However it is possible to show that there exists a unique set of
vi's, for which all polynomials of order 1 ≤
k ≤ n simultaneously obey the relation
b∫aω(q) pk(q) dq =
n∑i=1vi pk(qi) = 0.
(1.14)
This is done by solving the simultaneous set of equations
p0(q1)
p0(q2)
…
p0(qn)
p1(q1)
p1(q2)
…
p1(qn)
⋮
⋮
⋮
pn−1(q1)
pn−1(q2)
…
pn−1(qn)
v1
v2
⋮
vn
=
∫abω(q) p0(q) dq
0
⋮
0
(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 p0(q),
p1(q), …,
pn−1(q)
f(q) = ω(q)
n∑k=0ck pk(q)
(1.16)
where
ck =
b∫af(q) pk(q) dq/b∫aω(q) pk2(q) dq
(1.17)
is an expansion coefficient. Integrating both sides of eq. (1.16) over the
interval [a, b] yields:
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 wi and
abscissas qi 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=1vi g(qi)
(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.
Table 1.2. Gauss quadrature formulas for different
choices of a, b, and ω(q).
The qi's are the roots of the polynomial, and the
wi's are the corresponding weights.
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 Qn(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) = Qn(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 pn(q) of degree n,
is (2n − 1).
Theorem. Let {pk(q)} ∀ k
∈ [0, n] be a set of orthogonal polynomials on [a, b] with respect
to the inner product (1.12). Let qj, j = 1, 2,
…, n, be zeros of polynomial pn(q).
Then the quadrature rule given by eq. (1.1) has degree of exactness
(2n − 1).
Proof. Given two univariate polynomials f(q) and
pn(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) = pn(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(pn).
Multiplying both sides of eq. (1.22) by the weighting function ω(q) and integrating, yields:
The green integral in eq. (1.23) is zero because
the polynomial pn(q) is orthogonal to all polynomials
of degree ≤ (n − 1), according to eq. (1.12).
Since the rule (1.18) is interpolatory, it is exact for r(q).
Then
Here we have used the fact that qi's are the zeros of
pn(q), i.e. pn(qi) = 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) = Qn(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|ω|pn> ≠ 0. So, on one hand
as qi's are the zeros of pn(q),
i.e. pn(qi) = 0.
Hence, Qn(f) = Qn(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 ≈
ny∑j=1wy,j f(x,yj)
(1.25)
Next we approximate the integral by the product quadrature formula
∫∫f(x,y) dx dy ≈
nx∑i=1ny∑j=1wx,i wy,j f(xi,yj)
(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 ≈
nx∑i=1ny∑j=1nz∑k=1wx,i wy,j wz,k f(xi,yj,zk)
(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) =
x
y
z
=
f1(u1,u2,u3)
f2(u1,u2,u3)
f3(u1,u2,u3)
= f(u1,u2,u3)
(1.28)
where x, y, z are Cartesian coordinates, and
u1, u2, u3 are some
other coordinates. The functions f1, f2 and
f3 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 =
(u1,u2,u3), 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) =
(df1,df2,df3) =
df, which can be written more explicitly and with matrix
notation as
dr =
∂f1/∂u1
du1 +
∂f1/∂u2
du2 +
∂f1/∂u3
du3
∂f2/∂u1
du1 +
∂f2/∂u2
du2 +
∂f2/∂u3
du3
∂f3/∂u1
du1 +
∂f3/∂u2
du2 +
∂f3/∂u3
du3
=
∂ f/∂ u1 du1 +
∂ f/∂ u2 du2 +
∂ f/∂ u3 du3=
∂ f/∂ u du
(1.29)
where ∂ f/∂ uj
is a vector, whose components are the partial derivatives of x =
f1,
y =
f2, and z =
f3 with respect to uj,
∀ j = 1, 2, 3.
Left-multiplying the the vector dr by its transpose
yields the following equalities
dr† dr =
dr ⋅ dr = dr2 =
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 Jf and
∂(f1,f2,f3)/∂(u1,u2,u3),
is usually defined and arranged as follows:[N1Jacobian matrix]
J =
∂ f/∂ u =
(∂ f/∂ u1∂ f/∂ u2∂ f/∂ u3) =
∂f1/∂u1
∂f1/∂u2
∂f1/∂u3
∂f2/∂u1
∂f2/∂u2
∂f2/∂u3
∂f3/∂u1
∂f3/∂u2
∂f3/∂u3
(1.31)
or, component-wise:
Jij =
∂ fi/∂ uj.
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
gii =
|∂ f/∂ ui|2
and gij =
∂ f/∂ ui ⋅
∂ f/∂ uj
for i,j = 1, 2, 3.
Eq (1.30) gives the so-called first fundamental form or line
element
What is important for calculating integrals in 3D space, is the volume
element dV = dx dy dz, which can be
written as the absolute value of the mixed product of the three tangent vectors
in eq (1.29):
dV =
|∂ f/∂ u1 ⋅
(∂ f/∂ u2 ×
∂ f/∂ u3)|
du1 du2 du3= det( J )
du1 du2 du3=
√det( g )
du1 du2 du3
(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/2u +
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∫−1b − 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∫−1b − 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=1wiΦ(ui)
(1.41)
where the ui 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 − u2,
and the wi 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=1viφ(ri)
(1.42)
where the transformed weights, vi, are given by
vi =
b − a/2wi
(1.43)
in the case of Gauss-Legendre quadrature, and
vi =
(b −
a)2/4[(r − a)(b −
r)]1/2wi
(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=1wiΦ(ui)=
n∑i=1viφ(ri)
(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) = ln2u
is used, and yields the transformed weights:
vi = (b − a)
(ln
ri
− a/b − a)−2wi
(1.48)
Integrals over a semi-infinite interval [r0, ∞]
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 = r0 + Su ≡
f(u) ⇓ f'(u) = S,
f−1(r) =
r −
r0/S ≡ u,
f−1(r0) = 0,
f−1(∞) = ∞
(1.49)
which leads to integration rule:
∞∫r0φ(r) dr =
∞∫0Sφ(f(u)) du=
∞∫0
e−u e+uSφ(f(u)) du=
∞∫0
e−uΦ(u) du ≈
n∑i=1wiΦ(ui) =
n∑i=1vif(ri)
(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
vi = S
exp(ri
− r0/S)wi
(1.51)
2. - Volume integral
Evaluation of atomic and molecular properties requires calculation of a volume integral, which is an
integral over a three dimensional domain (i.e. a region V in ℝ3) of a function
f(r) = f(x,y,z) and is usually written as a triple integral:
I =
∫Vf(r) d3r =
∫Vf(r) dV =
∭Vf(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:
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,θ,φ) =
∂x/∂r
∂x/∂θ
∂x/∂φ
∂y/∂r
∂y/∂θ
∂y/∂φ
∂z/∂r
∂z/∂θ
∂z/∂φ
=
cosθ sinφ
−r sinθ sinφ
r cosθ cosφ
sinθ sinφ
r cosθ sinφ
r sinθ sinφ
cosθ
−r sinθ
0
(2.6)
and the volume element is:
dV = r2 sinθ dr dθ dφ
= dr dS
= r2 dr dΩ
(2.7)
where we have introduced the surface element
dS = r2 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/r2
= 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 Ω:
It can be calculated by the one-dimensional radial integration:
I =
+∞∫0r2g(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π∫0f(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
ci = ri, θi, φi, Ωi:
∫Vf(r) dV ≈
Nr∑i=1wr,iNθ∑j=1wθ,jNφ∑k=1wφ,k f(ri,θj,φk) ≈
Nr∑i=1wr,iNΩ∑j=1wΩ,j f(ri,Ωj)
(2.13)
2.1. - Partitioning into atomic subintegrals.
In a molecule, the density
ρ(r) has cusps at the
positions of the nuclei, so that a simple cartesian grid for the entire
molecule cannot account properly for the accumulation of density at that
positions. Then, following design principles put forward by Becke in 1988,
the molecular integration is broken up into separate, but overlapping atomic
contributions.[2Becke, A. D. J. Chem. Phys.1988, 88, 2547-2553.]
These single-center atomic subintegrals are then computed on
atomic grids. Finally, the results are summed together to give the molecular
integral (2.1).
Different partitioning schemes are available in PAMoC for experimental
and/or theoretical densities:
The great simplicity of the three-terms product
formula (Nr, 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.
2.3. -
Radial integral and transformations of the radial coordinate.
The general methodology to calculate the integral (2.11) is to map
the interval [0, ∞] to a finite interval [a, b] of the
variable q by a suitable variable transformation
r = r(q) ⇒
dr = r'(q) dq, and then to
use a quadrature formula defined on this interval.
The integral (2.11) takes then the form
where we have set G(q) =
ω−1(q)
r2(q) r'(q)
g[r(q)] and
vi =
wiω−1(qi)
r2(qi) r'(qi)
. 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.
The transformed coordinate qi 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 qi 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 [r0, ∞] or the finite
interval [r0, rmax], 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.
Table 2.1 − Mappings of the radial coordinate r, defined over the
interval [r0, rmax], onto a new coordinate q, whose
values are chosen as either equispaced points or roots of polynomials orthogonal in the
interval [a, b]. The option values of the associated
quadrature rules are also given.
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(qk) −
r0 = 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] −
r0
(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(qcenter) −
r0
(2.16)
As a consequence, the standardizing factor, σ, depends
only on the type of transformation, and not on how the points
qi are chosen nor on their number n. In addition,
PAMoC defines R as the product of σ and the
parameter Rcov, which eventually can be used to adjust the
grid to fit different atomic shapes:
R = σ Rcov
(2.17)
Typically, Rcov 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 ∈ [r0, ∞] into the finite interval [0, 1],
whose midpoint is qcenter =
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:
Becke and Ahlrichs transformations map the radial coordinate
into the finite interval [−1, +1], whose midpoint is
qcenter = 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
qcenter cannot be defined. Therefore its standardizing factor
is calculated by eq. (2.15) and depends on the value of the middle root and
n:
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.
Table 2.6 − Partitioning parameters used in the pruned
Lebedev (Nr,NΩ) and spherical
product (Nr,Nθ,Nφ) grids.
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
(Nr,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
(Nr,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 Nr-point
Euler-Maclaurin radial grid and the
NΩ-point Lebedev angular grid,
viz. EML(Nr,NΩ), or the
Nθ×Nφ-point
spherical product grid, viz.
EMSP(Nr,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 α1R,
α2R, α3R, and
α4R 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.
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 Nr
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
Ntot 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 xy
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 66, 183, 261,
381, 741, 1101, 1466,
861, 501, 381, 181
indicates that a six-point Lebedev angular grid is used at each of
the six innermost radial points, an 18-point angular grid at each of
the next three radial points, followed by a single 26-point grid, a
single 38-point grid, and so forth."
3. Lebedev angular quadrature
This section briefly discusses the evaluation of the angular integral (2.12)
using Lebedev quadrature:
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∫−1f(q) dq ≈
n∑i=1wif(qi)
(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 qi
and weights wi for i = 1, ..., n.
The evaluation points qi for quadrature order n are just the roots of the Legendre
polynomials Pn(q), which occur symmetrically about 0. In other words, with the n-th polynomial
normalized to give Pn(1) = 1, the i-th node, qi, is the i-th root of
Pn, with associated weight wi.
Table 3.1 − First six Legendre
polynomials with associated roots and weights.
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/2q +
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:
where n is the number of quadrature points used, the
qi are the roots of the Legendre polynomial
Pn(q) with associated weights wi,
and the νi =
b − a/2wi 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
Radial integral. The
Gauss-Legendre quadrature rule can also be used by PAMoC to estimate the
radial integral (2.11), in conjunction with a suitable variable transformation
to map the semi-infinite interval spanned by the radial coordinate r
into the finite interval [−1, +1] of the variable q. This is
accomplished by two transformations of the radial coordinate, due to Becke[2Becke, A. D. J. Chem. Phys.1988, 88, 2547-2553.]
and Ahlrichs,[7Treutler, O.; Ahlrichs, R. J. Chem. Phys.1995, 102, 346-354.]
which are described in section 7 with more details.
In addition, the linear transformation (3.2) is also available.
The appropriate weights for the Gauss-Legendre quadrature are given by
eqs. (3.6), (3.7), and (3.8), respectively,
where qi and wi are roots and weights of
the Gauss-Legendre quadrature rule, and ri 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, qi, combined with
Becke, Ahlrichs, and linear grids, ri ≡
r(qi).
Table 3.3 − Weights(a) of the
n-point Gauss-Legendre quadrature rule,
wi, 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 r0 = 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 r0 = 0.
(c)a = 0 and b = 10.
4. Gauss-Laguerre quadrature
Gauss-Laguerre quadrature provides an approximate solution to the definite integral of a function e−qf(q)
over the interval [0, +∞] as a weighted sum of function values at specified points within the domain of
integration:
+∞∫0
e−qf(q) dq ≈
n∑i=1wi f(qi)
(4.1)
It fits all polynomials of degree (2n − 1).
The evaluation points qi for quadrature order n are
just the roots of the Laguerre polynomials Ln(q),
with associated weight wi.
Table 4.1 − First six Laguerre
polynomials with associated roots and weights.
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,
Figure 4.1: Dependence of the value of the highest radial coordinate,
rn, on the number n of roots of the Laguerre
polynomial.
As the number of quadrature points n increases, the value of
rn increases accordingly (i.e.rn ≈ −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, Rcov,
times a parameter σ, the value of which is chosen so that the
the middle root of the quadrature is equal to Rcov. The
resulting values, for Rcov =
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
rn increases with n, but it seems to reach an
asymptotic value rasymp of about 5.7 units for
n > 80. This fact suggests to scale the radial coordinate furtherly
by the factor r∞/rasymp,
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∞/rasymp.
Gill and Chien pointed out that the Laguerre grid (4.4) is exact if
G(q) =
L2n−1(q),
which implies that
g(r) = R−1r−2
exp(−
r −
r0/R)L2n−1(r −
r0/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.]
As a test, the radial integrals of the two functions
ggauss(r) = exp(−r2) and
gexpon(r) = exp(−r), whose exact values
are √π/4 ≈ 0.44311346272637897 and 2,
respectively, were esimated by the Laguerre rule.
With non-standardized roots and weights, gexpon(r) was
integrated exactly by the 2-point formula, whereas
ggauss(r) needed the 49-point grid to reach an accuracy
of 7.0 (that is essentially the number of correct digits in the quadrature
result). In the latter case the summation converged at the 9th evaluation point,
so that the remaining 40 points were calculated unnecessarily.
On the contrary, using standardized roots and weights, the Gaussian function is
integrated by the 11-point quadrature formula with accuracy 6.8 and almost
exactly by the 15-point grid (accuracy 11.1), whereas the radial integral of
the exponential function is roughly approximated even by the 99-point quadrature
formula (accuracy 1.1). The additional correction of the standardizing
parameter by the integer factor 2 (corresponding to the value of the ratio
r∞/rasymp) does not facilitate the
integration of the Gaussian function, since 18 points are needed for an accuracy
of 10.0, but yields an accuracy of 6.7 for the 5-point quadrature of the
exponential function.
5. Generalized Gauss-Laguerre quadrature
GeneralizedGauss-Laguerre quadrature provides an approximate solution to the
definite integral of a function
qα e−qf(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:
+∞∫0qα e−qf(q) dq ≈
n∑i=1wi f(qi)
(5.1)
It fits all polynomials of degree (2n − 1).
The evaluation points qi for quadrature order n are
just the roots of the generalized Laguerre polynomial
Ln(α)(q),
with associated weight wi.
Table 5.1 − First three
generalized Laguerre polynomials with associated roots
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:
where G(q) = Rr2q−α
eqg(ri),
and ri and νi define the generalized
Laguerre grid:
ri =
r0 + Rqiνi =
Rri2qi−α
eqiwi
(5.3)
It is worth noting that, for r0 = 0 and
α = 2, the expression of the quadrature weights simplifies to
νi = R3
eqiwi.
In this case, the generalized Laguerre grid (5.3) is exact if
g(r) = R−3
exp(−r/R)
L2n−1(2)(r/R).
Table 5.2 − Roots(a) of the 11-point
Gauss-Laguerre and Generalized Gauss-Laguerre (α = 2)
quadrature rules, qi and standardized
ri ≡ r(qi).
Table 5.3 − Weights(a,b) of the 11-point
Gauss-Laguerre and Generalized Gauss-Laguerre (α = 2)
quadrature rules, wi and standardized
νi.
Gauss-Laguerre
Generalized Gauss-Laguerre
Gauss-Laguerre
Generalized Gauss-Laguerre
qi and ri (σ = 1)
σ = 1⁄q6
σ = 2⁄q6
qi and ri (σ = 1)
σ = 1⁄q6
σ = 2⁄q6
wi
σ = 1
σ =
1⁄q6
σ =
2⁄q6
wi
σ = 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 r0 = 0.
(a)R = 1 and r0 = 0.
(b) The weights are given in scientific notation, the power of ten being
shown in parentheses.
Figure 5.1: Comparison of Gauss-Laguerre (full symbols) and
Generalized Gauss-Laguerre (α = 2; open symbols) 11-point grids with no
standardization (blue circles, σ = 1),
with Gill-Chien standardization (red squares,
σ = 1/q6),[10Gill,
P. M. W.; Chien, S.-H. J. Comput. Chem.2003, 24,
732-740.] and with corrected standardization
(dark green triangles, σ =
2/q6 and σ = 3/q6).
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=1wi f(qi)
(6.1)
where n is the number of sample points used.
The qi are the roots of the Hermite polynomial
Hn(q) with associated weights wi.
Table 6.1 − First four Hermite
polynomials with associated roots and weights.
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:
rA = r − A and
rB = 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 N1 and N2 are normalization constants.
Boys[23Boys, S. F. Proc. Roy. Soc. A1950, 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 =
ζ1A +
ζ2B/ζ1 + ζ2
with exponent ζ3 =
ζ1 + ζ2:
for (t; k, k1, k2) =
(x; l, l1, l2),
(y; m, m1, m2), and
(z; n, n1, n2).
As the integrand GA(t) doesn't exactly correspond to the integrand of the quadrature formula, we need a change of variable
where qi are the roots of the Hermite polynomial of order
n and wi are the associated weights. The coordinates
ti are defined in terms of the roots qi,
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 FA(q) in eq. (6.9) is a
polynomial of degree (k +
k1 + k2), a number
n = (k + k1 +
k2)/2 is needed.
It is worth noting that in the case of monocentric integrals we have
where (k+k1+k2) =
(l+l1+l2),
(m+m1+m2),
(n+n1+n2) for t = x, y, z, respectively.
For even (k+k1+k2) an analytical solution is available for the integral
whereas for odd (k+k1+k2) the integral is equal to zero.
The final expression for the Mulliken unabridged cartesian multipole moments of atom A is
mA(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:
where G(q) =
1⁄2Rr2 eq2g(r),
and ri and νi define the Hermite grid:
ri =
r0 + Rqiνi =
1⁄2Rri2
eqi2wi
(6.14)
It is worth noting that, for r0 = 0, the expression
of the quadrature weights simplifies to
νi =
1⁄2R3
eqi2wi.
In this case, the Hermite grid (6.14) is exact if
g(r) =
2 R−3 exp(−r2/R)
H2n−1(r/R).
8. Gauss-Chebyshev quadrature of the second kind
Gauss-Chebyshev quadrature of the second kind provides an approximate solution
to the definite integral of a function
√1 − q2f(q)
over the interval [−1, +1] as a weighted sum of function values at
specified points within the domain of integration:
+1∫−1√1 − q2f(q)
dq ≈
n∑i=1wi f(qi)
(8.1)
It fits all polynomials of degree (2n − 1).
The evaluation points qi for quadrature order n are
just the roots of the Chebyshev polynomial of the second kind
Un(q),
with associated weight wi.
Table 8.1 − First nine Chebyshev
polynomials of the second kind with associated roots and weights.
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 −
q2 over the interval [−1, +1], i.e.:
+1∫−1Un(q) Um(q) √1 − q2 dq =
π/2δnm
(8.2)
As a consequence:
n∑i=1wi =
π/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].
r = r0 + R1 + q/1 − q
≡ r(q) ⇒
dr =
2 R/(1 − q)2
dq = 2
r − r0/
1 − q2
dq
(8.4)
The displacement of the origin from 0 to r0 was not
considered in the original paper, but has been added here for completeness.
At q = − 1, r(q) = r0,
and as q → 1, the sample points tend toward infinity.
Under this map, the radial integral becomes:
where G(q) =
2
r2
(r − r0)/(1
− q2)3/2g(r), and ri and νi
define the Becke grid:
ri =
r0 + R1 + qi/1 − qiνi = 2
ri2
(ri − r0)/(1 −
qi2)3/2wi =
2π/n + 1ri2 (ri −
r0)/√1 −
qi2qi = cos (i/n + 1π)
∀ i = 1, 2, ..., n
(8.6)
At q = 0, r(0) − r0 = 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)
= U2n−1(q), which implies
g(q) =
R3/4 r3/2 (r + R)3U2n−1(r − R/r + R)
(8.7)
In general, such g(r) possess a 1/r3/2
singularity at the origin and decay as 1/r9/2 at large
r.
which is normalized to map the center of the q interval onto
r(0) - r = R. At q = −1,
r(q) − r0 = 0,
and as q → 1, the sample points tend toward infinity.
Under this map, the radial integral becomes:
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:
ri =
b − a/2qi +
a + b/2
νi =
b − a/2
√1 − q2ri2wi =
π (b −
a)/2 (n + 1)ri2√1 − q2
∀ i = 1, 2, ..., n
(8.11)
with qi = 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, qi, combined with
Becke, Ahlrichs, and linear grids, ri ≡
r(qi).
Table 8.3 − Weights(a) of the
n-point Gauss-Chebyshev quadrature rule of the second kind,
wi, 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 r0 = 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 r0 = 0.
(c)a = 0 and b = 10.
Figure 8.1: Comparison of Becke (blue
circles), Ahlrichs (red squares), and
linear (dark green triangles) grids for
n = 11, evaluated using the roots of Legendre (open symbols) and 2nd
kind Chebyshev (full symbols) polynomials.
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.
9. Gauss-Gill quadrature
The Gauss-Gill quadrature provides an approximate solution to the definite
integral of a function ln2qf(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
ln2qf(q) dq ≈
n∑i=1wi f(qi)
(9.1)
It fits all polynomials of degree (2n − 1).
The evaluation points qi for quadrature order n are
the roots of the polynomial Qn(q) orthogonal on the
interval [0, 1] with respect to the weight function ln2q,
and wi are the associated weights.
Table 9.1 − First three Gill polynomials with associated
roots and weights
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 = r0 − R ln q
⇒
dr = −
R/q
dq
(9.2)
At q = 1, r(q) = r0, and
as q → 0, the sample points tend toward infinity.
Under this map, the radial integral becomes:
G(q) = R (r0 −
R ln q)2q−1 ln−2qg(r)
= R3r2
(r − r0)−2
e(r−r0)/Rg(r)
and ri and νi define the MultiExp
grid:
ri =
r0 − R ln qiνi =
R (r0 −
R ln qi)2qi−1
ln−2qiwi = R3r2
(r − r0)−2
e(r−r0)/Rwi
(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 =
RcovσGill-Chien,
with Rcov = 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 q2 = 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.qcenter =
1⁄2, can be used to compute the standardizing
factor as σcenter =
− ln−1(qcenter) =
ln−1(2) = 1.442695.
The MultiExp grid (9.4) is exact if
G(q) =
Q2n−1(q),
which implies
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,
where qi and wi are roots and weights of
the Gauss-Gill quadrature rule, and ri 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.
Table 9.2. − Evaluation Points (ri) of
Gauss-Gill quadrature rule combined with different Transformations of the
Radial Coordinate.(a)
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, r0 = 0, and
rmax = 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 qcenter = 1⁄2.
(c) The first column refers to the grid standardized by means of the
central point of the interval qcenter =
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.
Figure 9.1: MultiExp, Knowles' Log3, Handy, and linear transformations
of the radial coordinate r = r(q) using the roots,
qi, of the 11-point Gauss-Gill quadrature rule
(n = 11, r0 = 0, and rmax = 10).
Table 9.3. − Weights(a) (νi) of
Gauss-Gill quadrature rule combined with different Transformations of the
Radial Coordinate.
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 qcenter = 1⁄2.
(c) The first column refers to the grid standardized by means of the
central point of the interval qcenter =
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.
10. Evenly spaced quadrature: the extended or composite trapezoidal rule
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∫af(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:
where the evaluation points, qi, and the associated
weights, wi, are defined by
qi = a + (i − 1) h
(∀ i = 1, 2, ..., n)
w1 = wn =
h/2 and
w2 = w3 = … =
wn − 1 = h
(10.3)
h =
b − a/n − 1
and f(qi) ≡ fi,
f(a) ≡ f(q1)
≡ f1,
f(b) ≡ f(qn)
≡ fn.
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 [q2,
qn−1] and "semi-open" formulas refer to points
in the intervals [q1, qn−1]
and [q2, qn].
For instance, if f(q) has a singularity at the endpoint
qn, the last integral in the summation of Eq. (10.2)
can be approximated by:
qn∫qn−1f(q) dq ≈ h f(qn−1) ≡ h fn−1
(10.4)
and Eq. (10.2) yields the n-point semi-open trapezoidal formula:
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 w1 = 1⁄2h,
w2 = w3 = … =
wn−1 = h,
and wn = 3⁄2h.
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π∫0f(φ)dφ ≈
2π/nn∑i=1f(φ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 [r0, ∞] 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
+∞∫0r2g(r) dr =
b∫ar2(q) r'(q)
g(r) dq ≡
b∫aG(q) dq
(10.7)
where the integrand is now given by G(q) = r2(q) r'(q)
g(r). Applying the quadrature rule (10.2), yields
b∫aG(q) dq ≈
n∑i=1wiG(qi) =
n∑i=1νig(ri)
(10.8)
where the νi are the transformed weights
νi = r2(qi)
r'(qi) wi
(10.9)
and the qi and wi have expressions like
those of eq. (10.3).
r = r0 + Rqm/(1 − q)m
⇒
dr = m Rqm−1/(1 − q)m+1
dq
(10.10)
The displacement of the origin from 0 to r0 was not
considered in the original paper, but has been added here for completeness.
At the origin is q = 0, r(0) = r0,
r'(0) = 0, and G(0) =
r2(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
qi = i h =
i/n + 1
∀ i = 0, 1, 2, …,
n, n+1
(10.11)
which satisfies the conditions q0 = 0 and
qn+1 = 1. However the quadrature rule comprises
only n points qi (i = 1, 2, …,
n), because the first point does not contribute to the summation as
r'(q0) and G(q0) are
always zero (for m > 1), and the endpoint is discarded, since
G(qn+1) is singular. So, the n
evaluation points ri and the associated transformed weights
νi, which define the Handy grid, are
qi = i h,
w1 = w2 = … =
wn − 1 = h,
wn =
3/2h,
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 - r0/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'(q0) is not
null. Indeed in this case r'(q) = R/(1 − q),
r'(q0) = R, and ν0 is
proportional to Rr02. Then, if
r0 is not null, an (n+1)-point semi-open extended
trapezoid quadrature provides the n evaluation points
qi = (i − 1)h ∀
i = 1, 2, …, n, and the associated weights
w1 = 1⁄2h,
w2 = w3 = … =
wn − 1 = h,
wn = 3⁄2h, 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 R3q3m−1/
(1 − q)3m+1g[r(q)] equals
L(q), a continuous function that is linear
between the roots qi = i/(n + 1) and which
vanishes at q = 0 and q = 1. This implies, for
r0 = 0,
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.
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:
ri = r0 − R ln (1 −
qik)
νi = k Rqik−1/1 − qikwi ri2
∀ i = 1, 2, …, n
(10.15)
qi = i h,
w1 = w2 = … =
wn − 1 = h,
wn =
3/2h,
h =
1/n + 1,
k ≥ 2
For i = n + 1 (i.e. qn+1 = 1) both
ri 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. q0 = 0) it is
ri = r0 and ν0 = 0,
so that the first evaluation point, q0, doesn't contribute
to the integral and therefore has been excluded from the grid. However, if
k = 1, q0 can be excluded from the grid only if
r0 = 0, as ν0 is proportional to
R r02. So, as in the case of the Handy
grid for m = 1, if k = 1 and r0 is not null,
an (n+1)-point semi-open extended
trapezoid quadrature provides the n evaluation points
qi = (i − 1)h ∀
i = 1, 2, …, n, and the associated weights
w1 = 1⁄2h,
w2 = w3 = … =
wn − 1 = h,
wn = 3⁄2h, 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 -
r0/R = ln
2k/2k − 1 =
1/σk.
Here, σk =
ln−12k/2k − 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 r0 = 0 and R = 1, the Knowles grid is exact if
G(q) = kqk−1 ln2(1 −
qk)/
1 − qkg[r(q)] equals L(q),
which implies
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:
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. qn+1 = 0) both
ri 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. q1 = 1) it is
r1 = r0 and ν1 is
proportional to R r02, so that the first
evaluation point contributes to the integral only if
r0 ≠ 0. Otherwise, if r0 = 0, the
n evaluation points may be calculates as
qi = (n − i + 1)/(n + 1)
∀ i = 1, 2, …, n, and the associated
weights as
w1 = w2 = … =
wn−1 = h,
wn = 3⁄2h, with
h = 1/(n + 1).
At the center of the interval, it is q
= 1⁄2, i =
(n+2)⁄2, and r =
r0 + 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 = r0 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:
ri = r0 + R1 +
qi/1
− qi
= r0 + Ri − 1/n − i + 1
νi = 2 Rri2/(1
− qi)2wi
∀ i = 1, 2, …, n
(10.18)
qi =
2 i − n
− 2/n, h =
2/n and
w1 =
1/2h,
w2 = w3 = … =
wn − 1 = h,
wn =
3/2h,
r0 ≠ 0
For i = n + 1 (i.e. qn+1 = 1) both
ri 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. q1 = − 1) it is
r1 = r0 and ν1 is
proportional to R r02, so that the first
evaluation point contributes to the integral only if
r0 ≠ 0. Otherwise, if r0 = 0, the
n evaluation points may be calculates as
qi = (2 i − n − 1)/(n
+ 1) ∀ i = 1, 2, …, n, and the associated
weights as w1 = w2 = … =
wn−1 = h,
wn = 3⁄2h, with
h = 2/(n + 1).
At the center of the interval, it is q
= 0, i = (n+2)/2, and r =
r0 + 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 = r0 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:
qi =
2 i − n
− 1/n + 1, h =
2/n + 1 and
w1 = w2 = … =
wn − 1 = h,
wn =
3/2h
For i = n + 1 (i.e. qn+1 = 1) both
ri 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. q0 = − 1) it is
ri = r0 and ν0 = 0,
so that the first evaluation point, q0, 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 =
r0 + 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".
which provides a distribution of the sample points concentrated at the
origin r0 and with the furthest points at a distance from
the nucleus which does not exceed rmax.
At the origin and for m ≥ 2, it is q = 0, r(0) =
r0, r'(0) = 0, and G(0) =
r2(0) r'(0)
g[r(0)] = 0.
At the endpoint it is q = 1, r(1) = rmax,
r'(1) = m (rmax − r0),
and G(1) = r2(1)
r'(1) g[r(1)] =
mrmax2
(rmax − r0)
g(rmax).
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:
νi = m (rmax −
r0) wi ri2[1 +
(rmax − r0 −
2m)(1 −
qi)m−1]
qim−1/[1 + (rmax −
r0 − 2m) (1 −
qi)m]2
∀ i = 1, 2, ..., n
(10.22)
qi = (i − 1) h,
w1 = wn =
1/2h,
w2 = w3 = … =
wn − 1 = h,
h =
1/n − 1,
m ≥ 2
If m is equal to one, then r'(q0) is not
null and ν1 is proportional to
(rmax − r0)
r02 / (rmax −
r0 − 1).
Then, if r0 is not null, an n-point closed extended
trapezoid quadrature provides the n evaluation points
qi = (i − 1)h ∀
i = 1, 2, …, n, and the associated weights
w1 = wn =
1⁄2h,
w2 = w3 = … =
wn−1 = h,
h = 1/(n − 1), m = 1,
r0 ≠ 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
∈ [r0, rmax] is given by evenly spaced
points, then the corresponding grid is:
ri = r0 + qi
νi = wi ri2
∀ i = 1, 2, …, n
(10.23)
qi = (i − 1) h,
w1 = wn =
1/2h,
w2 = w3 = … =
wn − 1 = h, and
h =
rmax
− r0/n − 1
Such a grid provides sample points equally spaced from the origin
r0 to the furthest point at a distance rmax
from the origin. At the origin it is q1 = 0,
r1 = r0, r'1 = 1,
and G1 = r02g(r0). At the endpoint it is
qn = rmax − r0,
rn = rmax,
rn' = 1, and Gn =
rmax2g(rmax).
At the grid midpoint it is q(n+1)/2 =
rmax − r0, and
r(n+1)/2 = (r0 +
rmax)/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.
Table 10.2. − Evaluation Points (ri) of
Extended Trapezoid Quadrature over Equispaced Points (qi),
combined with different Transformations, r = r(q), of the
Radial Coordinate.(a)
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, r0 = 0, and
rmax = 10. (b) Both standardized (σ = 1.4427) and
non-standardized (σ = 1) grids are reported. (c) This grid is
standardized by construction. (d) Non-standardized grid.
Table 10.3. − Weights(a)
(νi) of Extended Trapezoid Quadrature over Equispaced
Points (qi), combined with different Transformations,
r = r(q), of the Radial
Coordinate.(b)
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, r0 = 0, and
rmax = 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
ri values which increase as qi values
decrease, whereas the latter gives ri values which increase
with qi 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 (r1 = 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.
Figure 10.4: Transformations of the
radial coordinate r = r(q) using a set of n
equally spaced points, qi, which are suitable for the
extended trapezoid quadrature rule (n = 11, r0 = 0,
and rmax = 10).
(a) Gill-Chien transformation (black circles,
eq. 10.15x, open circles refer to the non-standardized grid with σ = 1,
full circles refer to the standardized grid with σ = 1.4427).
(b) Knowles' Log1 grid (blue squares,
eq. 10.13 with k = 1, open squares refer to the non-standardized grid
with σ = 1, full squares refer to the standardized grid with
σ = 1.4427).
(c) Knowles' Log3 grid (orange squares,
eq. 10.13 with k = 3, open squares refer to the non-standardized grid
with σ = 1, full squares refer to the standardized grid with
σ = 7.4889).
(d) Handy grid (violet right triangles,
eq. 10.7 with m = 2, the grid is standardized by construction).
(e) modified Handy grid for finite intervals of r (violet left triangles, eq. 10.16 with m = 2 and
rmax = 10, the grid is standardized by construction).
(f) Becke's grid (dark green up
triangles, eq. 10.15y, the grid is standardized by construction).
(g) Ahlrichs' grid (magenta down triangles,
eq. 10.15z, the grid is standardized by construction).
(h) Linear grid for finite intervals of r (red stars, eq. 10.19 with rmax = 10, the grid is
not standardized).
where f(q) and its derivatives are continuous in
[a, b] and
qi = a + h
(i − 1).
The remainder term Em(f) has the expression:
Em(f) = (a − b)
h2m+2B2m+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 Bj that appear in eqs (10.24) and (10.25) are
Bernoulli numbers
(B0 = 1,
B1 = −1⁄2,
B2 = 1⁄6,
B4 = −1⁄30,
B6 = 1⁄42,
B8 = −1⁄30, ... ;
B3 = B5 = B7 = B9 = ... = 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 = q0 = 0 and tend
to zero at large distances, as
b = qn → ∞.
Under these conditions, eq. (10.24) is often written as
∞∫0f(q) dq = hn∑i=1fi
(10.26)
which is a simplified form of the trapezoidal rule.
11. Tests of radial grid accuracies
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
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.
When the grids of Tables 10.2 and 10.3 are used to estimate the slowly
converging radial integral of the function 1/(1 + r4) the
results shown in the final row of Table 11.2 are obtained, with the grid
ranking
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.
Table 11.3. − Accuracy(a) of radial
quadrature using the grids of Tables 9.2 and 9.3.
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.
Table 11.4. − Accuracy(a) of radial
quadrature using the grids of Tables 9.2 and 9.3.
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.
Concluding Remarks
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.
"An Investigation of the Performance of a Hybrid of
Hartree-Fock and Density Functional Theory"
Gill, P.M.W.; Johnson, B.J.; Pople, J.A.; Frisch, M.J. Intern. J.
Quantum Chem. Symp.1992, 26, 319-331.
"Efficient molecular numerical integration schemes"
Treutler, O.; Ahlrichs, R. J. Chem. Phys.1995, 102,
346-354.
"Improved radial grids for quadrature in molecular
density-functional calculations"
Mura, M. E.; Knowles, P. J. J. Chem. Phys.1996, 104,
9848-9858.
"Molecular integrals by numerical quadrature.
I. Radial integration"
Lindh, R.; Malmqvist, P.-A.; Gagliardi, L. Theor. Chem. Acc.2001, 106, 178-187.
"Radial Quadrature for Multiexponential Integrands"
Gill, P. M. W.; Chien, S.-H. J. Comput. Chem.2003,
24, 732-740.
"An Evaluation of the Radial Part of Numerical
Integration Commonly Used in DFT"
El-Sherbiny, A.; Poirier, R.A. J. Comp. Chem.2004,
25, 1378-1384.
(a) "A simple, reliable and efficient scheme
for automatic numerical integration"
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.
(a) "An adaptive numerical integrator for
molecular integrals"
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.
"Efficient density-functional theory integrations by
locally augmented radial grids"
Gräfenstein, J.; Cremer, D. J. Chem. Phys.2007127, 164113-7.
"A program to generate a basis set adaptive radial
quadrature grid for density functional theory"
Kakhiani, K.; Tsereteli, K.; Tsereteli, P. Computer Physics
Communications2009, 180, 256-268.
"The arrangement of atoms in crystals"
Bragg, W. L. Philos. Mag.1920, 40, 169-189.
"Atomic Radii in Crystals"
Slater, J. C. J. Chem. Phys.1964, 41, 3199-3204.
"Covalent radii revisited"
Cordero, 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.
"Atomic Shielding Constants"
Slater, J. C. Phys. Rev.1930, 36, 57-64.
"Atomic Screening Constants from SCF Functions. II.
Atoms with 37 to 86 Electrons"
Clementi, E.; Raimondi, D. L.; Reinhardt, W. P. J. Chem. Phys.1967, 47, 1300-1307.
"A local projection method for the linear combination of
atomic orbital implementation of density-functional theory"
Yang, W. J. Chem. Phys.1991, 94, 1208-1214.
"Distributed Multipole Analysis: Stability for Large
Basis Sets"
Stone, A. J.
J. Chem. Theory Comput.2005, 1, 1128-1132.
"Electronic Wave Functions. I. A General Method of
Calculation for the Stationary States of Any Molecular System"
Boys, S. F.
Proc. Roy. Soc. A1950, 200, 542-554.
"SG-0: A Small Standard Grid for DFT Quadrature on
Large Systems"
Chien, S.-H.; Gill, P. M. W. J. Comput. Chem.2006,
27, 730-739.
"Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights
of Interpolatory Quadrature"
Elhay, Sylvan; Kautsky, Jaroslav ACM Transactions on Mathematical
Software1987, 13, 399-415.
"Quadrature methods based on the Euler-Maclaurin formula
and on the Clenshaw-Curtis method of integration"
Smith, F. J. Numer. Math.1965, 7, 406-411.
Graham, R. L.; Knuth, D. E.; Patashnik, O.
"Concrete Mathematics", Second Edition, Addison-Wesley Publishing Company,
Inc., 1994, eq. (9.67), p. 469.
Notes
Jacobian matrix. In mathematics, the Jacobian matrix
(sometimes simply called "the Jacobian")
is the matrix of all first-order partial derivatives of a vector-valued
function. Suppose
f : ℝn →
ℝm
is a function which takes as input the vector
x ∈ ℝn
and produces as output the vector
f(x) ∈
ℝm.
Then the Jacobian matrix J of f is an m×n
matrix, usually defined and arranged as follows:
J =
∂ f/∂ x =
(∂ f/∂ x1
…
∂ f/∂ xn) =
∂f1/∂x1
…
∂f1/∂xn
⋮
⋱
⋮
∂fm/∂x1
…
∂fm/∂xn
(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(x1,…,xn), 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.
“Calculus: early transcendentals” by David Guichard and friends, Department of Mathematics, Whitman College, Walla Walla, Washington (accessed April 9, 2021).
