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

Question

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

##########################
# Quadrature weights/points
##########################

# for radial component
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 [
# Bragg-slater Radius for Hydrogen
rm = 0.35 # Angstrom
# variable substitution for the radial component 
R_quadrature = rm*(1+Mu_quadrature)/(1-Mu_quadrature)

# Chebyshev quadrature weights
Wr_quadrature = 16.0/3.0/(Nr+1)*sin_i**4

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

# file containing integration points (theta/phi in deg) and associated weights
fn = 'lebedenev_quadrature.txt'
Data = np.loadtxt(fn)
ThetaPhi_quadrature = Data[:,:2]

# convert to rad
ThetaPhi_quadrature = ThetaPhi_quadrature*pi/180.0
Wangular_quadrature = Data[:,2:3]

# make a tensor product of the Quadrature points
X_spherical_int = np.concatenate([
        np.tile(R_quadrature,(len(ThetaPhi_quadrature),1)),
        np.repeat(ThetaPhi_quadrature,len(R_quadrature),axis=0)   
    ],
    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([
        np.tile(Wr_quadrature,(len(Wangular_quadrature),1)),
        np.repeat(Wangular_quadrature,len(Wr_quadrature),axis=0)   
    ],
    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

# quadrature
I2 = np.einsum(
    'ij,ij',
    f(R_quadrature)*2.0*rm/(1-Mu_quadrature)**2,
    Wr_quadrature
)

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

# quadrature
I2 = 4.0*pi*np.einsum(
    'ij,ij',
    f(ThetaPhi_quadrature[:,1:2],None),
    Wangular_quadrature
)
 
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

##########################
# Quadrature weights/points
##########################

# for radial component
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 [
# Bragg-slater Radius for Hydrogen
rm = 0.35 # Angstrom
# variable substitution for the radial component 
R_quadrature = rm*(1+Mu_quadrature)/(1-Mu_quadrature)

# Chebyshev quadrature weights
Wr_quadrature = 16.0/3.0/(Nr+1)*sin_i**4

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

# file containing integration points (theta/phi in deg) and associated weights
fn = 'lebedenev_quadrature.txt'
Data = np.loadtxt(fn)
ThetaPhi_quadrature = Data[:,:2]

# convert to rad
ThetaPhi_quadrature = ThetaPhi_quadrature*pi/180.0
Wangular_quadrature = Data[:,2:3]

# make a tensor product of the Quadrature points
X_spherical_int = np.concatenate([
        np.tile(R_quadrature,(len(ThetaPhi_quadrature),1)),
        np.repeat(ThetaPhi_quadrature,len(R_quadrature),axis=0)   
    ],
    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([
        np.tile(Wr_quadrature,(len(Wangular_quadrature),1)),
        np.repeat(Wangular_quadrature,len(Wr_quadrature),axis=0)   
    ],
    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

# quadrature
I2 = np.einsum(
    'ij,ij',
    f(R_quadrature)*2.0*rm/(1-Mu_quadrature)**2,
    Wr_quadrature
)

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

# quadrature
I2 = 4.0*pi*np.einsum(
    'ij,ij',
    f(ThetaPhi_quadrature[:,1:2],None),
    Wangular_quadrature
)
 
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')

Answer 1

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.

Answer 2

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