Quantum Electrodynamics CCSD Example

Hamiltonian

In the QED-CC method for a single photon mode the Hamiltonian is of the form

$H = H_e + \sum_{pq}{d_{pq} E_{pq} (b + b^\dagger)} + \omega b^\dagger b$

with $H_e$ being the electronic Hamiltonian with adjusted integrals due to the dipole self-energy, and $b$ and $b^\dagger$ are boson annihilation and creation operators respectively for the cavity mode.

We define the electronic hamiltonian as before

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 * simplify(
           ∑(psym_tensor("g", 1:4...) * e(1:4...) * electron(1:4...), 1:4)
       );
julia> He = h + g∑_pq(F_pq E_pq)
- 2 ∑_pqi(g_pqii E_pq)
- 1/2 ∑_pqr(g_prrq E_pq)
+ ∑_pqi(g_piiq E_pq)
+ 1/2 ∑_pqrs(g_pqrs E_pq E_rs)

Then the bilinear and pure photon terms

julia> b = boson()b⁻
julia> H_bilinear = ∑(
           real_tensor("d", 1, 2) * E(1, 2) * (b' + b) * electron(1, 2), 1:2)∑_pq(d_pq E_pq b⁻)
+ ∑_pq(d_pq E_pq b†)
julia> H_photon = real_tensor("ω") * b'bω b† b⁻
julia> H = He + H_bilinear + H_photon∑_pq(F_pq E_pq)
- 2 ∑_pqi(g_pqii E_pq)
- 1/2 ∑_pqr(g_prrq E_pq)
+ ∑_pqi(g_piiq E_pq)
+ 1/2 ∑_pqrs(g_pqrs E_pq E_rs)
+ ∑_pq(d_pq E_pq b⁻)
+ ∑_pq(d_pq E_pq b†)
+ ω b† b⁻

Cluster Operator

For the QED-CCSD method the cluster operator is given by

$T = T1 + T2 + S0 + S1 + S2$

with the separate parts defined as

julia> T2 = 1//2 * ∑(psym_tensor("t", 1:4...) * E(1, 2) * E(3, 4) *
       occupied(2, 4) * virtual(1, 3), 1:4)1/2 ∑_aibj(t_aibj E_ai E_bj)
julia> S0 = real_tensor("s0") * b's0 b†
julia> S1 = ∑(real_tensor("s", 1, 2) * E(1, 2) * b' * occupied(2) * virtual(1), 1:2)∑_ai(s_ai E_ai b†)
julia> S2 = 1//2 * ∑(psym_tensor("s", 1:4...) * E(1, 2) * E(3, 4) * b' *
       occupied(2, 4) * virtual(1, 3), 1:4)1/2 ∑_aibj(s_aibj E_ai E_bj b†)

As for electronic CC, we disregard T1 as this is included in the integrals.

Energy

julia> Hbar_ccsd = simplify(bch(He, T2, 4));
julia> Hbar_ket_ccsd = simplify(act_on_ket(Hbar_ccsd, 2));
julia> E_ccsd = simplify_heavy(act_on_bra(Hbar_ket_ccsd))2 ∑_i(F_ii)
- 2 ∑_ij(g_iijj)
+ ∑_ij(g_ijji)
+ 2 ∑_iajb(g_iajb t_aibj)
- ∑_iajb(g_iajb t_ajbi)
julia> T = [T2, S0, S1, S2];
julia> Hbar = simplify(bch(H, T, 4));
julia> Hbar_ket = simplify(act_on_ket(Hbar, 3)); # Keeping 3 operators (E_ai E_bj b')
julia> E0 = simplify_heavy(act_on_bra(Hbar_ket))2 ∑_i(F_ii)
- 2 ∑_ij(g_iijj)
+ ∑_ij(g_ijji)
+ 2 ∑_i(s0 d_ii)
+ 2 ∑_ia(d_ia s_ai)
+ 2 ∑_iajb(g_iajb t_aibj)
- ∑_iajb(g_iajb t_ajbi)
julia> ΔE_QED_CCSD = E0 - E_ccsd # Additional terms from ccsd2 ∑_i(s0 d_ii)
+ 2 ∑_ia(d_ia s_ai)

Omega

Singles

julia> omega_ai_ccsd = project_biorthogonal(Hbar_ket_ccsd, E(1, 2));
julia> omega_0_ai = project_biorthogonal(Hbar_ket, E(1, 2));
julia> omega_0_ai - omega_ai_ccsd;
julia> look_for_tensor_replacements(ans, make_exchange_transformer("t", "u"));
julia> look_for_tensor_replacements(ans, make_exchange_transformer("s", "v"));
julia> Δomega_ai = anss0 d_ai
+ ∑_b(d_ab s_bi)
- ∑_j(d_ji s_aj)
+ 2 ∑_j(d_jj s_ai)
+ ∑_jb(d_jb v_aibj)
+ ∑_jb(s0 d_jb u_aibj)

Doubles

julia> omega_aibj_ccsd = project_biorthogonal(Hbar_ket_ccsd, E(1, 2) * E(3, 4));
julia> omega_0_aibj = project_biorthogonal(Hbar_ket, E(1, 2) * E(3, 4));
julia> omega_0_aibj - omega_aibj_ccsd;
julia> symmetrize(ans, make_permutation_mappings([(1, 2), (3, 4)]));
julia> simplify_heavy(ans);
julia> look_for_tensor_replacements(ans, make_exchange_transformer("t", "u"));
julia> Δomega_aibj_r, Δomega_aibj_ss, Δomega_aibj_ns =
           desymmetrize(ans, make_permutation_mappings([(1, 2), (3, 4)]));
julia> Δomega_aibj_rd_ai s_bj
+ ∑_c(d_ac s_bjci)
- ∑_k(d_ki s_akbj)
+ ∑_c(s0 d_ac t_bjci)
- ∑_k(s0 d_ki t_akbj)
- ∑_kc(d_kc s_ak t_bjci)
- ∑_kc(d_kc s_ci t_akbj)
+ ∑_kc(d_kc s_ai u_bjck)
julia> Δomega_aibj_ss2 ∑_k(d_kk s_aibj)
julia> Δomega_aibj_ns0

Photon

julia> project_biorthogonal(Hbar_ket, b');
julia> look_for_tensor_replacements(ans, make_exchange_transformer("s", "v"));
julia> omega_1 = ans2 ∑_i(d_ii)
+ s0 ω
+ 2 ∑_ia(F_ia s_ai)
+ ∑_iajb(g_iajb v_aibj)
+ 2 ∑_ia(s0 d_ia s_ai)

Singles Photon

julia> project_biorthogonal(Hbar_ket, E(1, 2) * b');
julia> simplify_heavy(ans);
julia> look_for_tensor_replacements(ans, make_exchange_transformer("t", "u"));
julia> look_for_tensor_replacements(ans, make_exchange_transformer("s", "v"));
julia> look_for_tensor_replacements(ans, make_exchange_transformer("g", "L"));
julia> omega_1_ai = ansd_ai
+ ω s_ai
+ ∑_b(F_ab s_bi)
- ∑_j(F_ji s_aj)
+ ∑_jb(F_jb v_aibj)
+ ∑_jb(d_jb u_aibj)
+ ∑_bj(s_bj L_aijb)
+ ∑_bjc(g_abjc v_bicj)
- ∑_jkb(g_jikb v_ajbk)
+ ∑_b(s0 d_ab s_bi)
- ∑_j(s0 d_ji s_aj)
+ ∑_jb(s0 d_jb v_aibj)
+ 2 ∑_jb(d_jb s_ai s_bj)
- 2 ∑_jb(d_jb s_aj s_bi)
+ ∑_bjkc(s_bj L_jbkc u_aick)
- ∑_jbkc(s_aj g_jbkc u_bick)
- ∑_bjkc(s_bi g_jbkc u_ajck)

Doubles Photon

julia> project_biorthogonal(Hbar_ket, E(1, 2) * E(3, 4) * b');
julia> symmetrize(ans, make_permutation_mappings([(1, 2), (3, 4)]));
julia> simplify_heavy(ans);
julia> look_for_tensor_replacements(ans, make_exchange_transformer("t", "u"));
julia> look_for_tensor_replacements(ans, make_exchange_transformer("s", "v"));
julia> look_for_tensor_replacements(ans, make_exchange_transformer("g", "L"));
julia> omega_1_aibj_r, omega_1_aibj_ss, omega_1_aibj_ns =
           desymmetrize(ans, make_permutation_mappings([(1, 2), (3, 4)]));
julia> omega_1_aibj_r∑_c(F_ac s_bjci)
- ∑_k(F_ki s_akbj)
+ ∑_c(d_ac t_bjci)
- ∑_k(d_ki t_akbj)
- ∑_k(s_ak g_bjki)
+ ∑_c(s_ci g_acbj)
- ∑_ck(g_acki s_bjck)
- ∑_ck(g_ackj s_bkci)
+ ∑_kc(g_aikc v_bjck)
+ ∑_c(s0 d_ac s_bjci)
- ∑_k(s0 d_ki s_akbj)
- ∑_kc(F_kc s_ak t_bjci)
- ∑_kc(F_kc s_ci t_akbj)
+ ∑_c(d_ac s_bj s_ci)
- ∑_k(d_ki s_ak s_bj)
- 2 ∑_kc(d_kc s_ak s_bjci)
- 2 ∑_kc(d_kc s_ci s_akbj)
+ ∑_kc(d_kc s_ai v_bjck)
+ ∑_ckd(s_ck L_adkc t_bjdi)
- ∑_ckl(s_ck L_kcli t_albj)
- ∑_kcd(s_ak g_bckd t_cjdi)
+ ∑_kcl(s_ak g_kcli t_bjcl)
+ ∑_kcl(s_ak g_kclj t_blci)
- ∑_cdk(s_ci g_adkc t_bjdk)
- ∑_cdk(s_ci g_bdkc t_akdj)
+ ∑_ckl(s_ci g_kjlc t_albk)
- ∑_klc(s_ak g_kilc u_bjcl)
+ ∑_ckd(s_ci g_ackd u_bjdk)
+ ∑_kcld(g_kcld s_akdi t_bjcl)
+ ∑_kcld(g_kcld s_akdj t_blci)
- ∑_kcld(g_kcld s_aibk u_cjdl)
- ∑_kcld(g_kcld s_aicj u_bkdl)
- ∑_kcld(g_kcld s_aicl u_bjdk)
- ∑_kcld(g_kcld t_aibk v_cjdl)
- ∑_kcld(g_kcld t_aicj v_bkdl)
+ ∑_kcld(g_kcld u_aick v_bjdl)
- ∑_kc(s0 d_kc s_ak t_bjci)
- ∑_kc(s0 d_kc s_ci t_akbj)
julia> omega_1_aibj_ssω s_aibj
+ ∑_cd(g_acbd s_cidj)
+ ∑_kl(g_kilj s_akbl)
+ 2 ∑_kc(d_kc s_ck s_aibj)
+ ∑_kcld(g_kcld s_akbl t_cidj)
+ ∑_kcld(g_kcld s_cidj t_akbl)
julia> omega_1_aibj_ns0