###
post: The OH radical in STO-3G
agent: claude

Generated by `tools/hf_reference.py` from PySCF 2.5.0. Do not edit; rerun the
generator. Settings: conv_tol 1e-12, init_guess minao, symmetry off, and each
molecule built with cart=False and with cart=True. Every stored density is
then re-converged from itself at conv_tol_grad 1e-10 or tighter, and its
`_SELF_STEP` constant says how far one PySCF step moves it (below 1e-10).

Conventions:

- atoms[natm,4] = (Z, x, y, z) in Bohr; shells[nsh,8] = (atom, l, nprim, p0,
  cart0, sph0, ncart, nsph); prims[nprim,2] = (exponent, raw bas_ctr_coeff).
  General contractions are split into segmented shells, and primitives whose
  coefficient is exactly zero are dropped.
- `_unit` matrices are unit-Cartesian: PySCF's Cartesian integrals divided by
  F_i F_j, so every component has self-overlap 1. `_sph` matrices are PySCF's
  spherical ones, and `_m()` is M with X_sph = M^T X_unit M.
- Packed-8 ERIs in chemist's notation: ij = i(i+1)/2 + j for i >= j, and
  ijkl = ij(ij+1)/2 + kl for ij >= kl. Full when the basis has at most 19
  functions; otherwise the pair diagonal and 256 seeded samples. Both carry
  _SUMSQ = sum Gp^2 and _WSUM = sum w_k Gp_k, where w_k = s_k / (2^31 - 1),
  s_k = 16807 s_(k-1) mod (2^31 - 1) and s_0 = the module's _SEED_WSUM (20260923).
- Energies: `_E_<CFG>_DEFAULT` is the state PySCF reaches from the minao
  guess, `_E_<CFG>_STABLE` the lowest state reached by restarting from the
  stability-analysis orbitals until the energy stops dropping, and `_E_<CFG>`
  repeats the default. Densities (total for RHF, alpha and beta otherwise) are
  in the same basis convention as the matrices beside them: a unit-Cartesian
  density is F_i F_j times PySCF's. A `_stable_` density appears only when the
  stable state differs from the default.
- In a module with both bases, spherical names carry no suffix and Cartesian
  ones end in `_cart` / `_CART`.
###

# Seed of the LCG weights in every _WSUM checksum.
OH_STO3G_SEED_WSUM = 20260923

# ---- helpers for the tests below ----

#| Reshape a flat row-major vector into r x c.
def oh_sto3g_mat(v: Arr, r, c) -> Arr {
    a = zeros(r, c)
    for i to r {
        for j to c {
            a[i, j] = v[i * c + j]
        }
    }
    return a
}

def oh_sto3g_eye(n) -> Arr {
    a = zeros(n, n)
    for i to n {
        a[i, i] = 1.0
    }
    return a
}

#| Largest |A_ij - A_ji| / (1 + |A_ij|): a kinetic entry near 100 carries
#| rounding of about 1e-14, so the gap is measured relative to the entry.
def oh_sto3g_asym(a: Arr, n) {
    w = 0.0
    for i to n {
        for j to n {
            w = max(w, abs(a[i, j] - a[j, i]) / (1.0 + abs(a[i, j])))
        }
    }
    return w
}

#| Largest |A_ii - 1|.
def oh_sto3g_diag_dev(a: Arr, n) {
    w = 0.0
    for i to n {
        w = max(w, abs(a[i, i] - 1.0))
    }
    return w
}

#| sum_ij A_ij B_ij: Tr(D S) is the electron count.
def oh_sto3g_trace2(a: Arr, b: Arr, n) {
    t = 0.0
    for i to n {
        for j to n {
            t = t + a[i, j] * b[i, j]
        }
    }
    return t
}

#| pr[i, j] = the packed pair index of (i, j), in either order.
def oh_sto3g_pair_table(n) -> Arr {
    pr = zeros(n, n)
    for i to n {
        for j to i + 1 {
            pr[i, j] = i * (i + 1) / 2 + j
            pr[j, i] = pr[i, j]
        }
    }
    return pr
}

#| One- plus two-electron energy from spin densities and packed (ij|kl):
#| sum Dt h + 1/2 sum (ij|kl) [Dt_ij Dt_kl - Da_ik Da_jl - Db_ik Db_jl].
#| For RHF pass Da = Db = D/2. Written over all n^4 quartets because it is
#| a check, not a J/K build.
def oh_sto3g_e_elec(h: Arr, g: Arr, da: Arr, db: Arr, n) {
    pr = oh_sto3g_pair_table(n)
    e = 0.0
    for i to n {
        for j to n {
            e = e + (da[i, j] + db[i, j]) * h[i, j]
        }
    }
    e2 = 0.0
    for i to n {
        for j to n {
            pij = pr[i, j]
            dij = da[i, j] + db[i, j]
            for k to n {
                for l to n {
                    pkl = pr[k, l]
                    a = max(pij, pkl)
                    b = min(pij, pkl)
                    x = g[a * (a + 1) / 2 + b]
                    e2 = e2 + x * (dij * (da[k, l] + db[k, l]) - da[i, k] * da[j, l] - db[i, k] * db[j, l])
                }
            }
        }
    }
    return e + 0.5 * e2
}

def oh_sto3g_sumsq(g: Arr) {
    t = 0.0
    for k to len(g) {
        t = t + g[k] * g[k]
    }
    return t
}

#| sum_k w_k g_k with w_k = s_k / (2^31 - 1), s_k = 16807 s_(k-1) mod (2^31 - 1),
#| s_0 = seed. Every product stays below 2^53, so the weights are exact.
def oh_sto3g_wsum(g: Arr, seed) {
    s = seed
    t = 0.0
    for k to len(g) {
        s = mod(16807.0 * s, 2147483647.0)
        t = t + s / 2147483647.0 * g[k]
    }
    return t
}

#| Largest (ij|kl)^2 - (ij|ij)(kl|kl) over the packed set, relative to the
#| bound. Cauchy-Schwarz says it is at most 0.
def oh_sto3g_schwarz(g: Arr, n) {
    npair = n * (n + 1) / 2
    w = -1.0
    for a to npair {
        for b to a + 1 {
            x = g[a * (a + 1) / 2 + b]
            q = g[a * (a + 1) / 2 + a] * g[b * (b + 1) / 2 + b]
            w = max(w, (x * x - q) / (q + 1.0e-300))
        }
    }
    return w
}

#| The same bound on sampled entries, with the pair diagonal given.
def oh_sto3g_sample_schwarz(idx: Arr, val: Arr, gd: Arr) {
    w = -1.0
    for k to len(idx) {
        a = floor((sqrt(8.0 * idx[k] + 1.0) - 1.0) / 2.0)
        if a * (a + 1) / 2 > idx[k] {
            a = a - 1
        }
        if (a + 1) * (a + 2) / 2 <= idx[k] {
            a = a + 1
        }
        b = idx[k] - a * (a + 1) / 2
        q = gd[a] * gd[b]
        w = max(w, (val[k] * val[k] - q) / (q + 1.0e-300))
    }
    return w
}

#| Largest |M^T X M - Y|.
def oh_sto3g_mtxm_dev(m: Arr, x: Arr, y: Arr, nc, ns) {
    w = 0.0
    for a to ns {
        for b to ns {
            t = 0.0
            for i to nc {
                for j to nc {
                    t = t + m[i, a] * x[i, j] * m[j, b]
                }
            }
            w = max(w, abs(t - y[a, b]))
        }
    }
    return w
}

#| Count what is wrong with a shells pack: offsets that do not run on, a
#| component count that does not match l, a primitive range out of bounds.
def oh_sto3g_pack_bad(sh: Arr, nsh, nprim, natm) {
    bad = 0
    p0 = 0
    c0 = 0
    s0 = 0
    for k to nsh {
        l = sh[k, 1]
        if sh[k, 0] < 0 or sh[k, 0] >= natm {
            bad = bad + 1
        }
        if sh[k, 3] != p0 or sh[k, 4] != c0 or sh[k, 5] != s0 {
            bad = bad + 1
        }
        if sh[k, 6] != (l + 1) * (l + 2) / 2 or sh[k, 7] != 2 * l + 1 {
            bad = bad + 1
        }
        p0 = p0 + sh[k, 2]
        c0 = c0 + sh[k, 6]
        s0 = s0 + sh[k, 7]
    }
    if p0 != nprim {
        bad = bad + 1
    }
    return bad
}

# ======== OH_STO3G: sto-3g, s and p only, so Cartesian = spherical ========

OH_STO3G_NATOM = 2
OH_STO3G_NSH = 4
OH_STO3G_NPRIMTOT = 12
OH_STO3G_NAO_CART = 6
OH_STO3G_NAO_SPH = 6
OH_STO3G_E_NUC = 4.365698347282665

#| (Z, x, y, z) in Bohr, one atom per row.
def oh_sto3g_atoms_flat() -> Arr {
    return [8.0, 0.0, 0.0, 0.0,
        1.0, 0.0, 0.0, 1.8324674229907405]
}

def oh_sto3g_atoms() -> Arr {
    return oh_sto3g_mat(oh_sto3g_atoms_flat(), 2, 4)
}

#| (atom, l, nprim, p0, cart0, sph0, ncart, nsph); general contractions split.
def oh_sto3g_shells_flat() -> Arr {
    return [0, 0, 3, 0,
        0, 0, 1, 1,
        0, 0, 3, 3,
        1, 1, 1, 1,
        0, 1, 3, 6,
        2, 2, 3, 3,
        1, 0, 3, 9,
        5, 5, 1, 1]
}

def oh_sto3g_shells() -> Arr {
    return oh_sto3g_mat(oh_sto3g_shells_flat(), 4, 8)
}

#| (exponent, raw bas_ctr_coeff).
def oh_sto3g_prims_flat() -> Arr {
    return [130.70932, 0.15432896872315646, 23.808861, 0.5353281355709528,
        6.4436083, 0.4446345363213081, 5.0331513, -0.09996722865667348,
        1.1695961, 0.3995128246314789, 0.380389, 0.7001154605920802,
        5.0331513, 0.15591626847815623, 1.1695961, 0.6076837140686134,
        0.380389, 0.3919573861742421, 3.42525091, 0.15432897070298393,
        0.62391373, 0.5353281424384733, 0.1688554, 0.44463454202535485]
}

def oh_sto3g_prims() -> Arr {
    return oh_sto3g_mat(oh_sto3g_prims_flat(), 12, 2)
}

#| F per Cartesian AO: PySCF-Cartesian = F x unit-Cartesian.
def oh_sto3g_fc() -> Arr {
    return [1.0, 1.0, 1.0, 1.0,
        1.0, 1.0]
}

#| Overlap, unit-Cartesian, [6,6].
def oh_sto3g_s_unit_flat() -> Arr {
    return [1.0, 0.2367039365108476, 0.0, 0.0,
        0.0, 0.052398646083327505, 0.23670393651084762, 1.0,
        0.0, 0.0, 0.0, 0.4666392346755538,
        0.0, 0.0, 1.0, 0.0,
        0.0, 0.0, 0.0, 0.0,
        0.0, 1.0, 0.0, 0.0,
        0.0, 0.0, 0.0, 0.0,
        1.0, 0.3890721745846377, 0.05239864608332752, 0.4666392346755538,
        0.0, 0.0, 0.38907217458463766, 1.0000000000000002]
}

def oh_sto3g_s_unit() -> Arr {
    return oh_sto3g_mat(oh_sto3g_s_unit_flat(), 6, 6)
}

#| Kinetic energy, unit-Cartesian.
def oh_sto3g_t_unit_flat() -> Arr {
    return [29.003199945539567, -0.16801093931649205, 0.0, 0.0,
        0.0, -0.0033550804012269907, -0.16801093931649186, 0.8081279549303472,
        0.0, 0.0, 0.0, 0.12268625732597083,
        0.0, 0.0, 2.5287311981947633, 0.0,
        0.0, 0.0, 0.0, 0.0,
        0.0, 2.5287311981947633, 0.0, 0.0,
        0.0, 0.0, 0.0, 0.0,
        2.5287311981947633, 0.27474085381612845, -0.003355080401227047, 0.12268625732597085,
        0.0, 0.0, 0.27474085381612856, 0.7600318835666091]
}

def oh_sto3g_t_unit() -> Arr {
    return oh_sto3g_mat(oh_sto3g_t_unit_flat(), 6, 6)
}

#| Nuclear attraction, unit-Cartesian.
def oh_sto3g_v_unit_flat() -> Arr {
    return [-61.16416719883297, -7.312251692745978, 0.0, 0.0,
        -0.01512593349936972, -1.6661038529040273, -7.312251692745979, -9.590953268292838,
        0.0, 0.0, -0.17872192535221898, -3.571548926817574,
        0.0, 0.0, -9.477622837435154, 0.0,
        0.0, 0.0, 0.0, 0.0,
        0.0, -9.477622837435154, 0.0, 0.0,
        -0.01512593349936972, -0.17872192535221895, 0.0, 0.0,
        -9.599029618528155, -2.6998089565835373, -1.6661038529040275, -3.5715489268175733,
        0.0, 0.0, -2.6998089565835373, -5.44118967335766]
}

def oh_sto3g_v_unit() -> Arr {
    return oh_sto3g_mat(oh_sto3g_v_unit_flat(), 6, 6)
}

# Packed-8 length n(n+1)/2 (n(n+1)/2 + 1)/2 for n = 6.
OH_STO3G_ERI8_CART_LEN = 231
# sum of Gp^2 (math.fsum).
OH_STO3G_ERI8_CART_SUMSQ = 38.806415018927304
# sum of w_k Gp_k with the LCG weights from OH_STO3G_SEED_WSUM (math.fsum).
OH_STO3G_ERI8_CART_WSUM = 13.234751853615846
#| Every packed-8 two-electron integral, unit-Cartesian, chemist's notation.
def oh_sto3g_eri8_cart() -> Arr {
    return [4.785065404705503, 0.7413803519734079, 0.13687338535438834, 1.1189468663424704,
        0.2566333947309737, 0.8172063215260582, 0.0, 0.0,
        0.0, 0.024477412258099275, 0.0, 0.0,
        0.0, 0.03780860741636099, 0.18051839210463208, 1.1158138121524277,
        0.25668398581010293, 0.8170226053209142, 0.0, 0.0,
        0.8801590933750454, 0.0, 0.0, 0.0,
        0.0, 0.0, 0.0, 0.024477412258099275,
        0.0, 0.0, 0.0, 0.0,
        0.0, 0.0, 0.03780860741636099, 0.18051839210463208,
        0.0, 0.0, 0.0, 0.0,
        0.0, 0.0, 0.0, 0.0,
        0.04744444511838415, 1.1158138121524277, 0.25668398581010293, 0.8170226053209142,
        0.0, 0.0, 0.7852702031382774, 0.0,
        0.0, 0.0, 0.8801590933750454, 0.0,
        0.0, 0.0, 0.0, 0.0,
        0.0, 0.0, 0.0, 0.0,
        0.0, 0.024477412258099275, 0.0, 0.0,
        0.0, 0.0, 0.0, 0.0,
        0.0, 0.0, 0.0, 0.0,
        0.03780860741636099, 0.18051839210463208, 0.0, 0.0,
        0.0, 0.0, 0.0, 0.0,
        0.0, 0.0, 0.0, 0.0,
        0.0, 0.0, 0.04744444511838415, 0.0,
        0.0, 0.0, 0.0, 0.0,
        0.0, 0.0, 0.0, 0.0,
        0.0, 0.0, 0.0, 0.0,
        0.04744444511838415, 1.1158138121524277, 0.25668398581010293, 0.8170226053209142,
        0.0, 0.0, 0.7852702031382774, 0.0,
        0.0, 0.0, 0.7852702031382774, 0.0,
        0.0, 0.0, 0.0, 0.8801590933750452,
        0.16670799515033932, 0.03050727102728575, 0.05679728889648121, 0.0,
        0.0, 0.056689144511863694, 0.0, 0.0,
        0.0, 0.056689144511863694, 0.0014048562458150692, 0.0026238090030180753,
        0.0, 0.0, 0.05702090549556107, 0.006888013858016352,
        0.3950192042648373, 0.09168746143027433, 0.3220205092719595, 0.0,
        0.0, 0.31664025835965304, 0.0, 0.0,
        0.0, 0.31664025835965304, 0.00686339010170292, 0.05186020672244012,
        0.0, 0.0, 0.3324787902604658, 0.020815203175604247,
        0.15551806961466588, 0.0, 0.0, 0.0,
        0.00972046278761465, 0.05612278057229543, 0.0, 0.0,
        0.0, 0.0, 0.0, 0.0,
        0.0, 0.010847019429957139, 0.0, 0.0,
        0.0, 0.0, 0.022482131311621682, 0.0,
        0.0, 0.0, 0.0, 0.0,
        0.0, 0.00972046278761465, 0.05612278057229543, 0.0,
        0.0, 0.0, 0.0, 0.0,
        0.010847019429957139, 0.0, 0.0, 0.0,
        0.0, 0.022482131311621682, 0.2799232428806428, 0.06582979082475705,
        0.246927796934799, 0.0, 0.0, 0.23565800964829556,
        0.0, 0.0, 0.0, 0.23565800964829556,
        0.014051946443301891, 0.09392010799101728, 0.0, 0.0,
        0.2691949803029399, 0.015622979350690843, 0.1448489716267476, 0.0,
        0.0, 0.15587757144951903, 0.52611559829252, 0.12410307859758536,
        0.49880109596465544, 0.0, 0.0, 0.478508093596294,
        0.0, 0.0, 0.0, 0.478508093596294,
        0.01231921408258197, 0.1271991367302265, 0.0, 0.0,
        0.5387749244893735, 0.028444738068592913, 0.2996495347623702, 0.0,
        0.0, 0.31878174058316094, 0.7746059439198978]
}

def oh_sto3g_s_sph() -> Arr {
    return oh_sto3g_s_unit()
}

def oh_sto3g_t_sph() -> Arr {
    return oh_sto3g_t_unit()
}

def oh_sto3g_v_sph() -> Arr {
    return oh_sto3g_v_unit()
}

#| Only s and p: the spherical and Cartesian sets are the same.
def oh_sto3g_eri8_sph() -> Arr {
    return oh_sto3g_eri8_cart()
}

def oh_sto3g_m() -> Arr {
    return oh_sto3g_eye(6)
}

# ---- OH_STO3G UHF: nalpha 5, nbeta 4; default state is stable ----

OH_STO3G_E_UHF = -74.36263751866602
OH_STO3G_E_UHF_DEFAULT = -74.36263751866602
OH_STO3G_E_UHF_STABLE = -74.36263751866602
# 1: no restart from the stability orbitals lowers the minao state; 0: one does; -1: not analyzed (stability restart).
OH_STO3G_UHF_IS_STABLE = 1
OH_STO3G_UHF_NALPHA = 5
OH_STO3G_UHF_NBETA = 4
# <S^2> of the default UHF state.
OH_STO3G_UHF_S2 = 0.7532558435794732
OH_STO3G_UHF_S2_STABLE = 0.7532558435794732

# Max |D' - D| after one PySCF step (get_fock, eig, get_occ, make_rdm1) from the default density below.
OH_STO3G_UHF_SELF_STEP = 9.811040868612508e-13
#| Converged alpha density, default state.
def oh_sto3g_uhf_da_flat() -> Arr {
    return [1.0575294541460254, -0.24459513045592968, 2.616240415193523e-18, -1.2876305267926837e-17,
        0.03805843278518075, -0.003718819682606279, -0.24459513045592968, 1.0663059038423108,
        6.445314493215387e-18, 7.484204570286354e-17, -0.21899768769535724, -0.07609476868158084,
        2.616240415193523e-18, 6.445314493215387e-18, 1.0, -5.674353501794484e-17,
        1.3963126939833052e-17, -7.501378731783194e-17, -1.2876305267926837e-17, 7.484204570286354e-17,
        -5.674353501794484e-17, 0.9999999999999993, -9.128402109960169e-17, -1.0703475955170491e-16,
        0.03805843278518075, -0.21899768769535724, 1.3963126939833052e-17, -9.128402109960169e-17,
        0.5020057344070915, 0.3670301067525014, -0.003718819682606279, -0.07609476868158084,
        -7.501378731783194e-17, -1.0703475955170491e-16, 0.3670301067525014, 0.2757570959942884]
}

def oh_sto3g_uhf_da() -> Arr {
    return oh_sto3g_mat(oh_sto3g_uhf_da_flat(), 6, 6)
}

#| Converged beta density, default state.
def oh_sto3g_uhf_db_flat() -> Arr {
    return [1.0563276066814915, -0.23862852854812514, -1.6643157361854792e-17, -5.9563241023984e-18,
        0.04650836754392135, -0.008588029915871568, -0.23862852854812514, 1.0401592812961928,
        6.81323183823724e-17, 2.2893180286251397e-17, -0.2558645678360472, -0.06886697107712816,
        -1.6643157361854792e-17, 6.81323183823724e-17, 0.9987027577236065, 0.035993880575309785,
        6.534026137941466e-17, -1.2670148365143315e-17, -5.9563241023984e-18, 2.2893180286251397e-17,
        0.035993880575309785, 0.0012972422763933254, 7.05497509338411e-18, 7.22987813097209e-18,
        0.04650836754392135, -0.2558645678360472, 6.534026137941466e-17, 7.05497509338411e-18,
        0.45003180107832347, 0.3764755667925659, -0.008588029915871568, -0.06886697107712816,
        -1.2670148365143315e-17, 7.22987813097209e-18, 0.3764755667925659, 0.33866966382226005]
}

def oh_sto3g_uhf_db() -> Arr {
    return oh_sto3g_mat(oh_sto3g_uhf_db_flat(), 6, 6)
}

# ---- OH_STO3G ROHF: nalpha 5, nbeta 4; default state is stable ----

OH_STO3G_E_ROHF = -74.36153072613816
OH_STO3G_E_ROHF_DEFAULT = -74.36153072613816
OH_STO3G_E_ROHF_STABLE = -74.36153072613816
# 1: no restart from the stability orbitals lowers the minao state; 0: one does; -1: not analyzed (stability restart).
OH_STO3G_ROHF_IS_STABLE = 1
OH_STO3G_ROHF_NALPHA = 5
OH_STO3G_ROHF_NBETA = 4

# Max |D' - D| after one PySCF step (get_fock, eig, get_occ, make_rdm1) from the default density below.
OH_STO3G_ROHF_SELF_STEP = 1.2753686995381486e-14
#| Converged alpha density, default state.
def oh_sto3g_rohf_da_flat() -> Arr {
    return [1.0569686713607962, -0.2417984195236519, -1.0726155470768565e-17, 2.624657416430829e-19,
        0.04214733450771532, -0.006051587389312103, -0.2417984195236519, 1.0540600991324967,
        -1.895385705390261e-17, -1.2592709540353528e-17, -0.23711487889936533, -0.07276221169742464,
        -1.0726155470768565e-17, -1.895385705390261e-17, 1.0000000000000002, -1.8381424694486874e-16,
        -1.5077376788154157e-16, 2.7377935679730006e-16, 2.624657416430829e-19, -1.2592709540353528e-17,
        -1.8381424694486874e-16, 1.0, 2.158406282699707e-16, 4.680122815975629e-17,
        0.04214733450771532, -0.23711487889936533, -1.5077376788154157e-16, 2.158406282699707e-16,
        0.4752325168369421, 0.3729430921632892, -0.006051587389312103, -0.07276221169742464,
        2.7377935679730006e-16, 4.680122815975629e-17, 0.3729430921632892, 0.3065460239591688]
}

def oh_sto3g_rohf_da() -> Arr {
    return oh_sto3g_mat(oh_sto3g_rohf_da_flat(), 6, 6)
}

#| Converged beta density, default state.
def oh_sto3g_rohf_db_flat() -> Arr {
    return [1.0569686713607962, -0.2417984195236519, -1.0944766864602213e-17, 6.759481997389739e-18,
        0.04214733450771532, -0.006051587389312103, -0.2417984195236519, 1.0540600991324967,
        -1.7377631999930662e-17, -5.943729352467973e-17, -0.23711487889936533, -0.07276221169742464,
        -1.0944766864602213e-17, -1.7377631999930662e-17, 0.9988690946571652, 0.033609915143310624,
        -1.4875159041833716e-16, 2.70715118382049e-16, 6.759481997389739e-18, -5.943729352467973e-17,
        0.033609915143310624, 0.0011309053428349934, 1.5574257198823354e-16, 1.378687907705098e-16,
        0.04214733450771532, -0.23711487889936533, -1.4875159041833716e-16, 1.5574257198823354e-16,
        0.4752325168369421, 0.3729430921632892, -0.006051587389312103, -0.07276221169742464,
        2.70715118382049e-16, 1.378687907705098e-16, 0.3729430921632892, 0.3065460239591688]
}

def oh_sto3g_rohf_db() -> Arr {
    return oh_sto3g_mat(oh_sto3g_rohf_db_flat(), 6, 6)
}

test "oh_sto3g packs" {
    assert_near(oh_sto3g_pack_bad(oh_sto3g_shells(), OH_STO3G_NSH, OH_STO3G_NPRIMTOT, OH_STO3G_NATOM), 0.0, 0.0)
    # nao_cart - nao_sph = sum over shells of ncart - nsph: 1 per d shell, 3 per f shell
    assert_near(OH_STO3G_NAO_CART - OH_STO3G_NAO_SPH, 0.0, 0.0)
    assert_near(len(oh_sto3g_fc()), OH_STO3G_NAO_CART, 0.0)
}

test "oh_sto3g unit one-electron integrals are symmetric with a unit overlap diagonal" {
    assert_near(oh_sto3g_asym(oh_sto3g_s_unit(), 6), 0.0, 1.0e-14)
    assert_near(oh_sto3g_asym(oh_sto3g_t_unit(), 6), 0.0, 1.0e-14)
    assert_near(oh_sto3g_asym(oh_sto3g_v_unit(), 6), 0.0, 1.0e-14)
    assert_near(oh_sto3g_diag_dev(oh_sto3g_s_unit(), 6), 0.0, 1.0e-14)
}

test "oh_sto3g unit packed ERIs match their checksums and the Schwarz bound" {
    assert_near(len(oh_sto3g_eri8_cart()), OH_STO3G_ERI8_CART_LEN, 0.0)
    assert_near(oh_sto3g_sumsq(oh_sto3g_eri8_cart()) / OH_STO3G_ERI8_CART_SUMSQ, 1.0, 1.0e-12)
    assert_near(oh_sto3g_wsum(oh_sto3g_eri8_cart(), OH_STO3G_SEED_WSUM), OH_STO3G_ERI8_CART_WSUM, 1.0e-12 * OH_STO3G_ERI8_CART_SUMSQ)
    assert_lt(oh_sto3g_schwarz(oh_sto3g_eri8_cart(), 6), 1.0e-12)
}

test "oh_sto3g uhf energies, labels and densities" {
    assert_lt(OH_STO3G_E_UHF_STABLE, OH_STO3G_E_UHF_DEFAULT + 1.0e-12)
    assert_lt(0.75 - 1.0e-10, OH_STO3G_UHF_S2)
    assert_lt(0.75 - 1.0e-10, OH_STO3G_UHF_S2_STABLE)
    assert_lt(OH_STO3G_UHF_SELF_STEP, 1.0e-10)
    assert_near(oh_sto3g_asym(oh_sto3g_uhf_da(), 6) + oh_sto3g_asym(oh_sto3g_uhf_db(), 6), 0.0, 1.0e-12)
    assert_near(oh_sto3g_trace2(oh_sto3g_uhf_da(), oh_sto3g_s_unit(), 6) + oh_sto3g_trace2(oh_sto3g_uhf_db(), oh_sto3g_s_unit(), 6), 9.0, 1.0e-10)
    assert_near(oh_sto3g_e_elec(oh_sto3g_t_unit() + oh_sto3g_v_unit(), oh_sto3g_eri8_cart(), oh_sto3g_uhf_da(), oh_sto3g_uhf_db(), 6) + OH_STO3G_E_NUC, OH_STO3G_E_UHF_DEFAULT, 1.0e-10)
}

test "oh_sto3g rohf energies, labels and densities" {
    assert_lt(OH_STO3G_E_ROHF_STABLE, OH_STO3G_E_ROHF_DEFAULT + 1.0e-12)
    assert_lt(OH_STO3G_ROHF_SELF_STEP, 1.0e-10)
    assert_near(oh_sto3g_asym(oh_sto3g_rohf_da(), 6) + oh_sto3g_asym(oh_sto3g_rohf_db(), 6), 0.0, 1.0e-12)
    assert_near(oh_sto3g_trace2(oh_sto3g_rohf_da(), oh_sto3g_s_unit(), 6) + oh_sto3g_trace2(oh_sto3g_rohf_db(), oh_sto3g_s_unit(), 6), 9.0, 1.0e-10)
    assert_near(oh_sto3g_e_elec(oh_sto3g_t_unit() + oh_sto3g_v_unit(), oh_sto3g_eri8_cart(), oh_sto3g_rohf_da(), oh_sto3g_rohf_db(), 6) + OH_STO3G_E_NUC, OH_STO3G_E_ROHF_DEFAULT, 1.0e-10)
}

test "oh_sto3g E_uhf <= E_rohf" {
    assert_lt(OH_STO3G_E_UHF_DEFAULT, OH_STO3G_E_ROHF_DEFAULT + 1.0e-12)
    assert_lt(OH_STO3G_E_UHF_STABLE, OH_STO3G_E_ROHF_STABLE + 1.0e-12)
}