# Becke partitioning : Computing integrals with Lebedev and Gauss-Chebyshev quadrature

## Preface:

This question is a follow up post from a previous question on MMSE: Usual way of computing exchange-correlation potential on a simple case H2 (RHF, LDA).

The radial variable substitution and the complete method are detailed in the paper: A multicenter numerical integration scheme for polyatomic molecules - The Journal of Chemical Physics.

## Question:

In the details sub-heading, I calculate the difference between the actual value of the integral (I1) and the result obtained with quadrature (I2). This should ideally be zero. That is I1 - I2 = 0. However, I get values that are not 0.

So how do I properly integrate with the help of the Lebedev + Gauss-Chebyshev quadrature?

## Details:

I will detail the different parts of the code and give the full version at the end.

#### Calculations to get the quadrature points :

import numpy as np
import matplotlib.pyplot as plt
from numpy.linalg import norm
pi = np.pi

##########################
##########################

Nr = 30
ii  = np.arange(1,Nr+1).reshape(-1,1)
sin_i = np.sin(ii*pi/(Nr+1))
cos_i = np.cos(ii*pi/(Nr+1))

# Chebyshev quadrature points on [-1,1]
Mu_quadrature = (Nr+1.0-2.0*ii)/(Nr+1.0) + 2.0/pi*(1.0 + 2.0/3.0*sin_i**2)*cos_i*sin_i

# Chebyshev quadrature points on [0, +infty [
rm = 0.35 # Angstrom
# variable substitution for the radial component

##########################
# Quadrature for solid angle / angular coordinates
##########################

# file containing integration points (theta/phi in deg) and associated weights

# make a tensor product of the Quadrature points
X_spherical_int = np.concatenate([
],
axis=1
)

R_int   = X_spherical_int[:,0:1]
Theta_int = X_spherical_int[:,1:2]
Phi_int = X_spherical_int[:,2:3]

Mu_int  = (R_int - rm)/(R_int + rm)

Wint = np.concatenate([
],
axis=1
)
Wint = np.prod(Wint,axis=1,keepdims=True)

X_cartesian_int = np.zeros_like(X_spherical_int)
R, Theta, Phi = R_int, Theta_int, Phi_int
Mu = (R - rm)/(R + rm)
X_cartesian_int[:,0:1] = R*np.cos(Theta)*np.sin(Phi)
X_cartesian_int[:,1:2] = R*np.sin(Theta)*np.sin(Phi)
X_cartesian_int[:,2:3] = R*np.cos(Phi)


where lebedenev_quadrature.txt is a text file containing Lebedev quadrature points obtained from the SPHERE_LEBEDEV_RULE dataset directory under the Sample Files section located at the bottom of the page. The file used in this code is lebedev_131.txt.

#### Testing Chebyshev quadrature on Gauss integral:

#################################
# Testing integration
#################################
def f(X):
return np.exp(-np.sum(X,axis=1,keepdims=True)**2)

# solution
I1 = np.sqrt(pi)/2

I2 = np.einsum(
'ij,ij',
)

print(f'Testing Gauss-Chebyshev quadrature with Gauss Integral :\n I1 = {I1} | I2 = {I2} | I1-I2 = {I1-I2}\n')


Output :

Testing Gauss-Chebyshev quadrature with Gauss Integral :
I1 = 0.8862269254527579 | I2 = 0.8862570079725882 | I1-I2 = -3.0082519830276766e-05


#### Testing Lebedev quadrature on unit sphere surface :

def f(Theta,Phi):
return np.ones_like(Theta)

# solution
I1 = 4*pi #np.exp(-1)*4*pi

I2 = 4.0*pi*np.einsum(
'ij,ij',
)

print(f'Testing Lebedev quadratude with Integral over unit sphere surface :\n I1 = {I1} | I2 = {I2} | I1-I2 = {I1-I2}\n')


Output:

Testing Lebedev quadrature with Integral over unit sphere surface :
I1 = 12.566370614359172 | I2 = 12.566370614362953 | I1-I2 = -3.780087354243733e-12


#### And testing both quadrature to compute the volume of unit sphere :

#""" Unit Sphere
def f(X):
R = norm(X,axis=1,keepdims=True)
return np.where(R<=1,1,0)
I1 = 4/3*pi
#"""

""" Testing with Gauss integral
def f(X):
return np.exp(-norm(X,axis=1,keepdims=True)**2)
I1 = 4*pi*np.sqrt(pi)/2
#"""

I2 = 4*pi*np.einsum(
'ij,ij',
f(X_cartesian_int)*2.0*rm/(1-Mu_int)**2,
Wint
)

print(f'Testing Lebedev + Chebyshev qudrature with Integral over unit sphere volume :\n I1 = {I1} | I2 = {I2} | I1-I2 = {I1-I2}\n')


Output:

Testing Lebedev + Chebyshev quadrature with Integral over unit sphere volume :
I1 = 4.1887902047863905 | I2 = 12.623540989301834 | I1-I2 = -8.434750784515444


#### Complete code :

import numpy as np
import matplotlib.pyplot as plt
from numpy.linalg import norm
pi = np.pi

##########################
##########################

Nr = 30
ii  = np.arange(1,Nr+1).reshape(-1,1)
sin_i = np.sin(ii*pi/(Nr+1))
cos_i = np.cos(ii*pi/(Nr+1))

# Chebyshev quadrature points on [-1,1]
Mu_quadrature = (Nr+1.0-2.0*ii)/(Nr+1.0) + 2.0/pi*(1.0 + 2.0/3.0*sin_i**2)*cos_i*sin_i

# Chebyshev quadrature points on [0, +infty [
rm = 0.35 # Angstrom
# variable substitution for the radial component

##########################
# Quadrature for solid angle / angular coordinates
##########################

# file containing integration points (theta/phi in deg) and associated weights

# make a tensor product of the Quadrature points
X_spherical_int = np.concatenate([
],
axis=1
)

R_int   = X_spherical_int[:,0:1]
Theta_int = X_spherical_int[:,1:2]
Phi_int = X_spherical_int[:,2:3]

Mu_int  = (R_int - rm)/(R_int + rm)

Wint = np.concatenate([
],
axis=1
)
Wint = np.prod(Wint,axis=1,keepdims=True)

X_cartesian_int = np.zeros_like(X_spherical_int)
R, Theta, Phi = R_int, Theta_int, Phi_int
Mu = (R - rm)/(R + rm)
X_cartesian_int[:,0:1] = R*np.cos(Theta)*np.sin(Phi)
X_cartesian_int[:,1:2] = R*np.sin(Theta)*np.sin(Phi)
X_cartesian_int[:,2:3] = R*np.cos(Phi)

#################################
# Testing integration
#################################
def f(X):
return np.exp(-np.sum(X,axis=1,keepdims=True)**2)

# solution
I1 = np.sqrt(pi)/2

I2 = np.einsum(
'ij,ij',
)

print(f'Testing Gauss-Chebyshev quadrature with Gauss Integral :\n I1 = {I1} | I2 = {I2} | I1-I2 = {I1-I2}\n')

def f(Theta,Phi):
return np.ones_like(Theta)

# solution
I1 = 4*pi #np.exp(-1)*4*pi

I2 = 4.0*pi*np.einsum(
'ij,ij',
)

print(f'Testing Lebedev quadrature with Integral over unit sphere surface :\n I1 = {I1} | I2 = {I2} | I1-I2 = {I1-I2}\n')

#""" Unit Sphere
def f(X):
R = norm(X,axis=1,keepdims=True)
return np.where(R<=1,1,0)
I1 = 4/3*pi
#"""

""" Testing with Gauss integral
def f(X):
return np.exp(-norm(X,axis=1,keepdims=True)**2)
I1 = 4*pi*np.sqrt(pi)/2
#"""

I2 = 4*pi*np.einsum(
'ij,ij',
f(X_cartesian_int)*2.0*rm/(1-Mu_int)**2,
Wint
)

print(f'Testing Lebedev + Chebyshev qudrature with Integral over unit sphere volume :\n I1 = {I1} | I2 = {I2} | I1-I2 = {I1-I2}\n')

• +1. Since the content was quite long, I thought it would be good to separate it into respective sub-headings. Please feel free to edit it further. Also the link to the radial variable substitution gives an timed out error: Your session has timed out. Please go back to the article page and click the PDF link again. Commented May 11, 2023 at 21:26
– mle
Commented May 12, 2023 at 5:26
• Could you please the following details by editing your question? You are testing Chebyshev, Lebedev, and both in the various parts of your code. You have listed the output you get. But what was the expected output? Commented May 12, 2023 at 8:50
• Also Are you trying to implement any of the Fortran scripts located at SPHERE_LEBEDEV_RULE dataset directory into Python? Commented May 12, 2023 at 8:57
• The difference between the actual value of the integral and the result obtained with quadrature should be I1 - I2 = 0.
– mle
Commented May 12, 2023 at 11:22

You are missing the key piece of the algorithm: employing Becke's multicenter integration scheme originally described in J. Chem. Phys. 88, 2547–2553 (1988).

An integral $$I =\int f({\bf r}) {\rm d}^3 r$$ that contains atomic cusps can be evaluated efficiently by introducing a resolution of the identity $$\sum_{A} w_{A}({\bf r})=1$$ with a set of atomic weight functions $$0 \leq w_A({\bf r}) \leq 1$$.

Inserting this resolution of the identity and rearranging the summation order, the original integral is rewritten as a sum of atomic integrals $$I = \sum_{A} I_{A}$$ which in turn are given by $$I_{A} = \int w_{A}({\bf r}) f({\bf r}) {\rm d}^3 r.$$ The weight functions $$w_A$$ are chosen such that these integrals can be efficiently computed using atom-centered coordinates $$I_{A} = \int w_{A}({\bf r}-{\bf r}_A) f({\bf r}-{\bf r}_A) r_A^2 {\rm d}r_A d\Omega_A,$$ where $$r^A$$ and $$\Omega^A$$ are the radius and solid angle measured from atom $$A$$. These integrals can now be computed by introducing quadrature grids for the radial $$r_A$$ and solid angle $$\Omega_A$$ coordinates. The key point here is that the additional $$r_A^2$$ factor from the Jacobian kills off any cusp on the nucleus $$A$$, while the weight factors are chosen in a way that off-center cusps are killed off by vanishing $$w_A$$.

There is a lot of flexibility in choosing the weight functions. Becke suggested one way to form $$w_A$$; however, more optimal schemes that converge more rapidly in the radial-angular quadrature have been suggested in later work. For recent work, see e.g. J. Chem. Phys. 149, 204111 (2018)

To compute integrals with atomic quadrature in practice, you can just use numgrid, for instance, which already implements all of the necessary equations and gives you the necessary quadature grids and weights when you feed it the atomic coordinates.

• Yes this is the integration method detailed in the article, but I am taking a step back. I want my quadrature integration to work before adding cell functions.
– mle
Commented May 12, 2023 at 11:19
• But still Ty, I will look into numgrid.
– mle
Commented May 12, 2023 at 11:32
• That is not what you ask above. Commented May 13, 2023 at 13:03

For the last integral including both quadratures, there are two problems :

1. the integrand : a missing factor of R²
2. the volume of a sphere was not an appropriate case I think

Testing with the new changes on the 3D Gauss integral gives good results, but I expected the error difference to be lower since with 30 Gauss-Chebyshev quadrature points I have a $$10^{-5}$$ error on the 1D Gauss integral and the Lebedev quadrature is of high order. At least the error decreases as the number of Gauss-Chebyshev quadrature points increases :

# testing full quadrature scheme with cell functions
# Gauss integral

def f(X):
return np.exp(-np.sum(np.square(X),axis=1,keepdims=True))

I1 = pi**(3/2)
I2 = 4*pi*np.einsum(
'ij,ij',
f(X_cartesian_int)*2.0*rm/(1-Mu)**2*R**2,
Wint
)

print(I1,I2,I1-I2)

Output : 5.568327996831708 5.570977047596403 -0.0026490507646954597