### Summary, reminder and hints

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 ∫a f(q) dq  =  n ∑i=1 qi+1 ∫qi f(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+1qi

 qi+1 ∫qi f(q) dq  ≈  f(ti) wi (2)

The integration (or quadrature) formula (2) represents the (signed) area of a rectangle with height f(ti) and width wi. It is usually referred to as the rectangle rule. [L1-3(L1) "Rectangle rule." Encyclopedia of Mathematics.
(L2) "Numerical Integration - Midpoint, Trapezoid, Simpson’s rule"
https://math.libretexts.org/@go/page/10269
(L3) Wikipedia contributors, "Riemann sum," Wikipedia, The Free Encyclopedia.
] Common choices for ti are the left-hand (ti = qi) and right-hand (ti = qi+1) endpoints of each subinterval as well as their average value (ti = qi+1 + qi2), which give the left-hand, right-hand and midpoint rules, respectively. In addition we consider the average of the left-hand and right-hand rules, which is known as the trapezoidal rule: [L2" Numerical Integration - Midpoint, Trapezoid, Simpson’s rule"
https://math.libretexts.org/@go/page/10269
, L4Wikipedia contributors, "Trapezoidal rule," Wikipedia, The Free Encyclopedia]

 qi+1 ∫qi f(q) dq  ≈    fi wi,   left-hand rule;   fi+1 wi,   right-hand rule;   f( qi + qi+1 2) wi,   midpoint rule;   fi + fi+1 2 wi,   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

 b ∫a f(q) dq  ≈   n ∑i=1 fi wi  ≡  Ln,   composite left-hand rule;  n ∑i=1 fi+1 wi  ≡  Rn,   composite right-hand rule;  n ∑i=1 f( qi + qi+1 2) wi  ≡  Mn,   composite midpoint rule;  1 2 n ∑i=1 (fi + fi+1) wi = 1 2 [ f1 w1 + fn+1 wn + n ∑i=2 fi (wi−1 + wi) ]  ≡  Tn,           composite trapezoidal rule. (4)

When the nodes of the partition are equally spaced, as is often the case, i.e. if wi = w = b − ani, the composite rules in (4) can be simplified for calculation efficiency

 b ∫a f(q) dq  ≈   w n ∑i=1 fi  ≡  Ln,   composite left-hand rule;  w n ∑i=1 fi+1  ≡  Rn,   composite right-hand rule;  w n ∑i=1 f( qi + qi+1 2)  ≡  Mn,   composite midpoint rule;  w [ f1 + fn+1 2  +  n ∑i=2 fi ]  ≡  Tn,   composite trapezoidal rule. (5)
Use the midpoint rule to estimate  ∫5

0
x2 dx  using five equal subintervals. Compare the result with the actual value of this integral.

Each subinterval has length w = 10 − 05 = 2. Therefore, the subintervals consist of

[0,2], [2,4], [4,6], [6,8], and [8,10].

The midpoints of these subintervals are { 1, 3, 5, 7, 9 }. Thus,

M5  =  2 ⋅ [ f(1) + f(3) + f(5) + f(7) + f(9) ] = 2 ⋅ [ 1 + 9 + 25 + 49 + 81 ] = 330.

Since

5 0 x2 dx = [ x33 ]10

0
= 10003,

the error in this approximation is | 10003 − 330 | = 103   or   10331000 ⋅ 100 = 1%.

The three composite rectangle rules are obtained directly from the definition of Riemann sums. [L3Wikipedia contributors, "Riemann sum," Wikipedia, The Free Encyclopedia.] The composite trapezoidal rule is not technically a Riemann sum, even if it is the average of the left and right Riemann sums, i.e. Tn = 12(Ln + Rn). Just as the trapezoidal rule is the average of the left-hand and right-hand rules for estimating definite integrals, a weighted average of the midpoint and trapezoidal rules may be used to obtain Simpson’s rule. It can be shown that S2n = ⅔Mn + ⅓Tn. [L2" Numerical Integration - Midpoint, Trapezoid, Simpson’s rule"
https://math.libretexts.org/@go/page/10269
, L5Wikipedia contributors, "Simpson's rule," Wikipedia, The Free Encyclopedia]

The trapezoid and Simpson's rules belong to a wide class of quadrature formulas based on the Lagrange interpolation [L6Wikipedia contributors, "Polynomial interpolation," Wikipedia, The Free Encyclopedia] of equally spaced points, known as Newton-Cotes formulas.[L7-8(L7) Weisstein, Eric W. "Newton-Cotes Formulas." From MathWorld − A Wolfram Web Resource.
(L8) Wikipedia contributors, "Newton–Cotes formulas," Wikipedia, The Free Encyclopedia.
] Recall that if f(q) is a function defined on a closed domain [a,b] of real numbers q, the interval [a,b] can be partitioned in n sub-intervals [qi, qi+1] of width wi = qi+1qi where the qi's are a finite sequence of n+1 numbers of the form a = q1  <  …  <  qn+1 = b. The Lagrange polynomial which interpolates the n+1 points (qi, fi) has order at most equal to n (it is important not to confuse the number of points of the quadrature rule with the order of the interpolating polynomial).

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

### 1. - Newton-Cotes formulas.

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=1 ln,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):

 b ∫a f(q) dq  ≈  b ∫a Ln(q) dq  =  n+1 ∑i=1 [ b ∫a ln,i(q) dq ] f(qi)  =  n+1 ∑i=1 wi f(qi) (1.2)

The expression of the Lagrange coefficients (cf Eq. (2.7) of the Introduction)

 ln,i(q)  =  n+1 ∏j=1; j≠i q − qj qi − qj (1.3)

can be considerably simplified in the case of equally spaced points, as we can write [1Apostol, T. M. Calculus; Volume II; Second Edition; John Wiley: New York: 1969. Section 15.7, p. 582.]

 qi  =  q1 + (i − 1) h     ∀   i = 1, 2, …, n+1 (1.4)

where h denotes the distance between adjacent points. Since   qiqj = (i − j)h,   eq (1.3) becomes

 ln,i(q)  =  n+1 ∏j=1; j≠i q − q1 − (j − 1) h (i − j) h  =  n+1 ∏j=1; j≠i s − j + 1 i − j, (1.5)

where

 s  =  q − q1 h. (1.6)

In the last term on the right of eq. (1.5) the product of the factors independent of s is

 n+1 ∏j=1; j≠i 1 i − j  =  ( i−1 ∏j=1 1 i − j ) ( n+1 ∏j=i+1 1 i − j )  =  1 (i − 1)! n+1 ∏j=i+1 (−1) j − i  =  (−1)n−i+1 (i − 1)! (n − i + 1)! (1.7)

Since q = q1 + s h, eq.(1.5) now becomes

 ln,i(s)  =  (−1)n−i+1 (i − 1)! (n − i + 1)! n+1 ∏j=1; j≠i (s − j + 1)  =  (−1)n−i+1 n! ( ni − 1 ) n+1 ∏j=1; j≠i (s − j + 1) (1.8)

and the integrals which in Eq. (1.2) are enclosed in square brackets become

 wi  =  b ∫a ln,i(q) dq  =  h n ∫0 ln,i(s) ds  =  (−1)n−i+1 h n! ( ni − 1 ) n ∫0 ( n+1 ∏j=1; j≠i (s − j + 1) ) ds (1.9)

for   i = 1, 2, …, n+1.  Note that for  q = q1 we have s = 0  and for  q = qn+1 = q1 + n h  we have  s = n,  and that  dq = h ds.

It should be noted that, based on our preferences, we are using 1-based indexing. However, it is very common practice to use zero-based numbering, [L9Wikipedia contributors, "Zero-based numbering," Wikipedia, The Free Encyclopedia.] which would have the advantage of writing eq. (1.9) in a somewhat simpler and more elegant way. [N1Eq. (1.9) using zero-based indexing.]

Table 1 lists some of the Newton–Cotes formulas of the closed type. Their derivation is illustrated in the exercises that follow the table.

Table 1.  Closed Newton-Cotes formulas.§
nb a f(q) dq  ≈  n+1 i=1 wi f(qi) common name
1 h 2 (f1 + f2) 2-point closed formula
or trapezoidal rule
2 1 3 h (f1 + 4 f2 + f3) 3-point closed formula
or Simpson's rule
3 3 8 h (f1 + 3 f2 + 3 f3 + f4) 4-point closed formula
or Simpson's second rule
4 2 45 h (7 f1 + 32 f2 + 12 f3 + 32 f4 + 7 f5) 5-point closed formula
or Boole's rule
5 5 288 h (19 f1 + 75 f2 + 50 f3 + 50 f4 + 75 f5 + 19 f6) 6-point closed formula
§ qi = a + (i−1)b−an = a + (i−1)h, for i = 1, 2, …, n+1. The notation fi is a shorthand for f(qi).
Find the (n+1)-point closed Newton-Cotes formulas for n = 1, 2, 3, and 4, using equations (1.1-1.3).

(♠)   n = 1:  2-point closed formula or trapezoidal rule. The Lagrange interpolating polynomial of eq. (1.1) becomes

L1(q)  =  l1,1(q) f1 + l1,2(q) f2  =  q − q2q1q2 f1 + q − q1q2q1 f2  =  f2f1 q2q1 q + q2 f1q1 f2 q2q1

where fi = f(qi)  and  h = q2q1 = b − a. Integrating over the domain [q1, q2] ≡ [a, b], yields

b a L1(q) dq  =  f1 b a l1,1(q) dq + f2 b a l1,2(q) dq  =  f2f1 q2q1 [ q2 2 ]q2

q1
+ q2 f1q1 f2 q2q1 [ q ]q2

q1
=  f2f1 q2q1 (q2q1) (q2 + q1) 2 + q2 f1q1 f2 q2q1 (q2q1)  =  1 2(q2q1) (f1 + f2)  =  1 2 h (f1 + f2)

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 = ba. 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

L2(q)  =  l2,1(q) f1 + l2,2(q) f2 + l2,3(q) f3  =  (q − q2) (q − q3) (q1q2) (q1q3) f1 + (q − q1) (q − q3) (q2q1) (q2q3) f2 + (q − q1) (q − q2) (q3q1) (q3q2) f3

We use integration by substitution to estimate the integral  ∫b

a
L2(q) dq. Since points q1 = a < q2 < q3 = b are equally spaced by h = q3q1 2, we have the following relations q2 = q1 + q3 2, q1 = q2h, and q3 = q2 + h, which lead to the expression

L2(q)  =  (q − q2) (q − q2h) 2 h2 f1 + h2 − (q − q2)2 h2 f2 + (q − q2) (q − q2 + h) 2 h2 f3  =  q − q2 2 h ( q − q2 h − 1 ) f1 + [1 − ( q − q2 h )2

] f2
+ q − q2 2 h ( q − q2 h + 1 ) f3

Introducing the new variable  u = q − q2 h,  we get

L2(u)  =  1 2 [ u (u − 1) f1  +  2 (1 − u2) f2  +  u (u + 1) f3 ]

Since  u = −1  for  q = q1 = q2h,  u = 1  for  q = q3 = q2 + h,  and  h du = dq,  we simply have to calculate the integral

b a L2(q) dq  =  h 1 −1 L2(u) du  =  h 2 ( f1 1 −1 u (u − 1) du  +  2 f2 1 −1 (1 − u2) du  +  f3 1 −1 u (u + 1) du )  =  h 2 ( f1 [ u3 3u2 2 ]1

−1
+  2 f2 [ uu3 3 ]1

−1
+  f3 [ u3 3 + u2 2 ]1

−1
)
=  h 3 (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

L3(q)  =  (q − q2) (q − q3) (q − q4) (q1q2) (q1q3) (q1q4) f1  +  (q − q1) (q − q3) (q − q4) (q2q1) (q2q3) (q2q4) f2  +  (q − q1) (q − q2) (q − q4) (q3q1) (q3q2) (q3q4) f3  +  (q − q1) (q − q2) (q − q3) (q4q1) (q4q2) (q4q3) f4

and, integrating, we get Simpson's second rule

b a L3(q) dq  =  3 8 h (f1 + 3 f2 + 3 f3 + f4)

which is exact if f is a polynomial up to cubic degree. Because of the 3⁄8 factor Simpson's second rule is also referred to as Simpson's 3⁄8 rule.

(♦)   n = 4:  5-point closed formula or Boole's rule. The Lagrange interpolating polynomial is

L4(q)  =  (q − q2) (q − q3) (q − q4) (q − q5) (q1q2) (q1q3) (q1q4) (q1q5) f1  +  (q − q1) (q − q3) (q − q4) (q − q5) (q2q1) (q2q3) (q2q4) (q2q5) f2  +  (q − q1) (q − q2) (q − q4) (q − q5) (q3q1) (q3q2) (q3q4) (q3q5) f3  +  (q − q1) (q − q2) (q − q3) (q − q5) (q4q1) (q4q2) (q4q3) (q4q5) f4  +  (q − q1) (q − q2) (q − q3) (q − q4) (q5q1) (q5q2) (q5q3) (q5q4) f5

and, integrating, we get Boole's rule

b a L4(q) dq  =  2 45 h (7 f1 + 32 f2 + 12 f3 + 32 f4 + 7 f5)

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) h n 0 n+1 j=1; j≠i (s − j + 1) ds wiαi h
11 −1s−1 1 2h 1 2h
21s 1 2h 1 2h
21 1 2 (s−1)(s−2) 2 3h 1 3h
2−1s(s−2) 4 3h 4 3h
31 2 s(s−1) 2 3h 1 3h
31 1 6 (s−1)(s−2)(s−3) 9 4h 3 8h
2 1 2 s(s−2)(s−3) 9 4h 9 8h
3 1 2 s(s−1)(s−3) 9 4h 9 8h
4 1 6 s(s−1)(s−2) 9 4h 3 8h
41 1 24 (s−1)(s−2)(s−3) (s−4) 112 15h 14 45h
2 1 6 s(s−2)(s−3)(s−4) 128 15h 64 45h
3 1 4 s(s−1)(s−3)(s−4) 32 15h 8 15h
4 1 6 s(s−1)(s−2)(s−4) 128 15h 64 45h
5 1 24 s(s−1)(s−2)(s−3) 112 15h 14 45h
51 1 120 (s−1)(s−2)(s−3) (s−4)(s−5) 475 12h 95 288h
2 1 24 s(s−2)(s−3)(s−4) (s−5) 125 4h 125 96h
3 1 12 s(s−1)(s−3)(s−4) (s−5) 125 12h 125 144h
4 1 12 s(s−1)(s−2)(s−4) (s−5) 125 12h 125 144h
5 1 24 s(s−1)(s−2)(s−3) (s−5) 125 4h 125 96h
6 1 120 s(s−1)(s−2)(s−3) (s−4) 475 12h 95 288h

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 \ i123 45678910 11
1 1 2 1 2
2 1 3 4 3 1 3
3 3 8 9 8 9 8 3 8
4 14 45 64 45 8 15 64 45 14 45
5 95 288 125 96 125 144 125 144 125 96 95 288
6 41 140 216 140 27 140 272 140 27 140 216 140 41 140
7 5257 17280 25039 17280 9621 17280 20923 17280 20923 17280 9621 17280 25039 17280 5257 17280
8 3956 14175 23252 14175 3712 14175 41984 14175 18160 14175 41984 14175 3712 14175 23252 14175 3956 14175
9 25713 89600 141669 89600 9720 89600 174096 89600 52002 89600 52002 89600 174096 89600 9720 89600 141669 89600 25713 89600
10 80335 299376 531500 299376 242625 299376 1362000 299376 1302750 299376 2136840 299376 1302750 299376 1362000 299376 242625 299376 531500 299376 80335 299376

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 ∫a f(q) dq  ≈  n ∑i=2 wi f(qi) (1.10)
and
 wi  =  (−1)n−i+1 h n! ( 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)]:

L2(q)  ≈  l2,2(q) f(q2)  =  (q − q1) (q − q3) (q2q1) (q2q3) f(q2)

Since the points are equally spaced by h = q3q1 2 = ba 2, we have the following relations q2 = q1 + q3 2 = a + b 2,   q1 = q2h,  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

f( a + b 2 ) b a l2,2(q) dq  =  f( a + b 2 ) 2 0 s (1 − s) ds  =

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 + b 2 ), 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 = q4q13 = b − a 3 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

L1(q)  =  l1,2(q) f2 + l1,3(q) f3  =  q − q3q2q3 f2 + q − q2q3q2 f3  =  1 h [(q3q) f2 + (q − q2) f3]

where the equality q3q2 = h has been used. Integrating over the domain [q1, q4] ≡ [q2h, q3+h], yields

b a L1(q) dq  =  f2 h q3+h q2h (q3q) dq  +  f3 h q3+h q2h (q − q2) dq  =  f2 2h (q3q2)(q3q2 + 2h)  +  f3 2h (q3q2)(q3q2 + 2h)  =  3 2 h (f2 + f3)  =  b − a 2 (f2 + f3)

The last term in the equation above is the 4-point open Newton-Cotes formula or trapezoidal rule, which is exact if f is a polynomial of first degree.

(♥)   n = 4:  5-point open formula or trapezoidal rule. The Lagrange interpolating polynomial of eq. (1.1) becomes

MISSING FORMULA

(♦)   n = 5:  6-point open formula or trapezoidal rule. The Lagrange interpolating polynomial of eq. (1.1) becomes

MISSING FORMULA

1. Apostol, T. M. Calculus; Volume II; Second Edition; John Wiley: New York: 1969. Section 15.7, p. 582. ISBN 0-471-00007-8
• Stewart, J. Calculus Early Transcendentals; Sixth Edition; Thomson Brooks/Cole: 2008. Section 5.2, p. 366. ISBN-13: 978-0-495-01166-8; ISBN-10: 0-495-01166-5.
1. Given a sequence of (n + 1) equally spaced points, indexed from 0 to n, so that qi = q0 + ih for each i = 0, 1, 2, …, n, eq. (1.9) becomes
 wi  =  (−1)n−i h n! ( ni ) n ∫0 n ∏j=0;j≠i (s − j) ds           ∀ i = 0, 1, 2, …, n.
Compare with Eq. (1.9) to see the difference between 0- and 1-based indexing.

1. "Rectangle rule." Encyclopedia of Mathematics. https://encyclopediaofmath.org/wiki/Rectangle_rule (accessed April 21, 2021).
2. "Numerical Integration - Midpoint, Trapezoid, Simpson’s rule". https://math.libretexts.org/@go/page/10269 (accessed April 21, 2021).
3. Wikipedia contributors, "Riemann sum," Wikipedia, The Free Encyclopedia, https://en.wikipedia.org/wiki/Riemann_sum (accessed April 11, 2021).
4. Wikipedia contributors, "Trapezoidal rule," Wikipedia, The Free Encyclopedia, https://en.wikipedia.org/wiki/Trapezoidal_rule (accessed April 21, 2021).
5. Wikipedia contributors, "Simpson's rule," Wikipedia, The Free Encyclopedia, https://en.wikipedia.org/wiki/Simpson%27s_rule (accessed April 22, 2021).
6. Wikipedia contributors, "Polynomial interpolation," Wikipedia, The Free Encyclopedia, https://en.wikipedia.org/wiki/Polynomial_interpolation (accessed April 11, 2021).
7. Weisstein, Eric W. "Newton-Cotes Formulas." From MathWorld − A Wolfram Web Resource. https://mathworld.wolfram.com/Newton-CotesFormulas.html (accessed May 7, 2021).
8. Wikipedia contributors, "Newton–Cotes formulas," Wikipedia, The Free Encyclopedia, https://en.wikipedia.org/wiki/Newton%E2%80%93Cotes_formulas (accessed May 7, 2021).
9. Wikipedia contributors, "Zero-based numbering," Wikipedia, The Free Encyclopedia, https://en.wikipedia.org/wiki/Zero-based_numbering (accessed May 10, 2021).

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

1.1. - Newton-Cotes quadrature formulas. The interpolatory procedure outlined above leads to an extremely useful and straightforward family of numerical integration techniques, which are known as Newton-Cotes formulas and apply to n equispaced quadrature points or nodes in the range [a, b], such that 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 Physik 1901, 46, 224-243.
] This implies that going to higher degrees does not always improve accuracy. Runge's phenomenon can be avoided by using a composite or extended rule.

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

 b = qn ∫ a = q1 f(q) dq = n−1 ∑ i=1 qi+1 ∫ qi f(q) dq ≈ 1 / 2 n−1 ∑ i=1 (qi+1 − qi) (fi + fi+1) = h ( 1 / 2 f1 + n−1 ∑ i=2 fi + 1 / 2 fn ) (1.10)

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

 b = qn ∫ a = q1 f(q) dq − h ( 1 / 2 f1 + n−1 ∑ i=2 fi + 1 / 2 fn ) = m ∑ k=1 B2k / (2k)! h2k [ f (2k−1)(a) − f (2k−1)(b) ] + Em(f) (1.11)

The numbers Bj that appear in eq. (1.11) are Bernoulli numbers (B0 = 1, B1 = −12, B2 = 16, B4 = −130, B6 = 142, B8 = −130, ... ; 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

1. pi(q) is a polynomial of degree i;
2. 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
 ≔ b ∫ a ω(q) pi(q) pj(q) dq = δij b ∫ a ω(q) pi2(q) dq ≕ δij (1.12)

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=1 vi 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 ≤ kn   simultaneously obey the relation

 b ∫ a ω(q) pk(q) dq = n ∑ i=1 vi 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=0 ck pk(q) (1.16)

where

 ck = b ∫ a f(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:

 b ∫ a f(q) dq = n ∑ k=0 ck b ∫ a ω(q) pk(q) dq eq. (1.14) = n ∑ k=0 ck n ∑ i=1 vi pk(qi) = n ∑ i=1 vi n ∑ k=0 ck pk(qi) eq. (1.16) = n ∑ i=1 wi f(qi) (1.18)

where

 wi = vi / ω(qi) (1.19)

Eq. (1.18) shows that the projection of the integrand f(q) into a basis of orthogonal polynomials occurs implicitly, with no need to explicitly evaluate the expansion coefficients. It states that, given an integer n, we can find a set of weigths 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=1 vi 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.
ω(q) Interval
[a,b]
Orthogonal polynomials weights radrule
1 [-1, 1] n-th Legendre polynomial
Pn(q) = 1 / 2n n/2⌋ k=0 (−1)k (2n − 2k)! / k! (nk)! (n − 2k)! qn−2k
wi = 2 / (1 − qi2) [P'n(qi)]2
(i = 1,2,…,n)
10
(1 − q2)½ [-1, 1] n-th Chebyshev polynomial of the second kind
Un(q) = n/2⌋ k=0 (-1)k ( n − k   k ) (2q)n−2k
with analytical roots:
qi = cos ( i π / n + 1 )     (i = 1,2,…, n)
wi = π / n + 1 (1 − qi2)
(i = 1,2,…, n)
20
ln2q [0, 1] Gill polynomials
Q0(q) = 1,      Q1(q) = 8q − 1,
Q2(q) = 7992q2 − 4104q + 217,      …
wi
(i = 1,2,…,n)
30
eq [0, ∞] n-th Laguerre polynomial
Ln(q) = n k=0 (−1)k / k! ( n / k ) qk
wi = qi / (n + 1)2 [Ln+1(qi) ]2
(i = 1,2,…, n)
40
qα eq,    α > −1 [0, ∞] n-th Generalized Laguerre polynomial
Ln(α)(q) = n k=0 ( n + α   n − k) (−1)k / k! qk
wi
(i = 1,2,…,n)
50
eq2 [−∞, ∞] n-th Hermite polynomial
Hn(q) = (−1)n eq²   dn / dqn eq²
wi = 2n−1 n! π / n2 [Hn−1(qi) ]2
(i = 1,2,…, n)
60

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

1.3. - Degree of exactness of a quadrature rule. In the following, we will use I(f) to denote the integral over a finite interval [a, b] that appears in the left hand side of eq. (1.1) or (1.18), and 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:

 b ∫ a f(q) ω(q) dq = b ∫ a pn(q) t(q) ω(q) dq + b ∫ a r(q) ω(q) dq = b ∫ a r(q) ω(q) dq (1.23)

The green integral in eq. (1.23) is zero because the polynomial 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

 b ∫ a r(q) ω(q) dq = n ∑ i=1 wi r(qi) = n ∑ i=1 wi pn(qi) t(qi) + n ∑ i=1 wi r(qi) = n ∑ i=1 wi [pn(qi) t(qi) + r(qi)] = n ∑ i=1 wi f(qi) (1.24)

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

 b ∫ a f(q) ω(q) dq = b ∫ a pn(q) t(q) ω(q) dq + b ∫ a r(q) ω(q) dq ≠ b ∫ a r(q) ω(q) dq

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

 n ∑ i=1 wi f(qi) = n ∑ i=1 wi [pn(qi) t(qi) + r(qi)] = n ∑ i=1 wi r(qi)

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=1 wy,j f(x,yj) (1.25)

Next we approximate the integral by the product quadrature formula

 ∫∫ f(x,y) dx dy ≈  nx ∑ i=1 ny ∑ j=1 wx,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=1 ny ∑ j=1 nz ∑ k=1 wx,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

 dr2 = dx2 + dy2 + dz2 = g11 du12 + g22 du22 + g33 du32 + 2 g12 du1 du2 + 2 g13 du1 du3 + 2 g23 du2 du3 (1.32)

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

 dV = | ∂ f / ∂ 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 / 2 u + a + b / 2 ≡ f(u) (1.35)

From this equation:

 f'(u) = df(u) / du = b − a / 2 (1.36)

and

 f−1(r) = 2r − a − b / b − a ≡ u      ⇒      f−1(a) = − 1,       f−1(b) = + 1 (1.37)

so that eq (1.34) can be rewritten as:

 b ∫ a φ(r) dr = +1 ∫ −1 b − a / 2 φ(f(u)) du (1.38)

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

 b ∫ a φ(r) dr = +1 ∫ −1 b − a / 2 ω(u) ω(u)−1 φ(f(u)) du = +1 ∫ −1 ω(u) Φ(u) du (1.39)

where it has been set:

 Φ(u) = b − a / 2 ω(u)−1 φ(f(u)) (1.40)

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

 b ∫ a φ(r) dr ≈ n ∑ i=1 wi Φ(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=1 vi φ(ri) (1.42)

where the transformed weights, vi, are given by

 vi = b − a / 2 wi (1.43)

in the case of Gauss-Legendre quadrature, and

 vi = (b − a)2 / 4[(r − a)(b − r)]1/2 wi (1.44)

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

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

 r = a + (b − a) u ≡ f(u) ⇓ f'(u) = b − a,             f−1(r) = r − a / b − a ≡ u,             f−1(a) = 0,             f−1(b) = 1 (1.45)

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

 b ∫ a φ(r) dr = 1 ∫ 0 (b − a) φ(f(u)) du = 1 ∫ 0 (b − a) ω(u) ω(u)−1 φ(f(u)) du = 1 ∫ 0 ω(u) Φ(u) du ≈ n ∑ i=1 wi Φ(ui) = n ∑ i=1 vi φ(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 )−2 wi (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 + S u ≡ 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 = ∞ ∫ 0 S φ(f(u)) du = ∞ ∫ 0 e−u e+u S φ(f(u)) du = ∞ ∫ 0 e−u Φ(u) du ≈ n ∑ i=1 wi Φ(ui) = n ∑ i=1 vi f(ri) (1.50)

where S is a scaling factor, ω(u) = eu 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 = ∫ V f(r) d3r = ∫ V f(r) dV = ∭ V f(x,y,z) dx dy dz (2.1)

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

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

 ∀  f(x,y,z) = gX(x) ⋅ gY(y) ⋅ gZ(z) ⇓ ∭ V f(x,y,z) dx dy dz = x2 ∫ x1 gX(x) dx ⋅ y2 ∫ y1 gY(y) dy ⋅ z2 ∫ z1 gZ(z) dz (2.2)

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

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

 x = r cosθ sinφ,          y = r sinθ sinφ,          z = r cosθ (2.5)

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

J(r,θ,φ) =
 ∂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 Ω:

 I = +∞ ∫ 0 r2 dr π ∫ 0 sinθ dθ 2π ∫ 0 f(r,θ,φ) dφ = +∞ ∫ 0 r2 dr ∫ Ω f(r,Ω) dΩ (2.10)

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

 I = +∞ ∫ 0 r2 g(r) dr (2.11)

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

 g(r) = π ∫ 0 sinθ dθ 2π ∫ 0 f(r,θ,φ) dφ = ∫ Ω f(r,Ω) dΩ (2.12)

The integral is numerically approximated as a sum of weighted function values, sampled over a given set of points ci = ri, θi, φi, Ωi:

 ∫ V f(r) dV ≈ Nr ∑ i=1 wr,i Nθ ∑ j=1 wθ,j Nφ ∑ k=1 wφ,k f(ri,θj,φk) ≈ Nr ∑ i=1 wr,i NΩ ∑ j=1 wΩ,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:

### 2.2. - Angular integral.

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.

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

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

 +∞ ∫ 0 r2 g(r) dr = b ∫ a ω(q) ω−1(q) r2(q) r'(q) g[r(q)] dq = b ∫ a ω(q) G(q) dq ≈ n ∑ i=1 wi G(qi) = n ∑ i=1 wi ω−1(qi) r2(qi) r'(qi) g[r(qi)] = n ∑ i=1 vi g(qi) (2.14)

where we have set G(q) = ω−1(q) r2(qr'(q) g[r(q)] and vi = wi ω−1(qi) r2(qir'(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.

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

The transformed coordinate 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.
{ri} grid rir0 qi {qi} radrule
[r0, ∞] MultiExpR ln qi exp [−(rir0)/R][0, 1]31, 71
[r0, ∞] KnowlesR ln (1 − qik) k1 − exp [−(rir0)/R][0, 1] 32, 72
[r0, ∞] Handy R qim / (1 − qi)m mrir0 / mRmrir0 [0, 1]33, 73
[r0, rmax] Handy (rmaxr0) qim / 1 + (rmaxr0 − 2m)(1 − qi)m    m (rir0) (rmaxr0 − 2m + 1) / (rir0)(rmaxr0 − 2m) + rmaxr0 [0, 1]34, 74
[r0, ∞] Becke R 1 + qi / 1 − qi rir0R / rir0 + R [−1, +1] 15, 25, 75
[r0, ∞] Ahlrichs R (1 + qi)α / ln 2 ln 1 − qi / 2  [−1, +1]16, 26, 76
[r0, rmax] linear rmaxr0 / 2(1 + qi) 2 ri − (rmax + r0) / rmaxr0 [−1, +1]17, 27, 77
[r0, rmax] linear (rmaxr0) qi rir0 / rmaxr0 [0, 1]37, 77
[r0, ∞] linear R qi rir0 / R[0, ∞] 47, 57, 77

2.3.4. - Standardization of quadrature grids. In order to compare the roots and weights from the various quadrature schemes, Gill and Chien[10Gill, P. M. W.; Chien, S.-H.
J. Comput. Chem. 2003, 24, 732-740.
] chose the R values so that the middle point of each quadrature is unity, i.e. r(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 = 12. Among the qudrature rules available in PAMoC, the Gauss-Gill and the extended trapezoidal rules have roots in the range [0, 1], and therefore can be used in combination with the MultiExp, Knowles and Handy grids, which are standardized by the following factors:

 σMutiExp = − 1 / ln qcenter = 1 / ln 2 = 1.442695 σKnowles = − 1 / ln (1 − qkcenter)      (e.g.: σLog1 = 1.442695  and  σLog3 = 7.488876) σHandy = ( 1 − qcenter / qcenter)m = 1

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:

 σLaguerre = 1 / q(n+1)/2

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

 H He 0.25 Li Be B C N O F Ne 1.45 1.05 0.85 0.70 0.65 0.60 0.50 Na Mg Al Si P S Cl Ar 1.80 1.50 1.25 1.10 1.00 1.00 1.00 K Ca Sc Ti V Cr Mn Fe Co Ni Cu Zn Ga Ge As Se Br Kr 2.20 1.80 1.60 1.40 1.35 1.40 1.40 1.40 1.35 1.35 1.35 1.35 1.30 1.25 1.15 1.15 1.15 Rb Sr Y Zr Nb Mo Tc Ru Rh Pd Ag Cd In Sn Sb Te I Xe 2.35 2.00 1.80 1.55 1.45 1.45 1.35 1.30 1.35 1.40 1.60 1.55 1.55 1.45 1.45 1.40 1.40 Cs Ba La Hf Ta W Re Os Ir Pt Au Hg Tl Pb Bi Po At Rn 2.60 2.15 1.95 1.55 1.45 1.35 1.35 1.30 1.35 1.35 1.35 1.50 1.90 1.80 1.60 1.90 ↓ Ce Pr Nd Pm Sm Eu Gd Tb Dy Ho Er Tm Yb Lu 1.85 1.85 1.85 1.85 1.85 1.85 1.80 1.75 1.75 1.75 1.75 1.75 1.75 1.75
 H He 0.31 0.28 Li Be B C N O F Ne 1.28 0.96 0.84 0.76 0.71 0.66 0.57 0.58 Na Mg Al Si P S Cl Ar 1.66 1.41 1.21 1.11 1.07 1.05 1.02 1.06 K Ca Sc Ti V Cr Mn Fe Co Ni Cu Zn Ga Ge As Se Br Kr 2.03 1.76 1.70 1.60 1.53 1.39 1.39 1.32 1.26 1.24 1.32 1.22 1.22 1.20 1.19 1.20 1.20 1.16 Rb Sr Y Zr Nb Mo Tc Ru Rh Pd Ag Cd In Sn Sb Te I Xe 2.20 1.95 1.90 1.75 1.64 1.54 1.47 1.46 1.42 1.39 1.45 1.44 1.42 1.39 1.39 1.38 1.39 1.40 Cs Ba La Hf Ta W Re Os Ir Pt Au Hg Tl Pb Bi Po At Rn 2.44 2.15 2.07 1.75 1.70 1.62 1.51 1.44 1.41 1.36 1.36 1.32 1.45 1.46 1.48 1.40 1.50 1.50 ↓ Ce Pr Nd Pm Sm Eu Gd Tb Dy Ho Er Tm Yb Lu 2.04 2.03 2.01 1.99 1.98 1.98 1.96 1.94 1.92 1.92 1.89 1.90 1.87 1.87

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

 H He 1.0000 0.5882 Li Be B C N O F Ne 3.0769 2.0513 1.5385 1.2308 1.0256 0.8791 0.7692 0.6838 Na Mg Al Si P S Cl Ar 4.0909 3.1579 2.5714 2.1687 1.8750 1.6514 1.4754 1.3333
 H He a0 0.31 Li Be B C N O F Ne 1.67 1.12 0.87 0.67 0.56 0.48 0.42 0.38 Na Mg Al Si P S Cl Ar 1.90 1.45 1.18 1.11 0.98 0.88 0.79 0.71 K Ca Sc Ti V Cr Mn Fe Co Ni Cu Zn Ga Ge As Se Br Kr 2.43 1.94 1.84 1.76 1.71 1.66 1.61 1.56 1.52 1.49 1.45 1.42 1.36 1.25 1.14 1.63 0.94 0.88 Rb Sr Y Zr Nb Mo Tc Ru Rh Pd Ag Cd In Sn Sb Te I Xe 2.65 2.19 2.12 2.06 1.98 1.90 1.83 1.78 1.73 1.69 1.65 1.61 1.56 1.45 1.33 1.23 1.15 1.08 Cs Ba La Hf Ta W Re Os Ir Pt Au Hg Tl Pb Bi Po At Rn 2.98 2.53 2.50 2.08 2.00 1.93 1.88 1.85 1.80 1.77 1.74 1.71 1.56 1.54 1.43 1.35 1.27 1.20 ↓ Ce Pr Nd Pm Sm Eu Gd Tb Dy Ho Er Tm Yb Lu 2.50 2.08 2.00 1.93 1.88 1.85 1.80 1.77 1.74 1.71 1.56 1.54 1.43 1.35 Figure 2.1: Covalent radius (bohr) as a function of the atomic number: (a) on mouse out, for atoms H to Ar, (b) on mouse over, for atoms H to Cm.
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.25000.50001.00004.5000
Li-Ne0.16670.50000.90003.5000
Na-Ar0.10000.40000.80002.5000
K-Kr 0.02000.10000.20003.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 86194 302 590 302
(110,770) 47 110230 350 770 350
(150,974) 53 146266 434 974 434
(185,1202)59 170350 5901202 590
(225,5720)131 434590120257201202
Spherical product grid
(Nr,Nθ,Nφ) L α1α2 α3α4 outer
region
(50,15,30)29 50128 200 450 200
(75,25,50)49 128200 4501250 450
(99,48,96)95 200450125046081250

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.

Table 2.7 − Lebedev partition, number Nr of radial points, scale factor R and total number Ntot of grid points in the SG-0 grid (from ref 26Chien, S.-H.; Gill, P. M. W.;
J. Comput. Chem. 2006, 27, 730-739.
).
Element Lebedev partition Nr R Ntot
H 66 183 261 381 741 1101 1466 861 501 381 181 231.301406
Li 66 183 261 381 741 1101 1466 861 501 381 181 231.951406
Be 64 182 261 382 741 861 1102 1465 501 381 181 62 232.201390
B 64 264 383 863 1466 381 62 231.451426
C 66 182 261 382 502 861 1101 1461 1702 1462 861 381 181 231.201390
N 66 183 261 382 742 1101 1702 1463 861 502 231.101414
O 65 181 262 381 504 861 1105 861 501 381 61 231.101154
F 64 382 504 742 1102 1462 1102 863 501 61 231.201494
Na 66 182 263 381 502 1108 742 62 262.301328
Mg 65 182 262 382 502 741 1102 1464 1101 861 382  181 61 262.201492
Al 66 182 261 382 502 741 861 1462 1702 1102 861 741 261 181 61 262.101496
Si 65 184 384 503 741 1102 1461 1703 861 501 61 261.301496
P 65 184 384 503 741 1102 1461 1703 861 501 61 261.301496
S 64 181 268 382 501 742 1101 1703 1461 1101 501 61 261.101456
Cl 64 187 262 382 501 741 1102 1703 1461 1101 861 61 261.451480

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:

 g(r) = π ∫ 0 sinθ dθ 2π ∫ 0 f(r,θ,φ) dφ = ∫ Ω f(r,Ω) dΩ  ≈  n ∑ i=1 wi f(r, Ωi) = n ∑ i=1 wi f(r, θi, φi) (9.1)
Table 9.1 − Degree L and points NΩ (L + 1)2 / 3 of Lebedev rules.
L NΩ L NΩ L NΩ L NΩ
3 6 19146 41590 89 2702
5 14 21170 47770 95 3074
7 26 23194 53974 101 3470
9 38 25230 591202 107 3890
11 50 27266 651454 113 4334
13 74 29 302 71 1730 119 4802
15 86 31 350 77 2030 125 5294
17 110 35 434 83 2354 131 5810

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

## 3. Gauss-Legendre quadrature

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

 +1 ∫ −1 f(q) dq ≈  n ∑ i=1 wi f(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.
Pn(q) = 1 / 2n n/2⌋ k=0 (−1)k (2n − 2k)! / k! (nk)! (n − 2k)! qn−2k qi wi = 2 / (1 − qi2) [P'n(qi)]2
P0(q) = 1
P1(q) = q q1 = 0w1 = 2
P2(q) = (3q2 − 1)/2 q1,2 = ∓ = ∓0.577350w1,2 = 1
P3(q) = (5q3 − 3q)/2 q1,3 = ∓ = ∓0.774597 w1,3 = 59 = 0.555556
q2 = 0w2 = 89 = 0.888889
P4(q) = (35q4 − 30q2 + 3)/8 q1,4 = ∓37 + 2765 = ∓0.861136 w1,4 = (18 − 30)/36 = 0.347855
q2,3 = ∓372765 = ∓0.339981 w2,3 = (18 + 30)/36 = 0.652145
P5(q) = (63q5 − 70q3 + 15q)/8 q1,5 = ∓⅓5 + 2107 = ∓0.906180 w1,5 = (322 − 1370)/900 = 0.236927
q2,4 = ∓⅓5 − 2107 = ∓0.538469 w2,4 = (322 + 1370)/900 = 0.478629
q3 = 0 w3 = 128/225 = 0.568889
P6(q) = (231q6 − 315q4 + 105q2 − 5)/16 q1,6 = ∓0.932470 w1,6 = 0.171324
q2,5 = ∓0.661209 w2,5 = 0.360762
q3,4 = ∓0.238619 w3,4 = 0.467914
The floor functionn/2⌋ specifies the largest integer less than or equal to n/2.

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

 r = b − a / 2 q + a + b / 2 ≡ r(q)   ⇒   dr = b − a / 2 dq (3.2)

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

 b ∫ a f(r) dr = +1 ∫ −1 b − a / 2 f(r) dq ≈  n ∑ i=1 b − a / 2 wi f(ri) = n ∑ i=1 νi f(ri) (3.3)

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 / 2 wi   are the transformed weights.

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

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

 q = cosθ   ⇒   dq = − sinθ dθ (3.4)

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

 π ∫ 0 f(θ) sinθ dθ = − −1 ∫ 1 f(q) dq = 1 ∫ −1 f(q) dq  ≈  n ∑ j=1 wj f(qj) (3.5)

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

 Becke:       νi = 2 / 1 − qi2 (ri − r0) ri2 wi (3.6) Ahlrichs:       νi = (ri − r0) ri2 ( α / 1 + qi − 1 / 1 − qi ln−1 1 − qi / 2) wi (3.7) Linear:       νi = b − a / 2 ri2 wi (3.8)

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.

i Gauss-Legendre Becke(a) Ahlrichs(a) Linear(b) i Gauss-Legendre Becke(b) Ahlrichs(b) Linear(c) 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. 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 eq f(q) over the interval [0, +∞] as a weighted sum of function values at specified points within the domain of integration:

 +∞ ∫ 0 e−q f(q) dq ≈  n ∑ i=1 wi 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.
Ln(q) = n k=0 ( n   k ) (−1)k / k! qk qi wi = qi / (n + 1)2 [ Ln+1(qi) ]2
L0(q) = 1
L1(q) = −q + 1 q1 = 1w1 = 1
L2(q) = ½(q2 − 4q + 2) q1 = 2 − 2 = 0.585786 w1 = ¼(2 + 2) = 0.853553
q2 = 2 + 2 = 3.41421 w2 = ¼(2 − 2) = 0.146447
L3(q) = ⅙(−q3 + 9q2 − 18q + 6) q1 = 0.415775 w1 = 0.711093
q2 = 2.29428 w2 = 0.278518
q3 = 6.28995 w3 = 0.010389
L4(q) = 1 / 24 (q4 − 16q3 + 72q2 − 96q + 24) q1 = 0.322548 w1 = 0.603154
q2 = 1.74576 w2 = 0.357419
q3 = 4.53662 w3 = 0.038888
q4 = 9.39507 w4 = 0.000539
L5(q) = 1 / 120 (−q5 + 25q4 − 200q3 + 600q2 − 600q + 120) q1 = 0.26356 w1 = 0.521756
q2 = 1.4134 w2 = 0.398667
q3 = 3.59643 w3 = 0.075942
q4 = 7.08581 w4 = 0.003612
q5 = 12.6408 w5 = 0.000023
L6(q) = 1 / 720 (q6 − 36q5 + 450q4 − 2400q3 + 5400q2 − 4320q + 720) q1 = 0.222847 w1 = 0.458965
q2 = 1.188932 w2 = 0.417001
q3 = 2.992736 w3 = 0.113373
q4 = 5.775144 w4 = 0.010399
q5 = 9.837467 w5 = 0.000261
q6 = 15.98287 w6 = 0.000001

Radial integral. PAMoC can use the Gauss-Laguerre rule to estimate the radial integral (2.11), provided that the keyword "radrule" is set to 47. Expressing the radial coordinate as a linear function of an adjustable scale factor R,

 r = r0 + R q ≡ r(q)    ⇒    dr = R dq, (4.2)

yields:

 +∞ ∫ 0 r2 g(r) dr = +∞ ∫ 0 R r2 e−q e+q g[r(q)] dq = +∞ ∫ 0 e−q G(q) dq ≈≈ n ∑ i=1 wi G(qi) = n ∑ i=1 wi R ri2 eqi g(ri) = n ∑ i=1 νi g(ri) (4.3)

where  G(q) = R r2 eq g(r),  and ri and vi define the Laguerre grid:

 ri = r0 + R qi   νi = R ri2 eqi wi = R ri2 qi  eqi / (n + 1)2 [Ln+1(qi)]2 (4.4) 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−1 r−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α eq f(q) over the interval [0, +∞] and for some real number α > −1 as a weighted sum of function values at specified points within the domain of integration:

 +∞ ∫ 0 qα e−q f(q) dq ≈  n ∑ i=1 wi 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
Ln(α)(q) = n k=0 ( n + α   n − k ) (−1)k / k! qk qi
L0(α)(q) = 1
L1(α)(q) = −q + α + 1 q1 = α + 1
L2(α)(q) = ½[q2 − (α + 2)(2qα − 1)] q1 = (α + 2) − α + 2
q2 = (α + 2) + α + 2
L3(α)(q) = ⅙{−q3 + (α + 3) [3q2 − 3(α + 2)q + (α + 1)(α + 2)]}

Radial integral.   PAMoC can use the Generalized Gauss-Laguerre rule to estimate the radial integral (2.11), provided that the keyword "radrule" is set to 57. Expressing the radial coordinate as a function of an adjustable scale factor R, as in eq. (4.2), yields:

 +∞ ∫ 0 r2 g(r) dr = +∞ ∫ 0 R r2 q−α qα e−q e+q g(r) dq = +∞ ∫ 0 qα e−q G(q) dq ≈≈ n ∑ i=1 wi G(qi) = n ∑ i=1 wi R ri2 qi−α eqi g(ri) = n ∑ i=1 νi g(ri) (5.2)

where  G(q) = R r2 q−α eq g(ri), and ri and νi define the generalized Laguerre grid:

 ri = r0 + R qi   νi = R ri2 qi−α eqi wi (5.3)

It is worth noting that, for r0 = 0  and  α = 2, the expression of the quadrature weights simplifies to νi = R3 eqi wi. In this case, the generalized Laguerre grid (5.3) is exact if g(r) = R−3 exp(−r/R) L2n−1(2)(r/R).

Gauss-Laguerre Generalized Gauss-Laguerre Gauss-Laguerre Generalized Gauss-Laguerre qi andri (σ = 1) σ = 1⁄q6 σ = 2⁄q6 qi andri (σ = 1) σ = 1⁄q6 σ = 2⁄q6 wi σ = 1 σ = σ = 1⁄q6 2⁄q6 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. 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).

## 6. Gauss-Hermite quadrature

Gauss-Hermite quadrature provides an approximate solution to the integral of a function e-q² f(q) over the interval [−∞, +∞] as a weighted sum of function values at specified points

 +∞ ∫ −∞ e-q² f(q) dq ≈  n ∑ i=1 wi 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.
Hn(q) = (−1)n eq²  dn / dqn eq² =
= n! n/2⌋ k=0 (−1)k (2 q)n−2k / k! (n − 2k)!
qi wi = 2n−1 n! π / n2 [ Hn−1(qi) ]2
H0(q) = 1
H1(q) = 2q q1 = 0 w1 = √π = 1.772454
H2(q) = 4q2 − 2 q1,2 = ∓√½ = ∓0.707107 w1,2 = ½√π = 0.886002
H3(q) = 8q3 − 12q q1,3 = ∓√(32) = ∓1.224745 w1,3 = 0.295409
q2 = 0 w2 = 1.181636
H4(q) = 16q4 − 48q2 + 12 q1,3 = ∓√[(3+√6)/2] = ∓1.650680 w1,3 = 0.081313
q2,4 = ∓√[(3-√6)/2] = ∓0.524648 w2,4 = 0.804914

Multipole moments of the overlap density between two GTO's. Stone pointed out that the unabridged Cartesian multipole moments of the overlap density between two cartesian Gaussian type orbitals (GTO's) can be efficiently and accurately estimated using Gauss-Hermite quadrature.[22Stone, A. J.
J. Chem. Theory Comput. 2005, 1, 1128-1132.
] The following integrals are needed:

 Sμ∊A,ν∊B(l,m,n)(A) = <χμA|xAl yAm zAn|χνB> = +∞ ∫ −∞ xAl yAm zAn χ*μA(r) χνB(r) d3r, (6.2))

where

 χμA(r) = N1 xAl1 yAm1 zAn1 exp(-ζ1 rA2)   χνB(r) = N2 xBl2 yBm2 zBn2 exp(-ζ2 rB2) (6.3)

rA = rA and rB = rB 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. A 1950, 200, 542-554.
] showed that the product of two 1s GTO's, one centered at A and the other centered at B, is another Gaussian centered at P = ζ1A + ζ2B / ζ1 + ζ2 with exponent ζ3 = ζ1 + ζ2:

 exp(−ζ1 rA2) ⋅ exp(−ζ2 rB2) = exp[− ζ1 ζ2 / ζ3 (A − B)2] ⋅ exp(− ζ3 rP2) ≡ KAB exp(− ζ3 rP2), (6.4)

where KAB = exp[ ζ1 ζ2 / ζ3 (AB)2] and rP = rP.  Since 1s GTO's are separable along the three Cartesian coordinates:

 exp(− ζ3 rP2) = exp(− ζ3 xP2) ⋅  exp(− ζ3 yP2) ⋅  exp(− ζ3 zP2), (6.5)

the integral  Sμ∊A,ν∊B(l,m,n)(A)  can be written as:

 Sμ∊A,ν∊B(l,m,n)(A) = N1 N2 KAB  ∏ (t; k, k1, k2) +∞ ∫ −∞ GA(t; k, k1, k2) dt, (6.6)

where

 GA(t; k, k1, k2) = (t − At)k+k1 (t − Bt)k2 exp[−ζ3(t − Pt)2 = tAk+k1 tBk2 exp[−ζ3tP2], (6.7)

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

 q = √ζ3 (t − Pt)   ⇔   t = q⁄√ζ3 + Pt   ⇒   dt = (1⁄√ζ3) dq, (6.8)

which, coupled with the integration by substitution, yields for each integral in eq. (6.6):

 +∞ ∫ −∞ GA(t) dt  =   +∞ ∫ −∞ (1⁄√ζ3) (q⁄√ζ3 + Pt - At)k+k1 (q⁄√ζ3 + Pt - Bt)k2 exp(−q2) dq   = =   +∞ ∫ −∞ FA(q) exp(−q2) dq   ≈   n ∑ i=1 wi FA(qi)   = = n ∑ i=1 (wi⁄√ζ3) (qi⁄√ζ3 + Pt - At)k+k1 (qi⁄√ζ3 + Pt - Bt)k2  =  n ∑ i=1 (wi⁄√ζ3) (ti − At)k+k1 (ti − Bt)k2 (6.9)

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

 Sμ∊A,ν∊A(l,m,n)(A) = <χμA|xAl yAm zAn|χνA> = N1 N2  ∏ t=x,y,z +∞ ∫ −∞ tAk+k1+k2 exp[-ζ3tA2] dt, (6.10)

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

 +∞ ∫ −∞ tAk+k1+k2 exp[-ζ3tA2] dt = ( π / ζ3 )½ (k+k1+k2−1)!! / (2 ζ3)(k+k1+k2)/2 , (6.11)

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:

 1⁄2 +∞ ∫ −∞ r2 g(r) dr = +∞ ∫ −∞ 1⁄2 R r2 e−q2 e+q2 g(r) dq = +∞ ∫ −∞ e−q2 G(q) dq  ≈≈ n ∑ i=1 wi G(qi) = n ∑ i=1 1⁄2  wi R ri2 eqi2 g(ri) = n ∑ i=1 νi g(ri) (6.13)

where  G(q) = 12 R r2 eq2 g(r), and ri and νi define the Hermite grid:

 ri = r0 + R qi   νi = 1⁄2 R ri2 eqi2 wi (6.14)

It is worth noting that, for r0 = 0,   the expression of the quadrature weights simplifies to νi = 12 R3 eqi2 wi. 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 − q2 f(q) over the interval [−1, +1] as a weighted sum of function values at specified points within the domain of integration:

 +1 ∫ −1 √1 − q2 f(q) dq  ≈  n ∑ i=1 wi 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.
Un(q) = n/2⌋ k=0(-1)k ( n − k   k ) (2q)n−2k qi = cos ( i / n + 1π ) wi = π / n + 1 sin2 ( i / n + 1π )
= π / n + 1 (1 − qi2)
U0(q) = 1
U1(q) = 2 q q1 = 0 w1 = π/2
U2(q) = 4 q2 − 1 q1,2 = ±12 w1,2 = π/4
U3(q) = 8 q3 − 4 q q2 = 0
q1,3 = ±1/2
w2 = π/4
w1,3 = π/8
U4(q) = 16 q4 − 12 q2 + 1 q1,2,3,4 = ±(3 ± 5)/8 w1,2 = π(5 − 5)/40
w3,4 = π(5 + 5)/40
U5(q) = 32 q5 − 32 q3 + 6 q q3 = 0
q1,5 = ±0.866025
q2,4 = ±0.5
w3 = 0.523599
w1,5 = 0.130900
w2,4 = 0.392699
U6(q) = 64 q6 − 80 q4 + 24 q2 − 1 q1,6 = ±0.900969
q2,5 = ±0.623490
q3,4 = ±0.222521
w1,6 = 0.084489
w2,5 = 0.274333
w3,4 = 0.426576
U7(q) = 128 q7 − 192 q5 + 80 q3 − 8 q q4 = 0
q1,7 = ±0.923880
q2,6 = ±0.707107
q3,5 = ±0.382683
w4 = 0.392699
w1,7 = 0.057509
w2,6 = 0.196350
w3,5 = 0.335190
U8(q) = 256 q8 − 448 q6 + 240 q4 − 40 q2 + 1 q1,8 = ±0.939693
q2,7 = ±0.766044
q3,6 = ±0.5
q4,5 = ±0.173648
w1,8 = 0.040833
w2,7 = 0.144226
w3,6 = 0.261799
w4,5 = 0.338540
U9(q) = 512 q9 − 1024 q7 + 672 q5 − 160 q3 + 10 q q5 = 0
q1,9 = ±0.951057
q2,8 = ±0.809017
q3,7 = ±0.587785
q4,6 = ±0.309017
w5 = 0.314159
w1,9 = 0.029999
w2,8 = 0.108539
w3,7 = 0.205620
w4,6 = 0.284160
The floor functionn/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 ∫ −1 Un(q) Um(q) √1 − q2 dq  =  π / 2 δnm (8.2)

As a consequence:

 n ∑ i=1 wi =  π / 2 (8.3)

Radial integral. The Gauss-Chebyshev second kind quadrature can be used to estimate the radial integral (2.11), provided that a suitable variable transformation is given to map the improper interval [0, ∞] to a finite interval [−1, +1].

Becke's grid. In 1988, Becke introduced the transformation:[2Becke, A. D.
J. Chem. Phys. 1988, 88, 2547-2553.
]

 r = r0 + R 1 + 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:

 +∞ ∫ 0 r2 g(r) dr = +1 ∫ −1 2  r2 (r − r0) / 1 − q2 g(r) dq = +1 ∫ −1 √1 − q2 G(q) dq ≈≈ n ∑ i=1 wi G(qi) = n ∑ i=1 2 wi ri2 (ri − r0) / (1 − qi2)3/2 g(r) = n ∑ i=1 νi g(ri) (8.5)

where  G(q) = 2  r2 (rr0) / (1 − q2)3/2 g(r), and ri and νi define the Becke grid:

 ri = r0 + R 1 + qi / 1 − qi   νi = 2 ri2 (ri − r0) / (1 − qi2)3/2 wi = 2π / n + 1 ri2 (ri − r0) / √1 − qi2   qi = 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)3 U2n−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.

Ahlrichs's grid. In 1995, Treutler and Ahlrichs introduced the transformation:[7Treutler, O.; Ahlrichs, R.
J. Chem. Phys. 1995, 102, 346-354.
]

 r = r0 − R / ln 2 (1 + q)α ln 1 − q / 2 ≡ r(q)   ⇒    dr = (r − r0) ( α / 1 + q − 1 / 1 − q ln−1 1 − q / 2 ) dq (8.8)

which is normalized to map the center of the q interval onto r(0) - r = R.   At q = −1,   r(q) − r0 = 0,   and as q → 1, the sample points tend toward infinity. Under this map, the radial integral becomes:

 +∞ ∫ 0 r2 g(r) dr = +1 ∫ −1 r2 (r − r0) ( α / 1 + q − 1 / 1 − q ln−1 1 − q / 2 ) √1 − q2 / √1 − q2 g(r) dq = +1 ∫ −1 √1 − q2 G(q) dq≈ n ∑ i=1 wi G(qi) = n ∑ i=1 π / n + 1 (1 − qi2) G(qi)= n ∑ i=1 π / n + 1 ri2 (ri − r0) ( α√ 1 − qi / 1 + qi − √ 1 + qi / 1 − qi ln−1 1 − qi / 2 ) g(r) = n ∑ i=1 νi g(ri) (8.9)

where ri and νi, with  α = 0.6,  define the Ahlrichs grid:

 ri = r0 − R / ln 2 (1 + qi)α ln 1 − qi / 2   νi = ri2 (ri − r0) ( α / 1 + qi − 1 / 1 − qi ln−1 1 − qi / 2 ) wi / √1 − qi2 = π / n + 1 ri2 (ri − r0) [ α √ 1 − qi / 1 + qi − √ 1 + qi / 1 − qi ln−1 1 − qi / 2 ]   qi = cos ( i / n + 1 π)          ∀ i = 1, 2, ..., n (8.10)

Linear grid. Although it is rather inappropriate for radial integration, the linear transformation (3.2) can also be used in combination with the 2nd kind Gauss-Chebyshev quadrature rule, yielding the following quadrature points and weights:

 ri = b − a / 2 qi + a + b / 2 νi = b − a / 2 √1 − q2 ri2 wi = π (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.

i 2nd kindGauss-Chebyshev Becke(a) Ahlrichs(a) Linear(b) i 2nd kindGauss-Chebyshev Becke(b) Ahlrichs(b) Linear(c) 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. 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 ln2 q f(q) over the interval [0, 1] as a weighted sum of function values at specified points within the domain of integration:[10Gill, P. M. W.; Chien, S.-H.
J. Comput. Chem. 2003, 24, 732-740.
]

 1 ∫ 0 ln2 q f(q) dq ≈  n ∑ i=1 wi 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 ln2 q, and wi are the associated weights.

Table 9.1 − First three Gill polynomials with associated roots and weights
nQn(q) qiwi
18q − 1 q1 = 1/8 w1 = 2.000000000000000
27992q2 − 4104q + 217 q1 = 0.059850992523974
q2 = 0.453662520989539
w1 = 1.669136108179106
w2 = 0.330863891820894
3 q1 = 0.036263311146964
q2 = 0.273148602374171
q3 = 0.653711089636059
w1 = 1.363830383647107
w2 = 0.565815459643824
w3 = 0.070354156709070

Roots and weights for n = 1-26, 30, 35, 40, 45, 50 can be obtained from Gill's website at the Australian National University (ANU). Their values for n = 1-3 are reported in Table 9.1

MultiExp grid for radial quadrature. The Gauss-Gill "log-squared" quadrature rule (9.1) was devised to estimate the radial integral (2.11), provided that a suitable variable transformation is given to map the improper interval [0, ∞] to a finite interval [0,1]. Gill and Chien proposed the transformation:[10Gill, P. M. W.; Chien, S.-H.
J. Comput. Chem. 2003, 24, 732-740.
]

 r = 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:

 +∞ ∫ 0 r2 g(r) dr = R 1 ∫ 0 r2 q−1 ln−2q ln2q g(r) dq = 1 ∫ 0 ln2q G(q) dq  ≈≈  n ∑ i=1 wi G(qi) = R n ∑ i=1 r2 qi−1 ln−2qi wi g(ri) = n ∑ i=1 νi g(ri) (9.3)

where

 G(q) = R (r0 − R ln q)2 q−1 ln−2q g(r) = R3 r2 (r − r0)−2 e(r−r0)/R g(r)

and ri and νi define the MultiExp grid:

 ri = r0 − R ln qi   νi = R (r0 − R ln qi)2 qi−1 ln−2 qi wi = R3 r2 (r − r0)−2 e(r−r0)/R wi (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 = 12, 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

 g(r) = R−3 r−2 (r − r0)2 e−(r−r0)/R Q2n−1(e−(r−r0)/R) (9.5)

For r0 = 0, eq. (9.5) can be written as

 g(r) = R−3 2n−1 ∑ k=0 ck e−(k+1)r/R (9.6)

Obviously, such g(r) are linear combinations of exponential functions, and for this reason Gill and Chien named "MultiExp" the grid defined by eq. (9.4). These authors concluded that "The accuracy with which a quadrature scheme estimates the value of an integral depends on how close the integrand is to the set of functions for which the quadrature is exact. Our principal interest lies in the development of efficient quadrature schemes for the types of integral that arise in atomic and molecular calculations using density functional theory. Because of shell structure, the radial electron density around an atom can usually be modeled well by a sum of exponential functions with a variety of decay constants. It has previously been noted[3Murray, C. W.; Handy, N. C.; Laming, G. J.
Mol. Phys. 1993, 78, 997-1014.
; 7Treutler, O.; Ahlrichs, R.
J. Chem. Phys. 1995, 102, 346-354.
] that such multiexponential character may limit the usefulness of the Laguerre grid, but it is clear that the MultiExp grid is well suited to such problems."[10Gill, P. M. W.; Chien, S.-H.
J. Comput. Chem. 2003, 24, 732-740.
]

Gauss-Gill quadrature with Knowles, Handy, and linear grids.  The Gauss-Gill quadrature can also be used in conjunction with any variable transformation which maps the semi-infinite interval spanned by the radial coordinate r into the finite interval [0, 1] of the variable q. PAMoC can use the transformations due to Knowles[8Mura, M. E.; Knowles, P. J.
J. Chem. Phys. 1996, 104, 9848-9858.
] and Handy,[3Murray, C. W.; Handy, N. C.; Laming, G. J.
Mol. Phys. 1993, 78, 997-1014.
] which are described in section 10 (eqs. 10.14 and 10.10, respectively) with more details, and of course some linear transformation like eq. (10.23). The appropriate weights for the Gauss-Gill quadrature are given by equations from (9.7) through (9.10), respectively,

 Knowles:       νi = k R ri2 qik−1 / (1 − qik) ln2 qi wi (9.7) Handy:       νi = m (ri − r0) ri2 / qi (1 − qi) ln2 qi wi (9.8) modified Handy:       νi = m (rmax − r0) ri2 [1 + (rmax − r0 − 2m)(1 − qi)m−1] qim−1 / [1 + (rmax − r0 − 2m) (1 − qi)m]2 ln2 qi wi (9.9) Linear:       νi = (rmax − r0) ri2 / ln2 qi wi (9.10)

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)
10.04550.05170.06560.0000 0.0000 0.0000 0.0000 0.00010.00000.0530
20.12370.14070.17850.0000 0.0007 0.0004 0.0015 0.00290.00210.3672
30.24020.27320.34650.0009 0.0124 0.0069 0.0116 0.02310.01610.9732
40.39950.45440.57630.0063 0.0846 0.0470 0.0511 0.10140.06801.8430
50.60880.69240.87830.0254 0.3416 0.1898 0.1710 0.33940.21372.9252
60.87921.00001.26850.0742 1.0000 0.5558 0.5037 1.00000.56454.1510
71.22921.39811.77340.1756 2.3655 1.3147 1.4234 2.82601.31685.4401
81.69121.92352.43990.3590 4.8379 2.6888 4.1468 8.23312.72476.7066
92.32972.64973.36110.6665 8.9806 4.9912 13.5684 26.93874.85707.8649
103.30443.75824.76721.171015.7793 8.7697 57.6650114.48827.22148.8364
115.24065.96047.56062.059327.747915.4216461.8325916.92359.02359.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 = 12. (c) The first column refers to the grid standardized by means of the central point of the interval qcenter = 12, 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.740 0 0 0 0 0 (-4)4.69
2(-3)1.48(-3)2.18(-1)4.450 0 0 (-8)1 (-8)6 (-8)2 (-2)6.24
3(-3)7.90(-2)1.16(-1)2.370 (-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 = 12. (c) The first column refers to the grid standardized by means of the central point of the interval qcenter = 12, 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 = ba:

 b ∫ a f(q)dq ≈ b − a / 2 [ f(a) + f(b) ] (10.1) Figure 10.1: Trapezoidal rule (on the left) and 4-point closed extended trapezoidal rule (on the right).

As shown in the right side of Figure 10.1, a better approxination is obtained by dividing the interval into a number of subintervals of equal length. Then the trapezoidal rule is used to estimate the area under the curve in each subinterval, and the results for all subintervals are summed together. For n equispaced points, using the trapezoidal rule  (n - 1)  times and adding the results gives the "extended" trapezoidal rule:

 b = qn ∫ a = q1 f(q) dq  =  n−1 ∑ i=1 qi+1 ∫ qi f(q) dq  ≈  1 / 2 n−1 ∑ i=1 (qi+1 − qi) (fi + fi+1)  =  h ( 1 / 2 f1 + n−1 ∑ i=2 fi + 1 / 2 fn )  =  n ∑ i=1 wi f(qi) (10.2)

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−1 f(q) dq ≈ h f(qn−1) ≡ h fn−1 (10.4)

and Eq. (10.2) yields the n-point semi-open trapezoidal formula:

 b = qn ∫ a = q1 f(q) dq ≈ h ( 1 / 2 f1 + n−2 ∑ i=2 fi + 3 / 2 fn−1 ) (10.5)

which is de facto reduced to a (n − 1)-point formula. In order that the semi-open trapezoidal formula has n evaluation points effectively, a (n + 1)-point formula is needed with n equal intervals of length   h = (b − a)/n,   and weights   w1 = 12h,   w2 = w3 = … = wn−1 = h,   and   wn = 32h.

The n-point trapezoid rule integrates exactly all spherical harmonics of degree L = n − 1 or less, so that, for a given L, a number of n = (L + 1) integration points is needed.

The following paragraphs describe two main applications of the extended trapezoidal rule in PAMoC.

φ quadrature over the interval [0, 2π]. In the case of explicit user request to estimate the angular integral by the two-term (Nθ, Nφ) product formula ("intgrid 8 | 9 | 10" and "intgrid -mmmlll"), PAMoC uses the closed extended trapezoidal rule for the φ quadrature over the interval [0, 2π]. In this case, using (n + 1) equispaced points, from φ1 = 0 to φn+1 = 2π with step-size h = 2π/n, and observing that f(0) = f(2π), we have:

 2π ∫ 0 f(φ)dφ ≈ 2π / n n ∑ i=1 f(φi), (10.6)

where    φi = (i − 1) (2π/n),   ∀ i = 1, 2, …, n.

Radial integral. The trapezoidal rule is also used by PAMoC to estimate the radial integral (2.11), in conjunction with a suitable variable transformation to map the improper interval [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.  rr(q),  so that  dr = r'(q) dq,  then the radial integral (2.11) can be solved by substitution, yielding

 +∞ ∫ 0 r2 g(r) dr  =  b ∫ a r2(q) r'(q) g(r) dq  ≡  b ∫ a G(q) dq (10.7)

where the integrand is now given by  G(q) = r2(qr'(q) g(r).  Applying the quadrature rule (10.2), yields

 b ∫ a G(q) dq  ≈  n ∑ i=1 wi G(qi)  =  n ∑ i=1 νi g(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).

Handy's grid. In 1993, Murray, Handy and Laming introduced the transformation:[3Murray, C. W.; Handy, N. C.; Laming, G. J.
Mol. Phys. 1993, 78, 997-1014.
]

 r = r0 + R qm / (1 − q)m    ⇒    dr = m R qm−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

 ri = r0 + R qim / (1 − qi)m = r0 + R im / (n − i + 1)m νi = wi m / qi (1 − qi) (ri − r0) ri2 = wi m (n + 1)2 / i (n − i + 1) (ri − r0) ri2          ∀ i = 1, 2, …, n (10.12) qi = i h,     w1 = w2 = … = wn − 1 = h,     wn = 3 / 2 h,     h = 1 / n + 1,     m ≥ 2

At the grid mid-point, if n is odd, we have i = (n + 1)/2, q(n + 1)/2 = 12, 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 R r02. 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 = 12 hw2 = w3 = … = wn − 1 = hwn = 32 h, with h = 1n.

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 R3 q3m−1 / (1 − q)3m+1 g[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,

 g(r) = m√R / m r(3m−1)/m ( m√R + m√r )2 L ( m√r / m√R + m√r) ) (10.13)

In general, such g(r) possess a 1/r(3m−2)/m singularity at the origin and decay as 1/r(3m+1)/m for large r.[10Gill, P. M. W.; Chien, S.-H.
J. Comput. Chem. 2003, 24, 732-740.
]

In PAMoC, radial integration by the semi-open extended trapezoidal rule, combined with the Handy grid with m = 2, can be selected by the keyword "radrule 73". This is the only option available in all versions of PAMoC before 2013.

Knowles' grid. In 1996, Mura and Knowles introduced the transformation:[8Mura, M. E.; Knowles, P. J.
J. Chem. Phys. 1996, 104, 9848-9858.
]

 r = r0 − R ln (1 − qk)    ⇒    dr = k R qk−1 / 1 − qk dq (10.14)

which, at both the origin and the endpoint, behaves as the Handy transformation. Therefore, sampling q over (n + 2) equispaced points according to eq. (10.11), leads to a (n + 2)-point semi-open quadrature formula, which de facto reduces to a n-point formula like (10.10), whose evaluation points and weights define the Knowles grid:

 ri = r0 − R ln (1 − qik) νi = k R qik−1 / 1 − qik wi ri2                            ∀ i = 1, 2, …, n (10.15) qi = i h,     w1 = w2 = … = wn − 1 = h,     wn = 3 / 2 h,     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 32.  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 = 12 hw2 = w3 = … = wn − 1 = hwn = 32 h, with h = 1n.

At the grid mid-point, if n is odd, we have i = (n + 1)/2, q(n + 1)/2 = 12, and r(n + 1)/2 - r0 / R = ln 2k / 2k − 1 = 1 / σk. Here, σk = ln−1 2k / 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(43) = 3.476059, σ3 = ln−1(87) = 7.488876, and σ4 = ln−1(1615) = 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) = k qk−1 ln2(1 − qk) / 1 − qk g[r(q)]  equals  L(q), which implies

 g(r) = exp(−r) / k r2 [1 − exp(−r)](k−1)/k L( k√1 − exp(−r) ) (10.16)

In general, such g(r) possess a 1/r(3k−2)/k singularity at the origin and decay exponentially.[10Gill, P. M. W.; Chien, S.-H.
J. Comput. Chem. 2003, 24, 732-740.
]

In PAMoC, radial integration by the semi-open extended trapezoidal rule, combined with Knowles' Log3 grid, is selected by the keyword "radrule 72".

Gill and Chien transformation of the radial coordinate. The (n+1)-point semi-open extended trapezoid rule can be used in combination with the transformation (9.2), introduced by Gill and Chien,[10Gill, P. M. W.; Chien, S.-H.
J. Comput. Chem. 2003,24, 732-740.
] to define the grid:

 ri = r0 − R ln qi νi = − R ri2 wi / qi                                                                 ∀ i = 1, 2, …, n (10.17) qi = (n − i + 1) h,     h = 1 / n     and w1 = 1 / 2 h,     w2 = w3 = … = wn − 1 = h,     wn = 3 / 2 h,     r0 ≠ 0.

Indeed, this transformation of the radial coordinate has been thought to be used in conjunction with the Gauss-Gill “log-squared” quadrature rule (9.1).[10Gill, P. M. W.; Chien, S.-H.
J. Comput. Chem. 2003,24, 732-740.
]   However, nothing prevents to apply it to a set of evenly spaced points.

For i = n + 1 (i.e.  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 32. 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 = (ni + 1)/(n + 1)   ∀ i = 1, 2, …, n, and the associated weights as   w1 = w2 = … = wn−1 = hwn = 32 h, with h = 1/(n + 1).

At the center of the interval, it is q = 12i = (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 + R 1 + qi / 1 − qi = r0 + R i − 1 / n − i + 1 νi = 2 R ri2 / (1 − qi)2 wi                                                 ∀ i = 1, 2, …, n (10.18) qi = 2 i − n − 2 / n,     h = 2 / n     and w1 = 1 / 2 h,     w2 = w3 = … = wn − 1 = h,     wn = 3 / 2 h,     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 32. 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 in − 1)/(n + 1)   ∀ i = 1, 2, …, n, and the associated weights as   w1 = w2 = … = wn−1 = hwn = 32 h, 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:

 ri = r0 − R / ln 2 (1 + qi)α ln 1 − qi / 2 νi = ri2 R / ln 2 (1 + qi)α−1 ( 1 + qi / 1 − qi − α ln 1 − qi / 2) wi                 ∀ i = 1, 2, …, n (10.19) qi = 2 i − n − 1 / n + 1,     h = 2 / n + 1     and w1 = w2 = … = wn − 1 = h,     wn = 3 / 2 h

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 32. 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".

Modified Handy's grid for finite intervals of the radial coordinate r. Figure 10.4 shows that Handy's grid comprises a number of points which are very far from the nucleus, and which probably contribute very little to the integral of typically exponentially (or Gaussian) decaying integrands. Such points can, of course, be systematically neglected within computer programs by dynamic pruning of the grid, which is reduced by about 20% of its original size.[3Murray, C. W.; Handy, N. C.; Laming, G. J.
Mol. Phys. 1993, 78, 997-1014.
, 4Gill, P. M. W.; Johnson, B. G.; Pople, J. A.
Chem. Phys. Lett. 1993, 209, 506-512.
, 10Gill, P. M. W.; Chien, S.-H.
J. Comput. Chem. 2003, 24, 732-740.
] On this line, PAMoC truncates the summation in the right hand side of eq. (10.10) at a value Nr < n such that qNrrmax, where rmax is large enough (usually 10 bohr) to be a good approximation to infinity. However. Mura and Knowles pointed out that it is preferable to use more appropriate grids from the outset[8Mura, M. E.; Knowles, P. J.
J. Chem. Phys. 1996, 104, 9848-9858.
] and indeed the Knowles grid is relatively compact, as shown in Figure 10.4. Therefore, a modified form of the Murray, Handy and Laming transformation[3Murray, C. W.; Handy, N. C.; Laming, G. J.
Mol. Phys. 1993, 78, 997-1014.
] is made available in PAMoC to map the interval [r0, rmax] to the interval [0, 1]:

 r(q) = r0 + (rmax − r0) qm / 1 + (rmax − r0 − 2m)(1 − q)m   ⇒   dr = m (rmax − r0) [1 + (rmax − r0 − 2m)(1 − q)m−1] qm−1 / [1 + (rmax − r0 − 2m) (1 − q)m]2  dq (10.20)

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 (rmaxr0), and G(1) = r2(1) r'(1) g[r(1)] = m rmax2 (rmaxr0) g(rmax). At the grid midpoint it is q = 12 and r(12) = 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:

 rmax ∫ r0 r2 g(r) dr = 1 ∫ 0 r2(q) r'(q) g[r(q)] dq ≡ 1 ∫ 0 G(q) dq ≈ 1 / n − 1 ( 1 / 2 G1 + n−1 ∑ i=2 Gi + 1 / 2 Gn ) = n ∑ i=1 wi Gi = n ∑ i=1 νi g(ri) (10.21)

The evalution points ri and their associated weights νi and wi define the modified Handy grid for the finite interval [r0, rmax]:

 ri = r0 + (rmax − r0) qim / 1 + (rmax − r0 − 2m)(1 − qi)m ν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 / 2 h,     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 (rmaxr0) r02 / (rmaxr0 − 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 = 12 hw2 = w3 = … = wn−1 = hh = 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 / 2 h,     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 = r02 g(r0).   At the endpoint it is qn = rmaxr0,   rn = rmax,   rn' = 1,  and  Gn = rmax2 g(rmax).   At the grid midpoint it is q(n+1)/2 = rmaxr0,  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)
10.08700.12550.0006 0.0043 0.0083 0.0139 0.09090.0428 0.9091
20.18230.26300.0046 0.0348 0.0400 0.0659 0.20000.1361 1.8182
30.28770.41500.0157 0.1179 0.1111 0.1782 0.33330.2738 2.7273
40.40550.58500.0377 0.2826 0.2500 0.3855 0.50000.4586 3.6364
50.53900.77760.0751 0.5623 0.5102 0.7418 0.71430.6970 4.5455
60.69311.00000.1335 1.0000 1.0000 1.3284 1.00001.0000 5.4545
70.87551.26300.2213 1.6570 1.9600 2.2581 1.40001.3854 6.3636
81.09861.58500.3514 2.6316 4.0000 3.6571 2.00001.8836 7.2727
91.38632.00000.5480 4.1036 9.0000 5.5862 3.00002.5508 8.1818
101.79182.58500.8644 6.4735 25.0000 7.8740 5.00003.5121 9.0909
112.48493.58501.470811.0145121.000010.000011.00005.157410.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.070 (-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).

Euler-Maclaurin quadrature. The extended trapezoid quadrature is closely related to the Euler-Maclaurin summation formula, which is an important tool in numerical analysis to estimate sums or integrals of functions. For an integral in the interval [a, b] and a step length h = (ba) / (n − 1), it takes the form:[16Smith, F. J.
Numer. Math. 1965, 7, 406-411.
, 17Graham, R. L.; Knuth, D. E.; Patashnik, O.
"Concrete Mathematics",
Second Edition, Addison-Wesley Publishing Company, Inc., 1994, eq. (9.67), p. 469.
]

 b = qn ∫ a = q1 f(q) dq = h ( 1 / 2 f1 + n−1 ∑ i=2 fi + 1 / 2 fn ) − m ∑ k=1 h2k B2k / (2k)! [ f (2k−1)(b) − f (2k−1)(a) ] + Em(f) (10.24)

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+2 B2m+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 = −12, B2 = 16, B4 = −130, B6 = 142, B8 = −130, ... ; 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

 ∞ ∫ 0 f(q) dq = h n ∑ i=1 fi (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

 Handy < modified Handy < Knowles' Log3 < Becke ≈ Linear < Gill-Chien < Ahlrichs

Adding a second, tighter Gaussian, worsens the linear grid accuracy, yielding the ranking

 Linear ≪ Handy < modified Handy < Becke ≈ Knowles' Log3 < Gill-Chien < Ahlrichs

Adding a third, even tighter Gaussian to g(r), leads to deterioration in all grid accuracies to yield

If the same test is done using an exponential function with unit exponent, the grids are ranked as

 Ahlrichs < Gill-Chien < linear < Becke ≈ Knowles' Log3 ≈ modified Handy < Handy

Adding one and two exponential functions with exponent increasingly higher makes the linear grid accuracy the worst of all and leaves the ranking of the remaining ones unchanged.

 Linear ≪ Knowles' Log3 ≈ Handy ≈ modified Handy < Ahlrichs ≈ Becke ≈ Gill-Chien

When the grids of Tables 10.2 and 10.3 are used to estimate the slowly converging radial integral of the function 1/(1 + r4) the results shown in the final row of Table 11.2 are obtained, with the grid ranking

 Gill-Chien < Ahlrichs ≈ modified Handy ≈ linear < Knowles' Log3 < Handy < Becke

favoring the grids which sample the integrand at large r.

Table 11.1. − Functions used to test radial quadrature rules.
ig(r)
1exp(−r2)
2j=0,1 10j exp(−10j r2)
3j=0,1,2 10j exp(−10j r2)
4exp(−r)
5j=0,1 102j exp(−10j r)
6j=0,1,2 102j exp(−10j r)
71/(1 + r4)
Table 11.2. − Accuracy(a) of radial quadrature using the grids of Tables 10.2 and 10.3.
g(r) Gill-Chien and
Knowles' Log1
σ = 1.4427
Knowles' Log3
σ = 7.4889
Handy modified
Handy
Becke Ahlrichs Linear
14.13.32.02.93.55.33.5
24.03.72.32.83.65.30.6
33.12.32.42.43.13.00.5
41.42.52.82.52.51.22.3
51.42.52.82.52.61.31.0
61.32.42.82.52.21.31.0
70.81.52.11.02.21.01.1
(a) As suggested by Gill and Chien,[10Gill, P. M. W.; Chien, S.-H.
J. Comput. Chem. 2003,24, 732-740.
] the quality of a quadrature estimation (“Approx”) of an integral is compared with the exact value (“Exact”) using the equation

Accuracy = −log10 | Approx / Exact − 1 |

The resulting values are essentially the number of correct digits in the quadrature result.

It is clear that, if r is null outside the interval [0,10], the Log1 grid should not be coupled with the Euler-Maclaurin quadrature because it covers less than half of the interval. On the other hand, the Handy grid should be discarded because two points (18%) fall outside the interval. Therefore the modified Handy grid and Knowles' Log3 grid remain as suitable choices for the Euler-Maclaurin radial quadrature.

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)
14.74.95.32.81.61.70.91.71.63.0
24.84.84.92.61.82.30.92.31.91.3
33.33.12.62.51.72.10.91.83.01.6
48.12.22.50.36.66.01.51.92.42.5
57.92.32.50.42.63.11.52.02.42.1
62.53.02.50.45.02.91.52.12.42.2
70.91.01.10.52.11.42.31.71.01.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.
i BeckeAhlrichsLinear BeckeAhlrichsLinear
12.23.42.62.33.73.6
22.33.51.32.43.21.2
32.34.01.32.52.41.4
42.83.92.62.53.52.5
52.93.92.02.53.62.4
63.52.91.72.52.82.3
73.71.21.02.61.11.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.

The quadrature methods examined so far use grid points and weights that are known beforehand. Other schemes have been designed, in which points and weights are chosen dynamically as the integrand is explored.[12(a) Pérez-Jordá, J. M.; San-Fabián, E.; Moscardó, F.
Comp. Phys. Commun. 1992, 70, 271-284.
(b) Pérez-Jordá, J. M.; San-Fabián, E.
Comp. Phys. Commun. 1993, 77, 46-56.
(c) Pérez-Jordá, J. M.; Becke A. D.; San-Fabián, E.
J. Chem. Phys. 1994, 100, 6520-6534.
; 13(a) Krack, M.; Köster, A. M.
J. Chem. Phys. 1998, 108, 3226-3234.
(b) Köster, A. M.; Flores-Moreno, R.; Reveles, J. U.
J. Chem. Phys. 2004, 121, 681-690.
; 14 Gräfenstein, J.; Cremer, D.
J. Chem. Phys. 2007 127, 164113-7.
; 15Kakhiani, K.; Tsereteli, K.; Tsereteli, P.
Computer Physics Communications 2009, 180, 256-268.
]

## References

1. "Über empirische Funktionen und die Interpolation zwischen äquidistanten Ordinaten"
Runge, C. Zeitschrift für Mathematik und Physik 1901, 46, 224-243.
2. "A multicenter numerical integration scheme for polyatomic molecules"
Becke, A. D. J. Chem. Phys. 1988, 88, 2547-2553.
3. "Quadrature schemes for integrals of density functional theory"
Murray, C. W.; Handy, N. C.; Laming, G. J. Mol. Phys. 1993, 78, 997-1014.
4. "A standard grid for density functional calculations"
Gill, P. M. W.; Johnson, B. G.; Pople, J. A. Chem. Phys. Lett. 1993, 209, 506-512.
5. (a) "Values of the nodes and weights of ninth to seventeenth order Gauss-Markov quadrature formulae invariant under

the octahedron group with inversion"

Lebedev, V.I. USSR Comp. Math. and Math. Phys. 1975, 15(1), 44-51. Zh. vychisl. Mat. mat. Fiz. 1975, 15(1), 48-54.

(b) "Quadratures on a sphere"

Lebedev, V.I. USSR Comp. Math. and Math. Phys. 1976, 16(2), 10-24. Zh. vychisl. Mat. mat. Fiz. 1976, 16(2), 293-306.

(c) "Spherical quadrature formulas exact to orders 25-29"

Lebedev, V.I. Siberian. Math. J. 1977, 18(1), 99-107. Sibirskii Matematicheskii Zhurnal 1977, 18(1), 132-142.

(d) "Quadrature formulas of orders 41, 47, and 53 for the sphere"

Lebedev, V. I.; Skorokhodov, A. L. Russian Acad. Sci. Dokl. Math. 1992, 45, 587-592.

(e) "A quadrature formula for the sphere of 59th algebraic order of accuracy"

Lebedev, V. I. Russian Acad. Sci. Dokl. Math. 1995, 50, 283-286.

(f) "A quadrature formula for the sphere of the 131-st algebraic order of accuracy"

Lebedev, V.I.; Laikov, D.N. Dokl. Math. 1999, 59, 477-481.

6. "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.
7. "Efficient molecular numerical integration schemes"
Treutler, O.; Ahlrichs, R. J. Chem. Phys. 1995, 102, 346-354.
8. "Improved radial grids for quadrature in molecular density-functional calculations"
Mura, M. E.; Knowles, P. J. J. Chem. Phys. 1996, 104, 9848-9858.
9. "Molecular integrals by numerical quadrature. I. Radial integration"
Lindh, R.; Malmqvist, P.-A.; Gagliardi, L. Theor. Chem. Acc. 2001, 106, 178-187.
Gill, P. M. W.; Chien, S.-H. J. Comput. Chem. 2003, 24, 732-740.
11. "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.
12. (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.

13. (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.

14. "Efficient density-functional theory integrations by locally augmented radial grids"
Gräfenstein, J.; Cremer, D. J. Chem. Phys. 2007 127, 164113-7.
15. "A program to generate a basis set adaptive radial quadrature grid for density functional theory"
Kakhiani, K.; Tsereteli, K.; Tsereteli, P. Computer Physics Communications 2009, 180, 256-268.
16. "The arrangement of atoms in crystals"
Bragg, W. L. Philos. Mag. 1920, 40, 169-189.
17. "Atomic Radii in Crystals"
Slater, J. C. J. Chem. Phys. 1964, 41, 3199-3204.
18. "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.
19. "Atomic Shielding Constants"
Slater, J. C. Phys. Rev. 1930, 36, 57-64.
20. "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.
21. "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.
22. "Distributed Multipole Analysis: Stability for Large Basis Sets"
Stone, A. J. J. Chem. Theory Comput. 2005, 1, 1128-1132.
23. "Electronic Wave Functions. I. A General Method of Calculation for the Stationary States of Any Molecular System"
Boys, S. F. Proc. Roy. Soc. A 1950, 200, 542-554.
24. "Gaussian quadrature formulae for ∫01−ln(x)f(x)dx"
Anderson, D.G. Math. Comp. 1965, 19, 477-481.
25. "Calculation of Gauss Quadrature Rules"
Golub, G.H.; Welsch, J.H. Math. Comp. 1969, 23, 221-230.
26. "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.
27. "Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of Interpolatory Quadrature"
Elhay, Sylvan; Kautsky, Jaroslav ACM Transactions on Mathematical Software 1987, 13, 399-415.
28. "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.
29. Graham, R. L.; Knuth, D. E.; Patashnik, O. "Concrete Mathematics", Second Edition, Addison-Wesley Publishing Company, Inc., 1994, eq. (9.67), p. 469.

## Notes

1. 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.

[Main source of this note: Wikipedia contributors. Jacobian matrix and determinant. Wikipedia, The Free Encyclopedia. Available at: https://en.wikipedia.org/w/index.php?title=Jacobian_matrix_and_determinant&oldid=834303707. Accessed April 13, 2018].