Improve Robustness of Poisson Solver

The grid library evaluates the Coulomb potential of molecules by numerically solving Poisson's equation. However, true electron densities possess sharp, infinite slope cusps at atomic nuclei. These physical singularities cause the existing numerical Initial Value Problem (IVP) and Boundary Value Problem (BVP) solvers to suffer from catastrophic instability and convergence failures. This project solves this by implementing a "Double-Split-Plus-Solve" architecture. Instead of forcing numerical solvers to handle the singularities, this strategy uses the BFit library and AtomDB to analytically fit and extract the sharp promolecular atomic cores and bulk bonding densities using Gaussians. Because Gaussians have exact analytical potentials, we calculate the potential of the cusps with zero numerical error. This leaves only a tiny, perfectly smooth residual density for grid's numerical BVP solver to integrate, allowing it to converge rapidly and accurately without crashing.