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(\phi) 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\).

[5]:

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)

[6]:

# 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, r_sectors=[], d_sectors=[degree], center=center)

[7]:

# 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:
[5.50204594e-11 4.97960712e-11 9.96328462e-11 1.50277207e-11
 3.25926654e-11 2.01591117e-11 3.26607444e-11 8.44455422e-11
 4.00808614e-11 5.06480266e-10 7.52775412e-11 3.48676058e-11
 4.13636728e-11 1.11618694e-11 4.14768677e-11 4.12894238e-11
 7.26674645e-11 8.43535809e-11 5.25728962e-11 6.84093102e-11]

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.

[8]:

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 - center, 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(f"The maximum error: {np.max(err)}")
print(f"The mean error: {np.mean(err)}")
print(f"The standard dev: {np.std(err)}")

The maximum error: 1.4653770128456663e-06
The mean error: 1.1070838082577226e-06
The standard dev: 4.600048234819135e-07