Coupled Cluster Singles and Doubles Example

Here is a detailed example of how to derive common expressions and equations for the closed shell CCSD method. This will also serve as a way of showcasing the various functionality of the package.

Hamiltonian

For the hamiltonian we will use the familiar expression for the non-relativistic electronic hamiltonian, however we will express it in terms of the fock matrix instead of the normal one-electron matrix

$F_{pq} = h_{pq} + \sum_{i}{\left(2 g_{pqii} - g_{piiq}\right)}$

$h_{pq} = F_{pq} - \sum_{i}{\left(2 g_{pqii} - g_{piiq}\right)}$

This is because every expression except the HF energy looks simpler expressed in terms of the fock matrix. In code we then define the hamiltonian as

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

We omit the nuclear repulsion term and we do not assume any symmetry in our integrals apart from the particle exchange symmetry in the two electron integrals. We are assuming to work with T1 transformed integrals, which allow us to ignore the T1 operator when deriving the equations.

We can derive the HF energy expression for the Hamiltonian in terms of the Fock matrix as

julia> E_HF = simplify_heavy(hf_expectation_value(H))2 ∑_i(F_ii)
- 2 ∑_ij(g_iijj)
+ ∑_ij(g_ijji)

Cluster Operator

The cluster operator we define then as only the T2 operator, since the effect of the T1 amplitudes is included in the T1 transformed integrals.

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)

Similarity transformed Hamiltonian

A central operator in coupled cluster theory is the similarity transformed Hamiltonian operator

$\bar H = e^{-T} H e^T$

This is usually computed using the Baker-Campbell-Hausdorff (BCH) expansion

$e^{-B} A e^B = A + [A, B] + \frac12 [[A, B], B] + \frac{1}{3!} [[[A, B], B], B] + ...$

A nice feature of this expansion for the electronic hamiltionian and the cluster operator is that it truncates after the 4th order term (5 terms). This is because the 5th order expression is zero

$[[[[[H, E_{ai}], E_{bj}], E_{ck}], E_{dl}], E_{em}] = 0$

This means that we can explicitly compute the full similarity transformed Hamiltonian with no fear of missing any terms. We supply the function bch(A, B, n) which computes the BCH expansion to nth order.

In code this looks like

julia> Hbar = simplify(bch(H, T2, 4)); # Avoid printing as this is quite long

Ground State Energy Expression and Equations

The coupled cluster ground state energy and equations are obtained by projecting the similarity constrained Schrödinger equation

$\bar H |\text{HF}\rangle = E |\text{HF}\rangle$

onto a set of bra states. The energy is obtained by projecting on the HF bra

$E_0 = \langle\text{HF}| \bar H |\text{HF}\rangle$

while the coupled cluster wave function is determined by projecting on a set of excited bra states determined from the cluster operator

$\langle\mu| \bar H |\text{HF}\rangle = 0$

In both cases we require the quantity $\bar H |\text{HF}\rangle$ which we can explicitly compute by the function act_on_ket(ex, [max_ops]).

julia> Hbar_ket = simplify(act_on_ket(Hbar, 2)); # Avoid printing long expression

Here we choose to throw away any terms with more than two operators as in CCSD we are never projecting on any bra that is more than doubly excited. This can speed up the projection and further operations by a bit.

Energy

To obtain the energy expression we project the Hbar_ket on HF from the left. This can be nicely done using the act_on_bra(ex, [max_ops]) function

julia> E0 = simplify_heavy(act_on_bra(Hbar_ket))2 ∑_i(F_ii)
- 2 ∑_ij(g_iijj)
+ ∑_ij(g_ijji)
+ 2 ∑_iajb(g_iajb t_aibj)
- ∑_iajb(g_iajb t_ajbi)

Since the interesting part here is that which differs from HF, we can get the correlation energy as

julia> E_corr = E0 - E_HF2 ∑_iajb(g_iajb t_aibj)
- ∑_iajb(g_iajb t_ajbi)

Here we see a very common pattern: the 2 * coulomb - 1 * exchange pattern where we have terms that differ by a prefactor of $-\frac12$ and an exchange of index 2 and 4 of a 4 index tensor, such as the two electron integrals or the T2 amplitudes. We usually define separate tensor names for these patterns such as

$L_{pqrs} = 2 g_{pqrs} - g_{psrq}$

and

$u_{aibj} = 2 t_{aibj} - t_{ajbi}$

We provide a pair of functions to look for this pattern in expressions. Firstly and mainly look_for_tensor_replacements function which is the actual engine that takes an expression and a tensor transformation function and looks for terms that are equal up to this transformation. We provide the utility function make_exchange_transformer which produces a transformer that looks for the 2 * coulomb - 1 * exchange pattern. We can show this in action on the correlation energy looking for both the patterns mentioned above

julia> look_for_tensor_replacements(E_corr, make_exchange_transformer("g", "L"))∑_iajb(L_iajb t_aibj)
julia> look_for_tensor_replacements(E_corr, make_exchange_transformer("t", "u"))∑_iajb(g_iajb u_aibj)

Omega Equations

Biorthogonal Bra

For the singles equations we want to compute the quantity

$\Omega_i^a = \bar{\left\langle_i^a\right|} \bar H |\text{HF}\rangle$

where the bra $\bar{\left\langle_i^a\right|}$ is the biorthogonal bra state for the singly excited determinant

$\left|_i^a\right\rangle = E_{ai} |\text{HF}\rangle$

defined such that

$\left\langle\bar{_i^a}|_j^b\right\rangle = \delta_{ij} \delta_{ab}$

For a singly excited determinant we can explicitly express this biorthogonal bra simply as

$\bar{\left\langle_i^a\right|} = \frac12 \langle\text{HF}| E_{ia}$

As a quick example we can show this

julia> bra = 1//2 * E(2, 1) * occupied(2) * virtual(1)1/2 E_ia
julia> ket = E(3, 4) * occupied(4) * virtual(3)E_ai
julia> hf_expectation_value(bra * ket)δ_ab δ_ij

However for higher excitations the biorthogonal bra states can not be as easily expressed explicitly. For doubly excited determinants we want a biorthogonal bra state $\bar{\left\langle_{ij}^{ab}\right|}$ such that

$\left\langle\bar{_{ij}^{ab}}|_{kl}^{cd}\right\rangle = P_{ai,bj} \delta_{ik} \delta_{jl} \delta_{ac} \delta_{bd} = \delta_{ik} \delta_{jl} \delta_{ac} \delta_{bd} + \delta_{jk} \delta_{il} \delta_{bc} \delta_{ad}$

with the permutation operator $P_{ai,bj}$ generating the equivalent permutations due to the particle symmetry in the doubly excited determinant

$\left|_{ij}^{ab}\right\rangle = \left|_{ji}^{ba}\right\rangle = E_{ai} E_{bj} |\text{HF}\rangle = E_{bj} E_{ai} |\text{HF}\rangle$

We can write the biorthogonal bra explicitly as

$\bar{\left\langle_{ij}^{ab}\right|} = \frac13 \left\langle_{ij}^{ab}\right| + \frac16 \left\langle_{ji}^{ab}\right| = \langle\text{HF}| \left( \frac13 E_{jb} E_{ia} + \frac16 E_{ib} E_{ja} \right)$

Which we can confirm using the code

julia> bra = (1//3 * E(4, 3) * E(2, 1) + 1//6 * E(2, 3) * E(4, 1)) *
           occupied(2, 4) * virtual(1, 3)1/6 E_ib E_ja
+ 1/3 E_jb E_ia
julia> ket = E(5, 6) * E(7, 8) * occupied(6, 8) * virtual(5, 7)E_ai E_bj
julia> hf_expectation_value(bra * ket)δ_ac δ_ik δ_bd δ_jl
+ δ_ad δ_il δ_bc δ_jk

Since we do not want to (and often can not, see note above) make an explicit expression for a biorthogonal bra, we provide the function project_biorthogonal which lets us insert Kronecker deltas according to the definition of the biorthogonality, and the symmetrize function to expand the permutations required by the biorthogonality. We can reproduce the same behaviour as for the biorthogonal doubles from above as

julia> ket = E(5, 6) * E(7, 8) * occupied(6, 8) * virtual(5, 7)E_ai E_bj
julia> project_biorthogonal(ket, E(1, 2) * E(3, 4))δ_ac δ_ik δ_bd δ_jl
julia> symmetrize(ans, make_permutation_mappings([(1, 2), (3, 4)]))δ_ac δ_ik δ_bd δ_jl
+ δ_ad δ_il δ_bc δ_jk

Here we supplied the "template" ket E(1, 2) * E(3, 4) which is that we want to project on the biorthogonal bra of.

Singles Equations

We can use the project_biorthogonal function to project Hbar_ket on the singles biorthogonal bra and get an expression for $\Omega_i^a$

julia> omega_ai = project_biorthogonal(Hbar_ket, E(1, 2))F_ai
+ 2 ∑_jb(F_jb t_aibj)
- ∑_jb(F_jb t_ajbi)
+ 2 ∑_bjc(g_abjc t_bicj)
- ∑_bjc(g_abjc t_bjci)
- 2 ∑_jkb(g_jikb t_ajbk)
+ ∑_jkb(g_jikb t_akbj)

Here we also see the 2 * coulomb - 1 * exchange pattern show up again, so we can simplify by.

julia> omega_ai = look_for_tensor_replacements(omega_ai,
           make_exchange_transformer("t", "u"))F_ai
+ ∑_jb(F_jb u_aibj)
+ ∑_bjc(g_abjc u_bicj)
- ∑_jkb(g_jikb u_ajbk)

Which produces a nice simplified expression for the singles part of the CCSD equations.

Doubles Equations

Here we want to derive the expression for the doubles $\Omega$

$\Omega_{ij}^{ab} = \bar{\left\langle_{ij}^{ab}\right|} \bar H |\text{HF}\rangle$

Just like for the singles equations we project on a biorthogonal bra, though now using a doubly excited ket as the template, and we need to symmetrize to include the permutations caused by the particle symmetry of the biorthogonal bra.

julia> project_biorthogonal(Hbar_ket, E(1, 2) * E(3, 4));
julia> omega_aibj = simplify_heavy(
           symmetrize(ans, make_permutation_mappings([(1, 2), (3, 4)])));
julia> omega_aibj_u = look_for_tensor_replacements(omega_aibj,
           make_exchange_transformer("t", "u"))g_aibj
+ ∑_c(F_ac t_bjci)
+ ∑_c(F_bc t_aicj)
- ∑_k(F_ki t_akbj)
- ∑_k(F_kj t_aibk)
+ ∑_cd(g_acbd t_cidj)
- ∑_ck(g_acki t_bjck)
- ∑_ck(g_ackj t_bkci)
- ∑_ck(g_bcki t_akcj)
- ∑_ck(g_bckj t_aick)
+ ∑_kl(g_kilj t_akbl)
+ ∑_kc(g_aikc u_bjck)
+ ∑_kc(g_bjkc u_aick)
+ ∑_kcld(g_kcld t_aicl t_bkdj)
+ ∑_kcld(g_kcld t_akbl t_cidj)
+ ∑_kcld(g_kcld t_akdj t_blci)
- ∑_kcld(g_kcld t_aibk u_cjdl)
- ∑_kcld(g_kcld t_aicj u_bkdl)
- ∑_kcld(g_kcld t_akbj u_cidl)
- ∑_kcld(g_kcld t_bjci u_akdl)
- ∑_kcld(g_kcld t_bjcl u_aidk)
+ ∑_kcld(g_kcld u_aick u_bjdl)

Here we needed to use the (sometimes) expensive simplify_heavy function to fully simplify, as well as recognizing the 2 * coulomb - 1 * exchange pattern for the T2 amplitudes. Now we have a correct and rather nice expression for the doubles equations, however, because of the particle symmetry of the bra that we introduced with the call to symmetrize some terms are equal after permuting the indices (a, i) <-> (b, j). Computing all of these redundant terms leads to a lot of unnecessary terms to code and compute, so it is much more efficient to only code the non-redundant ones, then symmetrize the final result numerically. To help with this, we provide the desymmetrize function which sort of acts like the reverse of symmetrize. We can call this function on our $\Omega_{ij}^{ab}$ expression

julia> omega_aibj_r, omega_aibj_ss, omega_aibj_ns = desymmetrize(
           omega_aibj_u, make_permutation_mappings([(1, 2), (3, 4)])
       );
julia> omega_aibj_r∑_c(F_ac t_bjci)
- ∑_k(F_ki t_akbj)
- ∑_ck(g_acki t_bjck)
- ∑_ck(g_ackj t_bkci)
+ ∑_kc(g_aikc u_bjck)
- ∑_kcld(g_kcld t_aibk u_cjdl)
- ∑_kcld(g_kcld t_aicj u_bkdl)
julia> omega_aibj_ssg_aibj
+ ∑_cd(g_acbd t_cidj)
+ ∑_kl(g_kilj t_akbl)
+ ∑_kcld(g_kcld t_akbl t_cidj)
+ ∑_kcld(g_kcld t_akdj t_blci)
+ ∑_kcld(g_kcld u_aick u_bjdl)
julia> omega_aibj_ns∑_kcld(g_kcld t_aicl t_bkdj)
- ∑_kcld(g_kcld t_bjcl u_aidk)

julia> using Serialization
julia> serialize("omega_aibj_r", omega_aibj_r)

The desymmetrize function returns three expressions. The first omega_aibj_r contains the terms for which redundant permutations exist elsewhere in the original expression. These needs to be symmetrized to obtain the correct expression. The second term omega_aibj_ss contains the "self-symmetric" terms. These are the terms which themselves carry the symmetry requested and can be evaluated as-is. The third and last term omega_aibj_ns are the "non-symmetric" terms. Interestingly, we have here gotten two terms in the non-symmetric expression, even though the result should be symmetric. Under closer inspection one can see that the two "non-symmetric" terms are indeed symmetric if one re-expands the $u_{aidk}$ tensor in terms of $t$. To avoid having non-symmetric terms we can instead directly desymmetrize the omega before looking for coulomb - exchange symmetry.

julia> omega_aibj_r, omega_aibj_ss, omega_aibj_ns = desymmetrize(
           omega_aibj, make_permutation_mappings([(1, 2), (3, 4)])
       );
julia> look_for_tensor_replacements(omega_aibj_r, make_exchange_transformer("g", "L"))∑_c(F_ac t_bjci)
- ∑_k(F_ki t_akbj)
+ ∑_kc(L_aikc t_bjck)
- ∑_kc(g_aikc t_bkcj)
- ∑_ck(g_ackj t_bkci)
- ∑_kcld(L_kcld t_aibk t_cjdl)
- ∑_kcld(L_kcld t_aicj t_bkdl)
- ∑_kcld(L_kcld t_aick t_bldj)
julia> look_for_tensor_replacements(omega_aibj_ss, make_exchange_transformer("g", "L"))g_aibj
+ ∑_cd(g_acbd t_cidj)
+ ∑_kl(g_kilj t_akbl)
+ 2 ∑_kcld(L_kcld t_aick t_bjdl)
+ ∑_kcld(g_kcld t_akbl t_cidj)
+ ∑_kcld(g_kcld t_akci t_bldj)
+ ∑_kcld(g_kcld t_akdj t_blci)
julia> omega_aibj_ns0

where we see that the non-symmetric part is indeed zero, but we have lost the ability to remove all prefactors. Counting total terms we can see that both strategies give 15 terms in total, so the cost should not be much different, but it can be nice to keep in mind. In any case, to evaluate the full doubles omega one needs to symmetrize the redundant expression, and add the direct evaluation of the self-symmetric and apparent non-symmetric expression which we can write as

$\Omega_{ij}^{ab} = (\Omega_{\text{ss}})_{ij}^{ab} + (\Omega_{\text{ns}})_{ij}^{ab} + P_{ij}^{ab} (\Omega_{\text{s}})_{ij}^{ab}$

Jacobian Transformation

The coupled cluster Jacobian is a central quantity for computing various properties of the coupled cluster wave function, such as density matrices and excitation energies. The Jacobian matrix $A_{\mu\nu}$ is given by

$A_{\mu\nu} = \langle\mu| e^{-T} [H, \tau_{\nu}] e^T |\text{HF}\rangle$

though one is rarely interested in the full matrix, but rather to compute eigenvalues and solve linear systems. The matrix is usually huge and very structured (though not very sparse), so it is way more efficient to code a direct linear transformation. The Jacobian is not symmetric so one needs to derive both the right and left transformations given by

$\rho_\mu = \sum_{\nu}{A_{\mu\nu} c_\nu}$

$\sigma_\nu = \sum_{\mu}{b_\mu A_{\mu\nu}}$

We can quite nicely derive expressions for both of these transformations using the code.

Right transformation

The transformed vector $\rho_\mu$ can be split into singles and doubles block like

$\rho_{ai} = \sum_{\nu}{A_{ai,\nu} c_\nu}$

$\rho_{aibj} = \sum_{\nu}{A_{aibj,\nu} c_\nu}$

Singles block

We then split the sum into blocks

$\rho_{ai} = \sum_{\nu}{A_{ai,\nu} c_\nu} = \sum_{ck}{A_{ai,ck} c_{ck}} + \sum_{ck \geq dl}{A_{ai,ckdl} c_{ckdl}}$

The singles sum is okay, but the triangular doubles sum $ck \geq dl$ is not something we can express using the sums in the code, so we have to square it up.

$\sum_{ck \geq dl}{A_{ai,ckdl} c_{ckdl}} = \frac12 \sum_{ckdl}{(1 + \delta_{ck,dl}) A_{ai,ckdl} c_{ckdl}} = \frac12 \sum_{ckdl}{A_{ai,ckdl} \tilde c_{ckdl}}$

Here we have to add a factor of 2 on the diagonal of $c_{ckdl}$ by defining the adjusted tensor

$\tilde c_{ckdl} = (1 + \delta_{ck,dl}) c_{ckdl}$

though we will not be using the $\tilde c$ name in the derivation, but rather remembering that we have to numerically scale the diagonal of $c_{ckdl}$ by 2 when coding the final expression.

We start by computing the projection

$\sum_{ck}{e^{-T} [H, E_{ck}] e^T |\text{HF}\rangle c_{ck}}$

julia> simplify(commutator(H, E(1, 2) * occupied(2) * virtual(1)));
julia> simplify(bch(ans, T2, 4));
julia> proj_ai = simplify(act_on_ket(ans, 2));
julia> proj_ck = simplify(∑(ans * real_tensor("c", 1, 2), [1, 2]));

This projection will also be useful when doing the doubles block. We also keep the proj_ai for the left transformation later. To get the singles block out we project on a biorthogonal singles bra

julia> ρ_ai_singles = project_biorthogonal(proj_ck, E(1, 2));
julia> ρ_ai_singles = look_for_tensor_replacements(ρ_ai_singles,
           make_exchange_transformer("t", "u"));
julia> ρ_ai_singles = look_for_tensor_replacements(ρ_ai_singles,
           make_exchange_transformer("g", "L"));
julia> ρ_ai_singles∑_b(F_ab c_bi)
- ∑_j(F_ji c_aj)
+ ∑_bj(c_bj L_aijb)
+ ∑_bjkc(c_bj L_jbkc u_aick)
- ∑_jbkc(c_aj g_jbkc u_bick)
- ∑_bjkc(c_bi g_jbkc u_ajck)

where we get a very concise expression for the singles-singles transformation. Next we compute the projection

$\frac12 \sum_{ckdl}{e^{-T} [H, E_{ck} E_{dl}] e^T |\text{HF}\rangle c_{ckdl}}$

julia> simplify(commutator(H, E(1, 2) * E(3, 4) * occupied(2, 4) * virtual(1, 3)));
julia> simplify(bch(ans, T2, 4));
julia> proj_aibj = simplify(act_on_ket(ans, 2));
julia> proj_ckdl = simplify(1//2 * ∑(ans * psym_tensor("c", 1:4...), 1:4));

This will also be useful for the doubles block later and we also keep the proj_aibj for the left transformation. Now we project on the singles bra to get

julia> ρ_ai_doubles = project_biorthogonal(proj_ckdl, E(1, 2));
julia> ρ_ai_doubles = look_for_tensor_replacements(ρ_ai_doubles,
           make_exchange_transformer("g", "L"));
julia> ρ_ai_doubles = look_for_tensor_replacements(ρ_ai_doubles,
           make_exchange_transformer("c", "cu"));
julia> ρ_ai_doubles∑_jb(F_jb cu_aibj)
+ ∑_bjc(L_abjc c_bicj)
- ∑_jkb(L_jikb c_ajbk)

Here we introduced another tensor pattern

$cu_{aibj} = 2 c_{aibj} - c_{ajbi}$

We can now combine the two parts of the $\rho_{ai}$ block

julia> ρ_ai = ρ_ai_singles + ρ_ai_doubles∑_b(F_ab c_bi)
- ∑_j(F_ji c_aj)
+ ∑_jb(F_jb cu_aibj)
+ ∑_bj(c_bj L_aijb)
+ ∑_bjc(L_abjc c_bicj)
- ∑_jkb(L_jikb c_ajbk)
+ ∑_bjkc(c_bj L_jbkc u_aick)
- ∑_jbkc(c_aj g_jbkc u_bick)
- ∑_bjkc(c_bi g_jbkc u_ajck)

Doubles block

We again split the sum into singles and doubles, squareup and include the factor of 2 on the diagonal

$\rho_{aibj} = \sum_{\nu}{A_{aibj,\nu} c_\nu} = \sum_{ck}{A_{aibj,ck} c_{ck}} + \frac12 \sum_{ckdl}{A_{aibj,ckdl} \tilde c_{ckdl}}$

Here we can reuse the right projection from above, projecting on the biorthogonal doubles bra instead.

First the singles sum

julia> project_biorthogonal(proj_ck, E(1, 2) * E(3, 4));
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("g", "L"));
julia> ρ_aibj_singles_r, ρ_aibj_singles_ss, ρ_aibj_singles_ns =
           desymmetrize(ans, make_permutation_mappings([(1, 2), (3, 4)]));
julia> ρ_aibj_singles_r- ∑_k(c_ak g_bjki)
+ ∑_c(c_ci g_acbj)
- ∑_kc(F_kc c_ak t_bjci)
- ∑_kc(F_kc c_ci t_akbj)
+ ∑_ckd(c_ck L_adkc t_bjdi)
- ∑_ckl(c_ck L_kcli t_albj)
- ∑_kcd(c_ak g_bckd t_cjdi)
+ ∑_kcl(c_ak g_kcli t_bjcl)
+ ∑_kcl(c_ak g_kclj t_blci)
- ∑_cdk(c_ci g_adkc t_bjdk)
- ∑_cdk(c_ci g_bdkc t_akdj)
+ ∑_ckl(c_ci g_kjlc t_albk)
- ∑_klc(c_ak g_kilc u_bjcl)
+ ∑_ckd(c_ci g_ackd u_bjdk)
julia> ρ_aibj_singles_ss0
julia> ρ_aibj_singles_ns0

Then the doubles sum

julia> project_biorthogonal(proj_ckdl, E(1, 2) * E(3, 4));
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("g", "L"));
julia> ρ_aibj_doubles_r, ρ_aibj_doubles_ss, ρ_aibj_doubles_ns =
           desymmetrize(ans, make_permutation_mappings([(1, 2), (3, 4)]));
julia> ρ_aibj_doubles_r∑_c(F_ac c_bjci)
- ∑_k(F_ki c_akbj)
+ ∑_kc(L_aikc c_bjck)
- ∑_kc(c_akci g_bjkc)
- ∑_kc(c_akcj g_bcki)
- ∑_kcld(L_kcld c_akdl t_bjci)
- ∑_kcld(L_kcld c_cidl t_akbj)
+ ∑_kcld(L_kcld c_aick u_bjdl)
+ ∑_kcdl(c_akci g_kdlc t_bjdl)
+ ∑_kcdl(c_akcj g_kdlc t_bldi)
- ∑_kcld(c_aibk g_kcld u_cjdl)
- ∑_ckld(c_aicj g_kcld u_bkdl)
- ∑_kcld(c_akci g_kcld u_bjdl)
julia> ρ_aibj_doubles_ss∑_kl(c_akbl g_kilj)
+ ∑_cd(c_cidj g_acbd)
+ ∑_klcd(c_akbl g_kcld t_cidj)
+ ∑_cdkl(c_cidj g_kcld t_akbl)
julia> ρ_aibj_doubles_ns0

Combining singles and doubles

julia> ρ_aibj_r = ρ_aibj_singles_r + ρ_aibj_doubles_r∑_c(F_ac c_bjci)
- ∑_k(F_ki c_akbj)
- ∑_k(c_ak g_bjki)
+ ∑_c(c_ci g_acbj)
+ ∑_kc(L_aikc c_bjck)
- ∑_kc(c_akci g_bjkc)
- ∑_kc(c_akcj g_bcki)
- ∑_kc(F_kc c_ak t_bjci)
- ∑_kc(F_kc c_ci t_akbj)
+ ∑_ckd(c_ck L_adkc t_bjdi)
- ∑_ckl(c_ck L_kcli t_albj)
- ∑_kcd(c_ak g_bckd t_cjdi)
+ ∑_kcl(c_ak g_kcli t_bjcl)
+ ∑_kcl(c_ak g_kclj t_blci)
- ∑_cdk(c_ci g_adkc t_bjdk)
- ∑_cdk(c_ci g_bdkc t_akdj)
+ ∑_ckl(c_ci g_kjlc t_albk)
- ∑_klc(c_ak g_kilc u_bjcl)
+ ∑_ckd(c_ci g_ackd u_bjdk)
- ∑_kcld(L_kcld c_akdl t_bjci)
- ∑_kcld(L_kcld c_cidl t_akbj)
+ ∑_kcld(L_kcld c_aick u_bjdl)
+ ∑_kcdl(c_akci g_kdlc t_bjdl)
+ ∑_kcdl(c_akcj g_kdlc t_bldi)
- ∑_kcld(c_aibk g_kcld u_cjdl)
- ∑_ckld(c_aicj g_kcld u_bkdl)
- ∑_kcld(c_akci g_kcld u_bjdl)
julia> ρ_aibj_ss = ρ_aibj_doubles_ss∑_kl(c_akbl g_kilj)
+ ∑_cd(c_cidj g_acbd)
+ ∑_klcd(c_akbl g_kcld t_cidj)
+ ∑_cdkl(c_cidj g_kcld t_akbl)

which concludes the jacobian right transformation.

Left transformation

The procedure for the left transformation is very similar to the right transformation with some key differences. We can still split the transformed vector into singles and doubles blocks

The transformed vector $\rho_\mu$ can be split into singles and doubles blocks like

$\sigma_{ai} = \sum_{\mu}{b_\mu A_{\mu,ai}}$

$\sigma_{aibj} = \sum_{\mu}{b_\mu A_{\mu,aibj}}$

Singles block

We again further split the sum into singles and doubles parts

$\sigma_{ai} = \sum_{\mu}{b_\mu A_{\mu,ai}} = \sum_{ck}{b_{ck} A_{ck,ai}} + \frac12 \sum_{ckdl}{(1 + δ_{ck,dl}) b_{ckdl} A_{ckdl,ai}}$

though now we get a nice cancellation due to the fact that the jacobian is defined in terms of the biorthonormal bra, so we can define an adjusted jacobian

$\tilde A_{ckdl,\mu} = \frac{1}{1 + δ_{ck,dl}} A_{ckdl,\mu}$

defined in terms of the biorthogonal bra. The factors of 2 on the diagonal thus cancel and we are left with the much nicer expression

$\sigma_{ai} = \sum_{ck}{b_{ck} A_{ck,ai}} + \frac12 \sum_{ckdl}{b_{ckdl} \tilde A_{ckdl,ai}}$

where we do not have to scale either the input vector $b$ or the output vector $\sigma$ on the diagonals, contrary to the right transformation.

We can reuse the projection from above and can directly compute the singles projection

julia> simplify(
           ∑(project_biorthogonal(proj_ai, E(5, 6)) * real_tensor("b", 5, 6), 5:6)
       );
julia> look_for_tensor_replacements(ans, make_exchange_transformer("t", "u"));
julia> look_for_tensor_replacements(ans, make_exchange_transformer("g", "L"));
julia> σ_ai_singles = ans- ∑_j(F_ij b_aj)
+ ∑_b(F_ba b_bi)
+ ∑_bj(b_bj L_iabj)
+ ∑_bjkc(b_bj L_iakc u_bjck)
- ∑_jbkc(b_aj g_ibkc u_bjck)
- ∑_bjkc(b_bi g_jakc u_bjck)

and similarly for the doubles

julia> project_biorthogonal(proj_ai, E(5, 6) * E(7, 8)) * psym_tensor("b", 5:8...);
julia> simplify_heavy(∑(ans, 5:8)); # No need to symmetrize because of sum
julia> look_for_tensor_replacements(ans, make_exchange_transformer("t", "u"));
julia> look_for_tensor_replacements(ans, make_exchange_transformer("g", "L"));
julia> σ_ai_doubles = ans- ∑_jbk(b_ajbk g_ijbk)
+ ∑_bcj(b_bicj g_bacj)
- ∑_bjck(F_ib b_ajck t_bjck)
- ∑_jbck(F_ja b_bick t_bjck)
+ ∑_bcjdk(L_iabc b_bjdk t_cjdk)
- ∑_jkbcl(L_iajk b_bkcl t_bjcl)
- ∑_jbkcd(b_ajbk g_icbd t_cjdk)
+ ∑_jbkcl(b_ajbk g_iclj t_bkcl)
+ ∑_jbkcl(b_ajbk g_iclk t_blcj)
- ∑_bcjdk(b_bicj g_bdka t_cjdk)
- ∑_bcjdk(b_bicj g_cdka t_bkdj)
+ ∑_bcjkl(b_bicj g_kalj t_bkcl)
- ∑_jbklc(b_ajbk g_ijlc u_bkcl)
+ ∑_bcjkd(b_bicj g_bakd u_cjdk)

And combined

julia> σ_ai = σ_ai_singles + σ_ai_doubles- ∑_j(F_ij b_aj)
+ ∑_b(F_ba b_bi)
+ ∑_bj(b_bj L_iabj)
- ∑_jbk(b_ajbk g_ijbk)
+ ∑_bcj(b_bicj g_bacj)
- ∑_bjck(F_ib b_ajck t_bjck)
- ∑_jbck(F_ja b_bick t_bjck)
+ ∑_bjkc(b_bj L_iakc u_bjck)
- ∑_jbkc(b_aj g_ibkc u_bjck)
- ∑_bjkc(b_bi g_jakc u_bjck)
+ ∑_bcjdk(L_iabc b_bjdk t_cjdk)
- ∑_jkbcl(L_iajk b_bkcl t_bjcl)
- ∑_jbkcd(b_ajbk g_icbd t_cjdk)
+ ∑_jbkcl(b_ajbk g_iclj t_bkcl)
+ ∑_jbkcl(b_ajbk g_iclk t_blcj)
- ∑_bcjdk(b_bicj g_bdka t_cjdk)
- ∑_bcjdk(b_bicj g_cdka t_bkdj)
+ ∑_bcjkl(b_bicj g_kalj t_bkcl)
- ∑_jbklc(b_ajbk g_ijlc u_bkcl)
+ ∑_bcjkd(b_bicj g_bakd u_cjdk)

Doubles block

The same cancellations hold as above, so the procedure is almost identical, but using the other ket projection.

Singles:

julia> simplify(
           ∑(project_biorthogonal(proj_aibj, E(5, 6)) * real_tensor("b", 5, 6), 5:6)
       );
julia> look_for_tensor_replacements(ans, make_exchange_transformer("g", "L"));
julia> σ_aibj_singles_r, σ_aibj_singles_ss, σ_aibj_singles_ns =
           desymmetrize(ans, make_permutation_mappings([(1, 2), (3, 4)]));
julia> σ_aibj_singles_r2 F_ia b_bj
- F_ib b_aj
- ∑_k(b_ak L_ikjb)
+ ∑_c(b_ci L_jbca)
julia> σ_aibj_singles_ss0
julia> σ_aibj_singles_ns0

Doubles:

julia> project_biorthogonal(proj_aibj, E(5, 6) * E(7, 8)) * psym_tensor("b", 5:8...);
julia> simplify_heavy(∑(ans, 5:8)); # No need to symmetrize because of sum
julia> look_for_tensor_replacements(ans, make_exchange_transformer("t", "u"));
julia> look_for_tensor_replacements(ans, make_exchange_transformer("g", "L"));
julia> σ_aibj_doubles_r, σ_aibj_doubles_ss, σ_aibj_doubles_ns =
           desymmetrize(ans, make_permutation_mappings([(1, 2), (3, 4)]));
julia> σ_aibj_doubles_r- ∑_k(F_ik b_akbj)
+ ∑_c(F_ca b_bjci)
+ ∑_ck(L_iack b_bjck)
- ∑_ck(b_ajck g_ibck)
- ∑_kc(b_akcj g_ikcb)
- ∑_ckdl(L_iajc b_bkdl t_ckdl)
- ∑_kcdl(L_iakb b_cjdl t_ckdl)
+ ∑_kcdl(L_iakc b_bjdl u_ckdl)
+ ∑_ckdl(b_ajck g_idlb t_ckdl)
+ ∑_kcdl(b_akcj g_idlb t_cldk)
- ∑_kcld(b_aibk g_jcld u_ckdl)
- ∑_ckld(b_aicj g_kbld u_ckdl)
- ∑_ckld(b_ajck g_ibld u_ckdl)
julia> σ_aibj_doubles_ss∑_kl(b_akbl g_ikjl)
+ ∑_cd(b_cidj g_cadb)
+ ∑_klcd(b_akbl g_icjd t_ckdl)
+ ∑_cdkl(b_cidj g_kalb t_ckdl)
julia> σ_aibj_doubles_ns0

Combined:

julia> σ_aibj_r = σ_aibj_singles_r + σ_aibj_doubles_r2 F_ia b_bj
- F_ib b_aj
- ∑_k(F_ik b_akbj)
+ ∑_c(F_ca b_bjci)
- ∑_k(b_ak L_ikjb)
+ ∑_c(b_ci L_jbca)
+ ∑_ck(L_iack b_bjck)
- ∑_ck(b_ajck g_ibck)
- ∑_kc(b_akcj g_ikcb)
- ∑_ckdl(L_iajc b_bkdl t_ckdl)
- ∑_kcdl(L_iakb b_cjdl t_ckdl)
+ ∑_kcdl(L_iakc b_bjdl u_ckdl)
+ ∑_ckdl(b_ajck g_idlb t_ckdl)
+ ∑_kcdl(b_akcj g_idlb t_cldk)
- ∑_kcld(b_aibk g_jcld u_ckdl)
- ∑_ckld(b_aicj g_kbld u_ckdl)
- ∑_ckld(b_ajck g_ibld u_ckdl)
julia> σ_aibj_ss = σ_aibj_doubles_ss∑_kl(b_akbl g_ikjl)
+ ∑_cd(b_cidj g_cadb)
+ ∑_klcd(b_akbl g_icjd t_ckdl)
+ ∑_cdkl(b_cidj g_kalb t_ckdl)

Which concludes the left transformation.

$\eta_\nu$ expression

Another quantity related to the Jacobian is the $\eta_\nu$ expression given by

$\eta_\nu = \langle\text{HF}| e^{-T} [H, \tau_\nu] e^T |\text{HF}\rangle$

Since we already have the projections

$e^{-T} [H, \tau_\nu] e^T |\text{HF}\rangle$

from the jacobian transformation, the expressions for $\eta$ are simply obtained by projecting those on a HF bra.

julia> η_ai = act_on_bra(proj_ai)2 F_ia
julia> act_on_bra(proj_aibj);
julia> look_for_tensor_replacements(ans, make_exchange_transformer("g", "L"));
julia> η_aibj = ans2 L_iajb

One-electron density matrices

One frequently wants the one-electron density matrix to get properties like dipole, quadrupole, magnetic moments, etc. both for a state and transition moments. In coupled cluster theory we have asymmetric density matrices which we can express as

$D_{pq} = \langle L| E_{pq} |R\rangle$

where the left and right states are given in general by

$\langle L| = \left( L_0 \langle\text{HF}| + \sum_{\mu}{L_\mu \langle\mu|} \right) e^{-T}$

$|R\rangle = e^T \left( |\text{HF}\rangle R_0 + \sum_{\nu}{|\nu\rangle R_\nu} \right)$

The left and right coefficient vectors can represent either ground or excited states giving us the possibility to derive general expressions for ground and excited states, as well as transition densities.

Firstly we will need the similarity transformed operator for all parts

julia> bch_Epq = bch(E(1, 2) * electron(1, 2), T2, 4)E_pq
- 1/2 ∑_ajbk(δ_ij t_ajbk E_bk E_ap)
- 1/2 ∑_ajbk(δ_ik t_ajbk E_bp E_aj)
+ 1/2 ∑_bicj(δ_ab t_bicj E_cj E_pi)
+ 1/2 ∑_bicj(δ_ac t_bicj E_pj E_bi)
- 1/8 ∑_bjckdlem(δ_ij δ_ad t_bjck t_dlem E_ck E_bl E_em)
- 1/8 ∑_bjckdlem(δ_ij δ_ae t_bjck t_dlem E_ck E_dl E_bm)
- 1/8 ∑_bjckdlem(δ_ik δ_ad t_bjck t_dlem E_cl E_em E_bj)
- 1/8 ∑_bjckdlem(δ_ik δ_ae t_bjck t_dlem E_dl E_cm E_bj)
- 1/8 ∑_bjckdlem(δ_il δ_ab t_bjck t_dlem E_ck E_dj E_em)
- 1/8 ∑_bjckdlem(δ_il δ_ac t_bjck t_dlem E_dk E_em E_bj)
- 1/8 ∑_bjckdlem(δ_im δ_ab t_bjck t_dlem E_ck E_dl E_ej)
- 1/8 ∑_bjckdlem(δ_im δ_ac t_bjck t_dlem E_dl E_ek E_bj)

Ref-Ref contributions

The terms arising from the $L_0$ and $R_0$ (reference-reference) are quite straightforward to derive:

$D_{pq}^\text{ref-ref} = L_0 \langle\text{HF}| e^{-T} E_{pq} e^T |\text{HF}\rangle R_0$

julia> real_tensor("L0") * real_tensor("R0") * bch_Epq;
julia> hf_expectation_value(ans)2 δ_ij L0 R0

where we see that only the HF density survives.

$\mu$-Ref contributions

The terms arising from $L_\mu$ and $R_0$ are the only additional terms needed for the ground state density as $R_\nu = 0$ for the ground state. We will omit the $R_0$ in the derivation to reduce clutter. We can split the term into singles and doubles parts

$D_{pq}^{\mu\text{-ref}} = \sum_{\mu}{L_\mu \langle\mu|} e^{-T} E_{pq} e^T |\text{HF}\rangle \\= \sum_{ai}{L_{ai} \tilde{\langle_i^a|}} e^{-T} E_{pq} e^T |\text{HF}\rangle + \sum_{ai \geq bj}{L_{aibj} \tilde{\langle_{ij}^{ab}|}} e^{-T} E_{pq} e^T |\text{HF}\rangle$

in both cases we want the projection

$e^{-T} E_{pq} e^T |\text{HF}\rangle$

julia> ref_ket = simplify_heavy(act_on_ket(bch_Epq, 2));

Singles:

julia> ∑(project_biorthogonal(ref_ket, E(3, 4)) * real_tensor("L", 3, 4), 3:4);
julia> simplify(ans);
julia> look_for_tensor_replacements(ans, make_exchange_transformer("t", "u"))L_ai
+ ∑_bj(L_bj u_aibj)

For doubles we have the same cancellation as for the Jacobian left transformation

$\sum_{ai \geq bj}{L_{aibj} \tilde{\langle_{ij}^{ab}|}} e^{-T} E_{pq} e^T |\text{HF}\rangle = \frac12 \sum_{aibj}{L_{aibj} \bar{\langle_{ij}^{ab}|}} e^{-T} E_{pq} e^T |\text{HF}\rangle$

julia> project_biorthogonal(ref_ket, E(3, 4) * E(5, 6));
julia> ∑(ans * psym_tensor("L", 3:6...), 3:6);
julia> D_mu_ref = simplify_heavy(ans)∑_icj(L_aicj t_bicj)
- ∑_abk(L_ajbk t_aibk)

julia> D_ij = D_mu_ref * occupied(1, 2)- ∑_abk(L_ajbk t_aibk)
julia> D_ab = D_mu_ref * virtual(1, 2)∑_icj(L_aicj t_bicj)

julia> disable_external_index_translation()
julia> D_mu_ref # Here we see index ₁ and ₂ as we specified in E(1, 2)∑_iaj(L_₁iaj t_₂iaj C(₁∈v, ₂∈v))
- ∑_abi(L_a₂bi t_a₁bi C(₁∈o, ₂∈o))
julia> enable_external_index_translation()

Ref-$\nu$ contributions

$D_{pq}^{\text{ref-}\nu} = L_0 \sum_{\nu}{\langle\text{HF}| e^{-T} E_{pq} e^T |\nu\rangle R_\nu} \\= L_0 \sum_{ai}{\langle\text{HF}| e^{-T} E_{pq} e^T |_i^a\rangle R_{ai}} + L_0 \sum_{ai \geq bj}{ \langle\text{HF}| e^{-T} E_{pq} e^T |_{ij}^{ab}\rangle R_{aibj} }$

Singles:

julia> ∑(bch_Epq * E(3, 4) * occupied(4) * virtual(3) * real_tensor("R", 3, 4), 3:4);
julia> singles_ket = simplify(act_on_ket(ans, 2));
julia> simplify_heavy(act_on_bra(singles_ket))2 R_ai

For doubles we scale the $R_{aibj}$ by 2 on the diagonal

$\sum_{ai \geq bj}{ \langle\text{HF}| e^{-T} E_{pq} e^T |_{ij}^{ab}\rangle R_{aibj} } = \frac12 \sum_{aibj}{ (1 + \delta_{ai,bj}) \langle\text{HF}| e^{-T} E_{pq} e^T |_{ij}^{ab}\rangle R_{aibj} } = \frac12 \sum_{aibj}{ \langle\text{HF}| e^{-T} E_{pq} e^T |_{ij}^{ab}\rangle \tilde R_{aibj} }$

julia> bch_Epq * E(3, 4) * E(5, 6) * occupied(4, 6) * virtual(3, 5);
julia> simplify(1//2 * ∑(ans * psym_tensor("R", 3:6...), 3:6));
julia> doubles_ket = simplify(act_on_ket(ans, 2));
julia> simplify_heavy(act_on_bra(doubles_ket))0

$\mu$-$\nu$ contributions

$D_{pq}^{\mu\text{-}\nu} = \sum_{\mu\nu}{ L_\mu \langle\mu| e^{-T} E_{pq} e^T |\nu\rangle R_\nu } \\= \sum_{aick}{ L_{ai} \tilde{\langle_i^a|} e^{-T} E_{pq} e^T |_k^c\rangle R_{ck} } \\+ \sum_{aibjck}{ L_{aibj} \tilde{\langle_{ij}^{ab}|} e^{-T} E_{pq} e^T |_k^c\rangle R_{ck} } \\+ \sum_{aickdl}{ L_{ai} \tilde{\langle_i^a|} e^{-T} E_{pq} e^T |_{kl}^{cd}\rangle R_{ckdl} } \\+ \sum_{aibjckdl}{ L_{aibj} \tilde{\langle_{ij}^{ab}|} e^{-T} E_{pq} e^T |_{kl}^{cd}\rangle R_{ckdl} }$

Singles-Singles:

julia> project_biorthogonal(singles_ket, E(3, 4));
julia> simplify(∑(ans * real_tensor("L", 3, 4), 3:4));
julia> D_ss = ans∑_i(L_ai R_bi)
- ∑_a(L_aj R_ai)
+ 2 ∑_ak(δ_ij L_ak R_ak)

Doubles-Singles:

julia> project_biorthogonal(singles_ket, E(3, 4) * E(5, 6));
julia> simplify_heavy(∑(ans * psym_tensor("L", 3:6...), 3:6));
julia> look_for_tensor_replacements(ans, make_exchange_transformer("t", "u"));
julia> D_ds = ans∑_bj(R_bj L_aibj)
- ∑_jbck(R_aj L_bjck t_bick)
- ∑_bjck(R_bi L_bjck t_ajck)
+ ∑_bjck(R_bj L_bjck u_aick)

Singles-Doubles:

julia> project_biorthogonal(doubles_ket, E(3, 4));
julia> simplify(∑(ans * real_tensor("L", 3, 4), 3:4));
julia> look_for_tensor_replacements(ans, make_exchange_transformer("R", "Ru"));
julia> D_sd = ans∑_bj(L_bj Ru_aibj)

Doubles-Doubles:

julia> project_biorthogonal(doubles_ket, E(3, 4) * E(5, 6));
julia> simplify_heavy(∑(ans * psym_tensor("L", 3:6...), 3:6));
julia> D_dd = ans∑_icj(L_aicj R_bicj)
- ∑_abk(L_ajbk R_aibk)
+ ∑_akbl(δ_ij L_akbl R_akbl)

Combining and separating into blocks:

julia> D = D_ss + D_sd + D_ds + D_dd∑_i(L_ai R_bi)
- ∑_a(L_aj R_ai)
+ ∑_bj(L_bj Ru_aibj)
+ ∑_bj(R_bj L_aibj)
+ ∑_icj(L_aicj R_bicj)
- ∑_abk(L_ajbk R_aibk)
- ∑_jbck(R_aj L_bjck t_bick)
- ∑_bjck(R_bi L_bjck t_ajck)
+ ∑_bjck(R_bj L_bjck u_aick)
+ 2 ∑_ak(δ_ij L_ak R_ak)
+ ∑_akbl(δ_ij L_akbl R_akbl)
julia> D_ij = D * occupied(1, 2)- ∑_a(L_aj R_ai)
- ∑_abk(L_ajbk R_aibk)
+ 2 ∑_ak(δ_ij L_ak R_ak)
+ ∑_akbl(δ_ij L_akbl R_akbl)
julia> D_ia = D * occupied(1) * virtual(2)∑_bj(L_bj Ru_aibj)
- ∑_jbck(R_aj L_bjck t_bick)
- ∑_bjck(R_bi L_bjck t_ajck)
+ ∑_bjck(R_bj L_bjck u_aick)
julia> D_ai = D * occupied(2) * virtual(1)∑_bj(R_bj L_aibj)
julia> D_ab = D * virtual(1, 2)∑_i(L_ai R_bi)
+ ∑_icj(L_aicj R_bicj)