# grid.poisson module¶

Poisson solver module.

This module solves the following Poisson equation:

$\nabla^2 V(r) = -4\pi \rho(r),$

for some Coulomb potential $$V(r)$$ and charge density $$\rho(r)$$ over a centered atomic grid. It is recommended to use the boundary value problem for handing singularities near the origin of the atomic grid.

interpolate_laplacian(molgrid, func_vals)

Return a function that interpolates the Laplacian of a function.

$\nabla^2 f = \frac{1}{r}\frac{\partial^2 rf}{\partial r^2} - \frac{\hat{L}}{r^2},$

such that the angular momentum operator satisfies $$\hat{L}(Y_l^m) = l (l + 1) Y_l^m$$. Expanding f in terms of spherical harmonic expansion, we get that

$\nabla^2 f = \sum_l \sum_m \bigg[ \frac{\partial^2 \rho_{lm}(r)}{\partial r^2} + \frac{2}{r} \frac{\partial \rho_{lm}(r)}{\partial r} - \frac{l(l+1)}{r^2}\rho_{lm}(r) \bigg] Y_l^m,$

where $$\rho_{lm}^f$$ is the lth, mth radial component of function f.

Parameters
• molgrid (Union[MolGrid, AtomGrid]) – Atomic grid that can integrate spherical functions and interpolate radial components.

• func_vals (ndarray(N,)) – The function values evaluated on all $$N$$ points on the atomic grid.

Returns

Function that interpolates the Laplacian of a function whose input is Cartesian points. The float value is the cutoff where radial points smaller than the cutoff are replaced with the cutoff. Computing the Laplacian at r=0 can cause problems depending on the function provided.

Return type

callable[ndarray(M,3), float -> ndarray(M,)]

Warning

• Since $$\rho_{lm}$$ and its derivatives are being interpolated and due to division by powers of $$r$$, it is recommended to be very careful of having values near zero.

solve_poisson_bvp(molgrid, func_vals, transform, boundary=None, include_origin=True, remove_large_pts=1000000.0, ode_params=None)

Return interpolation of the solution to the Poisson equation solved as a boundary value problem.

The Poisson equation solves for function $$g$$ of the following:

$\nabla^2 g = (-4\pi) f,$

for a fixed function $$f$$, where $$\nabla^2$$ is the Laplacian. This is transformed to an set of ODE problems as a boundary value problem.

If boundary is not provided, then the boundary of $$g$$ for large r is set to $$\int \int \int f(r, \theta, \phi) / r$$. The solution $$g$$ is assumed to be zero at the origin $$g(0, \theta, \phi) = 0$$. Use solve_poisson_ivp if this assumption isn’t needed.

Parameters
• molgrid (Union[MolGrid, AtomGrid]) – Molecular or atomic grid that is used for integration and expanding func into real spherical harmonic basis.

• func_vals (ndarray(N,)) – The function values evaluated on all $$N$$ points on the molecular grid.

• transform (BaseTransform, optional) – Transformation from infinite domain $$r \in [0, \infty)$$ to another domain that is a finite.

• boundary (float, optional) – The boundary value of $$g$$ in the limit of r to infinity.

• include_origin (bool, optional) – If true, will add r=0 point when solving for the ode only. If false, it is recommended to have many radial points near the origin.

• remove_large_pts (float, optional) – If true, will remove any points larger than remove_large_pts when solving for the ode only.

• ode_params (dict, optional) – The parameters for the ode solver. See grid.ode.solve_ode_bvp for all options.

Returns

The solution to Poisson equaiton/potential $$g : \mathbb{R}^3 \rightarrow \mathbb{R}$$.

Return type

callable(ndarray(N, 3) -> float)

Examples

>>> # Set up of the radial grid
>>> # Set up the atomic grid with degree 10 at each radial point. Molecular grid works as well.
>>> degree = 10
>>> # Set the charge distribution to be unit-charge density and evaluate on atomic grid points.
>>> def charge_distribution(x, alpha=0.1):
>>>    r = np.linalg.norm(x, axis=1)
>>>    return (alpha / np.pi)**(3.0 / 2.0) * np.exp(-alpha * r**2.0)
>>> func_vals = charge_distribution(atomic_grid.points)
>>> # Solve the Poisson equation with Becke transformation
>>> transform = BeckeRTransform(1e-6, 1.5, trim_inf=True)
>>> potential = solve_poisson_bvp(
>>>      atgrid, func_vals, InverseRTransform(tf), include_origin=True,
>>>      remove_large_pts=1e6, ode_params={"tol" : 1e-6, "max_nodes": 20000},
>>> )
>>> actual = potential(atgrid.points)


References

1

Becke, A. D., & Dickson, R. M. (1988). Numerical solution of Poissons equation in polyatomic molecules. The Journal of chemical physics, 89(5), 2993-2997.

solve_poisson_ivp(molgrid, func_vals, transform, r_interval=(1000, 1e-05), ode_params=None)

Return interpolation of the solution to the Poisson equation solved as an initial value problem.

The Poisson equation solves for function $$g$$ of the following:

$\nabla^2 g = (-4\pi) f,$

for a fixed function $$f$$, where $$\nabla^2$$ is the Laplacian. This is transformed to a set of ODE problems as an initial value problem.

Ihe initial value problem is chosen so that the boundary of $$g$$ for large r is set to $$\int \int \int f(r, \theta, \phi) / r$$. Depending on $$f$$, this function has difficulty in capturing the origin $$r=0$$ region, and is recommended to keep the final interval $$a$$ close to zero.

Parameters
• molgrid (Union[MolGrid, AtomGrid]) – Molecular or atomic grid that is used for integration and expanding func into real spherical harmonic basis.

• func_vals (ndarray(N,)) – The function values evaluated on all $$N$$ points on the molecular grid.

• transform (BaseTransform, optional) – Transformation from infinite domain $$r \in [0, \infty)$$ to another domain that is a finite.

• r_interval (tuple, optional) – The interval $$(b, a)$$ of $$r$$ for which the ODE solver will start from and end, where $$b>a$$. The value $$b$$ should be large as it determines the asymptotic region of $$g$$ and value $$a$$ is recommended to be small but not zero depending on $$f$$.

• ode_params (dict, optional) – The parameters for the ode solver. See grid.ode.solve_ode_ivp for all options.

Returns

The solution to Poisson equaiton/potential $$g : \mathbb{R}^3 \rightarrow \mathbb{R}$$.

Return type

callable(ndarray(N, 3) -> float)

Examples

>>> # Set up of the radial grid
>>> oned_grid = Trapezoidal(10000)
>>> tf = LinearFiniteRTransform(0.0, 1000)
>>> # Set up the atomic grid with degree 10 at each radial point. Molecular grid works as well.
`