Deriving Box 13.2 from "Molecular Electronic Structure Theory"

Here we show how to reproduce box 13.2 from the book Molecular Electronic Structure Theory

Setup:

julia> using SpinAdaptedSecondQuantization
julia> h = ∑((
               real_tensor("F", 1, 2) +
               ∑((-2 * psym_tensor("g", 1, 2, 3, 3) +
                   psym_tensor("g", 1, 3, 3, 2)) * occupied(3), [3])
           ) * E(1, 2) * electron(1, 2), 1:2);
julia> g = 1//2 * ∑(psym_tensor("g", 1:4...) * e(1:4...) * electron(1:4...), 1:4);
julia> H = h + g∑_pq(F_pq E_pq)
- 2 ∑_pqi(g_pqii E_pq)
+ ∑_pqi(g_piiq E_pq)
+ 1/2 ∑_pqrs(g_pqrs E_pq E_rs)
- 1/2 ∑_pqrs(δ_qr g_pqrs E_ps)
julia> Eai(a, i) = E(a, i) * virtual(a) * occupied(i)Eai (generic function with 1 method)

Entry 1:

julia> simplify_heavy(act_on_ket(H))2 ∑_i(F_ii)
- 2 ∑_ij(g_iijj)
+ ∑_ij(g_ijji)
+ ∑_ai(F_ai E_ai)
+ 1/2 ∑_aibj(g_aibj E_ai E_bj)

Since we do not have both the one electron matrix $h_{pq}$ and the Fock matrix $F_{pq}$ we do not perfectly reproduce the form of the HF energy from box 13.2, but the rest of the box is reproduced nicely.

Entry 2:

julia> simplify_heavy(act_on_ket(commutator(H, Eai(1, 2))));
julia> look_for_tensor_replacements(ans, make_exchange_transformer("g", "L"))2 F_ia
- ∑_j(F_ij E_aj)
+ ∑_b(F_ba E_bi)
+ ∑_bj(L_iabj E_bj)
- ∑_jbk(g_ijbk E_aj E_bk)
+ ∑_bcj(g_bacj E_bi E_cj)

Entry 3:

julia> simplify_heavy(act_on_ket(commutator(H, Eai(1, 2), Eai(3, 4))));
julia> look_for_tensor_replacements(ans, make_exchange_transformer("g", "L"));
julia> r, ss, ns = desymmetrize(ans, make_permutation_mappings([(1, 2), (3, 4)]));
julia> r # Terms that have P_ij^ab in front- F_ib E_aj
- ∑_k(L_ikjb E_ak)
+ ∑_c(L_jbca E_ci)
- ∑_ck(g_ibck E_aj E_ck)
- ∑_kc(g_ikcb E_ak E_cj)
julia> ss # Symmetric terms2 L_iajb
+ ∑_kl(g_ikjl E_ak E_bl)
+ ∑_cd(g_cadb E_ci E_dj)
julia> ns # Non symmetric terms0

Entry 4:

julia> simplify_heavy(act_on_ket(commutator(H, Eai(1, 2), Eai(3, 4), Eai(5, 6))));
julia> look_for_tensor_replacements(ans, make_exchange_transformer("g", "L"));
julia> r, ss, ns = desymmetrize(ans,
           make_permutation_mappings([(1, 2), (3, 4), (5, 6)]));
julia> r # Terms that have P_ijk^abc in front- L_ibkc E_aj
+ ∑_l(g_ibkl E_aj E_cl)
- ∑_d(g_ibdc E_aj E_dk)
julia> ss # Symmetric terms0
julia> ns # Non symmetric terms0

Entry 5:

julia> x = simplify_heavy(
           act_on_ket(commutator(H, Eai(1, 2), Eai(3, 4), Eai(5, 6), Eai(7, 8))));
julia> r, ss, ns = desymmetrize(x,
           make_permutation_mappings([(1, 2), (3, 4), (5, 6)]));
julia> r # The first part with P_ijk^abcg_ibkd E_aj E_cl
+ g_iblc E_aj E_dk
julia> r, ss, ns = desymmetrize(x,
           make_permutation_mappings([(1, 2), (3, 4), (5, 6), (7, 8)]));
julia> r # The second part with P_ijkl^abcd1/2 g_ibkd E_aj E_cl