GSoC 2025: Improve Robustness of Poisson Solver

[PaulWAyers]

Description

Improve the robustness of the Poisson solvers that grid uses to evaluate the Coulomb potential due to a charge distribution,

$$ \Phi(\mathbf{r}) = \int \frac{\rho(\mathbf{r}')}{|\mathbf{r} - \mathbf{r}'|} d\mathbf{r}' $$

As pointed out by Becke, this is more efficiently treated, numerically, as a solution to the corresponding Poisson equation.

$$ \nabla^2 \Phi(\mathbf{r}) = - 4 \pi \rho(\mathbf{r}) $$

📚 Package Description and Impact

Grid is a pure Python library for numerical integration, interpolation and differentiation of interest for the quantum chemistry community.

👷 What will you do?

Grid has two different Poisson solvers, one using a boundary-value method (which seems more robust, but is not perfect #215) and one using an initial-value method which is less robust (see #162). You will attempt to find robust default parameters for these methods. It is conceivable that one needs (much) larger grids for solving the Poisson equation than we routinely use in numerical integration. (Our Gaussian-quadrature-ish grids are especially adapted to integration; for differentiation their accuracy is diminished. Solving the Poisson equation is essentially a differentiation.)

One approach, which would potentially make the method more robust, would be to use a screened Coulomb kernel instead; that approach can remove the singularity in the equation and, by taking the limit to zero screening, give robust results.

Another approach, which might be much easier, would be to use the strategy proposed in #16. That is, estimate the density with a linear combination of Gaussians using BFit, so that one only needs to solve the Poisson equation for the error in the fit.

🏁 Expected Outcomes

1. Explore the boundary value solver. Find good values for the default parameters.
2. See if the initial value problem can be made more robust or, failing that, if circumstances where it is robust can be identified.
3. Develop a protocol where a sufficient grid, algorithm, and algorithm parameters can be selected, allowing high-level porcelain to be designed for the function.
4. Explore whether using the "promolecular" trick is beneficial.
5. Explore whether screened Coulomb potentials behave better.
6. Write comprehensive tests and documentation for all new functionality.
7. Write tutorial Jupyter notebooks that show how to use the new functionality.

Required skills	Python, OOP
Preferred skills	Be comfortable with math and numerical algorithms. Experience with scientific programming can help
Project size	350 hours, Large
Difficulty	Medium 🤔

🙋 Mentors

Marco Martínez-González	mmg870630_at_gmail_dot_com	@marco-2023
Derrick Yang	yxt1991_at_gmail_dot_com	@tczorro
Esteban Vöhringer-Martinez	estebanvohringer_at_qcmmlab_dot_com	@evohringer
Ali Tehrani	19at27_at_queensu_dot_ca	@Ali-Tehrani

[majilacodes]

Respected mentors @marco-2023 @tczorro @evohringer @Ali-Tehrani,
I wanted to express my keen interest in contributing to the Grid library through GSoC 2025, particularly in improving the robustness of the Poisson solvers. I have a background in Python, object-oriented programming, and scientific computing, and I am comfortable with numerical algorithms and mathematical modeling.

I have been reviewing the existing issues (#215 and #162) and have some thoughts:

1. Have any preliminary experiments been conducted on screened Coulomb potentials, and if so, what challenges were encountered?
2. Regarding the promolecular trick, is the primary goal error reduction in charge density estimation or computational efficiency?
3. I also have experience in data analysis and visualization, which I believe will be useful for analyzing solver performance and identifying optimal parameters. Would there be specific metrics or benchmarking techniques you'd recommend focusing on?

I would love to discuss potential approaches and gain insights into the best way to contribute meaningfully. Any guidance or recommended readings would be greatly appreciated, looking forward to your response!

[PaulWAyers]

1. I don't know of any tests with the current version of Grid. It was something we looked at with a previous version. It's a good idea to see if that changes things.
2. The main goal is that you are learning a much smaller quantity with faster decaying asymptotics; numerical precision should be higher and convergence should be more robust.
3. I think our main metric will be the (screened) Coulomb metric, so just the electrostatic potential times the charge distribution, integrated. Having a pointwise-accurate electrostatic potential (small maximum absolute error) is also important.

The main lead reference on this is the one listed by Becke in the project description. Looking at papers that cite that one is a good way to proceed (though it is very highly cited, so you will need to filter the citing references).

[majilacodes]

Thank you for your response Sir @PaulWAyers,

Over the past days, I spent a good amount of time researching and working with the grid module. Here are the issues I encountered:

1. The most significant challenge is handling singularities near the origin, the BVP solver assumes Φ(0)=0, which may not be valid for all physical scenarios (such as when nucleus is positioned at the origin, the potential will be strongly negative at that point, not zero), whereas without much surprise, the IVP solver documentation explicitly mentions: "difficulty in capturing the origin region".

2. The transform parameter must map from [0,∞) to a finite domain. Improper choice of transform can lead to domain incompatibility errors (I ran into one such error with BeckeRTransform since it produces a (-1,1) domain, no validation is performed to ensure the transform has appropriate properties)

3. The solvers may become unstable for highly irregular charge distributions as I supsect the Becke & Dickson approach might have limitations for complex molecular systems.

Having said that, here are the solutions I am currently working on:

1. I've been exploring FEniCS as a complementary approach for solving the Poisson equation. I'm particularly interested in how FEniCS handles the same singularity issues that I'm addressing in Grid, especially near r→0 where the l(l+1)/r² terms dominate.
   Its approach of expressing the PDE as a variational problem (finding u∈V such that a(u,v)=L(v) for all v∈V̂) provides an elegant mathematical framework that might offer insights for improving Grid's solvers. I'm investigating whether some of these techniques could be adapted to enhance the robustness of Grid's boundary value solver. What are your suggestions regarding this?

2. For the regularization of Singularities, I am currently working to implement screened Coulomb potentials (Yukawa-type potentials) to regularize the r⁻¹ singularity. I suggest we could add automatic refinement of grid points near the origin when needed or include proper error handling for division by zero and small r values.

3. I am also really interested in building what I would call a 'Parameter Optimization Framework' which is basically a robust parameter selection heuristic based on charge distribution characteristics, that can help us in automatic parameter selection based on error estimates. Here, we could also add validation checks for transform domain compatibility, an issue I mentioned above.

That would be all from my side for now as I explore and research more options, looking forward to hear your thoughts on this.

Best Regards,
Akshat Majila