GSoC 2026 : Improve Robustness of Poisson Solver

[Ao-chuba]

GSoC 2026 Final Report

Contributor: Ao-chuba
Organization: theochem
Project: GRID — Improve Robustness of Poisson Solver
Mentors: @marco-2023, @PaulWAyers, @FarnazH
Project Discussion Thread: orgs/theochem/discussions/11

---

The Problem

GRID evaluates Coulomb (electrostatic) potentials of molecules by numerically solving Poisson's equation. The existing solvers -solve_poisson_ivp and solve_poisson_bvp - struggle with nuclear cusp singularities: the physical infinite-slope spike in the electron density at each atomic nucleus.

This causes:

• Numerical instability and convergence failures for sharp densities (large Gaussian exponents).
• The IVP solver diverging entirely for cusp-like densities (α ≥ 500).
• Accuracy degrading by orders of magnitude as densities become more localized.

The benchmark study (#315) showed the full picture:

α (sharpness)	IVP Relative L2 Error	BVP Relative L2 Error	IVP/BVP Ratio
0.5	1.31e-03	5.41e-07	2,420× worse
10.0	1.47e-02	6.56e-07	22,400× worse
50.0	1.64e-01	1.41e-06	116,000× worse
500.0	Diverged (>10.0)	7.85e-06	—

Even the better BVP solver has no mechanism to handle the nuclear cusp analytically, it simply gets harder as the density sharpens. The goal of this GSoC project was to fix this fundamentally, not just for a specific parameter range.

---

The Solution: Double-Split-Plus-Solve Architecture

Instead of forcing numerical solvers to handle the singularities, My mentor and I designed a "Double-Split-Plus-Solve" architecture that analytically eliminates the nuclear cusps before any numerical integration happens.

ρ_total = ρ_core (analytical) + ρ_bonding (fitted) + ρ_residual (smooth)
V_total = V_core (exact)      + V_bonding (exact)   + V_residual (BVP, now trivially easy)

The key insight: Gaussians have exact analytical Coulomb potentials. If we can express the sharp parts of the density as Gaussians, their potentials can be computed with zero numerical error, leaving only a smooth, well-behaved residual for the BVP solver.

---

Work Delivered

PR 1 — Analytical Coulomb Potential Engine

#304 ·

Implemented the core mathematical engine: coulomb_potential() in src/grid/coulomb.py.

• coulomb_gaussian_s: exact potential of s-type Gaussians $e^{-\alpha r^2}$, with correct $r \to 0$ limit of $2\sqrt{\alpha/\pi}$.
• coulomb_gaussian_p: exact potential of normalized p-type Gaussians $\frac{2}{3}\frac{\alpha^{5/2}}{\pi^{3/2}} r^2 e^{-\alpha r^2}$, with $r \to 0$ limit $\frac{10}{3}\sqrt{\alpha/\pi}$.
• coulomb_potential: multi-center evaluator accumulating contributions from multiple fitted Gaussians at arbitrary 3D points.

This gives zero numerical error for the analytical core subtraction step.

---

PR 2 — Pre-fitted Atomic Gaussian Parameters

#306 ·

Added load_atomic_gaussian_params(element) and the data file src/grid/data/atomic_gauss_params.json.

Key design decision: The initial plan was to use BFit's GreedyKLSCF. However, BFit's output coefficients were in an internal unnormalized convention incompatible with coulomb_potential. The approach was switched to Non-Negative Least Squares (NNLS) fitting against a fixed UGBS exponent grid, which directly produces coefficients in the correct normalized Gaussian convention. Pre-fitted parameters for H, C, N, O, and Cl are included.

This keeps GRID completely free of BFit and AtomDB at runtime, the fitting is done offline, and only the output JSON is shipped with GRID.

---

PR 3 — Split 1: Analytical Core Subtraction

#319 ·

Implemented solve_poisson_robust() in src/grid/robust_poisson.py with the first split stage.

For each atomic center, the pre-fitted Gaussian parameters are loaded, the analytical core density is subtracted from the input, and its exact Coulomb potential is computed and added directly to the final result. The smooth residual is then passed to solve_poisson_bvp.

# Split 1 formula
V_total = coulomb_potential(core_params, points)  # zero numerical error
         + solve_poisson_bvp(ρ - ρ_core, ...)    # now trivially smooth

---

PR 4 — Split 2: NNLS Bonding Density Fitting

#320 ·

Extended solve_poisson_robust() with the split2=True flag. After the analytical core subtraction (Split 1), the post-Split-1 residual is additionally fitted with a geometric Gaussian basis (np.geomspace(0.05, 5000, 20)) using scipy.optimize.nnls.

This step introduces a sequential greedy fitting strategy:

• For each atom, an NNLS problem is solved using the geometric basis placed at that center.
• Non-Negative Least Squares (NNLS) ensures we only subtract physically meaningful (positive) density contributions.
• The fitted density is immediately subtracted before moving to the next atom, meaning subsequent atoms fit a partially reduced residual.

This captures the remaining bonding and polarization density analytically, leaving an incredibly flat and smooth final residual for the BVP solver to handle with maximum accuracy.

---

PR 5 — Tutorial Notebook

#327 ·

Added examples/Robust_Poisson_Solver.ipynb — a complete, executable tutorial demonstrating the full pipeline:

1. The problem: plain BVP on a sharp Gaussian.
2. Split 1: analytical core subtraction.
3. Split 1+2: NNLS bonding density fitting.
4. Convergence benchmark sweeping α from 0.5 to 200.
5. Real atomic density test: Hydrogen pro-atom from grid/data/proatoms/.

Accuracy results from the notebook (α = 50.0):

Method	Relative L2 Error
Plain solve_poisson_bvp	1.48e-06
solve_poisson_robust (Split 1)	1.35e-06
solve_poisson_robust (Split 1+2)	1.66e-07

Split 1+2 achieves a ~9× improvement over the plain BVP solver on a sharp Gaussian density. More importantly, the Hydrogen pro-atom test in step 5 proves that the Double-Split architecture successfully handles true Slater-type physical cusps without divergence, producing a stable, finite, and non-negative potential everywhere.

---

Study — Screened Coulomb vs. Double-Split

#328

An independent mathematical study benchmarking the Screened Coulomb (Yukawa) potential $V_\mu(r) = e^{-\mu r}/r$ as an alternative approach. The conclusion is definitive:

μ (screening)	Relative L2 Error vs. exact 1/r
0.001	2.70e-04
0.100	2.28e-02
0.500	9.41e-02
1.000	16.6%

Screening exponentially suppresses outer tails but does not resolve the $r \to 0$ nuclear cusp singularity and introduces massive physical distortion across the entire grid. The Double-Split architecture is the only approach that preserves exact $1/r$ physics while analytically eliminating the cusp.

---

Links

Resource	Link
Project Discussion	orgs/theochem/discussions/11
PR 1: Coulomb Potential Engine	theochem/grid#304
PR 2: Atomic Gaussian Parameters	theochem/grid#306
Study: IVP vs BVP Benchmarks	theochem/grid#315
PR 3: Split 1 Solver	theochem/grid#319
PR 4: Split 2 Solver	theochem/grid#320
PR 5: Tutorial Notebook	theochem/grid#327
Study: Screened Coulomb vs Double-Split	theochem/grid#328

[PaulWAyers]

Looks great @Ao-chuba . Thanks for doing GSoC with QC-Devs!

[Ao-chuba]

Thank you @PaulWAyers! Working with the QC-Devs and theochem team has been an amazing experience. I really appreciate all the mentorship and support throughout the summer!