GitHub

Defining new index spaces

Defining your own index spaces can be a powerful tool when working with methods where you are dealing with other types of indices than the default provided spaces GeneralOrbital, OccupiedOrbital and VirtualOrbital. A few examples could be

  • Auxiliary index for integral decompositions (RI, Cholesky)
  • Multiple photon modes for QED methods
  • Subdividing GeneralOrbital to have different classes of fermions, for example:
    • Dividing into active and inactive orbitals for active space methods
    • Dividing into orbitals for different kinds of physical fermions such as protons and positrons (NEO-CCSD)

RI/Cholesky index

The first example is adding a simple index space for use as an auxiliary basis for integral decompositions. We define a new index space using the new_space function

julia> using SpinAdaptedSecondQuantization
julia> AuxiliaryIndex = new_space(:AuxiliaryIndex, "aux", "JKLMN")5

The variable AuxiliaryIndex now contains the id for the new index space which are ordered indices starting at 1. This can also be obtained from the global dictionary index_ids

julia> SASQ.index_ids[:AuxiliaryIndex]5

The second parameter "aux" is the "short name" for the space, used when printing terms without index translation. The last parameter "JKLMN" are the names to give indices in the given space when translated.

Note

Running new_space for the same key multiple times in the same session will not do anything, meaning it is completely safe to keep index space definitions in the top of a script that gets run multiple times.

We can now use our index space

julia> x = real_tensor("L", 1, 2, 3) * constrain(1 => AuxiliaryIndex) * electron(2, 3)L_Jpq
julia> disable_external_index_translation()
julia> x # Here the short name "aux" gets usedL_₁₂₃ C(₁∈aux, ₂∈g, ₃∈g)
julia> enable_external_index_translation()

It can be nice to define a convenience function to constrain many indices to the new space

julia> auxiliary(indices...) = constrain(J => AuxiliaryIndex for J in indices)auxiliary (generic function with 1 method)
julia> real_tensor("S", 1, 2) * auxiliary(1, 2)S_JK

We can now, for example, define the two electron operator in terms of auxiliary indices

julia> g = 1//2 * ∑(real_tensor("L", 5, 1, 2) * real_tensor("L", 5, 3, 4) *
           e(1, 2, 3, 4) * electron(1, 2, 3, 4) * auxiliary(5), 1:5)1/2 ∑_pqrsJ(L_Jpq L_Jrs E_pq E_rs)
- 1/2 ∑_pqrsJ(δ_qr L_Jpq L_Jrs E_ps)
julia> simplify_heavy(hf_expectation_value(g))2 ∑_Jij(L_Jii L_Jjj)
- ∑_Jij(L_Jij L_Jji)

Multilevel indices

Say we want to divide our electronic indices into an active and an inactive space. Then we can define the following spaces

julia> ActiveOrbital = new_space(:ActiveOrbital, "a", "pqrstuv")6
julia> InactiveOrbital = new_space(:InactiveOrbital, "i", "pqrstuv")7

Note that the same names have been given to these indices as the general electronic indices, so they will be printed the same when translated. To be able to distinguish them we can add a color to be printed for a certain index space (see set_color for available colors)

julia> set_color(ActiveOrbital, :cyan)
julia> set_color(InactiveOrbital, :red)

now if we use the indices we see that they are no longer identical

julia> ∑(real_tensor("g", 1, 2, 3, 4) *
           electron(1, 2) * constrain(3 => ActiveOrbital, 4 => InactiveOrbital), 1:4)∑_pqpp(g_pqpp)

However, we want these two new spaces to act in such a way that

ActiveOrbital ⊆ GeneralOrbital
InctiveOrbital ⊆ GeneralOrbital
ActiveOrbital ∪ InctiveOrbital == GeneralOrbital
ActiveOrbital ∩ InctiveOrbital == {}

To setup these relations we call the function add_space_sum

julia> add_space_sum(ActiveOrbital, InactiveOrbital, GeneralOrbital)
julia> ActiveOrbital ⊆ GeneralOrbitaltrue
julia> InactiveOrbital ⊆ GeneralOrbitaltrue

This for example allows the following simplification to occur

julia> ∑(real_tensor("h", 1, 2) * E(1, 2) *
           electron(1) * constrain(2 => ActiveOrbital), 1:2)∑_pp(h_pp E_pp)
julia> ans + ∑(real_tensor("h", 1, 2) * E(1, 2) *
           electron(1) * constrain(2 => InactiveOrbital), 1:2)∑_pp(h_pp E_pp)
+ ∑_pp(h_pp E_pp)
julia> simplify(ans)∑_pq(h_pq E_pq)

Now we also want to subdivide the VirtualOrbital and OccupiedOrbital index spaces into active and inactive subspaces, so we define these spaces

julia> begin
           ActiveVirtualOrbital = new_space(:ActiveVirtualOrbital, "av", "abcdefg");
           ActiveOccupiedOrbital = new_space(:ActiveOccupiedOrbital, "ao", "ijklmno");
           InactiveVirtualOrbital = new_space(:InactiveVirtualOrbital, "iv", "abcdefg");
           InactiveOccupiedOrbital = new_space(:InactiveOccupiedOrbital, "io", "ijklmno");
       
           set_color(ActiveVirtualOrbital, :cyan)
           set_color(ActiveOccupiedOrbital, :cyan)
       
           set_color(InactiveVirtualOrbital, :red)
           set_color(InactiveOccupiedOrbital, :red)
       end

Then we call add_space_sum to signal how each space can split into two smaller spaces

julia> add_space_sum(ActiveOrbital, InactiveOrbital, GeneralOrbital)
julia> add_space_sum(ActiveVirtualOrbital, InactiveVirtualOrbital, VirtualOrbital)
julia> add_space_sum(ActiveOccupiedOrbital, InactiveOccupiedOrbital, OccupiedOrbital)
julia> add_space_sum(ActiveOccupiedOrbital, ActiveVirtualOrbital, ActiveOrbital)
julia> add_space_sum(InactiveOccupiedOrbital, InactiveVirtualOrbital, InactiveOrbital)

In addition to this we now have a new relation we want to specify, namely the following

OccupiedOrbital ∩ ActiveOrbital == ActiveOccupiedOrbital
OccupiedOrbital ∩ InactiveOrbital == InactiveOccupiedOrbital
...

To achieve this we call the add_space_intersection function

julia> add_space_intersection(ActiveOrbital, OccupiedOrbital, ActiveOccupiedOrbital)
julia> add_space_intersection(ActiveOrbital, VirtualOrbital, ActiveVirtualOrbital)
julia> add_space_intersection(InactiveOrbital, OccupiedOrbital, InactiveOccupiedOrbital)
julia> add_space_intersection(InactiveOrbital, VirtualOrbital, InactiveVirtualOrbital)

This is what allows the following simplification to happen

julia> t = real_tensor("h", 1, 2) * constrain(1 => ActiveOrbital, 2 => InactiveOrbital)h_pp
julia> o = E(1, 2) * occupied(1, 2)E_ij
julia> t * oh_ii E_ii

where we see that the two indices got constrained to the smaller ActiveOccupiedOrbital and InactiveOccupiedOrbital spaces respectively.

We can make some convenience constructor functions

julia> active(indices...) = constrain(p => ActiveOrbital for p in indices)active (generic function with 1 method)
julia> inactive(indices...) = constrain(p => InactiveOrbital for p in indices)inactive (generic function with 1 method)
julia> aocc(indices...) = constrain(p => ActiveOccupiedOrbital for p in indices)aocc (generic function with 1 method)
julia> iocc(indices...) = constrain(p => InactiveOccupiedOrbital for p in indices)iocc (generic function with 1 method)
julia> avir(indices...) = constrain(p => ActiveVirtualOrbital for p in indices)avir (generic function with 1 method)
julia> ivir(indices...) = constrain(p => InactiveVirtualOrbital for p in indices)ivir (generic function with 1 method)

HF interaction energy

We can now use this to, for example, get an expression for the interaction energy between the active and inactive spaces at the HF level.

We define the full hamiltonian as

julia> h = ∑(rsym_tensor("h", 1, 2) * E(1, 2) * electron(1, 2), 1:2)∑_pq(h_pq E_pq)
julia> g = 1//2 * simplify(
           ∑(rsym_tensor("g", 1:4...) * e(1:4...) * electron(1:4...), 1:4)
       )- 1/2 ∑_pqr(g_prqr E_pq)
+ 1/2 ∑_pqrs(g_pqrs E_pq E_rs)
julia> H = h + g∑_pq(h_pq E_pq)
- 1/2 ∑_pqr(g_prqr E_pq)
+ 1/2 ∑_pqrs(g_pqrs E_pq E_rs)

which can give us the full HF energy

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

We can now also define Hamiltonians acting on only the active and inactive spaces as

julia> h_active = ∑(rsym_tensor("h", 1, 2) * E(1, 2) * active(1, 2), 1:2)∑_pq(h_pq E_pq)
julia> g_active = 1//2 * simplify(
           ∑(rsym_tensor("g", 1:4...) * e(1:4...) * active(1:4...), 1:4)
       )- 1/2 ∑_pqr(g_prqr E_pq)
+ 1/2 ∑_pqrs(g_pqrs E_pq E_rs)
julia> h_inactive = ∑(rsym_tensor("h", 1, 2) * E(1, 2) * inactive(1, 2), 1:2)∑_pq(h_pq E_pq)
julia> g_inactive = 1//2 * simplify(
           ∑(rsym_tensor("g", 1:4...) * e(1:4...) * inactive(1:4...), 1:4)
       )- 1/2 ∑_pqr(g_prqr E_pq)
+ 1/2 ∑_pqrs(g_pqrs E_pq E_rs)
julia> H_active = h_active + g_active∑_pq(h_pq E_pq)
- 1/2 ∑_pqr(g_prqr E_pq)
+ 1/2 ∑_pqrs(g_pqrs E_pq E_rs)
julia> H_inactive = h_inactive + g_inactive∑_pq(h_pq E_pq)
- 1/2 ∑_pqr(g_prqr E_pq)
+ 1/2 ∑_pqrs(g_pqrs E_pq E_rs)

We can now define the HF energies of the subspaces

julia> E_active = simplify_heavy(hf_expectation_value(H_active))2 ∑_i(h_ii)
+ 2 ∑_ij(g_iijj)
- ∑_ij(g_ijij)
julia> E_inactive = simplify_heavy(hf_expectation_value(H_inactive))2 ∑_i(h_ii)
+ 2 ∑_ij(g_iijj)
- ∑_ij(g_ijij)

Then finally we can define an interaction energy by subtracting these from the full energy

julia> E_int = simplify_heavy(E_HF - E_active - E_inactive,
           (OccupiedOrbital => (ActiveOccupiedOrbital, InactiveOccupiedOrbital),))4 ∑_ii(g_iiii)
- 2 ∑_ii(g_iiii)

Note here the use of the extra argument to the simplify_heavy function. This defaults to splitting general indices into occupied and virtual, but it was now more useful to split occupied orbitals into active and inactive parts.

Multilevel coupled cluster

If we now define the Hamiltonian in terms of the fock matrix as usual for coupled cluster derivations

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)

Then as an interesting example we can define the T2 cluster operator acting on active orbitals for one excitation, then for all orbitals for the other

julia> T2 = 1//2 * ∑(psym_tensor("t", 1:4...) * E(1, 2) * E(3, 4) *
           occupied(2, 4) * virtual(1, 3) * active(1, 2, 3, 4), 1:4) +
           ∑(psym_tensor("t", 1:4...) * E(1, 2) * E(3, 4) *
           occupied(2, 4) * virtual(1, 3) * active(1, 2) * inactive(3, 4), 1:4)1/2 ∑_aibj(t_aibj E_ai E_bj)
+ ∑_aiai(t_aiai E_ai E_ai)

We then define the similarity transformed Hamiltonian and its right projection as usual

julia> Hbar = simplify(bch(H, T2, 4));
julia> Hbar_ket = simplify(act_on_ket(Hbar, 2));

Energy:

julia> simplify_heavy(act_on_bra(Hbar_ket));
julia> look_for_tensor_replacements(ans,
           make_exchange_transformer("g", "L"));
julia> E0 = ans2 ∑_i(F_ii)
- ∑_ij(L_iijj)
+ ∑_iajb(L_iajb t_aibj)
+ 2 ∑_iaia(L_iaia t_aiai)

Singles:

julia> project_biorthogonal(Hbar_ket, E(1, 2));
julia> simplify_heavy(ans, [
           GeneralOrbital => (OccupiedOrbital, VirtualOrbital),
           OccupiedOrbital => (ActiveOccupiedOrbital, InactiveOccupiedOrbital),
           VirtualOrbital => (ActiveVirtualOrbital, InactiveVirtualOrbital),
           GeneralOrbital => (ActiveOrbital, InactiveOrbital),
           ActiveOrbital => (ActiveOccupiedOrbital, ActiveVirtualOrbital),
           InactiveOrbital => (InactiveOccupiedOrbital, InactiveVirtualOrbital),
       ]);
julia> look_for_tensor_replacements(ans,
           make_exchange_transformer("t", "u"));
julia> look_for_tensor_replacements(ans,
           make_exchange_transformer("g", "L"));
julia> omega_ai = ansF_ai
+ 2 ∑_ia(F_ia t_aiai)
+ 2 ∑_ia(F_ia t_aiai)
- ∑_ia(F_ia t_aiai)
- ∑_ia(F_ia t_aiai)
+ ∑_jb(F_jb u_aibj)
+ ∑_aia(L_aaia t_aiai)
+ ∑_aia(L_aaia t_aiai)
- ∑_iia(L_iiia t_aiai)
- ∑_iia(L_iiia t_aiai)
+ ∑_ajb(g_aajb u_aibj)
- ∑_ijb(g_iijb u_aibj)

Since the different terms here have different constraints for the external indices, it can be useful to separate this into active/inactive blocks:

julia> omega_ai_AA = omega_ai * active(1, 2)F_ai
+ 2 ∑_ia(F_ia t_aiai)
+ ∑_jb(F_jb u_aibj)
+ ∑_bia(L_abia t_biai)
- ∑_jia(L_jiia t_ajai)
+ ∑_bjc(g_abjc u_bicj)
- ∑_jkb(g_jikb u_ajbk)
julia> omega_ai_AI = omega_ai * active(1) * inactive(2)F_ai
- ∑_ia(F_ia t_aiai)
+ ∑_aib(L_aaib t_aibi)
- ∑_ija(L_iija t_aiaj)
- ∑_ijb(g_iijb u_aibj)
julia> omega_ai_IA = omega_ai * active(2) * inactive(1)F_ai
- ∑_ia(F_ia t_aiai)
+ ∑_aib(L_aaib t_aibi)
- ∑_ija(L_iija t_aiaj)
+ ∑_ajb(g_aajb u_aibj)
julia> omega_ai_II = omega_ai * inactive(1, 2)F_ai
+ 2 ∑_ia(F_ia t_aiai)
+ ∑_bia(L_abia t_biai)
- ∑_jia(L_jiia t_ajai)