18 point Lebedev rule

[PaulWAyers]

The 18 point degree 5 Lebedev rule is nonstandard.

Procedure:
It should be generated by the 6 points:
(+/- 1, 0, 0)
(0,+/-1,0)
(0,0,+/-1)
(+/-1/sqrt(2),+/-1/sqrt(2),0) (plus the 6 combinations)

Then evaluate the spherical harmonics up to degree 5 and determine the integration weights such that the errors are zero. (Solve linear equations). If the solution is not unique, add degree 6 functions as needed and do least squares.

[PaulWAyers]

$$ \sum_{k=0}^{17} w_k Y_l^m(\theta_k,\phi_k) = 0$$

then solve for $w_k$. (Note that for $\ell = m = 0$ the integral is $\sqrt{4 \pi}$, depending on the normalization convention. This gives 1+3+5+7+9+11 equations with 18 unknowns.

In the unlikely event that the equations are undetermined, either minimize the error for $\ell=6$ or choose the least-norm solution. (I'd choose the latter; it's easier.)

[Ali-Tehrani]

This is addressed in the pull-request above

The code to generate the Lebedev grid is given below. Note that the integral of spherical harmonics at l=0 needs to be divided by 4pi since Lebedev weights are assumed to be multiplied by 4pi.

import numpy as np
from grid.utils import (
    convert_cart_to_sph,
    generate_real_spherical_harmonics,
    generate_real_spherical_harmonics_scipy
)


factor = 1.0 / np.sqrt(2.0)
pts = np.array([
    [1.0, 0.0, 0.0],
    [-1.0, 0.0, 0.0],
    [0.0, 1.0, 0.0],
    [0.0, -1.0, 0.0],
    [0.0, 0.0, 1.0],
    [0.0, 0.0, -1.0],
    [factor, factor, 0.0],
    [factor, -factor, 0.0],
    [-factor, factor, 0.0],
    [-factor, -factor, 0.0],
    [factor, 0.0, factor],
    [factor, 0.0, -factor],
    [-factor, 0.0, factor],
    [-factor, 0.0, -factor],
    [0.0, factor, factor],
    [0.0, factor, -factor],
    [0.0, -factor, factor],
    [0.0, -factor, -factor],
])
print(f"Number of points {len(pts)}")
print(f"Norm of pts {np.linalg.norm(pts, axis=1)}")

# Generate matrix A, and vector B,  Ax = B
r, theta, phi = convert_cart_to_sph(pts).T
A = generate_real_spherical_harmonics_scipy(5, theta, phi) # Rows are degrees and columns are points
b = np.array([np.sqrt(4.0 * np.pi) / (4.0 * np.pi)] + [0.0] * 35)

# Solve for x in Ax=B
sol, residues, rank, _ = np.linalg.lstsq(A, b)

# Print Results
print(f"Weights: {sol}")
print(f"Residues: {residues}")
print(f"Rank: {rank}")


save_dict = {
    "points": pts,
    "weights": sol,
    "degree": 5,
    "size": len(pts)
}
np.savez("lebedev_5_18.npz", **save_dict)

The results are

Number of points 18
Norm of pts [1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
Weights: [0.03333333 0.03333333 0.03333333 0.03333333 0.03333333 0.03333333
 0.06666667 0.06666667 0.06666667 0.06666667 0.06666667 0.06666667
 0.06666667 0.06666667 0.06666667 0.06666667 0.06666667 0.06666667]
Residues: [1.52913967e-32]
Rank: 18