###
post: What NIST says an atom's energy is
date: 2026-09-23

Generated by `tools/nist_atomic_reference.py` from NIST Standard Reference
Database 141, Atomic Reference Data for Electronic Structure Calculations:
<https://math.nist.gov/DFTdata/atomdata/tables/ptable.html>

Total energies, their decomposition into kinetic, Hartree, electron-nucleus
and exchange-correlation parts, and the orbital eigenvalues, for the ground
configuration of each atom. NIST's stated target is one microHartree, and
the values were reproduced across several independent codes and machines --
which is what makes this an oracle rather than another code's opinion.

Only the first column, the plain local density approximation, is taken here.
The relativistic variants are what a Dirac solver would be held to.

Do not edit this file.
###

# ---- H (Z = 1), 1s ----
H_Z = 1
H_LDA_ETOT = -0.445671
H_LDA_EKIN = 0.425027
H_LDA_ECOUL = 0.282827
H_LDA_EENUC = -0.920999
H_LDA_EXC = -0.232525
#| The occupied orbital eigenvalues in LDA, in the order 1s.
def h_lda_eigs() {
    return [-0.233471]
}

# ---- He (Z = 2), 1s^2 ----
HE_Z = 2
HE_LDA_ETOT = -2.834836
HE_LDA_EKIN = 2.767922
HE_LDA_ECOUL = 1.99612
HE_LDA_EENUC = -6.625564
HE_LDA_EXC = -0.973314
#| The occupied orbital eigenvalues in LDA, in the order 1s.
def he_lda_eigs() {
    return [-0.570425]
}

# ---- Li (Z = 3), [He] 2s^1 ----
LI_Z = 3
LI_LDA_ETOT = -7.335195
LI_LDA_EKIN = 7.237484
LI_LDA_ECOUL = 3.990755
LI_LDA_EENUC = -16.911041
LI_LDA_EXC = -1.652393
#| The occupied orbital eigenvalues in LDA, in the order 1s 2s.
def li_lda_eigs() {
    return [-1.878564, -0.10554]
}

# ---- Be (Z = 4), [He] 2s^2 ----
BE_Z = 4
BE_LDA_ETOT = -14.447209
BE_LDA_EKIN = 14.309424
BE_LDA_ECOUL = 7.115257
BE_LDA_EENUC = -33.357034
BE_LDA_EXC = -2.514856
#| The occupied orbital eigenvalues in LDA, in the order 1s 2s.
def be_lda_eigs() {
    return [-3.856411, -0.205744]
}

# ---- B (Z = 5), [He] 2s^2 2p^1 ----
B_Z = 5
B_LDA_ETOT = -24.344198
B_LDA_EKIN = 24.161047
B_LDA_ECOUL = 11.503002
B_LDA_EENUC = -56.484587
B_LDA_EXC = -3.52366
#| The occupied orbital eigenvalues in LDA, in the order 1s 2s 2p.
def b_lda_eigs() {
    return [-6.564347, -0.344701, -0.136603]
}

# ---- C (Z = 6), [He] 2s^2 2p^2 ----
C_Z = 6
C_LDA_ETOT = -37.425749
C_LDA_EKIN = 37.190391
C_LDA_ECOUL = 17.627997
C_LDA_EENUC = -87.515412
C_LDA_EXC = -4.728724
#| The occupied orbital eigenvalues in LDA, in the order 1s 2s 2p.
def c_lda_eigs() {
    return [-9.947718, -0.500866, -0.199186]
}

# ---- N (Z = 7), [He] 2s^2 2p^3 ----
N_Z = 7
N_LDA_ETOT = -54.025016
N_LDA_EKIN = 53.731727
N_LDA_ECOUL = 25.799845
N_LDA_EENUC = -127.416229
N_LDA_EXC = -6.140359
#| The occupied orbital eigenvalues in LDA, in the order 1s 2s 2p.
def n_lda_eigs() {
    return [-14.011501, -0.676151, -0.266297]
}

# ---- O (Z = 8), [He] 2s^2 2p^4 ----
O_Z = 8
O_LDA_ETOT = -74.473077
O_LDA_EKIN = 74.116881
O_LDA_ECOUL = 36.331102
O_LDA_EENUC = -177.152578
O_LDA_EXC = -7.768482
#| The occupied orbital eigenvalues in LDA, in the order 1s 2s 2p.
def o_lda_eigs() {
    return [-18.758245, -0.871362, -0.338381]
}

# ---- F (Z = 9), [He] 2s^2 2p^5 ----
F_Z = 9
F_LDA_ETOT = -99.099648
F_LDA_EKIN = 98.676131
F_LDA_ECOUL = 49.535385
F_LDA_EENUC = -237.688846
F_LDA_EXC = -9.622319
#| The occupied orbital eigenvalues in LDA, in the order 1s 2s 2p.
def f_lda_eigs() {
    return [-24.189391, -1.086859, -0.415606]
}

# ---- Ne (Z = 10), [He] 2s^2 2p^6 ----
NE_Z = 10
NE_LDA_ETOT = -128.233481
NE_LDA_EKIN = 127.738667
NE_LDA_ECOUL = 65.726488
NE_LDA_EENUC = -309.988206
NE_LDA_EXC = -11.71043
#| The occupied orbital eigenvalues in LDA, in the order 1s 2s 2p.
def ne_lda_eigs() {
    return [-30.305855, -1.322809, -0.498034]
}

# ---- Na (Z = 11), [Ne] 3s^1 ----
NA_Z = 11
NA_LDA_ETOT = -161.44006
NA_LDA_EKIN = 160.896571
NA_LDA_ECOUL = 79.732497
NA_LDA_EENUC = -388.546372
NA_LDA_EXC = -13.522757
#| The occupied orbital eigenvalues in LDA, in the order 1s 2s 2p 3s.
def na_lda_eigs() {
    return [-37.719975, -2.063401, -1.060636, -0.103415]
}

# ---- Mg (Z = 12), [Ne] 3s^2 ----
MG_Z = 12
MG_LDA_ETOT = -199.139406
MG_LDA_EKIN = 198.541505
MG_LDA_ECOUL = 95.67329
MG_LDA_EENUC = -477.899149
MG_LDA_EXC = -15.455051
#| The occupied orbital eigenvalues in LDA, in the order 1s 2s 2p 3s.
def mg_lda_eigs() {
    return [-45.973167, -2.903746, -1.71897, -0.175427]
}

# ---- Al (Z = 13), [Ne] 3s^2 3p^1 ----
AL_Z = 13
AL_LDA_ETOT = -241.315573
AL_LDA_EKIN = 240.663489
AL_LDA_ECOUL = 112.670733
AL_LDA_EENUC = -577.205757
AL_LDA_EXC = -17.444038
#| The occupied orbital eigenvalues in LDA, in the order 1s 2s 2p 3s 3p.
def al_lda_eigs() {
    return [-55.156044, -3.934827, -2.564018, -0.286883, -0.102545]
}

# ---- Si (Z = 14), [Ne] 3s^2 3p^2 ----
SI_Z = 14
SI_LDA_ETOT = -288.198397
SI_LDA_EKIN = 287.48774
SI_LDA_ECOUL = 131.767814
SI_LDA_EENUC = -687.900675
SI_LDA_EXC = -19.553276
#| The occupied orbital eigenvalues in LDA, in the order 1s 2s 2p 3s 3p.
def si_lda_eigs() {
    return [-65.184426, -5.075056, -3.514938, -0.398139, -0.153293]
}

# ---- P (Z = 15), [Ne] 3s^2 3p^3 ----
P_Z = 15
P_LDA_ETOT = -339.946219
P_LDA_EKIN = 339.173459
P_LDA_ECOUL = 153.065007
P_LDA_EENUC = -810.4001
P_LDA_EXC = -21.784585
#| The occupied orbital eigenvalues in LDA, in the order 1s 2s 2p 3s 3p.
def p_lda_eigs() {
    return [-76.061897, -6.329346, -4.576617, -0.512364, -0.20608]
}

# ---- S (Z = 16), [Ne] 3s^2 3p^4 ----
S_Z = 16
S_LDA_ETOT = -396.716081
S_LDA_EKIN = 395.878143
S_LDA_ECOUL = 176.68795
S_LDA_EENUC = -945.141155
S_LDA_EXC = -24.141019
#| The occupied orbital eigenvalues in LDA, in the order 1s 2s 2p 3s 3p.
def s_lda_eigs() {
    return [-87.789937, -7.69994, -5.751257, -0.630912, -0.261676]
}

# ---- Cl (Z = 17), [Ne] 3s^2 3p^5 ----
CL_Z = 17
CL_LDA_ETOT = -458.664179
CL_LDA_EKIN = 457.758293
CL_LDA_ECOUL = 202.772337
CL_LDA_EENUC = -1092.568993
CL_LDA_EXC = -26.625816
#| The occupied orbital eigenvalues in LDA, in the order 1s 2s 2p 3s 3p.
def cl_lda_eigs() {
    return [-100.369229, -9.187993, -7.039982, -0.754458, -0.32038]
}

# ---- Ar (Z = 18), [Ne] 3s^2 3p^6 ----
AR_Z = 18
AR_LDA_ETOT = -525.946195
AR_LDA_EKIN = 524.969812
AR_LDA_ECOUL = 231.458124
AR_LDA_EENUC = -1253.131982
AR_LDA_EXC = -29.242149
#| The occupied orbital eigenvalues in LDA, in the order 1s 2s 2p 3s 3p.
def ar_lda_eigs() {
    return [-113.800134, -10.794172, -8.443439, -0.883384, -0.38233]
}

test "the NIST energies add up the way NIST says they do" {
    # Etot = Ekin + Ecoul + Eenuc + Exc, exactly, for every atom taken
    assert_near(H_LDA_ETOT, H_LDA_EKIN + H_LDA_ECOUL + H_LDA_EENUC + H_LDA_EXC, 2.0e-6)
    assert_near(HE_LDA_ETOT, HE_LDA_EKIN + HE_LDA_ECOUL + HE_LDA_EENUC + HE_LDA_EXC, 2.0e-6)
    assert_near(LI_LDA_ETOT, LI_LDA_EKIN + LI_LDA_ECOUL + LI_LDA_EENUC + LI_LDA_EXC, 2.0e-6)
    assert_near(BE_LDA_ETOT, BE_LDA_EKIN + BE_LDA_ECOUL + BE_LDA_EENUC + BE_LDA_EXC, 2.0e-6)
    assert_near(B_LDA_ETOT, B_LDA_EKIN + B_LDA_ECOUL + B_LDA_EENUC + B_LDA_EXC, 2.0e-6)
    assert_near(C_LDA_ETOT, C_LDA_EKIN + C_LDA_ECOUL + C_LDA_EENUC + C_LDA_EXC, 2.0e-6)
    assert_near(N_LDA_ETOT, N_LDA_EKIN + N_LDA_ECOUL + N_LDA_EENUC + N_LDA_EXC, 2.0e-6)
    assert_near(O_LDA_ETOT, O_LDA_EKIN + O_LDA_ECOUL + O_LDA_EENUC + O_LDA_EXC, 2.0e-6)
    assert_near(F_LDA_ETOT, F_LDA_EKIN + F_LDA_ECOUL + F_LDA_EENUC + F_LDA_EXC, 2.0e-6)
    assert_near(NE_LDA_ETOT, NE_LDA_EKIN + NE_LDA_ECOUL + NE_LDA_EENUC + NE_LDA_EXC, 2.0e-6)
    assert_near(NA_LDA_ETOT, NA_LDA_EKIN + NA_LDA_ECOUL + NA_LDA_EENUC + NA_LDA_EXC, 2.0e-6)
    assert_near(MG_LDA_ETOT, MG_LDA_EKIN + MG_LDA_ECOUL + MG_LDA_EENUC + MG_LDA_EXC, 2.0e-6)
    assert_near(AL_LDA_ETOT, AL_LDA_EKIN + AL_LDA_ECOUL + AL_LDA_EENUC + AL_LDA_EXC, 2.0e-6)
    assert_near(SI_LDA_ETOT, SI_LDA_EKIN + SI_LDA_ECOUL + SI_LDA_EENUC + SI_LDA_EXC, 2.0e-6)
    assert_near(P_LDA_ETOT, P_LDA_EKIN + P_LDA_ECOUL + P_LDA_EENUC + P_LDA_EXC, 2.0e-6)
    assert_near(S_LDA_ETOT, S_LDA_EKIN + S_LDA_ECOUL + S_LDA_EENUC + S_LDA_EXC, 2.0e-6)
    assert_near(CL_LDA_ETOT, CL_LDA_EKIN + CL_LDA_ECOUL + CL_LDA_EENUC + CL_LDA_EXC, 2.0e-6)
    assert_near(AR_LDA_ETOT, AR_LDA_EKIN + AR_LDA_ECOUL + AR_LDA_EENUC + AR_LDA_EXC, 2.0e-6)
}

test "the atoms are bound, their shells are ordered, and heavier is deeper" {
    # Every occupied level is below zero and they rise with n and l, which
    # is the shell structure the table is a table of.
    e_h = h_lda_eigs()
    for i to len(e_h) {
        assert_lt(e_h[i], 0.0)
    }
    for i to len(e_h) - 1 {
        assert_lt(e_h[i], e_h[i + 1])
    }
    e_he = he_lda_eigs()
    for i to len(e_he) {
        assert_lt(e_he[i], 0.0)
    }
    for i to len(e_he) - 1 {
        assert_lt(e_he[i], e_he[i + 1])
    }
    e_li = li_lda_eigs()
    for i to len(e_li) {
        assert_lt(e_li[i], 0.0)
    }
    for i to len(e_li) - 1 {
        assert_lt(e_li[i], e_li[i + 1])
    }
    e_be = be_lda_eigs()
    for i to len(e_be) {
        assert_lt(e_be[i], 0.0)
    }
    for i to len(e_be) - 1 {
        assert_lt(e_be[i], e_be[i + 1])
    }
    e_b = b_lda_eigs()
    for i to len(e_b) {
        assert_lt(e_b[i], 0.0)
    }
    for i to len(e_b) - 1 {
        assert_lt(e_b[i], e_b[i + 1])
    }
    e_c = c_lda_eigs()
    for i to len(e_c) {
        assert_lt(e_c[i], 0.0)
    }
    for i to len(e_c) - 1 {
        assert_lt(e_c[i], e_c[i + 1])
    }
    e_n = n_lda_eigs()
    for i to len(e_n) {
        assert_lt(e_n[i], 0.0)
    }
    for i to len(e_n) - 1 {
        assert_lt(e_n[i], e_n[i + 1])
    }
    e_o = o_lda_eigs()
    for i to len(e_o) {
        assert_lt(e_o[i], 0.0)
    }
    for i to len(e_o) - 1 {
        assert_lt(e_o[i], e_o[i + 1])
    }
    e_f = f_lda_eigs()
    for i to len(e_f) {
        assert_lt(e_f[i], 0.0)
    }
    for i to len(e_f) - 1 {
        assert_lt(e_f[i], e_f[i + 1])
    }
    e_ne = ne_lda_eigs()
    for i to len(e_ne) {
        assert_lt(e_ne[i], 0.0)
    }
    for i to len(e_ne) - 1 {
        assert_lt(e_ne[i], e_ne[i + 1])
    }
    e_na = na_lda_eigs()
    for i to len(e_na) {
        assert_lt(e_na[i], 0.0)
    }
    for i to len(e_na) - 1 {
        assert_lt(e_na[i], e_na[i + 1])
    }
    e_mg = mg_lda_eigs()
    for i to len(e_mg) {
        assert_lt(e_mg[i], 0.0)
    }
    for i to len(e_mg) - 1 {
        assert_lt(e_mg[i], e_mg[i + 1])
    }
    e_al = al_lda_eigs()
    for i to len(e_al) {
        assert_lt(e_al[i], 0.0)
    }
    for i to len(e_al) - 1 {
        assert_lt(e_al[i], e_al[i + 1])
    }
    e_si = si_lda_eigs()
    for i to len(e_si) {
        assert_lt(e_si[i], 0.0)
    }
    for i to len(e_si) - 1 {
        assert_lt(e_si[i], e_si[i + 1])
    }
    e_p = p_lda_eigs()
    for i to len(e_p) {
        assert_lt(e_p[i], 0.0)
    }
    for i to len(e_p) - 1 {
        assert_lt(e_p[i], e_p[i + 1])
    }
    e_s = s_lda_eigs()
    for i to len(e_s) {
        assert_lt(e_s[i], 0.0)
    }
    for i to len(e_s) - 1 {
        assert_lt(e_s[i], e_s[i + 1])
    }
    e_cl = cl_lda_eigs()
    for i to len(e_cl) {
        assert_lt(e_cl[i], 0.0)
    }
    for i to len(e_cl) - 1 {
        assert_lt(e_cl[i], e_cl[i + 1])
    }
    e_ar = ar_lda_eigs()
    for i to len(e_ar) {
        assert_lt(e_ar[i], 0.0)
    }
    for i to len(e_ar) - 1 {
        assert_lt(e_ar[i], e_ar[i + 1])
    }
    # and the total energy falls as the nucleus grows
    assert_lt(HE_LDA_ETOT, H_LDA_ETOT)
    assert_lt(LI_LDA_ETOT, HE_LDA_ETOT)
    assert_lt(BE_LDA_ETOT, LI_LDA_ETOT)
    assert_lt(B_LDA_ETOT, BE_LDA_ETOT)
    assert_lt(C_LDA_ETOT, B_LDA_ETOT)
    assert_lt(N_LDA_ETOT, C_LDA_ETOT)
    assert_lt(O_LDA_ETOT, N_LDA_ETOT)
    assert_lt(F_LDA_ETOT, O_LDA_ETOT)
    assert_lt(NE_LDA_ETOT, F_LDA_ETOT)
    assert_lt(NA_LDA_ETOT, NE_LDA_ETOT)
    assert_lt(MG_LDA_ETOT, NA_LDA_ETOT)
    assert_lt(AL_LDA_ETOT, MG_LDA_ETOT)
    assert_lt(SI_LDA_ETOT, AL_LDA_ETOT)
    assert_lt(P_LDA_ETOT, SI_LDA_ETOT)
    assert_lt(S_LDA_ETOT, P_LDA_ETOT)
    assert_lt(CL_LDA_ETOT, S_LDA_ETOT)
    assert_lt(AR_LDA_ETOT, CL_LDA_ETOT)
}

test "the virial theorem does NOT hold for these, and that is the point" {
    # For the exact solution of a Coulomb Hamiltonian, 2T + V = 0, so the
    # kinetic energy would be minus the total. It is not, and asserting that
    # it was is how this test started. Ekin here is the NON-interacting
    # Kohn-Sham kinetic energy, which is not the interacting one, and the
    # Levy-Perdew relation that does hold in density functional theory
    # carries a scaling term in the exchange-correlation energy.
    #
    # NIST's own numbers say so: for oxygen T is 74.116881 against -Etot of
    # 74.473077, half a percent apart. Left here as an assertion of the gap
    # rather than deleted, so nobody puts it back.
    assert_lt(1.0e-4, abs(H_LDA_EKIN + H_LDA_ETOT) / abs(H_LDA_ETOT))
    assert_lt(1.0e-4, abs(HE_LDA_EKIN + HE_LDA_ETOT) / abs(HE_LDA_ETOT))
    assert_lt(1.0e-4, abs(LI_LDA_EKIN + LI_LDA_ETOT) / abs(LI_LDA_ETOT))
    assert_lt(1.0e-4, abs(BE_LDA_EKIN + BE_LDA_ETOT) / abs(BE_LDA_ETOT))
    assert_lt(1.0e-4, abs(B_LDA_EKIN + B_LDA_ETOT) / abs(B_LDA_ETOT))
    assert_lt(1.0e-4, abs(C_LDA_EKIN + C_LDA_ETOT) / abs(C_LDA_ETOT))
    assert_lt(1.0e-4, abs(N_LDA_EKIN + N_LDA_ETOT) / abs(N_LDA_ETOT))
    assert_lt(1.0e-4, abs(O_LDA_EKIN + O_LDA_ETOT) / abs(O_LDA_ETOT))
    assert_lt(1.0e-4, abs(F_LDA_EKIN + F_LDA_ETOT) / abs(F_LDA_ETOT))
    assert_lt(1.0e-4, abs(NE_LDA_EKIN + NE_LDA_ETOT) / abs(NE_LDA_ETOT))
    assert_lt(1.0e-4, abs(NA_LDA_EKIN + NA_LDA_ETOT) / abs(NA_LDA_ETOT))
    assert_lt(1.0e-4, abs(MG_LDA_EKIN + MG_LDA_ETOT) / abs(MG_LDA_ETOT))
    assert_lt(1.0e-4, abs(AL_LDA_EKIN + AL_LDA_ETOT) / abs(AL_LDA_ETOT))
    assert_lt(1.0e-4, abs(SI_LDA_EKIN + SI_LDA_ETOT) / abs(SI_LDA_ETOT))
    assert_lt(1.0e-4, abs(P_LDA_EKIN + P_LDA_ETOT) / abs(P_LDA_ETOT))
    assert_lt(1.0e-4, abs(S_LDA_EKIN + S_LDA_ETOT) / abs(S_LDA_ETOT))
    assert_lt(1.0e-4, abs(CL_LDA_EKIN + CL_LDA_ETOT) / abs(CL_LDA_ETOT))
    assert_lt(1.0e-4, abs(AR_LDA_EKIN + AR_LDA_ETOT) / abs(AR_LDA_ETOT))
}