###
post: Cartesian to spherical, one shell at a time
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).

M_l maps the unit-Cartesian components of one shell (PySCF order: lx from l
down, then ly) to PySCF's spherical ones (p as x, y, z; d as xy, yz, z^2,
xz, x^2 - y^2; f and g as m = -l..l): M_l = diag(F) cart2sph(l, 'sp'), with
F^2 = 4 pi (2i-1)!! (2j-1)!! (2k-1)!! / (2l+1)!! for l >= 2 and F = 1 below.

Also the two constants pyscf/gto/test/test_mole.py:520 and :529 assert: the
s and p factors of the unnormalized cart2sph, 1/(2 sqrt(pi)) and
sqrt(3/(4 pi)).
###

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

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

# gto.cart2sph(0): test_mole.py:520 asserts 0.282094791773878143.
C2S_S_FACTOR = 0.28209479177387814
# gto.cart2sph(1) diagonal: test_mole.py:529 asserts 0.488602511902919921.
C2S_P_FACTOR = 0.4886025119029199
C2S_LMAX = 4
C2S_PI = 3.141592653589793

#| M_0, [1,1].
def c2s_m0_flat() -> Arr {
    return [1.0]
}

def c2s_m0() -> Arr {
    return c2s_mat(c2s_m0_flat(), 1, 1)
}

#| Powers (i, j, k) of each component.
def c2s_powers0_flat() -> Arr {
    return [0, 0, 0]
}

def c2s_powers0() -> Arr {
    return c2s_mat(c2s_powers0_flat(), 1, 3)
}

#| F^2 per component.
def c2s_fc2_0() -> Arr {
    return [1.0]
}

#| M_1, [3,3].
def c2s_m1_flat() -> Arr {
    return [1.0, 0.0, 0.0, 0.0,
        1.0, 0.0, 0.0, 0.0,
        1.0]
}

def c2s_m1() -> Arr {
    return c2s_mat(c2s_m1_flat(), 3, 3)
}

#| Powers (i, j, k) of each component.
def c2s_powers1_flat() -> Arr {
    return [1, 0, 0, 0,
        1, 0, 0, 0,
        1]
}

def c2s_powers1() -> Arr {
    return c2s_mat(c2s_powers1_flat(), 3, 3)
}

#| F^2 per component.
def c2s_fc2_1() -> Arr {
    return [1.0, 1.0, 1.0]
}

#| M_2, [6,5].
def c2s_m2_flat() -> Arr {
    return [0.0, 0.0, -0.49999999999999994, 0.0,
        0.8660254037844387, 1.0, 0.0, 0.0,
        0.0, 0.0, 0.0, 0.0,
        0.0, 1.0, 0.0, 0.0,
        0.0, -0.49999999999999994, 0.0, -0.8660254037844387,
        0.0, 1.0, 0.0, 0.0,
        0.0, 0.0, 0.0, 0.9999999999999999,
        0.0, 0.0]
}

def c2s_m2() -> Arr {
    return c2s_mat(c2s_m2_flat(), 6, 5)
}

#| Powers (i, j, k) of each component.
def c2s_powers2_flat() -> Arr {
    return [2, 0, 0, 1,
        1, 0, 1, 0,
        1, 0, 2, 0,
        0, 1, 1, 0,
        0, 2]
}

def c2s_powers2() -> Arr {
    return c2s_mat(c2s_powers2_flat(), 6, 3)
}

#| F^2 per component.
def c2s_fc2_2() -> Arr {
    return [2.5132741228718345, 0.8377580409572781, 0.8377580409572781, 2.5132741228718345,
        0.8377580409572781, 2.5132741228718345]
}

#| M_3, [10,7].
def c2s_m3_flat() -> Arr {
    return [0.0, 0.0, 0.0, 0.0,
        -0.6123724356957945, 0.0, 0.7905694150420948, 1.0606601717798212,
        0.0, -0.27386127875258304, 0.0, 0.0,
        0.0, 0.0, 0.0, 0.0,
        0.0, -0.6708203932499369, 0.0, 0.8660254037844387,
        0.0, 0.0, 0.0, 0.0,
        0.0, -0.27386127875258304, 0.0, -1.0606601717798212,
        0.0, 1.0000000000000002, 0.0, 0.0,
        0.0, 0.0, 0.0, 0.0,
        0.0, 0.0, 0.0, 1.0954451150103321,
        0.0, 0.0, -0.7905694150420948, 0.0,
        -0.6123724356957945, 0.0, 0.0, 0.0,
        0.0, 0.0, 0.0, 0.0,
        -0.6708203932499369, 0.0, -0.8660254037844387, 0.0,
        0.0, 0.0, 1.0954451150103321, 0.0,
        0.0, 0.0, 0.0, 0.0,
        0.0, 0.0, 0.9999999999999999, 0.0,
        0.0, 0.0]
}

def c2s_m3() -> Arr {
    return c2s_mat(c2s_m3_flat(), 10, 7)
}

#| Powers (i, j, k) of each component.
def c2s_powers3_flat() -> Arr {
    return [3, 0, 0, 2,
        1, 0, 2, 0,
        1, 1, 2, 0,
        1, 1, 1, 1,
        0, 2, 0, 3,
        0, 0, 2, 1,
        0, 1, 2, 0,
        0, 3]
}

def c2s_powers3() -> Arr {
    return c2s_mat(c2s_powers3_flat(), 10, 3)
}

#| F^2 per component.
def c2s_fc2_3() -> Arr {
    return [1.7951958020513101, 0.35903916041026207, 0.35903916041026207, 0.35903916041026207,
        0.11967972013675403, 0.35903916041026207, 1.7951958020513101, 0.35903916041026207,
        0.35903916041026207, 1.7951958020513101]
}

#| M_4, [15,9].
def c2s_m4_flat() -> Arr {
    return [0.0, 0.0, 0.0, 0.0,
        0.37500000000000006, 0.0, -0.5590169943749475, 0.0,
        0.739509972887452, 1.1180339887498947, 0.0, -0.42257712736425823,
        0.0, 0.0, 0.0, 0.0,
        0.0, 0.0, 0.0, 0.0,
        0.0, 0.0, 0.0, -0.8964214570007952,
        0.0, 0.7905694150420948, 0.0, 0.0,
        0.0, 0.0, 0.0, 0.21957751641342,
        0.0, 0.0, 0.0, -1.299038105676658,
        0.0, 1.0606601717798212, 0.0, -0.4008918628686366,
        0.0, 0.0, 0.0, 0.0,
        0.0, 0.0, 0.0, 0.0,
        0.0, -0.87831006565368, 0.0, 0.9819805060619659,
        0.0, 0.0, -1.1180339887498947, 0.0,
        -0.42257712736425823, 0.0, 0.0, 0.0,
        0.0, 0.0, 0.0, 0.0,
        0.0, 0.0, 0.0, 0.0,
        -0.4008918628686366, 0.0, -1.0606601717798212, 0.0,
        0.0, 0.0, 1.1338934190276817, 0.0,
        0.0, 0.0, 0.0, 0.0,
        0.0, 0.0, 0.0, 0.0,
        0.0, 0.0, 1.1952286093343936, 0.0,
        0.0, 0.0, 0.0, 0.0,
        0.0, 0.0, 0.37500000000000006, 0.0,
        0.5590169943749475, 0.0, 0.739509972887452, 0.0,
        -0.7905694150420948, 0.0, -0.8964214570007952, 0.0,
        0.0, 0.0, 0.0, 0.0,
        0.0, 0.0, 0.0, 0.0,
        -0.87831006565368, 0.0, -0.9819805060619659, 0.0,
        0.0, 0.0, 0.0, 0.0,
        1.1952286093343936, 0.0, 0.0, 0.0,
        0.0, 0.0, 0.0, 0.0,
        0.0, 0.0, 1.0, 0.0,
        0.0, 0.0, 0.0]
}

def c2s_m4() -> Arr {
    return c2s_mat(c2s_m4_flat(), 15, 9)
}

#| Powers (i, j, k) of each component.
def c2s_powers4_flat() -> Arr {
    return [4, 0, 0, 3,
        1, 0, 3, 0,
        1, 2, 2, 0,
        2, 1, 1, 2,
        0, 2, 1, 3,
        0, 1, 2, 1,
        1, 1, 2, 1,
        0, 3, 0, 4,
        0, 0, 3, 1,
        0, 2, 2, 0,
        1, 3, 0, 0,
        4]
}

def c2s_powers4() -> Arr {
    return c2s_mat(c2s_powers4_flat(), 15, 3)
}

#| F^2 per component.
def c2s_fc2_4() -> Arr {
    return [1.3962634015954636, 0.19946620022792336, 0.19946620022792336, 0.11967972013675403,
        0.03989324004558467, 0.11967972013675403, 0.19946620022792336, 0.03989324004558467,
        0.03989324004558467, 0.19946620022792336, 1.3962634015954636, 0.19946620022792336,
        0.11967972013675403, 0.19946620022792336, 1.3962634015954636]
}

#| (n)!! with (-1)!! = 1.
def c2s_dfact(n) {
    r = 1.0
    k = n
    while k > 1 {
        r = r * k
        k = k - 2
    }
    return r
}

#| Overlap of two unit-normalized components on one center with one exponent:
#| per axis (a + b - 1)!! / sqrt((2a - 1)!! (2b - 1)!!) when a + b is even, else 0.
def c2s_shell_overlap(pw: Arr, nc) -> Arr {
    s = zeros(nc, nc)
    for i to nc {
        for j to nc {
            t = 1.0
            for d to shape(pw)[1] {
                a = pw[i, d]
                b = pw[j, d]
                if mod(a + b, 2.0) == 1.0 {
                    t = 0.0
                } else {
                    t = t * c2s_dfact(a + b - 1) / sqrt(c2s_dfact(2 * a - 1) * c2s_dfact(2 * b - 1))
                }
            }
            s[i, j] = t
        }
    }
    return s
}

#| Largest relative gap between the tabulated F^2 and the closed form.
def c2s_fc2_dev(pw: Arr, f2: Arr, nc) {
    w = 0.0
    for i to nc {
        a = pw[i, 0]
        b = pw[i, 1]
        c = pw[i, 2]
        l = a + b + c
        x = 1.0
        if l >= 2 {
            x = 4.0 * C2S_PI * c2s_dfact(2 * a - 1) * c2s_dfact(2 * b - 1) * c2s_dfact(2 * c - 1) / c2s_dfact(2 * l + 1)
        }
        w = max(w, abs(f2[i] - x) / x)
    }
    return w
}

test "c2s spherical functions are orthonormal and F matches its closed form" {
    assert_near(C2S_S_FACTOR, 1.0 / (2.0 * sqrt(C2S_PI)), 1.0e-16)
    assert_near(C2S_P_FACTOR, sqrt(3.0 / (4.0 * C2S_PI)), 1.0e-16)
    assert_lt(c2s_mtxm_dev(c2s_m0(), c2s_shell_overlap(c2s_powers0(), 1), c2s_eye(1), 1, 1), 1.0e-14)
    assert_lt(c2s_fc2_dev(c2s_powers0(), c2s_fc2_0(), 1), 1.0e-15)
    assert_lt(c2s_mtxm_dev(c2s_m1(), c2s_shell_overlap(c2s_powers1(), 3), c2s_eye(3), 3, 3), 1.0e-14)
    assert_lt(c2s_fc2_dev(c2s_powers1(), c2s_fc2_1(), 3), 1.0e-15)
    assert_lt(c2s_mtxm_dev(c2s_m2(), c2s_shell_overlap(c2s_powers2(), 6), c2s_eye(5), 6, 5), 1.0e-14)
    assert_lt(c2s_fc2_dev(c2s_powers2(), c2s_fc2_2(), 6), 1.0e-15)
    assert_lt(c2s_mtxm_dev(c2s_m3(), c2s_shell_overlap(c2s_powers3(), 10), c2s_eye(7), 10, 7), 1.0e-14)
    assert_lt(c2s_fc2_dev(c2s_powers3(), c2s_fc2_3(), 10), 1.0e-15)
    assert_lt(c2s_mtxm_dev(c2s_m4(), c2s_shell_overlap(c2s_powers4(), 15), c2s_eye(9), 15, 9), 1.0e-14)
    assert_lt(c2s_fc2_dev(c2s_powers4(), c2s_fc2_4(), 15), 1.0e-15)
}