Interpolation and Solving Poisson’s Equation

Details

It is known based on the properties of the real spherical harmonics \(Y_{lm}\) that any \(L_2\) function \(f : \mathbb{R}^3 \rightarrow \mathbb{R}\) can be decomposed as

\[f(r, \theta, \phi) = \sum_{l=0}^\infty \sum_{m=-l}^l w_{lm}(r) Y_{lm}(\theta, \phi),\]

where the unknown \(w_{lm}\) is a function of the radial component and is found as follows. For a fixed \(r\), the radial component \(w_{lm}\) is computed as the integral

\[w_{lm}(r) = \int \int f(r, \theta, \phi) Y_{lm}(\theta, \phi) \sin(\theta) d\theta d\phi\]

over various values of \(r\) and then interpolated \(\tilde{w}_{lm}\) using a cubic spline. This integral can be done using the angular grid module. The interpolation of \(f\) is simply then

\[f(r, \theta, \phi) \approx \sum_{l=0}^{L_{max}} \sum_{m=-l}^l \tilde{w}_{lm}(r) Y_{lm}(\theta, \phi),\]

where \(L_{max}\) is the maximum chosen degree \(l\) of the real spherical harmonics.

Example: Unit-charge distribution

This example tries to interpolate the unit-charge distribution given by \(f(r, \theta, \phi) = \bigg( \frac{\alpha}{\pi} \bigg)^{1.5} e^{-\alpha r_A^2}\), where \(\alpha = 0.25\) and \(r_A\) is the radius centered at the coordinate \([0, 1, 0]\). Since, this function does not depend on the angles, then \(w_{lm} = 0\) when \(l > 0\).

import numpy as np

center = [0, 1, 0]

def charge_distribution(x, alpha=0.25):
    r = np.linalg.norm(x - center, axis=1)
    return (alpha / np.pi) ** (3.0 / 2.0) * np.exp(-alpha * r**2.0)

# Construct Atomic Grid
from grid.onedgrid import GaussLegendre
from grid.rtransform import BeckeRTransform
from grid.atomgrid import AtomGrid

oned = GaussLegendre(npoints=100)
btf = BeckeRTransform(rmin=1e-30, R=1.5)
radial = btf.transform_1d_grid(oned)
degree = 29
atgrid = AtomGrid.from_pruned(radial, 1, sectors_r=[], sectors_degree=[degree], center=center)

# Construct the interpolation of the charge distribution.
density = charge_distribution(atgrid.points)   # evaluate the charge distribution on the atomic grid points.
interpolate_func = atgrid.interpolate(density)

# Compute the density form the interpolation on random points
random_pts = np.vstack(
    (
        np.random.uniform(-1., 1.0, size=20),
        np.random.uniform(0.5, 1.5, size=20),
        np.random.uniform(-1.0, 1.0, size=20)
    )
).T
interpolate = interpolate_func(random_pts)
true = charge_distribution(random_pts)

print("Error difference between interpolation and true:")
print(np.abs(interpolate - true))
assert np.all(np.abs(interpolate - true) < 1e-6)


# First derivative wrt to Cartesian coordinates (x, y, z) can also be calculated
derivative = interpolate_func(random_pts, deriv=1)
# First derivative wrt to Spherical coordinates (r, \theta, \phi)
derivative = interpolate_func(random_pts, deriv=1, deriv_spherical=True)
# Higher order derivative wrt to r can be calculated
sec_deriv = interpolate_func(random_pts, deriv=2, only_radial_deriv=True)

Error difference between interpolation and true:
[2.03885450e-11 5.64345203e-11 1.09432265e-10 6.12496338e-11
 9.15886811e-12 5.53942031e-11 1.82077388e-10 7.62214354e-11
 1.08020218e-10 1.26516332e-11 6.38427089e-11 1.36835335e-12
 2.29083766e-12 1.07915340e-10 1.61463135e-11 8.81852612e-11
 4.39340925e-12 6.39232084e-10 7.09602863e-12 7.00513605e-12]

Solving Poisson Equation

The Poisson equation is defined as \(\nabla^2 V = - 4\pi f(\vec{\textbf{r}})\), where \(V\) is the unknown potential, \(\nabla^2\) is the Laplacian and \(f\) is the charge distribution. It is well-known that the solution to this is given by \(V(\vec{\textbf{r}}) = \int \frac{f(\vec{\textbf{r}}^\prime)}{|\vec{\textbf{r}} - \vec{\textbf{r}}^\prime| } d\vec{\textbf{r}}^\prime\). For the unit-charge distribution given above, the solution is given by \(V(r, \theta, \phi) = \frac{\text{erf} \bigg[ \sqrt{\alpha} r_A \bigg]}{r_A},\) where ERF is the error function. Grid offers two methods of solving the Poisson equation over an atomic grid. The recommended method is the solve_poisson_bvp and is used in this example.

from grid.rtransform import InverseRTransform
from grid.poisson import solve_poisson_bvp
from scipy.special import erf


def actual_potential(x, alpha=0.25):
    r_PC = np.linalg.norm(x, axis=1)
    desired = erf(np.sqrt(alpha) * r_PC) / r_PC
    desired[r_PC == 0.0] = 0.0
    return desired

#Set the charge distribution to be unit-charge density and evaluate on atomic grid points.
func_vals = charge_distribution(atgrid.points)

# Solve for the potential as an initial value problem and evaluate it over the atomic grid.
potential = solve_poisson_bvp(
        atgrid,
        func_vals,
        InverseRTransform(btf),
        remove_large_pts=10.0,
        include_origin=True,
)
potential_values = potential(atgrid.points)


err = np.abs(actual_potential(atgrid.points) - potential_values)
print(np.max(err), np.mean(err), np.std(err))

0.10495702872626828 0.030589794409081324 0.02550959949900786