added support for triplet excitation operator

src/operator.jl

@@ -8,6 +8,7 @@ must extend.
 abstract type Operator end
 
 include("operators/singlet_excitation_operator.jl")
+include("operators/triplet_excitation_operator.jl")
 include("operators/fermion_operator.jl")
 include("operators/boson_operator.jl")
 

src/operators/commutation_relations.jl

@@ -54,3 +54,43 @@ end
 function reductive_commutator(::FermionOperator, ::BosonOperator)
     return (1, zero(Expression{Int64}))
 end
+
+function reductive_commutator(::TripletExcitationOperator, ::BosonOperator)
+    return (1, zero(Expression{Int64}))
+end
+
+function reductive_commutator(t::TripletExcitationOperator,
+    e::SingletExcitationOperator)
+
+    p = t.p
+    q = t.q
+
+    r = e.p
+    s = e.q
+
+    (1, δ(r, q) * τ(p, s) - δ(p, s) * τ(r, q))
+end
+
+function reductive_commutator(t1::TripletExcitationOperator,
+    t2::TripletExcitationOperator)
+
+    p = t1.p
+    q = t1.q
+
+    r = t2.p
+    s = t2.q
+
+    (1, δ(r, q) * E(p, s) - δ(p, s) * E(r, q))
+end
+
+function reductive_commutator(e::TripletExcitationOperator, a::FermionOperator)
+    p = e.p
+    q = e.q
+    r = a.p
+
+    (1, if a.dag
+        δ(q, r) * fermiondag(p, a.spin)
+    else
+        -δ(p, r) * fermion(q, a.spin)
+    end * (a.spin == α ? 1 : -1))
+end

src/operators/sorting.jl

@@ -1,5 +1,7 @@
 # Implement ordering of new operator types here:
 
+# current sorting is: T, E, a, b
+
 function Base.isless(::FermionOperator, ::SingletExcitationOperator)
     false
 end
@@ -11,3 +13,15 @@ end
 function Base.isless(::BosonOperator, ::FermionOperator)
     false
 end
+
+function Base.isless(::SingletExcitationOperator, ::TripletExcitationOperator)
+    false
+end
+
+function Base.isless(::FermionOperator, ::TripletExcitationOperator)
+    false
+end
+
+function Base.isless(::BosonOperator, ::TripletExcitationOperator)
+    false
+end

src/operators/triplet_excitation_operator.jl

@@ -0,0 +1,64 @@
+export τ
+
+"""
+    SingletExcitationOperator
+
+The T_pq type operator.
+"""
+struct TripletExcitationOperator <: Operator
+    p::Int
+    q::Int
+end
+
+function Base.show(io::IO,
+    (
+        e, constraints, translation
+    )::Tuple{TripletExcitationOperator,Constraints,IndexTranslation})
+    print(io, "T_")
+    print_mo_index(io, constraints, translation, e.p, e.q)
+end
+
+function exchange_indices(e::TripletExcitationOperator, mapping)
+    TripletExcitationOperator(
+        exchange_index(e.p, mapping),
+        exchange_index(e.q, mapping)
+    )
+end
+
+function get_all_indices(e::TripletExcitationOperator)
+    (e.p, e.q)
+end
+
+function Base.:(==)(a::TripletExcitationOperator, b::TripletExcitationOperator)
+    (a.p, a.q) == (b.p, b.q)
+end
+
+function Base.isless(a::TripletExcitationOperator, b::TripletExcitationOperator)
+    (a.p, a.q) < (b.p, b.q)
+end
+
+"""
+    τ(p, q)
+
+Constructs an expression containing a single triplet excitation operator.
+"""
+τ(p, q) = Expression(TripletExcitationOperator(p, q))
+
+function convert_to_elementary_operators(o::TripletExcitationOperator)
+    Expression(
+        [
+            (fermiondag(o.p, α)*fermion(o.q, α))[1]
+            -(fermiondag(o.p, β)*fermion(o.q, β))[1]
+        ]
+    )
+end
+
+function act_on_ket(op::TripletExcitationOperator)
+    p = op.p
+    q = op.q
+    τ(p, q) * virtual(p) * occupied(q)
+end
+
+function Base.adjoint(op::TripletExcitationOperator)
+    TripletExcitationOperator(op.q, op.p)
+end