###
post: Helium: STO-3G, the bare dication, and cc-pVDZ
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.
HE_SEED_WSUM = 20260923

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

#| Reshape a flat row-major vector into r x c.
def he_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 he_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 he_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 he_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 he_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 he_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 he_e_elec(h: Arr, g: Arr, da: Arr, db: Arr, n) {
    pr = he_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 he_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 he_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 he_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 he_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 he_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 he_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
}

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

HE_STO3G_NATOM = 1
HE_STO3G_NSH = 1
HE_STO3G_NPRIMTOT = 3
HE_STO3G_NAO_CART = 1
HE_STO3G_NAO_SPH = 1
HE_STO3G_E_NUC = 0.0

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

def he_sto3g_atoms() -> Arr {
    return he_mat(he_sto3g_atoms_flat(), 1, 4)
}

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

def he_sto3g_shells() -> Arr {
    return he_mat(he_sto3g_shells_flat(), 1, 8)
}

#| (exponent, raw bas_ctr_coeff).
def he_sto3g_prims_flat() -> Arr {
    return [6.36242139, 0.15432897044824204, 1.158923, 0.5353281415548381,
        0.31364979, 0.4446345412914224]
}

def he_sto3g_prims() -> Arr {
    return he_mat(he_sto3g_prims_flat(), 3, 2)
}

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

#| Overlap, unit-Cartesian, [1,1].
def he_sto3g_s_unit_flat() -> Arr {
    return [1.0000000000000002]
}

def he_sto3g_s_unit() -> Arr {
    return he_mat(he_sto3g_s_unit_flat(), 1, 1)
}

#| Kinetic energy, unit-Cartesian.
def he_sto3g_t_unit_flat() -> Arr {
    return [1.411763182618484]
}

def he_sto3g_t_unit() -> Arr {
    return he_mat(he_sto3g_t_unit_flat(), 1, 1)
}

#| Nuclear attraction, unit-Cartesian.
def he_sto3g_v_unit_flat() -> Arr {
    return [-3.343511632756007]
}

def he_sto3g_v_unit() -> Arr {
    return he_mat(he_sto3g_v_unit_flat(), 1, 1)
}

# Packed-8 length n(n+1)/2 (n(n+1)/2 + 1)/2 for n = 1.
HE_STO3G_ERI8_CART_LEN = 1
# sum of Gp^2 (math.fsum).
HE_STO3G_ERI8_CART_SUMSQ = 1.1145298174583458
# sum of w_k Gp_k with the LCG weights from HE_SEED_WSUM (math.fsum).
HE_STO3G_ERI8_CART_WSUM = 0.6011915020908758
#| Every packed-8 two-electron integral, unit-Cartesian, chemist's notation.
def he_sto3g_eri8_cart() -> Arr {
    return [1.0557129427350722]
}

def he_sto3g_s_sph() -> Arr {
    return he_sto3g_s_unit()
}

def he_sto3g_t_sph() -> Arr {
    return he_sto3g_t_unit()
}

def he_sto3g_v_sph() -> Arr {
    return he_sto3g_v_unit()
}

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

def he_sto3g_m() -> Arr {
    return he_eye(1)
}

# ---- HE_STO3G RHF: nalpha 1, nbeta 1; default state is stable ----

HE_STO3G_E_RHF = -2.807783957539974
HE_STO3G_E_RHF_DEFAULT = -2.807783957539974
HE_STO3G_E_RHF_STABLE = -2.807783957539974
# 1: no restart from the stability orbitals lowers the minao state; 0: one does; -1: not analyzed (no virtual orbitals: nothing to rotate).
HE_STO3G_RHF_IS_STABLE = 1
HE_STO3G_RHF_NALPHA = 1
HE_STO3G_RHF_NBETA = 1

# Max |D' - D| after one PySCF step (get_fock, eig, get_occ, make_rdm1) from the default density below.
HE_STO3G_RHF_SELF_STEP = 0.0
#| Converged total density (alpha + beta), default state.
def he_sto3g_rhf_d_flat() -> Arr {
    return [2.0]
}

def he_sto3g_rhf_d() -> Arr {
    return he_mat(he_sto3g_rhf_d_flat(), 1, 1)
}

# ---- HE_STO3G UHF: nalpha 1, nbeta 1; default state is stable ----

HE_STO3G_E_UHF = -2.807783957539974
HE_STO3G_E_UHF_DEFAULT = -2.807783957539974
HE_STO3G_E_UHF_STABLE = -2.807783957539974
# 1: no restart from the stability orbitals lowers the minao state; 0: one does; -1: not analyzed (no virtual orbitals: nothing to rotate).
HE_STO3G_UHF_IS_STABLE = 1
HE_STO3G_UHF_NALPHA = 1
HE_STO3G_UHF_NBETA = 1
# <S^2> of the default UHF state.
HE_STO3G_UHF_S2 = -2.220446049250313e-16
HE_STO3G_UHF_S2_STABLE = -2.220446049250313e-16

# Max |D' - D| after one PySCF step (get_fock, eig, get_occ, make_rdm1) from the default density below.
HE_STO3G_UHF_SELF_STEP = 0.0
#| Converged alpha density, default state.
def he_sto3g_uhf_da_flat() -> Arr {
    return [1.0]
}

def he_sto3g_uhf_da() -> Arr {
    return he_mat(he_sto3g_uhf_da_flat(), 1, 1)
}

#| Converged beta density, default state.
def he_sto3g_uhf_db_flat() -> Arr {
    return [1.0]
}

def he_sto3g_uhf_db() -> Arr {
    return he_mat(he_sto3g_uhf_db_flat(), 1, 1)
}

# ---- HE_STO3G DICATION_RHF: nalpha 0, nbeta 0; default state is stable ----

HE_STO3G_E_DICATION_RHF = 0.0
HE_STO3G_E_DICATION_RHF_DEFAULT = 0.0
HE_STO3G_E_DICATION_RHF_STABLE = 0.0
# 1: no restart from the stability orbitals lowers the minao state; 0: one does; -1: not analyzed (one electron: E is the lowest H_core level).
HE_STO3G_DICATION_RHF_IS_STABLE = 1
HE_STO3G_DICATION_RHF_NALPHA = 0
HE_STO3G_DICATION_RHF_NBETA = 0

# Max |D' - D| after one PySCF step (get_fock, eig, get_occ, make_rdm1) from the default density below.
HE_STO3G_DICATION_RHF_SELF_STEP = 0.0
#| Converged total density (alpha + beta), default state.
def he_sto3g_dication_rhf_d_flat() -> Arr {
    return [0.0]
}

def he_sto3g_dication_rhf_d() -> Arr {
    return he_mat(he_sto3g_dication_rhf_d_flat(), 1, 1)
}

# Lowest eigenvalue of H_core C = S C e.
HE_STO3G_EPS1E = -1.931748450137523
# (11|11), pyscf/gto/test/test_mole.py:942.
HE_STO3G_ERI_1111 = 1.0557129427350722
# ======== HE_CCPVDZ: cc-pvdz, s and p only, so Cartesian = spherical ========

HE_CCPVDZ_NATOM = 1
HE_CCPVDZ_NSH = 3
HE_CCPVDZ_NPRIMTOT = 5
HE_CCPVDZ_NAO_CART = 5
HE_CCPVDZ_NAO_SPH = 5
HE_CCPVDZ_E_NUC = 0.0

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

def he_ccpvdz_atoms() -> Arr {
    return he_mat(he_ccpvdz_atoms_flat(), 1, 4)
}

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

def he_ccpvdz_shells() -> Arr {
    return he_mat(he_ccpvdz_shells_flat(), 3, 8)
}

#| (exponent, raw bas_ctr_coeff).
def he_ccpvdz_prims_flat() -> Arr {
    return [38.36, 0.04017607538324854, 5.77, 0.2613680747694884,
        1.24, 0.7930712394954358, 0.2976, 1.0,
        1.275, 1.0]
}

def he_ccpvdz_prims() -> Arr {
    return he_mat(he_ccpvdz_prims_flat(), 5, 2)
}

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

#| Overlap, unit-Cartesian, [5,5].
def he_ccpvdz_s_unit_flat() -> Arr {
    return [1.0, 0.6341677316935693, 0.0, 0.0,
        0.0, 0.6341677316935693, 1.0, 0.0,
        0.0, 0.0, 0.0, 0.0,
        1.0000000000000004, 0.0, 0.0, 0.0,
        0.0, 0.0, 1.0000000000000004, 0.0,
        0.0, 0.0, 0.0, 0.0,
        1.0000000000000004]
}

def he_ccpvdz_s_unit() -> Arr {
    return he_mat(he_ccpvdz_s_unit_flat(), 5, 5)
}

#| Kinetic energy, unit-Cartesian.
def he_ccpvdz_t_unit_flat() -> Arr {
    return [2.9224349544537818, 0.46665999220234144, 0.0, 0.0,
        0.0, 0.4666599922023412, 0.44639999999999996, 0.0,
        0.0, 0.0, 0.0, 0.0,
        3.187500000000001, 0.0, 0.0, 0.0,
        0.0, 0.0, 3.187500000000001, 0.0,
        0.0, 0.0, 0.0, 0.0,
        3.187500000000001]
}

def he_ccpvdz_t_unit() -> Arr {
    return he_mat(he_ccpvdz_t_unit_flat(), 5, 5)
}

#| Nuclear attraction, unit-Cartesian.
def he_ccpvdz_v_unit_flat() -> Arr {
    return [-4.802687652973955, -2.0124487537790783, 0.0, 0.0,
        0.0, -2.0124487537790783, -1.7410711381518775, 0.0,
        0.0, 0.0, 0.0, 0.0,
        -2.402502709564773, 0.0, 0.0, 0.0,
        0.0, 0.0, -2.402502709564773, 0.0,
        0.0, 0.0, 0.0, 0.0,
        -2.402502709564773]
}

def he_ccpvdz_v_unit() -> Arr {
    return he_mat(he_ccpvdz_v_unit_flat(), 5, 5)
}

# Packed-8 length n(n+1)/2 (n(n+1)/2 + 1)/2 for n = 5.
HE_CCPVDZ_ERI8_CART_LEN = 120
# sum of Gp^2 (math.fsum).
HE_CCPVDZ_ERI8_CART_SUMSQ = 16.826897023586064
# sum of w_k Gp_k with the LCG weights from HE_SEED_WSUM (math.fsum).
HE_CCPVDZ_ERI8_CART_WSUM = 9.025911808873676
#| Every packed-8 two-electron integral, unit-Cartesian, chemist's notation.
def he_ccpvdz_eri8_cart() -> Arr {
    return [1.5430357804430757, 0.7849792460524572, 0.4267744902797824, 0.8034603752199505,
        0.4759077297853284, 0.6155616041576865, 0.0, 0.0,
        0.0, 0.1773074541945639, 0.0, 0.0,
        0.0, 0.1425235475592445, 0.13000788595743887, 1.099798622537393,
        0.627782301689996, 0.7344037258323985, 0.0, 0.0,
        1.0405308991232871, 0.0, 0.0, 0.0,
        0.0, 0.0, 0.0, 0.1773074541945639,
        0.0, 0.0, 0.0, 0.0,
        0.0, 0.0, 0.1425235475592445, 0.13000788595743887,
        0.0, 0.0, 0.0, 0.0,
        0.0, 0.0, 0.0, 0.0,
        0.06370597341571146, 1.099798622537393, 0.627782301689996, 0.7344037258323985,
        0.0, 0.0, 0.9131189522918641, 0.0,
        0.0, 0.0, 1.0405308991232871, 0.0,
        0.0, 0.0, 0.0, 0.0,
        0.0, 0.0, 0.0, 0.0,
        0.0, 0.1773074541945639, 0.0, 0.0,
        0.0, 0.0, 0.0, 0.0,
        0.0, 0.0, 0.0, 0.0,
        0.1425235475592445, 0.13000788595743887, 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.06370597341571146, 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.06370597341571146, 1.099798622537393, 0.627782301689996, 0.7344037258323985,
        0.0, 0.0, 0.9131189522918641, 0.0,
        0.0, 0.0, 0.9131189522918641, 0.0,
        0.0, 0.0, 0.0, 1.0405308991232871]
}

def he_ccpvdz_s_sph() -> Arr {
    return he_ccpvdz_s_unit()
}

def he_ccpvdz_t_sph() -> Arr {
    return he_ccpvdz_t_unit()
}

def he_ccpvdz_v_sph() -> Arr {
    return he_ccpvdz_v_unit()
}

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

def he_ccpvdz_m() -> Arr {
    return he_eye(5)
}

# ---- HE_CCPVDZ RHF: nalpha 1, nbeta 1; default state is stable ----

HE_CCPVDZ_E_RHF = -2.8551604772427406
HE_CCPVDZ_E_RHF_DEFAULT = -2.8551604772427406
HE_CCPVDZ_E_RHF_STABLE = -2.8551604772427406
# 1: no restart from the stability orbitals lowers the minao state; 0: one does; -1: not analyzed (stability restart).
HE_CCPVDZ_RHF_IS_STABLE = 1
HE_CCPVDZ_RHF_NALPHA = 1
HE_CCPVDZ_RHF_NBETA = 1

# Max |D' - D| after one PySCF step (get_fock, eig, get_occ, make_rdm1) from the default density below.
HE_CCPVDZ_RHF_SELF_STEP = 4.773959005888173e-15
#| Converged total density (alpha + beta), default state.
def he_ccpvdz_rhf_d_flat() -> Arr {
    return [0.7023880950276224, 0.608056008867343, 0.0, 0.0,
        0.0, 0.608056008867343, 0.5263929052002826, 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.0, 0.0, 0.0, 0.0,
        0.0]
}

def he_ccpvdz_rhf_d() -> Arr {
    return he_mat(he_ccpvdz_rhf_d_flat(), 5, 5)
}

# ---- HE_CCPVDZ UHF: nalpha 1, nbeta 1; default state is stable ----

HE_CCPVDZ_E_UHF = -2.85516047724274
HE_CCPVDZ_E_UHF_DEFAULT = -2.85516047724274
HE_CCPVDZ_E_UHF_STABLE = -2.85516047724274
# 1: no restart from the stability orbitals lowers the minao state; 0: one does; -1: not analyzed (stability restart).
HE_CCPVDZ_UHF_IS_STABLE = 1
HE_CCPVDZ_UHF_NALPHA = 1
HE_CCPVDZ_UHF_NBETA = 1
# <S^2> of the default UHF state.
HE_CCPVDZ_UHF_S2 = -2.220446049250313e-16
HE_CCPVDZ_UHF_S2_STABLE = -2.220446049250313e-16

# Max |D' - D| after one PySCF step (get_fock, eig, get_occ, make_rdm1) from the default density below.
HE_CCPVDZ_UHF_SELF_STEP = 2.3869795029440866e-15
#| Converged alpha density, default state.
def he_ccpvdz_uhf_da_flat() -> Arr {
    return [0.3511940475138112, 0.3040280044336715, 0.0, 0.0,
        0.0, 0.3040280044336715, 0.2631964526001413, 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.0, 0.0, 0.0, 0.0,
        0.0]
}

def he_ccpvdz_uhf_da() -> Arr {
    return he_mat(he_ccpvdz_uhf_da_flat(), 5, 5)
}

#| Converged beta density, default state.
def he_ccpvdz_uhf_db_flat() -> Arr {
    return [0.3511940475138112, 0.3040280044336715, 0.0, 0.0,
        0.0, 0.3040280044336715, 0.2631964526001413, 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.0, 0.0, 0.0, 0.0,
        0.0]
}

def he_ccpvdz_uhf_db() -> Arr {
    return he_mat(he_ccpvdz_uhf_db_flat(), 5, 5)
}

# Lowest eigenvalue of H_core C = S C e.
HE_CCPVDZ_EPS1E = -1.9936233377269412
test "he_sto3g packs" {
    assert_near(he_pack_bad(he_sto3g_shells(), HE_STO3G_NSH, HE_STO3G_NPRIMTOT, HE_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(HE_STO3G_NAO_CART - HE_STO3G_NAO_SPH, 0.0, 0.0)
    assert_near(len(he_sto3g_fc()), HE_STO3G_NAO_CART, 0.0)
}

test "he_sto3g unit one-electron integrals are symmetric with a unit overlap diagonal" {
    assert_near(he_asym(he_sto3g_s_unit(), 1), 0.0, 1.0e-14)
    assert_near(he_asym(he_sto3g_t_unit(), 1), 0.0, 1.0e-14)
    assert_near(he_asym(he_sto3g_v_unit(), 1), 0.0, 1.0e-14)
    assert_near(he_diag_dev(he_sto3g_s_unit(), 1), 0.0, 1.0e-14)
}

test "he_sto3g unit packed ERIs match their checksums and the Schwarz bound" {
    assert_near(len(he_sto3g_eri8_cart()), HE_STO3G_ERI8_CART_LEN, 0.0)
    assert_near(he_sumsq(he_sto3g_eri8_cart()) / HE_STO3G_ERI8_CART_SUMSQ, 1.0, 1.0e-12)
    assert_near(he_wsum(he_sto3g_eri8_cart(), HE_SEED_WSUM), HE_STO3G_ERI8_CART_WSUM, 1.0e-12 * HE_STO3G_ERI8_CART_SUMSQ)
    assert_lt(he_schwarz(he_sto3g_eri8_cart(), 1), 1.0e-12)
}

test "he_sto3g rhf energies, labels and densities" {
    assert_lt(HE_STO3G_E_RHF_STABLE, HE_STO3G_E_RHF_DEFAULT + 1.0e-12)
    assert_lt(HE_STO3G_RHF_SELF_STEP, 1.0e-10)
    assert_near(he_asym(he_sto3g_rhf_d(), 1), 0.0, 1.0e-12)
    assert_near(he_trace2(he_sto3g_rhf_d(), he_sto3g_s_unit(), 1), 2.0, 1.0e-10)
    assert_near(he_e_elec(he_sto3g_t_unit() + he_sto3g_v_unit(), he_sto3g_eri8_cart(), he_sto3g_rhf_d() * 0.5, he_sto3g_rhf_d() * 0.5, 1) + HE_STO3G_E_NUC, HE_STO3G_E_RHF_DEFAULT, 1.0e-10)
}

test "he_sto3g uhf energies, labels and densities" {
    assert_lt(HE_STO3G_E_UHF_STABLE, HE_STO3G_E_UHF_DEFAULT + 1.0e-12)
    assert_lt(0.0 - 1.0e-10, HE_STO3G_UHF_S2)
    assert_lt(0.0 - 1.0e-10, HE_STO3G_UHF_S2_STABLE)
    assert_lt(HE_STO3G_UHF_SELF_STEP, 1.0e-10)
    assert_near(he_asym(he_sto3g_uhf_da(), 1) + he_asym(he_sto3g_uhf_db(), 1), 0.0, 1.0e-12)
    assert_near(he_trace2(he_sto3g_uhf_da(), he_sto3g_s_unit(), 1) + he_trace2(he_sto3g_uhf_db(), he_sto3g_s_unit(), 1), 2.0, 1.0e-10)
    assert_near(he_e_elec(he_sto3g_t_unit() + he_sto3g_v_unit(), he_sto3g_eri8_cart(), he_sto3g_uhf_da(), he_sto3g_uhf_db(), 1) + HE_STO3G_E_NUC, HE_STO3G_E_UHF_DEFAULT, 1.0e-10)
}

test "he_sto3g dication_rhf energies, labels and densities" {
    assert_lt(HE_STO3G_E_DICATION_RHF_STABLE, HE_STO3G_E_DICATION_RHF_DEFAULT + 1.0e-12)
    assert_lt(HE_STO3G_DICATION_RHF_SELF_STEP, 1.0e-10)
    assert_near(he_asym(he_sto3g_dication_rhf_d(), 1), 0.0, 1.0e-12)
    assert_near(he_trace2(he_sto3g_dication_rhf_d(), he_sto3g_s_unit(), 1), 0.0, 1.0e-10)
    assert_near(he_e_elec(he_sto3g_t_unit() + he_sto3g_v_unit(), he_sto3g_eri8_cart(), he_sto3g_dication_rhf_d() * 0.5, he_sto3g_dication_rhf_d() * 0.5, 1) + HE_STO3G_E_NUC, HE_STO3G_E_DICATION_RHF_DEFAULT, 1.0e-10)
}

test "he_ccpvdz packs" {
    assert_near(he_pack_bad(he_ccpvdz_shells(), HE_CCPVDZ_NSH, HE_CCPVDZ_NPRIMTOT, HE_CCPVDZ_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(HE_CCPVDZ_NAO_CART - HE_CCPVDZ_NAO_SPH, 0.0, 0.0)
    assert_near(len(he_ccpvdz_fc()), HE_CCPVDZ_NAO_CART, 0.0)
}

test "he_ccpvdz unit one-electron integrals are symmetric with a unit overlap diagonal" {
    assert_near(he_asym(he_ccpvdz_s_unit(), 5), 0.0, 1.0e-14)
    assert_near(he_asym(he_ccpvdz_t_unit(), 5), 0.0, 1.0e-14)
    assert_near(he_asym(he_ccpvdz_v_unit(), 5), 0.0, 1.0e-14)
    assert_near(he_diag_dev(he_ccpvdz_s_unit(), 5), 0.0, 1.0e-14)
}

test "he_ccpvdz unit packed ERIs match their checksums and the Schwarz bound" {
    assert_near(len(he_ccpvdz_eri8_cart()), HE_CCPVDZ_ERI8_CART_LEN, 0.0)
    assert_near(he_sumsq(he_ccpvdz_eri8_cart()) / HE_CCPVDZ_ERI8_CART_SUMSQ, 1.0, 1.0e-12)
    assert_near(he_wsum(he_ccpvdz_eri8_cart(), HE_SEED_WSUM), HE_CCPVDZ_ERI8_CART_WSUM, 1.0e-12 * HE_CCPVDZ_ERI8_CART_SUMSQ)
    assert_lt(he_schwarz(he_ccpvdz_eri8_cart(), 5), 1.0e-12)
}

test "he_ccpvdz rhf energies, labels and densities" {
    assert_lt(HE_CCPVDZ_E_RHF_STABLE, HE_CCPVDZ_E_RHF_DEFAULT + 1.0e-12)
    assert_lt(HE_CCPVDZ_RHF_SELF_STEP, 1.0e-10)
    assert_near(he_asym(he_ccpvdz_rhf_d(), 5), 0.0, 1.0e-12)
    assert_near(he_trace2(he_ccpvdz_rhf_d(), he_ccpvdz_s_unit(), 5), 2.0, 1.0e-10)
    assert_near(he_e_elec(he_ccpvdz_t_unit() + he_ccpvdz_v_unit(), he_ccpvdz_eri8_cart(), he_ccpvdz_rhf_d() * 0.5, he_ccpvdz_rhf_d() * 0.5, 5) + HE_CCPVDZ_E_NUC, HE_CCPVDZ_E_RHF_DEFAULT, 1.0e-10)
}

test "he_ccpvdz uhf energies, labels and densities" {
    assert_lt(HE_CCPVDZ_E_UHF_STABLE, HE_CCPVDZ_E_UHF_DEFAULT + 1.0e-12)
    assert_lt(0.0 - 1.0e-10, HE_CCPVDZ_UHF_S2)
    assert_lt(0.0 - 1.0e-10, HE_CCPVDZ_UHF_S2_STABLE)
    assert_lt(HE_CCPVDZ_UHF_SELF_STEP, 1.0e-10)
    assert_near(he_asym(he_ccpvdz_uhf_da(), 5) + he_asym(he_ccpvdz_uhf_db(), 5), 0.0, 1.0e-12)
    assert_near(he_trace2(he_ccpvdz_uhf_da(), he_ccpvdz_s_unit(), 5) + he_trace2(he_ccpvdz_uhf_db(), he_ccpvdz_s_unit(), 5), 2.0, 1.0e-10)
    assert_near(he_e_elec(he_ccpvdz_t_unit() + he_ccpvdz_v_unit(), he_ccpvdz_eri8_cart(), he_ccpvdz_uhf_da(), he_ccpvdz_uhf_db(), 5) + HE_CCPVDZ_E_NUC, HE_CCPVDZ_E_UHF_DEFAULT, 1.0e-10)
}