grid.poisson module¶
Poisson solver module.
This module solves the following Poisson equation:
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
- 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 >>> radial_grid = Trapezoidal(10000) >>> # Set up the atomic grid with degree 10 at each radial point. Molecular grid works as well. >>> degree = 10 >>> atomic_grid = AtomGrid(radial, degrees=[degree]) >>> # 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 Poisson`s 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) >>> radial_grid = tf.transform_1d_grid(oned) >>> # Set up the atomic grid with degree 10 at each radial point. Molecular grid works as well. >>> atomic_grid = AtomGrid(radial_grid, degrees=[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 for the potential as an initial value problem and evaluate it over the atomic grid. >>> potential = solve_poisson_ivp( >>> atgrid, func_vals, InverseRTransform(tf), r_interval=(1000, 1e-3), >>> ode_params={"method" : "DOP853", "atol": 1e-8}, >>> ) >>> potential_values = potential(atgrid.points)