# 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$$.

[1]:

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)

[5]:

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

[3]:

# 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:
[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.

[7]:

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