-- a confined basis for Fe: all orbitals up to 5s, 4p and 4d, the empty 4p and 4d contracted
G, F, x, Atom = CalculateAtomicBasisFunctions(26, "[Ar] 3d6 4s2", {{"Method", "LDA"},
  {"NMax", {s = 5, p = 4, d = 4}},
  {"Confinement", {Rho = 4, N = 8}},
  {"Q", {["4p"] = 1, ["4d"] = 2}},
  {"P", {["4d"] = 0.8}}})
for i, o in ipairs(Atom.Orbitals) do
  print(string.format("%-6s occupation %5.3f  Q %3.1f  rc %5.2f  energy %11.5f  <r> %7.4f", o.Label, o.Occupation, o.Q, o.rc, o.Energy, o.RExpectation[1]))
end
print(string.format("smallest overlap eigenvalue %.3f, orthonormality error %.1e", Atom.MinOverlapEigenvalue, Atom.OrthonormalityError))