Implementing your own types

In specific cases you might want to implement your own operator and tensor types. This is a simple matter of implementing the correct set of functions on the new types. Here is a short tutorial for how to implement your own type.

Implementing your own tensor type

As an example we will implement a tensor where every single index permutation is equivalent. This might not be very useful, but it is simple to implement.

We start by defining the new tensor type with convenient member variables

julia> using SpinAdaptedSecondQuantization
julia> struct FullySymmetricTensor <: SASQ.Tensor # Must be subtype to Tensor
           symbol::String
           indices::Vector{Int}
       end

Next we must implement the following set of functions:

• get_symbol(t), should return the name of the tensor
• get_indices(t), should return ordered list of indices
• exchange_indices(t, mapping), should return tensor after indices are changed according to the given mapping
• reorder_indices(t, permutation) should return tensor after reordering indices according to the given permutation

Optionally we can implement the functions, but they have reasonable default implementations

• Base.show(io::IO, t), if one wants to print the tensor in a custom way.
• get_permutations(t), needed for print_eT_function_generator to take advantage of available permutations.

Here we will only implement the required functions.

julia> SASQ.get_symbol(t::FullySymmetricTensor) = t.symbol
julia> SASQ.get_indices(t::FullySymmetricTensor) = t.indices
julia> function SASQ.exchange_indices(t::FullySymmetricTensor, mapping)
           # Use internal function exchange_index on each index
           new_inds = [SASQ.exchange_index(p, mapping) for p in t.indices]
           sort!(new_inds) # Sort indices as any permutation is equivalent
           FullySymmetricTensor(t.symbol, new_inds)
       end
julia> function SASQ.reorder_indices(t::FullySymmetricTensor, permutation)
           new_inds = t.indices[permutation]
           sort!(ind)
           FullySymmetricTensor(t.symbol, new_inds)
       end

It is also useful to make a constructor that produces an Expression containing that tensor:

julia> function fsym_tensor(symbol, indices...)
           SASQ.Expression(FullySymmetricTensor(symbol, sort!(collect(indices))))
       endfsym_tensor (generic function with 1 method)

Now we can use our new tensor type

julia> fsym_tensor("x", 2, 3, 2, 7, 5) # Sorts indicesx_₂₂₃₅₇
julia> ∑(fsym_tensor("x", 1, 2, 3) * E(1, 2) * E(1, 3) * electron(1, 2, 3), 1:3)∑_pqr(x_pqr E_pq E_pr)
julia> act_on_ket(ans)∑_aij(x_aij E_aj E_ai)
+ 4 ∑_ijk(δ_ijk x_ijk)
+ ∑_abi(δ_ab x_abi E_ai)
julia> simplify(ans)4 ∑_i(x_iii)
+ ∑_ai(x_aai E_ai)
+ ∑_aij(x_aij E_ai E_aj)

Implementing your own operator type

Here as an example we will implement the "hole" excitation operator defined as

$E^V_{pq} = a_{p\alpha} a^\dagger_{q\alpha} + a_{p\beta} a^\dagger_{q\beta} = 2 \delta_{pq} - E_{qp}$

julia> using SpinAdaptedSecondQuantization
julia> Ev(p, q) = 2 * δ(p, q) - E(q, p)Ev (generic function with 1 method)
julia> Ev(1, 2) * electron(1, 2)- E_qp
+ 2 δ_pq

We start by defining our new type

julia> struct HoleExcitationOperator <: SASQ.Operator # Must be subtype of Operator
           p::Int
           q::Int
       end

The set of functions to implement is similar to those needed for implementing a tensor above. Required functions are:

• Base.show for printing.
• exchange_indices
• get_all_indices
• act_on_ket
• Base.:(==)
• Base.isless
• Base.adjoint

The Base.show function takes in a tuple containing the information to correctly print the names of indices.

julia> function Base.show(io::IO,
       (
           o, constraints, translation
       )::Tuple{HoleExcitationOperator,SASQ.Constraints,SASQ.IndexTranslation})
           print(io, "Ev_")
           SASQ.print_mo_index(io, constraints, translation, o.p, o.q)
       end

julia> function SASQ.exchange_indices(o::HoleExcitationOperator, mapping)
           HoleExcitationOperator(
               SASQ.exchange_index(o.p, mapping),
               SASQ.exchange_index(o.q, mapping)
           )
       end
julia> function SASQ.get_all_indices(o::HoleExcitationOperator)
           (o.p, o.q)
       end
julia> function Base.:(==)(a::HoleExcitationOperator, b::HoleExcitationOperator)
           (a.p, a.q) == (b.p, b.q)
       end
julia> function Base.isless(a::HoleExcitationOperator, b::HoleExcitationOperator)
           (a.p, a.q) < (b.p, b.q)
       end
julia> function Base.adjoint(o::HoleExcitationOperator)
           HoleExcitationOperator(o.q, o.p)
       end
julia> function SASQ.act_on_ket(o::HoleExcitationOperator)
           p, q = o.p, o.q
       
           2 * δ(p, q) * virtual(p, q) -
           E(q, p) * occupied(p) * virtual(q)
       end

Then there are functions that are implemented for each pair of operator types. These only needs to be implemented for the pairs you plan on using in the same expressions. Here we will only be implementing these for the normal SingletExcitationOperator.

• Base.isless To define the lexicographic ordering of different operator types
• reductive_commutator To define commutation relations

The first function is simple as we only need to decide on an order of the operator. Say we want to have E_pq operators come before Ev_pq operators, then we define the function

julia> function Base.isless(
           ::Type{HoleExcitationOperator},
           ::Type{SASQ.SingletExcitationOperator})
           false
       end

next is the reductive_commutator which gives the commutation relation between different operators, giving both the expression for the commutator and whether it is a commutator or anticommutator. Here both the commutators we want to implement are commutators, so we return a positive 1 to indicate no sign change when commuting.

julia> function SASQ.reductive_commutator(
           a::HoleExcitationOperator,
           b::HoleExcitationOperator)
       
           p, q = a.p, a.q
           r, s = b.p, b.q
       
           (1, δ(p, s) * E(q, r) - δ(q, r) * E(s, p))
       end
julia> function SASQ.reductive_commutator(
           a::SASQ.SingletExcitationOperator,
           b::HoleExcitationOperator)
       
           p, q = a.p, a.q
           r, s = b.p, b.q
       
           (1, δ(p, r) * E(s, q) - δ(q, s) * E(p, r))
       end

Then finally we make a constructor that produces an Expression

julia> Ev(p, q) = SASQ.Expression(HoleExcitationOperator(p, q))Ev (generic function with 1 method)

Then we can start using the new operator

julia> Ev(1, 2) * electron(1, 2)Ev_pq
julia> commutator(Ev(1, 2), E(3, 4)) * electron(1:4...)- δ_pr E_qs
+ δ_qs E_rp
julia> ∑(fsym_tensor("x", 1, 2, 3, 4) * E(1, 2) * Ev(3, 4) * electron(1:4...), 1:4)∑_pqrs(x_pqrs E_pq Ev_rs)
julia> simplify_heavy(act_on_ket(ans))2 ∑_ai(x_aaii)
+ ∑_aib(x_aibb E_ai)
- ∑_aij(x_aijj E_ai)
- ∑_aibj(x_aibj E_ai E_bj)