Public API

E(p, q)

Constructs an expression containing a single excitation operator.

Example

julia> using SpinAdaptedSecondQuantization

julia> E(1, 2) * electron(1, 2)
E_pq

julia> E(1, 2) * E(3, 4) * occupied(2, 4) * virtual(1, 3)
E_ai E_bj

act_eT_on_bra(braop::Expression, T::Expression; [max_n=Inf], [max_ops=Inf])

computes

$\langle\text{HF}| \text{braop}\ e^T$

This can be useful for evaluating certain coupled cluster expressions from left to right. Since the cluster operator T is purely exciting acting $e^T$ to the left will terminate after only a few terms depending on the left state being projected on.

act_on_bra(x::Expression, [max_ops=Inf])

Projects x on a HF bra by projecting the adjoint on a ket and ajointing the result.

act_on_bra(x, max_ops=Inf) = act_on_ket(x', max_ops)'

See act_on_ket

act_on_ket(::Operator)

All typed extending Operator must overload this function. It should return the expression you end up with after acting the operator on a normal RHF ket.

Example:

$E_{p q} |\text{HF}⟩ = \left(2 δ_{p q} C(p∈O,q∈O) + E_{p q} C(p∈V,q∈O)\right) |\text{HF}⟩$

Which we can also do in code:

julia> using SpinAdaptedSecondQuantization

julia> E(1, 2) * electron(1, 2)
E_pq

julia> ket = act_on_ket(ans)
E_ai
+ 2 δ_ij

julia> disable_external_index_translation()

julia> ket
E_₁₂ C(₁∈v, ₂∈o)
+ 2 δ_₁₂ C(₁∈o, ₂∈o)

julia> enable_external_index_translation()

act_on_ket(ex::Expression, [max_ops=Inf])

Project ex onto a $|\text{HF}\rangle$ ket to the right. This is done termwise and for each term the folliwing recursive formula is used to project an arbitrary string of operators on a ket.

The base case is that there is an implementation of act_on_ket(o) for each operator type in the string such that we have

$A_i |\text{HF}\rangle = \tilde A_i |\text{HF}\rangle$

where $\tilde A_i$ is constrained and simplified from the projection. The recursive step is then

$A_1 A_2 ... A_n |\text{HF}\rangle = A_1 A_2 ... \tilde A_n |\text{HF}\rangle = [A_1 A_2 ..., \tilde A_n]_\Gamma |\text{HF}\rangle + \Gamma \tilde A_n (A_1 A_2 ... A_{n-1}) |\text{HF}\rangle$

where the $[A, B]_\Gamma$ notation is the reductive_commutator which chooses whether to compute a commutator or anticommutator depending on the input to reduce the number of operators in the output.

The optional parameter max_ops specifies the maximum number of operators to keep in any terms. If only a few excitations are needed, specifying this can greatly speed up the projection.

Example

julia> using SpinAdaptedSecondQuantization

julia> E(1, 2) * E(3, 4) * electron(1:4...)
E_pq E_rs

julia> ket = simplify(act_on_ket(ans))
E_ai E_bj
+ 2 δ_ij E_ak
- δ_ik E_aj
+ δ_bc E_ai
+ 2 δ_jk E_ai
+ 4 δ_ij δ_kl
+ 2 δ_ij δ_ab

julia> disable_external_index_translation()

julia> ket
E_₁₂ E_₃₄ C(₁∈v, ₂∈o, ₃∈v, ₄∈o)
+ 2 δ_₁₂ E_₃₄ C(₁∈o, ₂∈o, ₃∈v, ₄∈o)
- δ_₁₄ E_₃₂ C(₁∈o, ₂∈o, ₃∈v, ₄∈o)
+ δ_₂₃ E_₁₄ C(₁∈v, ₂∈v, ₃∈v, ₄∈o)
+ 2 δ_₃₄ E_₁₂ C(₁∈v, ₂∈o, ₃∈o, ₄∈o)
+ 4 δ_₁₂ δ_₃₄ C(₁∈o, ₂∈o, ₃∈o, ₄∈o)
+ 2 δ_₁₄ δ_₂₃ C(₁∈o, ₂∈v, ₃∈v, ₄∈o)

julia> enable_external_index_translation()

add_space_intersection(a::Int, b::Int, c::Int)

Specifies that the intersection of spaces a and b is space c. See Defining new index spaces for example.

add_space_sum(a::Int, b::Int, c::Int)

Specifies that the union of spaces a and b is space c. See Defining new index spaces for example.

anticommutator(a::Expression, b::Expression)

Computes anticommutator between a and b. This is computed termwise using the reductive_commutator function, then adding/subtracting the swapped result if the reductive commutator for a given pair of terms turned out to be a commutator.

bch(A, Bs::AbstractArray, n)

Equivalent to bch(A, sum(Bs), n) if all Bs commute with eachother, but evalutating in a smart way in order to avoid computing equivalent terms many times, for example instead of computing both $rac12 [[A, B_1], B_2]$ and $rac12 [[A, B_2], B_1]$, it will only compute the first one and multiply by 2. This can be much faster if one has many commuting terms in Bs.

bch(A, B, n)

Compute Baker-Campbell-Haussdorff expansion

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

up to nth order. Terminates if commutator becomes zero.

boson()

Constructs a boson annihilation operator.

Example

julia> using SpinAdaptedSecondQuantization

julia> boson()
b⁻

bosondag()

Constructs a boson creation operator.

Example

julia> using SpinAdaptedSecondQuantization

julia> bosondag()
b†

julia> boson()'
b†

commutator(A, B, C, ...)

Computes the nested commutator $[...[[A, B], C], ...]$

commutator(A::Expression, B::Expression, n::Integer)

Compute n nested commutators $[...[[[A], B], B]..., B]$.

commutator(a::Expression, b::Expression)

Computes commutator between a and b. This is computed termwise using the reductive_commutator function, then adding/subtracting the swapped result if the reductive commutator for a given pair of terms turned out to be an anti-commutator.

convert_to_elementary_operators(ex::Expression)

Converts all operators to their representation using "elementary" annihilation and creation operators if implemented for all operator types in ex.

delta(indices...)

Constructs an expression consisting of a single Kronecker delta of the given indices.

desymmetrize(ex::Expression, mappings)

Function to look for terms which are equal under the index permutations given by the mappings. See symmetrize and make_permutation_mappings.

Returns three expressions:

1. The first expression contains the terms where redundant copies were found which are equal under the given symmetry. These needs to be symmetrized to obtain the original expression, usually after numerically evaluating the non-redundant terms.
2. The second expression contains the "self-symmetric" terms, which on their own are invariant under the given symmetries and do not need to be symmetrized.
3. The last expression contains any left over "non-symmetric" terms. This is expected to be zero when desymmetrizing an expression with a known symmetry such as CCSD omega equations.

Examples

Setting up an example expression:

julia> using SpinAdaptedSecondQuantization

julia> ∑(real_tensor("F", 1, 5) * psym_tensor("t", 3, 4, 5, 2) * occupied(2, 4) * virtual(1, 3, 5), [5])
∑_c(F_ac t_bjci)

julia> symmetrize(ans, make_permutation_mappings([(1, 2), (3, 4)]))
∑_c(F_ac t_bjci)
+ ∑_c(F_bc t_aicj)

julia> ans + ∑(psym_tensor("g", 1, 5, 3, 6) * psym_tensor("t", 5, 2, 6, 4) * occupied(2, 4) * virtual(1, 3, 5, 6), 5:6)
∑_c(F_ac t_bjci)
+ ∑_c(F_bc t_aicj)
+ ∑_cd(g_acbd t_cidj)

julia> ans + real_tensor("x", 1, 2, 3, 4) * occupied(2, 4) * virtual(1, 3)
x_aibj
+ ∑_c(F_ac t_bjci)
+ ∑_c(F_bc t_aicj)
+ ∑_cd(g_acbd t_cidj)

Desymmetrizing:

julia> r, ss, ns = desymmetrize(ans * 1//1, make_permutation_mappings([(1, 2), (3, 4)]));

julia> r # Only one of the two redundant terms survives
∑_c(F_ac t_bjci)

julia> ss # The self-symmetric term
∑_cd(g_acbd t_cidj)

julia> ns # The non-symmetric term
x_aibj

disable_color(s)

Disables the coloring for index space s.

disable_external_index_translation()

Disables external indices to be printed with names according to their constraints. This can be useful as a debugging tool when dealing with terms with different constraints on external indices in the same expression.

See enable_external_index_translation for example.

Enabled by default. Re-enable using enable_external_index_translation

disable_internal_index_translation()

Disables internal indices to be printed with names according to their constraints. This can be useful as a debugging tool if summation indices are getting reordered/changed in weird ways.

See enable_internal_index_translation for example.

Enabled by default. Re-enable using enable_internal_index_translation

e(p, q, r, s) = E(p, q) * E(r, s) - δ(r, q) * E(p, s)

Alias for the two electron singlet excitation operator.

Example:

julia> using SpinAdaptedSecondQuantization

julia> e(1, 2, 3, 4) * electron(1, 2, 3, 4)
E_pq E_rs
- δ_qr E_ps

electron(indices...)

Constraining the given indices to the GeneralOrbital index space.

Example

julia> using SpinAdaptedSecondQuantization

julia> electron(1, 2)
C(p∈g, q∈g)

julia> E(1, 2) * ans
E_pq

enable_external_index_translation()

Enables external indices to be printed with names according to their constraints.

julia> using SpinAdaptedSecondQuantization

julia> ex = E(1, 2) * occupied(1, 2) + E(1, 2) * virtual(1, 2)
E_ab
+ E_ij
julia> disable_external_index_translation()

julia> ex # Now it is clearer that they are the same indices
E_₁₂ C(₁∈v, ₂∈v)
+ E_₁₂ C(₁∈o, ₂∈o)
julia> enable_external_index_translation()

julia> ex # back to default behavour
E_ab
+ E_ij

Enabled by default. Disabled using disable_external_index_translation.

enable_internal_index_translation()

Enables internal indices to be printed with names according to their constraints.

Example

julia> using SpinAdaptedSecondQuantization

julia> ex = ∑(real_tensor("h", 1, 2) * E(1, 2) * electron(1, 2), 1:2)
∑_pq(h_pq E_pq)

julia> disable_internal_index_translation()

julia> ex
∑_₁₂(h_₁₂ E_₁₂ C(₁∈g, ₂∈g))

julia> ex *= E(1, 2) * occupied(2) * virtual(1) # Indices shifted up
∑_₃₄(h_₃₄ E_₃₄ E_ai C(₃∈g, ₄∈g))

julia> enable_internal_index_translation()

julia> ex # Shift in indices is invisible when translated
∑_pq(h_pq E_pq E_ai)

Enabled by default.

fermion(p, σ)

Constructs a fermion annihilation operator for orbital p with spin σ.

Example

julia> using SpinAdaptedSecondQuantization

julia> fermion(1, α) * electron(1)
a⁻_pα

fermiondag(p, σ)

Constructs a fermion creation operator for orbital p with spin σ.

Example

julia> using SpinAdaptedSecondQuantization

julia> fermiondag(1, α) * electron(1)
a†_pα

julia> fermion(1, α)' * electron(1)
a†_pα

hf_expectation_value(ex::Expression) = act_on_ket(ex, 0)

Alias for calling act_on_ket keeping no operators, which gives the HF expectation value.

look_for_tensor_replacements(ex::Expression, transformer)

Looks for pairs of terms in ex that are equal up to a tensor transformation given by transformer. The most common use case is to recognize patterns such as $L_{pqrs} = 2 g_{pqrs} - g_{psrq}$ where a convencience function to produce the relevant transformer provided (see make_exchange_transformer).

DocTestSetup = quote
    using SpinAdaptedSecondQuantization
    h = ∑(real_tensor("h", 1, 2) * E(1, 2) * electron(1, 2), 1:2)
    g = 1//2 * simplify(
        ∑(psym_tensor("g", 1:4...) * e(1:4...) * electron(1:4...), 1:4)
    )
    H = h + g + real_tensor("h_nuc")
end

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

julia> look_for_tensor_replacements(E_HF, make_exchange_transformer("g", "L"))
h_nuc
+ 2 ∑_i(h_ii)
+ ∑_ij(L_iijj)

make_exchange_transformer(name_from, name_to)

Convenience function to make tensor transformer for use in look_for_tensor_replacements.

make_permutation_mappings(index_groups)

Convenience function to make mapping for all permutations of groups of indices. See symmetrize and desymmetrize.

new_space(id, shortname::String, names::String)

Adds a new index space. See Defining new index spaces for example.

occupied(indices...)

Constraining the given indices to the OccupiedOrbital index space.

Example

julia> using SpinAdaptedSecondQuantization

julia> occupied(1, 2)
C(i∈o, j∈o)

julia> E(1, 2) * ans
E_ij

Print term to numpy's einsum routine

Used with https://github.com/alexancp/einsumpath-to-eT to generate Fortran subroutines for eT.

print_eT_function_generator(name, ex, symbol, indices, translation,
                [wf_type], [tensor_translation], [noinput_tensors], [outperms])

This function generates a script to be run with another script for making eT compatible fortran code.

print_julia_function(name, ex::Expression, [symbol="X"];
                [tensor_translation], [explicit_tensor_blocks])

Returns String with a runnable julia function to evaluate the given expression ex. name specifies the name of the function, and the optional parameter symbol specifies the name of the accumulator variable used inside the function. If not specified it will use the generic name X.

The optional keyword argument tensor_translation is a dictionary for replacing the name of certain tensors.

The optional keyword argument explicit_tensor_blocks is a vector of names of tensors you want subblocks to be taken as separate inputs to the function rather than using @view to get subblocks.

The returned function uses the TensorOperations.jl package for tensor contractions.

Example

julia> using SpinAdaptedSecondQuantization

julia> L = psym_tensor("L", 1, 3, 4, 5);

julia> R = psym_tensor("R", 2, 3, 4, 5);

julia> ex = ∑(L * R * occupied(3, 5) * virtual(1, 2, 4), 3:5)
∑_icj(L_aicj R_bicj)

julia> print_julia_function("density", ex) |> print
function density(no, nv, L, R)
    nao = no + nv
    o = 1:no
    v = no+1:nao
    
    X = zeros(nv, nv)
    L_vovo = @view L[v,o,v,o]
    R_vovo = @view R[v,o,v,o]
    @tensoropt (a=>10χ,b=>10χ,c=>10χ,i=>χ,j=>χ) begin
        X[a,b] += L_vovo[a,i,c,j]*R_vovo[b,i,c,j]
    end
    X
end

julia> print_julia_function("density", ex; explicit_tensor_blocks=["L", "R"]) |> print
function density(no, nv, L_vovo, R_vovo)
    nao = no + nv
    o = 1:no
    v = no+1:nao
    
    X = zeros(nv, nv)
    @tensoropt (a=>10χ,b=>10χ,c=>10χ,i=>χ,j=>χ) begin
        X[a,b] += L_vovo[a,i,c,j]*R_vovo[b,i,c,j]
    end
    X
end

project_biorthogonal(ex::Expression, template::Expression)

This has the effect of projecting the expression ex on the biorthogonal bra state of the template expression.

For a singlet double excitation, the following will be equivalent:

bra = (
        1//3 * E(1, 2) * E(3, 4) +
        1//6 * E(1, 4) * E(3, 2)
    )' * occupied(2, 4) * virtual(1, 3)

result = hf_expectation_value(bra * ex)

ket = E(1, 2) * E(3, 4) * occupied(2, 4) * virtual(1, 3)
result = project_biorthogonal(act_on_ket(ex), ket)
result = symmetrize(result, make_permutation_mappings([(1, 2), (3, 4)]))

set_color(s, color)

Sets the color of a index space s to the specified color. Allowed colors are given by the Base.text_colors dict. This includes integers in the range 0-255 which the corresponding colors can be found in this table on Wikipedia.

simplify(ex::Expression)

Applies various simplification strategies to the given expression. The most important are:

• Remove Kronecker deltas when indices are summed over.
• Perform various operations to try to lower the lexicographic ordering of the each term including:
  • Lowering indices showing up in a Kronecker delta to the lowest of the indices in said delta.
  • Lowering summation indices to lowest available indices.
  • Reorder summation indices to lower the lexicographic ordering of the order they show up within the term.
  • Bubble sort operators swapping neighbouring operators if:
    • They commute
    • The swap would lower the terms lexicographic ordering.
  • Look for chains of operators that are zero.
• Look for pairs of terms differing only in index constraints and scalar prefactors, and tries to combine them.

Examples:

Collapsing deltas:

julia> using SpinAdaptedSecondQuantization

julia> ∑(real_tensor("h", 1, 2) * electron(1, 2) * δ(1, 2), [2])
∑_q(δ_pq h_pq)

julia> simplify(ans)
h_pp

Combining terms differing only by index constraints and scalars:

julia> real_tensor("h", 1, 2) * (occupied(1, 2) + occupied(1) * virtual(2))
h_ia
+ h_ij

julia> simplify(ans)
h_ip

julia> ∑(real_tensor("h", 1, 2) * E(1, 2) * electron(1, 2), [1, 2])
∑_pq(h_pq E_pq)

julia> ans - ∑(real_tensor("h", 1, 2) * E(1, 2) * electron(1) * occupied(2), [1, 2])
∑_pq(h_pq E_pq)
- ∑_pi(h_pi E_pi)

julia> simplify(ans)
∑_pa(h_pa E_pa)

Sorting commuting operators:

julia> E(3, 4) * E(1, 2) * electron(1:4...)
E_rs E_pq

julia> simplify(ans) # not able to swap as commutator is not zero
E_rs E_pq

julia> ans * occupied(2, 4) * virtual(1, 3)
E_bj E_ai

julia> simplify(ans) # now able to swap
E_ai E_bj

julia> E(3, 4) * E(5, 6) * E(1, 2) * electron(5, 6) * occupied(2, 4) * virtual(1, 3)
E_bj E_pq E_ai

julia> simplify(ans) # Can not swap first and last because of non commuting between
E_bj E_pq E_ai

julia> E(5, 6) * E(3, 4) * E(1, 2) * electron(5, 6) * occupied(2, 4) * virtual(1, 3)
E_pq E_bj E_ai

julia> simplify(ans) # Can swap adjacent last two, but not first
E_pq E_ai E_bj

Illegal operator chains:

julia> fermion(1, α) * fermion(1, α) * electron(1)
a⁻_pα a⁻_pα

julia> simplify(ans)
0

julia> E(1, 2) * E(1, 3) * E(1, 4) * occupied(2, 3, 4) * virtual(1)
E_ai E_aj E_ak

julia> simplify(ans)
0

simplify_heavy(ex::Expression,
    [mapping=GeneralOrbital => (OccupiedOrbital, VirtualOrbital)])

Similar to simplify, but enabling a few extra simplification strategies that can be quite expensive, especially as the number of summation indices in a term grows large.

The most expensive step is to iterate over every permutation of the summation indices of each term in order to find the order that has the least lexiographic ordering, which naturally has a scaling of O(n!) where n is the larget number of summation indices in a single term. This is quite managable for "small" numbers of indices (less than ~8), but quickly becomes unfeasible.

The optional mapping argument defaults to splitting terms with general indices into a term with the index being occupied plus a term where the index is virtual, then performing simplification when all terms are split and recombining terms according to the same splitting if possible. This often allows for a few additional cancellations to occur than using only simplify If dealing with other index spaces, this splitting can be specified can be specified to be some other set of index spaces. See split_indices.

Examples:

HF energy with full real orbital symmetry:

julia> using SpinAdaptedSecondQuantization

julia> h = ∑(rsym_tensor("h", 1, 2) * E(1, 2) * electron(1, 2), 1:2)
∑_pq(h_pq E_pq)

julia> g = 1//2 * ∑(rsym_tensor("g", 1:4...) * e(1:4...) * electron(1:4...), 1:4)
1/2 ∑_pqrs(g_pqrs E_pq E_rs)
- 1/2 ∑_pqrs(δ_qr g_pqrs E_ps)

julia> H = simplify(h + g + real_tensor("h_nuc"));

julia> E_HF = simplify(hf_expectation_value(H)) # Not nice
h_nuc
+ 2 ∑_i(h_ii)
+ 2 ∑_ij(g_iijj)
- ∑_pi(g_pipi)
+ ∑_ia(g_iaia)

julia> simplify_heavy(E_HF) # Much nicer!
h_nuc
+ 2 ∑_i(h_ii)
+ 2 ∑_ij(g_iijj)
- ∑_ij(g_ijij)

Simplification by combining terms differing by index constraints

julia> ∑(real_tensor("h", 1, 2) * E(1, 2) * electron(1, 2), [1, 2])
∑_pq(h_pq E_pq)

julia> ans - ∑(real_tensor("h", 1, 2) * E(1, 2) * occupied(1, 2), [1, 2])
∑_pq(h_pq E_pq)
- ∑_ij(h_ij E_ij)

julia> ans - ∑(real_tensor("h", 1, 2) * E(1, 2) * virtual(1, 2), [1, 2])
∑_pq(h_pq E_pq)
- ∑_ab(h_ab E_ab)
- ∑_ij(h_ij E_ij)

julia> simplify(ans)
∑_pq(h_pq E_pq)
- ∑_ab(h_ab E_ab)
- ∑_ij(h_ij E_ij)

julia> simplify_heavy(ans)
∑_ai(h_ai E_ai)
+ ∑_ia(h_ia E_ia)

symmetrize(ex::Expression, mappings)

Function to expand the expression ex including all permutations of indices given by mappings. This is common to use to expand all permutations among neighbouring pairs of indices, such as to compute $P_{ijk}^{abc} \Omega_{ijk}^{abc}$. For this specific purpose we provide the convenience function make_permutation_mappings which produces all permutations between groups of indices.

DocTestSetup = :(using SpinAdaptedSecondQuantization)

julia> x = real_tensor("t", 1, 2, 3, 4, 5, 6) *
           occupied(2, 4, 6) * virtual(1, 3, 5)
t_aibjck

julia> symmetrize(x, make_permutation_mappings([(1, 2), (3, 4), (5, 6)]))
t_aibjck
+ t_aickbj
+ t_bjaick
+ t_bjckai
+ t_ckaibj
+ t_ckbjai

virtual(indices...)

Constraining the given indices to the VirtualOrbital index space.

Example

julia> using SpinAdaptedSecondQuantization

julia> virtual(1, 2)
C(a∈v, b∈v)

julia> E(1, 2) * ans
E_ab

δ(indices...)

Unicode alias for delta.

τ(p, q)

Constructs an expression containing a single triplet excitation operator.

Example

julia> using SpinAdaptedSecondQuantization

julia> τ(1, 2) * electron(1, 2)
T_pq

julia> convert_to_elementary_operators(ans)
a†_pα a⁻_qα
- a†_pβ a⁻_qβ