###
post: Water in STO-3G, as PySCF has it
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.
WATERSTO3G_SEED_WSUM = 20260923

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

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

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

WATERSTO3G_NATOM = 3
WATERSTO3G_NSH = 5
WATERSTO3G_NPRIMTOT = 15
WATERSTO3G_NAO_CART = 7
WATERSTO3G_NAO_SPH = 7
WATERSTO3G_E_NUC = 9.188258417746113

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

def watersto3g_atoms() -> Arr {
    return watersto3g_mat(watersto3g_atoms_flat(), 3, 4)
}

#| (atom, l, nprim, p0, cart0, sph0, ncart, nsph); general contractions split.
def watersto3g_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,
        2, 0, 3, 12,
        6, 6, 1, 1]
}

def watersto3g_shells() -> Arr {
    return watersto3g_mat(watersto3g_shells_flat(), 5, 8)
}

#| (exponent, raw bas_ctr_coeff).
def watersto3g_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,
        3.42525091, 0.15432897070298393, 0.62391373, 0.5353281424384733,
        0.1688554, 0.44463454202535485]
}

def watersto3g_prims() -> Arr {
    return watersto3g_mat(watersto3g_prims_flat(), 15, 2)
}

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

#| Overlap, unit-Cartesian, [7,7].
def watersto3g_s_unit_flat() -> Arr {
    return [1.0, 0.2367039365108476, 0.0, 0.0,
        0.0, 0.05388499213036836, 0.05388499213036836, 0.23670393651084762,
        1.0, 0.0, 0.0, 0.0,
        0.47433176700934565, 0.47433176700934565, 0.0, 0.0,
        1.0, 0.0, 0.0, 0.0,
        0.0, 0.0, 0.0, 0.0,
        1.0, 0.0, -0.31072668893058497, 0.31072668893058497,
        0.0, 0.0, 0.0, 0.0,
        1.0, 0.24094658705713787, 0.24094658705713787, 0.05388499213036836,
        0.47433176700934565, 0.0, -0.310726688930585, 0.2409465870571379,
        1.0000000000000002, 0.2516808424829575, 0.05388499213036836, 0.47433176700934565,
        0.0, 0.310726688930585, 0.2409465870571379, 0.2516808424829575,
        1.0000000000000002]
}

def watersto3g_s_unit() -> Arr {
    return watersto3g_mat(watersto3g_s_unit_flat(), 7, 7)
}

#| Kinetic energy, unit-Cartesian.
def watersto3g_t_unit_flat() -> Arr {
    return [29.003199945539567, -0.16801093931649205, 0.0, 0.0,
        0.0, -0.0025622436863846736, -0.0025622436863846736, -0.16801093931649186,
        0.8081279549303472, 0.0, 0.0, 0.0,
        0.12826283563609547, 0.12826283563609547, 0.0, 0.0,
        2.5287311981947633, 0.0, 0.0, 0.0,
        0.0, 0.0, 0.0, 0.0,
        2.5287311981947633, 0.0, -0.22420987735782355, 0.22420987735782355,
        0.0, 0.0, 0.0, 0.0,
        2.5287311981947633, 0.17385891414668747, 0.17385891414668747, -0.0025622436863848106,
        0.1282628356360955, 0.0, -0.22420987735782355, 0.17385891414668747,
        0.7600318835666091, 0.008476840173892185, -0.0025622436863848106, 0.1282628356360955,
        0.0, 0.22420987735782355, 0.17385891414668747, 0.008476840173892185,
        0.7600318835666091]
}

def watersto3g_t_unit() -> Arr {
    return watersto3g_mat(watersto3g_t_unit_flat(), 7, 7)
}

#| Nuclear attraction, unit-Cartesian.
def watersto3g_v_unit_flat() -> Arr {
    return [-61.72329688309335, -7.444599890028586, 0.0, 0.0,
        -0.018996438617907625, -1.7439691717269703, -1.7439691717269703, -7.444599890028587,
        -10.142196980858547, 0.0, 0.0, -0.2235045224420784,
        -3.8650803747486626, -3.8650803747486626, 0.0, 0.0,
        -9.98577456465457, 0.0, 0.0, 0.0,
        0.0, 0.0, 0.0, 0.0,
        -10.141625718175106, 0.0, 2.2512619716810565, -2.2512619716810565,
        -0.018996438617907625, -0.2235045224420784, 0.0, 0.0,
        -10.079486404510178, -1.8186761733468688, -1.8186761733468688, -1.7439691717269699,
        -3.8650803747486635, 0.0, 2.2512619716810565, -1.8186761733468693,
        -5.833953772754734, -1.6156328893155159, -1.7439691717269699, -3.8650803747486635,
        0.0, -2.2512619716810565, -1.8186761733468693, -1.6156328893155159,
        -5.833953772754734]
}

def watersto3g_v_unit() -> Arr {
    return watersto3g_mat(watersto3g_v_unit_flat(), 7, 7)
}

# Packed-8 length n(n+1)/2 (n(n+1)/2 + 1)/2 for n = 7.
WATERSTO3G_ERI8_CART_LEN = 406
# sum of Gp^2 (math.fsum).
WATERSTO3G_ERI8_CART_SUMSQ = 42.53442033100022
# sum of w_k Gp_k with the LCG weights from WATERSTO3G_SEED_WSUM (math.fsum).
WATERSTO3G_ERI8_CART_WSUM = 16.599013267525155
#| Every packed-8 two-electron integral, unit-Cartesian, chemist's notation.
def watersto3g_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.17149760215075188, 0.03138026252454724, 0.058411020607098886, 0.0,
        0.0, 0.05830033803863549, -0.0011404458740247871, -0.0021287514894520864,
        -0.0, 0.05851231603108477, 0.0008843351757629452, 0.0016506963333003625,
        0.0, -0.00016437395187282464, 0.058427798421395526, 0.007287667504468778,
        0.40314276722614695, 0.09355840257683681, 0.3282054625686653, 0.0,
        0.0, 0.32275150936706104, -0.005528234084780467, -0.04163336715672436,
        0.0, 0.3327779203648954, 0.004286754831923558, 0.03228373384543884,
        0.0, -0.007774773125137019, 0.3287802964535281, 0.021841218083353386,
        0.1610574050156167, 0.0, 0.0, 0.0,
        0.00995838348187528, 0.057352550822734084, 0.0, -0.0,
        0.0, -0.0087251345308759, 0.0, 0.0,
        0.0, 0.006765725191049079, -0.0, 0.0,
        0.0, 0.0, 0.023384406798298535, -0.2246162474029574,
        -0.05281917251202047, -0.19790493314831978, -0.0, 0.0,
        -0.18885696593222556, 0.01270466589522253, 0.08124319475632823, 0.0,
        -0.2122250269491666, -0.0021295479215783786, -0.018525505930012905, -0.0,
        0.011354554619973863, -0.19241527884474716, -0.012894396033370066, -0.1179623905290038,
        -0.0, 0.10882134008503214, 0.17417402539700913, 0.04095753535608456,
        0.1534612889802691, 0.0, 0.0, 0.14644522985761746,
        -0.002129547921578379, -0.018525505930012905, -0.0, 0.15103405928654226,
        0.01160969739200277, 0.07171777140518795, 0.0, -0.012283447443397495,
        0.16273590042419403, 0.009998692829046535, 0.09147149701522486, 0.0,
        -0.06625030361864281, 0.0747568350996768, 0.5316455102911697, 0.12538713207499905,
        0.503017128197289, 0.0, 0.0, 0.48252986910623535,
        -0.009893347990075012, -0.10153005969362508, -0.0, 0.5205293644644061,
        0.007671592166676399, 0.07872938578620595, 0.0, -0.02946592308486956,
        0.5053786099923891, 0.02955808822799496, 0.3059654949183104, 0.0,
        -0.25601150702596864, 0.19851883041511711, 0.7746059439198978, 0.17149760215075188,
        0.03138026252454724, 0.058411020607098886, 0.0, 0.0,
        0.05830033803863549, 0.0011404458740247871, 0.0021287514894520864, 0.0,
        0.05851231603108477, 0.0008843351757629452, 0.0016506963333003625, 0.0,
        0.00016437395187282464, 0.058427798421395526, 0.007174976730941921, 0.021157976896380663,
        0.0, -0.011379448659054677, 0.009691905057181172, 0.028284618344378603,
        0.007287667504468778, 0.40314276722614695, 0.09355840257683681, 0.3282054625686653,
        0.0, 0.0, 0.32275150936706104, 0.005528234084780467,
        0.04163336715672436, 0.0, 0.3327779203648954, 0.004286754831923558,
        0.03228373384543884, 0.0, 0.007774773125137019, 0.3287802964535281,
        0.021157976896380666, 0.13176845795280362, 0.0, -0.06625682657974014,
        0.07279964032663618, 0.20712426341320792, 0.021841218083353386, 0.1610574050156167,
        0.0, 0.0, 0.0, 0.00995838348187528,
        0.057352550822734084, 0.0, 0.0, 0.0,
        0.0087251345308759, 0.0, 0.0, 0.0,
        0.006765725191049079, 0.0, 0.0, 0.0,
        0.0, 0.018535644365079516, 0.0, 0.0,
        0.0, 0.0, 0.0, 0.023384406798298535,
        0.2246162474029574, 0.05281917251202047, 0.19790493314831978, 0.0,
        0.0, 0.18885696593222556, 0.01270466589522253, 0.08124319475632823,
        0.0, 0.2122250269491666, 0.0021295479215783786, 0.018525505930012905,
        0.0, 0.011354554619973863, 0.19241527884474716, 0.01137944865905468,
        0.06625682657974015, 0.0, -0.015211875272395026, 0.035749386434473854,
        0.09134914223130425, 0.012894396033370066, 0.1179623905290038, 0.0,
        0.10882134008503214, 0.17417402539700913, 0.04095753535608456, 0.1534612889802691,
        0.0, 0.0, 0.14644522985761746, 0.002129547921578379,
        0.018525505930012905, 0.0, 0.15103405928654226, 0.01160969739200277,
        0.07171777140518795, 0.0, 0.012283447443397495, 0.16273590042419403,
        0.009691905057181174, 0.07279964032663618, 0.0, -0.035749386434473854,
        0.05814038739590601, 0.12561198241700153, 0.009998692829046535, 0.09147149701522486,
        0.0, 0.06625030361864281, 0.0747568350996768, 0.16062464969163207,
        0.037603131804085574, 0.14256179987993847, 0.0, 0.0,
        0.1392986937878778, -6.098637220230962e-20, -4.336808689942018e-19, 0.0,
        0.14209840853997588, 0.0027387266582381965, 0.024242314322331773, 0.0,
        -2.168404344971009e-19, 0.14613657007671566, 0.008689101144829238, 0.06908077040502726,
        0.0, -0.04356088219013673, 0.04410270131771886, 0.13441103605110244,
        0.008689101144829238, 0.06908077040502725, 0.0, 0.04356088219013673,
        0.04410270131771886, 0.035842000821326006, 0.5316455102911697, 0.12538713207499905,
        0.503017128197289, 0.0, 0.0, 0.48252986910623535,
        0.009893347990075012, 0.10153005969362508, 0.0, 0.5205293644644061,
        0.007671592166676399, 0.07872938578620595, 0.0, 0.02946592308486956,
        0.5053786099923891, 0.028284618344378603, 0.20712426341320792, 0.0,
        -0.09134914223130423, 0.1256119824170015, 0.3429248398814874, 0.02955808822799496,
        0.3059654949183104, 0.0, 0.25601150702596864, 0.19851883041511711,
        0.13441103605110247, 0.7746059439198978]
}

def watersto3g_s_sph() -> Arr {
    return watersto3g_s_unit()
}

def watersto3g_t_sph() -> Arr {
    return watersto3g_t_unit()
}

def watersto3g_v_sph() -> Arr {
    return watersto3g_v_unit()
}

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

def watersto3g_m() -> Arr {
    return watersto3g_eye(7)
}

# ---- WATERSTO3G RHF: nalpha 5, nbeta 5; default state is stable ----

WATERSTO3G_E_RHF = -74.9630631297292
WATERSTO3G_E_RHF_DEFAULT = -74.9630631297292
WATERSTO3G_E_RHF_STABLE = -74.9630631297292
# 1: no restart from the stability orbitals lowers the minao state; 0: one does; -1: not analyzed (stability restart).
WATERSTO3G_RHF_IS_STABLE = 1
WATERSTO3G_RHF_NALPHA = 5
WATERSTO3G_RHF_NBETA = 5

# Max |D' - D| after one PySCF step (get_fock, eig, get_occ, make_rdm1) from the default density below.
WATERSTO3G_RHF_SELF_STEP = 1.9912960169676808e-12
#| Converged total density (alpha + beta), default state.
def watersto3g_rhf_d_flat() -> Arr {
    return [2.1062938996739544, -0.44624534025120155, -2.278239885142663e-17, -2.783706566767363e-18,
        0.10858551289154836, -0.028298324308520625, -0.028298324308520802, -0.44624534025120155,
        1.9682426718396986, -4.925649979552284e-17, 1.962921605668791e-16, -0.6175183083251482,
        -0.034551732473745975, -0.034551732473744996, -2.278239885142663e-17, -4.925649979552284e-17,
        2.000000000000001, 3.684679138317929e-17, -1.1530729233643215e-16, 3.299259305880717e-16,
        2.714482613392324e-16, -2.783706566767363e-18, 1.962921605668791e-16, 3.684679138317929e-17,
        0.7358325126501221, 4.331459062822807e-16, -0.5399712976911439, 0.539971297691143,
        0.10858551289154836, -0.6175183083251482, -1.1530729233643215e-16, 4.331459062822807e-16,
        1.2387277266608012, 0.47298803951197305, 0.47298803951197366, -0.028298324308520625,
        -0.034551732473745975, 3.299259305880717e-16, -0.5399712976911439, 0.47298803951197305,
        0.6014828074476598, -0.1910046059724499, -0.028298324308520802, -0.034551732473744996,
        2.714482613392324e-16, 0.539971297691143, 0.47298803951197366, -0.1910046059724499,
        0.6014828074476587]
}

def watersto3g_rhf_d() -> Arr {
    return watersto3g_mat(watersto3g_rhf_d_flat(), 7, 7)
}

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

WATERSTO3G_E_UHF = -74.96306312972918
WATERSTO3G_E_UHF_DEFAULT = -74.96306312972918
WATERSTO3G_E_UHF_STABLE = -74.96306312972918
# 1: no restart from the stability orbitals lowers the minao state; 0: one does; -1: not analyzed (stability restart).
WATERSTO3G_UHF_IS_STABLE = 1
WATERSTO3G_UHF_NALPHA = 5
WATERSTO3G_UHF_NBETA = 5
# <S^2> of the default UHF state.
WATERSTO3G_UHF_S2 = 8.881784197001252e-16
WATERSTO3G_UHF_S2_STABLE = 8.881784197001252e-16

# Max |D' - D| after one PySCF step (get_fock, eig, get_occ, make_rdm1) from the default density below.
WATERSTO3G_UHF_SELF_STEP = 1.4954704141700859e-13
#| Converged alpha density, default state.
def watersto3g_uhf_da_flat() -> Arr {
    return [1.0531469498370278, -0.22312267012588524, 1.1523061623645346e-18, 3.881812870396987e-17,
        0.05429275644574441, -0.014149162154072336, -0.01414916215407241, -0.22312267012588524,
        0.984121335921395, 5.545266185202522e-18, -1.7820437208028039e-16, -0.30875915416267397,
        -0.01727586623780832, -0.01727586623780794, 1.1523061623645346e-18, 5.545266185202522e-18,
        1.0, 1.5214359607583777e-16, 3.824413402972076e-17, 4.264552603502837e-17,
        -8.887398033093843e-17, 3.881812870396987e-17, -1.7820437208028039e-16, 1.5214359607583777e-16,
        0.3679162563262958, 1.876944693267598e-16, -0.26998564884550225, 0.26998564884550236,
        0.05429275644574441, -0.30875915416267397, 3.824413402972076e-17, 1.876944693267598e-16,
        0.6193638633291378, 0.23649401975646317, 0.23649401975646284, -0.014149162154072336,
        -0.01727586623780832, 4.264552603502837e-17, -0.26998564884550225, 0.23649401975646317,
        0.30074140372349223, -0.09550230298502883, -0.01414916215407241, -0.01727586623780794,
        -8.887398033093843e-17, 0.26998564884550236, 0.23649401975646284, -0.09550230298502883,
        0.30074140372349195]
}

def watersto3g_uhf_da() -> Arr {
    return watersto3g_mat(watersto3g_uhf_da_flat(), 7, 7)
}

#| Converged beta density, default state.
def watersto3g_uhf_db_flat() -> Arr {
    return [1.053146949837031, -0.22312267012592005, -2.3354890396317193e-17, 2.1753305146439885e-17,
        0.05429275644582206, -0.014149162154025806, -0.014149162154025868, -0.22312267012592005,
        0.9841213359217038, 1.5362721760992175e-16, -5.883770572360201e-17, -0.30875915416312827,
        -0.01727586623818764, -0.017275866238187363, -2.3354890396317193e-17, 1.5362721760992175e-16,
        0.9999999999999998, 8.946033325334457e-17, 1.1445416562734433e-16, -9.193405592150583e-17,
        -1.767066753043936e-16, 2.1753305146439885e-17, -5.883770572360201e-17, 8.946033325334457e-17,
        0.3679162563255405, 1.7172875604643448e-16, -0.26998564884554493, 0.26998564884554477,
        0.05429275644582206, -0.30875915416312827, 1.1445416562734433e-16, 1.7172875604643448e-16,
        0.6193638633281551, 0.23649401975678822, 0.23649401975678822, -0.014149162154025806,
        -0.01727586623818764, -9.193405592150583e-17, -0.26998564884554493, 0.23649401975678822,
        0.30074140372439373, -0.09550230298506573, -0.014149162154025868, -0.017275866238187363,
        -1.767066753043936e-16, 0.26998564884554477, 0.23649401975678822, -0.09550230298506573,
        0.3007414037243934]
}

def watersto3g_uhf_db() -> Arr {
    return watersto3g_mat(watersto3g_uhf_db_flat(), 7, 7)
}

test "watersto3g packs" {
    assert_near(watersto3g_pack_bad(watersto3g_shells(), WATERSTO3G_NSH, WATERSTO3G_NPRIMTOT, WATERSTO3G_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(WATERSTO3G_NAO_CART - WATERSTO3G_NAO_SPH, 0.0, 0.0)
    assert_near(len(watersto3g_fc()), WATERSTO3G_NAO_CART, 0.0)
}

test "watersto3g unit one-electron integrals are symmetric with a unit overlap diagonal" {
    assert_near(watersto3g_asym(watersto3g_s_unit(), 7), 0.0, 1.0e-14)
    assert_near(watersto3g_asym(watersto3g_t_unit(), 7), 0.0, 1.0e-14)
    assert_near(watersto3g_asym(watersto3g_v_unit(), 7), 0.0, 1.0e-14)
    assert_near(watersto3g_diag_dev(watersto3g_s_unit(), 7), 0.0, 1.0e-14)
}

test "watersto3g unit packed ERIs match their checksums and the Schwarz bound" {
    assert_near(len(watersto3g_eri8_cart()), WATERSTO3G_ERI8_CART_LEN, 0.0)
    assert_near(watersto3g_sumsq(watersto3g_eri8_cart()) / WATERSTO3G_ERI8_CART_SUMSQ, 1.0, 1.0e-12)
    assert_near(watersto3g_wsum(watersto3g_eri8_cart(), WATERSTO3G_SEED_WSUM), WATERSTO3G_ERI8_CART_WSUM, 1.0e-12 * WATERSTO3G_ERI8_CART_SUMSQ)
    assert_lt(watersto3g_schwarz(watersto3g_eri8_cart(), 7), 1.0e-12)
}

test "watersto3g rhf energies, labels and densities" {
    assert_lt(WATERSTO3G_E_RHF_STABLE, WATERSTO3G_E_RHF_DEFAULT + 1.0e-12)
    assert_lt(WATERSTO3G_RHF_SELF_STEP, 1.0e-10)
    assert_near(watersto3g_asym(watersto3g_rhf_d(), 7), 0.0, 1.0e-12)
    assert_near(watersto3g_trace2(watersto3g_rhf_d(), watersto3g_s_unit(), 7), 10.0, 1.0e-10)
    assert_near(watersto3g_e_elec(watersto3g_t_unit() + watersto3g_v_unit(), watersto3g_eri8_cart(), watersto3g_rhf_d() * 0.5, watersto3g_rhf_d() * 0.5, 7) + WATERSTO3G_E_NUC, WATERSTO3G_E_RHF_DEFAULT, 1.0e-10)
}

test "watersto3g uhf energies, labels and densities" {
    assert_lt(WATERSTO3G_E_UHF_STABLE, WATERSTO3G_E_UHF_DEFAULT + 1.0e-12)
    assert_lt(0.0 - 1.0e-10, WATERSTO3G_UHF_S2)
    assert_lt(0.0 - 1.0e-10, WATERSTO3G_UHF_S2_STABLE)
    assert_lt(WATERSTO3G_UHF_SELF_STEP, 1.0e-10)
    assert_near(watersto3g_asym(watersto3g_uhf_da(), 7) + watersto3g_asym(watersto3g_uhf_db(), 7), 0.0, 1.0e-12)
    assert_near(watersto3g_trace2(watersto3g_uhf_da(), watersto3g_s_unit(), 7) + watersto3g_trace2(watersto3g_uhf_db(), watersto3g_s_unit(), 7), 10.0, 1.0e-10)
    assert_near(watersto3g_e_elec(watersto3g_t_unit() + watersto3g_v_unit(), watersto3g_eri8_cart(), watersto3g_uhf_da(), watersto3g_uhf_db(), 7) + WATERSTO3G_E_NUC, WATERSTO3G_E_UHF_DEFAULT, 1.0e-10)
}