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')
Your session has timed out. Please go back to the article page and click the PDF link again.
$\endgroup$