[BUG] : Incorrect normalization in screening.py upper bound calculations

[Ao-chuba]

@FarnazH @PaulWAyers @Ali-Tehrani
There is a discrepancy in the normalization constants used in gbasis/screening.py compared to the rest of the library. Functions like compute_primitive_upper_bound and compute_primitive_cutoff_radius utilize a 1D normalization factor, whereas gbasis primarily operates with 3D Gaussian basis functions. Specifically, in
screening.py (lines 115-116/208-212):
fix was mentioned as TODO: Fix this, it fails for some reason

161: # TODO: Fix this, it fails for some reason, it is needed for screening of 1rdms

[Ao-chuba]

While looking into the TODO and profiling the formulas, I noticed something I'm a bit confused, perhaps I'm misunderstanding the design intent here?
When evaluating a 3D carbon p-orbital contraction at r=1.0 Bohr using _eval_deriv_contractions, the maximum true value is roughly 0.392. However, the 1D formulas in compute_primitive_upper_bound (using the 1/4th power normalization) predict a bound of 0.364, which doesn't seem to bound the 3D function. It looks like the test_compute_contraction_upper_bound currently passes because compute_contraction_upper_bound multiplies the primitive bound by max_norm_cart (the 3D normalization tensor). This multiplies the 1D norm by the 3D norm, inflating the bound to ~5620.18. Since the test simply checks if 0.392 <= 5620.18, it passes but this bound seems orders of magnitude too loose to be useful for 1RDMS screening.
Was this 1D $\times$ 3D normalization multiplication done deliberately to ensure the bound was safely above any possible value for the time being, or is this the mathematical discrepancy the TODO was referring to?