# 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)
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

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