Fair enough. Is there any way to get the normalization constants for the radial and angular parts separately?

The Schlegel paper does not spell out the real form, but assuming the normalization is the same as the complex form, then the “orthonormal” normalization is what is being used there.

I am still confused, though, because under the “orthonormal” convention, Y_00 = sqrt(1/(4 Pi)), but I have verified that in psi4, Y_00 = 1. The test for this is to print out the GaussianShell.coef() for an s-function. For this, the value returned = (2 alpha/Pi)^0.75, where alpha is the Gaussian width. This is the formula for the normalization of a Gaussian, so it must mean that Y_00 = 1.