-- the full relativistic Hamiltonian of Ne from its own Dirac-Fock orbitals
G, F, x, Atom = CalculateAtomicBasisFunctions(10)
indices, NF = CreateAtomicIndicesList(Atom.Labels, {{"Kappas", true}})
H = NewOperator("Dirac", NF, 10, indices, G, F, Atom.Nucleus) + NewOperator("AU", NF, indices, G, F)
psi = NewWavefunction(NF, 0, {{string.rep("1", NF), 1}})
print(string.format("Etot = %.10f Hartree, <psi|H|psi> = %.10f Hartree", Atom.Etot, psi*H*psi))