Mathematica implementation of a primitive (and wrong) algorithm for computing the first ionization energy of a given atom

gist.m

eclist[z_] := ElementData[z, "ElectronConfigurationString"][[1]] /.
  Superscript[a_, b_] -> {a, b}

orbitals[z_] := Partition[Flatten[If[StringQ[eclist[z][[1]]],
    Join[
     orbitals[StringReplace[eclist[z][[1]], ("[" ~~ base__ ~~ "]") -> base]],
     Drop[eclist[z], 1]],
    eclist[z]]], 3]

lastOrbital[z_] := Last[orbitals[z]]

doesShellExist[z_, n_] := Count[orbitals[z], {n, _, _}] > 0

countShellElectrons[z_, n_] := 
 Total[#[[3]] & /@ Cases[orbitals[z], {n, _, _}]]

countOrbitalElectrons[z_, n_, l_] := 
 Total[#[[3]] & /@ Cases[orbitals[z], {n, l, _}]]

countShellSPElectrons[z_, n_] := 
 countOrbitalElectrons[z, n, "s"] + countOrbitalElectrons[z, n, "p"]

countShellSPDElectrons[z_, n_] := 
 countShellSPElectrons[z, n] + countOrbitalElectrons[z, n, "d"]

calculateInnerContribs[z_, n_, use085_] := If[doesShellExist[z, n],
  If[use085, 0.85, 1]*countShellElectrons[z, n] + 
   calculateInnerContribs[z, n - 1, False],
  0]

shieldConst[z_] := With[{l = lastOrbital[z]}, If[l[[1]] == 1, (l[[3]] - 1)*0.3,
   If[l[[2]] == "d" \[Or] l[[2]] == "f",
    If[l[[2]] == "d",
     ((l[[3]] - 1)*0.35) + calculateInnerContribs[z, l[[1]] - 1, False] + 
      countShellSPElectrons[z, l[[1]]],
     ((l[[3]] - 1)*0.35) + calculateInnerContribs[z, l[[1]] - 1, False] + 
      countShellSPDElectrons[z, l[[1]]]],
    ((l[[3]] - 1)*0.35) + calculateInnerContribs[z, l[[1]] - 1, True]]]]

zeff[z_] := ElementData[z, "AtomicNumber"] - shieldConst[z]

neff[z_] := With[{n = lastOrbital[z][[1]]}, Piecewise[{
    {1, n == 1}, {2, n == 2}, {3, n == 3},
    {3.7, n == 4}, {4, n == 5}, {4.2, n == 6}}]]

ei[z_] := (zeff[z]^2/neff[z]^2)*13.6*Quantity["eV"]