PR1 : implement analytical Coulomb potential for Gaussian densities

[Ao-chuba]

First inital PR for Improvement of the robustness of poisson equation .
This PR implements the core mathematical engine for evaluating the exact analytical electrostatic (Coulomb) potential of $s$-type and $p$-type radial Gaussian charge densities. The math is derived from the exact formulations and $r \to 0$ limits outlined in the course project reference (https://qchem.qc-edu.org/problems/HF_dft.html). This PR contains the standalone coulomb.py module and an isolated unit test suite validating math correctness across various distances and normalization schemes, establishing the foundation for the subsequent robust Poisson solver scaffolding.
related to issue -> #288

[Copilot]

Pull request overview

Adds a new standalone coulomb module providing exact analytical Coulomb potential evaluations for s-type and p-type radial Gaussian charge densities, together with a higher-level coulomb_potential aggregator that sums contributions from a set of centered Gaussians. This is the first building block towards a more robust Poisson solver (issue #288).

Changes:

• New src/grid/coulomb.py with coulomb_gaussian_s, coulomb_gaussian_p, and coulomb_potential, including correct handling of the r → 0 limit for both normalized and unnormalized conventions.
• New src/grid/tests/test_coulomb.py exercising the analytical formulas at r=0 and at non-zero radii for both normalizations, plus an aggregate s-type-only consistency test against a hand-rolled reference.

Reviewed changes

Copilot reviewed 2 out of 2 changed files in this pull request and generated 4 comments.

File	Description
src/grid/coulomb.py	New module implementing analytical Coulomb potentials of s- and p-type Gaussian densities and a multi-Gaussian aggregator.
src/grid/tests/test_coulomb.py	New unit tests validating the s/p formulas at the origin and at non-zero radii, and an aggregate s-only check.

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[Ao-chuba]

as per copilot suggestions i took into consideration for some of its inputs:

• Added validation inside coulomb_potential ensuring that coeffs_p, expons_p, and centers_p are either all provided or all omitted to prevent silent dropping of p-type density.
• Added checks to ensure that the coefficients, exponents, and centers have the same length. This prevents silent truncation of Gaussians due to Python's zip behavior.
• Added checks to ensure alpha > 0 and beta > 0. This prevents invalid math operations (like taking np.sqrt of negative numbers or division by zero).
• changed number 1e-12 to a named module-level constant (_R_ZERO_THRESHOLD) to avoid duplication and improve readability.

[marco-2023]

The code can and should be shorter and simpler (read):

out = erf(sqrt_a * r) / r
out[r < _R_ZERO_THRESHOLD] = 2.0 * sqrt_a / np.sqrt(np.pi)

Code is repeated inside the if-else. It is better to compute out and if not normalizes multiply by the prefactor at the end (you can even do it with return)

# Apply normalization factor only if needed
if normalized:
    return out

return prefactor * out

[Ao-chuba]

prefactor = 1.0 if normalized else (np.pi / alpha) ** 1.5
out[mask_nonzero] = prefactor * erf(sqrt_a * r_nonzero) / r_nonzero
out[mask_zero]    = prefactor * 2.0 * sqrt_a / np.sqrt(np.pi)

Would be better if its structured like this instead?

[marco-2023]

I don't think we should have two masks just for one or two points. My idea was something like:

out = erf(sqrt_a * r) / r
out[r < _R_ZERO_THRESHOLD] = 2.0 * sqrt_a / np.sqrt(np.pi)

if normalized:
    return out

return prefactor * out

[marco-2023]

This pattern is even better

np.divide(erf(sqrt_a * r), r, out = out, where=r >= _R_ZERO_THRESHOLD)
out[r < _R_ZERO_THRESHOLD] = 2.0 * sqrt_a / np.sqrt(np.pi)

if normalized:
    return out

return prefactor * out

[marco-2023]

This function can also be simplified, similar to the previous one.
Compute term 1 and term 2, and use if-else for prefactors at return.
We can also get rid of masks.

[marco-2023]

It is a better practice to store the reference value in one variable, the result in another, and compare them at the end. And being consistent with this across tests, similar to what you did for test_coulomb_gaussian_p

[marco-2023]

Nice work! I think this is a good progress.
Please simplify the functions.
Homogenize the tests

I think this module should not be in the root of the package, but we can leave that for another moment.