From 6453a6bf640c007f9a0c47291283de95a3549622 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Fri, 10 Jan 2025 17:45:24 -0500 Subject: [PATCH] Start adding back missing functionality --- Project.toml | 6 + src/ITensors.jl | 53 ++ src/LazyApply/LazyApply.jl | 386 +++++++++ src/Ops/Ops.jl | 4 + src/Ops/op.jl | 335 ++++++++ src/Ops/ops_itensor.jl | 73 ++ src/Ops/trotter.jl | 35 + src/SiteTypes/SiteTypes.jl | 14 + src/SiteTypes/SiteTypesChainRulesCoreExt.jl | 7 + src/SiteTypes/sitetype.jl | 827 ++++++++++++++++++++ src/SiteTypes/sitetypes/aliases.jl | 40 + src/SiteTypes/sitetypes/boson.jl | 32 + src/SiteTypes/sitetypes/electron.jl | 336 ++++++++ src/SiteTypes/sitetypes/fermion.jl | 111 +++ src/SiteTypes/sitetypes/generic_sites.jl | 48 ++ src/SiteTypes/sitetypes/qubit.jl | 502 ++++++++++++ src/SiteTypes/sitetypes/qudit.jl | 110 +++ src/SiteTypes/sitetypes/spinhalf.jl | 66 ++ src/SiteTypes/sitetypes/spinone.jl | 143 ++++ src/SiteTypes/sitetypes/tj.jl | 234 ++++++ 20 files changed, 3362 insertions(+) create mode 100644 src/LazyApply/LazyApply.jl create mode 100644 src/Ops/Ops.jl create mode 100644 src/Ops/op.jl create mode 100644 src/Ops/ops_itensor.jl create mode 100644 src/Ops/trotter.jl create mode 100644 src/SiteTypes/SiteTypes.jl create mode 100644 src/SiteTypes/SiteTypesChainRulesCoreExt.jl create mode 100644 src/SiteTypes/sitetype.jl create mode 100644 src/SiteTypes/sitetypes/aliases.jl create mode 100644 src/SiteTypes/sitetypes/boson.jl create mode 100644 src/SiteTypes/sitetypes/electron.jl create mode 100644 src/SiteTypes/sitetypes/fermion.jl create mode 100644 src/SiteTypes/sitetypes/generic_sites.jl create mode 100644 src/SiteTypes/sitetypes/qubit.jl create mode 100644 src/SiteTypes/sitetypes/qudit.jl create mode 100644 src/SiteTypes/sitetypes/spinhalf.jl create mode 100644 src/SiteTypes/sitetypes/spinone.jl create mode 100644 src/SiteTypes/sitetypes/tj.jl diff --git a/Project.toml b/Project.toml index 9cac7b7c30..0a3118ce1b 100644 --- a/Project.toml +++ b/Project.toml @@ -4,8 +4,14 @@ authors = ["ITensor developers and contributors"] version = "0.8.0" [deps] +ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" ITensorBase = "4795dd04-0d67-49bb-8f44-b89c448a1dc7" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +TensorAlgebra = "68bd88dc-f39d-4e12-b2ca-f046b68fcc6a" [compat] +ChainRulesCore = "1.25.1" ITensorBase = "0.1.2" +LinearAlgebra = "1.10.0" +TensorAlgebra = "0.1.1" julia = "1.10" diff --git a/src/ITensors.jl b/src/ITensors.jl index 9d8e2c28ee..12a0ee25b4 100644 --- a/src/ITensors.jl +++ b/src/ITensors.jl @@ -1,5 +1,58 @@ module ITensors +using TensorAlgebra: contract using ITensorBase: ITensor, Index +include("SiteTypes/SiteTypes.jl") +using .SiteTypes: SiteTypes +include("LazyApply/LazyApply.jl") +using .LazyApply: LazyApply +include("Ops/Ops.jl") +using .Ops: Ops + +# TODO: Used in `ITensorMPS.jl`, define in `BackendSelection.jl`. +struct Algorithm{algname} end +macro Algorithm_str(algname) + return :(Algorithm{$(Expr(:quote, Symbol(algname)))}) +end + +# TODO: Used in `ITensorMPS.jl`, decide where or if to define it. +function outer end + +# TODO: Used in `ITensorMPS.jl`, decide where or if to define it. +struct Apply{Args} + args::Args +end + +# TODO: Used in `ITensorMPS.jl`, decide where or if to define it. +macro ts_str(tags) end + +# TODO: Used in `ITensorMPS.jl`, decide where or if to define it. +struct OneITensor end + +# TODO: Used in `ITensorMPS.jl`, define in `ITensorBase.jl`. +function addtags end +function addtags! end +function commonind end +function commoninds end +function dag end +function noprime end +function noprime! end +function prime end +function prime! end +function removetags end +function removetags! end +function replaceprime end +function replaceprime! end +function replacetags end +function replacetags! end +function setprime end +function setprime! end +function settags end +function settags! end +function sim end +function swapprime end +function swapprime! end +function uniqueind end +function uniqueinds end end diff --git a/src/LazyApply/LazyApply.jl b/src/LazyApply/LazyApply.jl new file mode 100644 index 0000000000..f0e69fa71f --- /dev/null +++ b/src/LazyApply/LazyApply.jl @@ -0,0 +1,386 @@ +module LazyApply + +import Base: + ==, + +, + -, + *, + /, + ^, + exp, + adjoint, + copy, + show, + getindex, + length, + isless, + iterate, + firstindex, + lastindex, + keys, + reverse, + size + +export Applied, Scaled, Sum, Prod, Exp, coefficient, argument, expand, materialize, terms + +struct Applied{F,Args<:Tuple,Kwargs<:NamedTuple} + f::F + args::Args + kwargs::Kwargs +end +Applied(f, args::Tuple) = Applied(f, args, (;)) + +materialize(x) = x +function materialize(a::Applied) + return a.f(materialize.(a.args)...; a.kwargs...) +end + +function (a1::Applied == a2::Applied) + return a1.f == a2.f && a1.args == a2.args && a1.kwargs == a2.kwargs +end + +# +# Applied algebra +# + +# Used for dispatch +const Scaled{C<:Number,A} = Applied{typeof(*),Tuple{C,A},NamedTuple{(),Tuple{}}} +const Sum{A} = Applied{typeof(sum),Tuple{Vector{A}},NamedTuple{(),Tuple{}}} +const Prod{A} = Applied{typeof(prod),Tuple{Vector{A}},NamedTuple{(),Tuple{}}} + +# Some convenient empty constructors +Sum{A}() where {A} = Applied(sum, (A[],)) +Prod{A}() where {A} = Applied(prod, (A[],)) + +coefficient(co::Scaled{C}) where {C} = co.args[1] +argument(co::Scaled{C}) where {C} = co.args[2] + +# +# Generic algebra +# + +# 1.3 * Op("X", 1) + 1.3 * Op("X", 2) +# 1.3 * Op("X", 1) * Op("X", 2) + 1.3 * Op("X", 3) * Op("X", 4) +function (a1::Scaled{C,A} + a2::Scaled{C,A}) where {C,A} + return Sum{Scaled{C,A}}() + a1 + a2 +end + +function (a1::Prod{A} + a2::Prod{A}) where {A} + return Sum{Prod{A}}() + a1 + a2 +end + +(c::Number * a::Scaled{C}) where {C} = (c * coefficient(a)) * argument(a) +(a::Scaled{C} * c::Number) where {C} = (coefficient(a) * c) * argument(a) + +-(a::Scaled{C}) where {C} = (-one(C) * a) +-(a::Sum) = (-1 * a) +-(a::Prod) = (-1 * a) + +(os::Sum{A} + o::A) where {A} = Applied(sum, (vcat(os.args[1], [o]),)) +(o::A + os::Sum{A}) where {A} = Applied(sum, (vcat([o], os.args[1]),)) + +(a1::Sum{A} - a2::A) where {A} = a1 + (-a2) +(a1::A - a2::Sum{A}) where {A} = a1 + (-a2) + +(a1::Sum{A} - a2::Prod{A}) where {A} = a1 + (-a2) +(a1::Sum{A} - a2::Scaled{C,Prod{A}}) where {C,A} = a1 + (-a2) +(a1::Sum{A} - a2::Sum{Scaled{C,Prod{A}}}) where {C,A} = a1 + (-a2) + +(a1::Prod{A} * a2::A) where {A} = Applied(prod, (vcat(only(a1.args), [a2]),)) +(a1::A * a2::Prod{A}) where {A} = Applied(prod, (vcat([a1], only(a2.args)),)) + +# Fixes ambiguity error with: +# *(a1::Applied, a2::Sum) +# *(os::Prod{A}, o::A) +(a1::Prod{Sum{A}} * a2::Sum{A}) where {A} = Applied(prod, (vcat(only(a1.args), [a2]),)) + +# 1.3 * Op("X", 1) + 1 * Op("X", 2) +# 1.3 * Op("X", 1) * Op("X", 2) + 1 * Op("X", 3) +# 1.3 * Op("X", 1) * Op("X", 2) + 1 * Op("X", 3) * Op("X", 4) +function (co1::Scaled{C1,A} + co2::Scaled{C2,A}) where {C1,C2,A} + c1, c2 = promote(coefficient(co1), coefficient(co2)) + return c1 * argument(co1) + c2 * argument(co2) +end + +# (1.3 * Op("X", 1)) * (1.3 * Op("X", 2)) +function (co1::Scaled{C1} * co2::Scaled{C2}) where {C1,C2} + c = coefficient(co1) * coefficient(co2) + o = argument(co1) * argument(co2) + return c * o +end + +function (a1::Prod{A} * a2::Scaled{C,A}) where {C,A} + return coefficient(a2) * (a1 * argument(a2)) +end + +function (a1::Prod{A} + a2::Scaled{C,A}) where {C,A} + return one(C) * a1 + Prod{A}() * a2 +end + +# (Op("X", 1) + Op("X", 2)) + (Op("X", 3) + Op("X", 4)) +# (Op("X", 1) * Op("X", 2) + Op("X", 3) * Op("X", 4)) + (Op("X", 5) * Op("X", 6) + Op("X", 7) * Op("X", 8)) +(a1::Sum{A} + a2::Sum{A}) where {A} = Applied(sum, (vcat(a1.args[1], a2.args[1]),)) +(a1::Sum{A} - a2::Sum{A}) where {A} = a1 + (-a2) + +(a1::Prod{A} * a2::Prod{A}) where {A} = Applied(prod, (vcat(only(a1.args), only(a2.args)),)) + +(os::Sum{Scaled{C,A}} + o::A) where {C,A} = os + one(C) * o +(o::A + os::Sum{Scaled{C,A}}) where {C,A} = one(C) * o + os + +# Op("X", 1) + Op("X", 2) + 1.3 * Op("X", 3) +(os::Sum{A} + co::Scaled{C,A}) where {C,A} = one(C) * os + co + +# 1.3 * Op("X", 1) + (Op("X", 2) + Op("X", 3)) +(co::Scaled{C,A} + os::Sum{A}) where {C,A} = co + one(C) * os + +# 1.3 * (Op("X", 1) + Op("X", 2)) +(c::Number * os::Sum) = Applied(sum, (c * os.args[1],)) + +(a1::Applied * a2::Sum) = Applied(sum, (map(a -> a1 * a, only(a2.args)),)) +(a1::Sum * a2::Applied) = Applied(sum, (map(a -> a * a2, only(a1.args)),)) +(a1::Sum * a2::Sum) = Applied(prod, ([a1, a2],)) + +function _expand(a1::Sum, a2::Sum) + return Applied(sum, (vec([a1[i] * a2[j] for i in 1:length(a1), j in 1:length(a2)]),)) +end + +function expand(a::Prod) + if length(a) == 1 + return a[1] + elseif length(a) ≥ 2 + a12 = _expand(a[1], a[2]) + return expand(Applied(prod, (vcat([a12], a[3:end]),))) + end +end + +# (Op("X", 1) + Op("X", 2)) * 1.3 +(os::Sum * c::Number) = c * os + +# (Op("X", 1) + Op("X", 2)) / 1.3 +(os::Sum / c::Number) = inv(c) * os + +# Promotions +function (co1::Scaled{C,Prod{A}} + co2::Scaled{C,A}) where {C,A} + return co1 + coefficient(co2) * Applied(prod, ([argument(co2)],)) +end + +function (a1::Scaled - a2::Scaled) + return a1 + (-a2) +end + +function (a1::Prod{A} + a2::A) where {A} + return a1 + Applied(prod, ([a2],)) +end + +function (a1::Sum{A} + a2::Prod{A}) where {A} + return Prod{A}() * a1 + a2 +end + +function (a1::Sum{A} + a2::Sum{Scaled{C,Prod{A}}}) where {C,A} + return (one(C) * Prod{A}() * a1) + a2 +end + +function (a1::Prod{A} - a2::A) where {A} + return a1 + (-a2) +end + +function (co1::Sum{Scaled{C,Prod{A}}} + co2::Scaled{C,A}) where {C,A} + return co1 + coefficient(co2) * Applied(prod, ([argument(co2)],)) +end + +function (a1::Sum{Scaled{C1,Prod{A}}} - a2::Scaled{C2,A}) where {C1,C2,A} + return a1 + (-a2) +end + +function (a1::Sum{Scaled{C,Prod{A}}} - a2::Prod{A}) where {C,A} + return a1 + (-a2) +end + +function (a1::Sum{Scaled{C1,Prod{A}}} - a2::Scaled{C2,Prod{A}}) where {C1,C2,A} + return a1 + (-a2) +end + +function (a1::Sum{A} + a2::Scaled{C,Prod{A}}) where {C,A} + return Sum{Scaled{C,Prod{A}}}() + a1 + a2 +end + +function (a1::Sum{Scaled{C1,Prod{A}}} + a2::Scaled{C2,A}) where {C1,C2,A} + C = promote_type(C1, C2) + return one(C) * a1 + one(C) * a2 +end + +# (::Sum{Scaled{Bool,Prod{Op}}} + ::Scaled{Float64,Prod{Op}}) +function (a1::Sum{Scaled{C1,A}} + a2::Scaled{C2,A}) where {C1,C2,A} + C = promote_type(C1, C2) + return one(C) * a1 + one(C) * a2 +end + +# TODO: Is this needed? It seems like: +# +# (a1::Sum{A} + a2::A) +# +# is not being called. +function (a1::Sum{Scaled{C,A}} + a2::Scaled{C,A}) where {C,A} + return Applied(sum, (vcat(only(a1.args), [a2]),)) +end + +function (a1::Sum{Scaled{C,Prod{A}}} + a2::Sum{A}) where {C,A} + a2 = one(C) * a2 + a2 = Prod{A}() * a2 + return a1 + one(C) * Prod{A}() * a2 +end + +function (a1::Sum{Prod{A}} + a2::A) where {A} + return a1 + (Prod{A}() * a2) +end + +function (a1::Sum{Prod{A}} + a2::Scaled{C,A}) where {C,A} + return a1 + (Prod{A}() * a2) +end + +function (a1::Sum{Scaled{C,Prod{A}}} + a2::A) where {C,A} + return a1 + one(C) * a2 +end +(a1::Sum{Scaled{C,Prod{A}}} - a2::A) where {C,A} = a1 + (-a2) + +function (a1::Sum{Scaled{C,Prod{A}}} + a2::Sum{Scaled{C,A}}) where {C,A} + return a1 + (Prod{A}() * a2) +end + +function (o::A + os::Sum{Scaled{C,Prod{A}}}) where {C,A} + return one(C) * o + os +end + +function (a::Sum^n::Int) + r = a + for _ in 2:n + r *= a + end + return r +end + +function (a::Prod^n::Int) + r = a + for _ in 2:n + r *= a + end + return r +end + +exp(a::Applied) = Applied(exp, (a,)) + +const Exp{A} = Applied{typeof(exp),Tuple{A},NamedTuple{(),Tuple{}}} +const Adjoint{A} = Applied{typeof(adjoint),Tuple{A},NamedTuple{(),Tuple{}}} + +argument(a::Exp) = a.args[1] + +(c::Number * e::Exp) = Applied(*, (c, e)) +(e::Exp * c::Number) = c * e +(e1::Exp * e2::Exp) = Applied(prod, ([e1, e2],)) +(e1::Applied * e2::Exp) = Applied(prod, ([e1, e2],)) +(e1::Exp * e2::Applied) = Applied(prod, ([e1, e2],)) + +function reverse(a::Prod) + return Applied(prod, (reverse(only(a.args)),)) +end + +adjoint(a::Prod) = Applied(prod, (map(adjoint, reverse(only(a.args))),)) + +# +# Convenient indexing +# + +getindex(a::Union{Sum,Prod}, I...) = only(a.args)[I...] +iterate(a::Union{Sum,Prod}, args...) = iterate(only(a.args), args...) +size(a::Union{Sum,Prod}) = size(only(a.args)) +length(a::Union{Sum,Prod}) = length(only(a.args)) +firstindex(a::Union{Sum,Prod}) = 1 +lastindex(a::Union{Sum,Prod}) = length(a) +keys(a::Union{Sum,Prod}) = 1:length(a) + +length(a::Scaled{C,<:Sum}) where {C} = length(argument(a)) +length(a::Scaled{C,<:Prod}) where {C} = length(argument(a)) +getindex(a::Scaled{C,<:Sum}, I...) where {C} = getindex(argument(a), I...) +getindex(a::Scaled{C,<:Prod}, I...) where {C} = getindex(argument(a), I...) +lastindex(a::Scaled{C,<:Sum}) where {C} = lastindex(argument(a)) +lastindex(a::Scaled{C,<:Prod}) where {C} = lastindex(argument(a)) + +# +# Functions convenient for OpSum code +# + +terms(a::Union{Sum,Prod}) = only(a.args) +terms(a::Scaled{C,<:Union{Sum,Prod}}) where {C} = terms(argument(a)) +copy(a::Applied) = Applied(deepcopy(a.f), deepcopy(a.args), deepcopy(a.kwargs)) +Sum(a::Vector) = Applied(sum, (a,)) +Prod(a::Vector) = Applied(prod, (a,)) +function isless(a1::Applied{F}, a2::Applied{F}) where {F} + return (isless(a1.args, a2.args) && isless(a1.kwargs, a2.kwargs)) +end + +# +# Printing +# + +function show(io::IO, ::MIME"text/plain", a::Sum) + print(io, "sum(\n") + for n in eachindex(a) + print(io, " ", a[n]) + if n ≠ lastindex(a) + print(io, "\n") + end + end + print(io, "\n)") + return nothing +end +show(io::IO, a::Sum) = show(io, MIME("text/plain"), a) + +function show(io::IO, ::MIME"text/plain", a::Prod) + print(io, "prod(\n") + for n in eachindex(a) + print(io, " ", a[n]) + if n ≠ lastindex(a) + print(io, "\n") + end + end + print(io, "\n)") + return nothing +end +show(io::IO, a::Prod) = show(io, MIME("text/plain"), a) + +function show(io::IO, m::MIME"text/plain", a::Exp) + print(io, a.f, "(") + for n in 1:length(a.args) + print(io, a.args[n]) + if n < length(a.args) + print(io, ", ") + end + end + print(io, ")") + return nothing +end +show(io::IO, a::Exp) = show(io, MIME("text/plain"), a) + +function show(io::IO, m::MIME"text/plain", a::Applied) + print(io, a.f, "(") + for n in eachindex(a.args) + print(io, a.args[n]) + if n < length(a.args) + print(io, ", ") + end + end + if !isempty(a.kwargs) + print(io, "; ") + for n in 1:length(a.kwargs) + print(io, keys(a.kwargs)[n], "=", a.kwargs[n]) + if n < length(a.kwargs) + print(io, ", ") + end + end + end + print(io, ")") + return nothing +end +show(io::IO, a::Applied) = show(io, MIME("text/plain"), a) + +end diff --git a/src/Ops/Ops.jl b/src/Ops/Ops.jl new file mode 100644 index 0000000000..fcbff091ea --- /dev/null +++ b/src/Ops/Ops.jl @@ -0,0 +1,4 @@ +module Ops +include("op.jl") +include("trotter.jl") +end diff --git a/src/Ops/op.jl b/src/Ops/op.jl new file mode 100644 index 0000000000..886bb969b9 --- /dev/null +++ b/src/Ops/op.jl @@ -0,0 +1,335 @@ +using ..LazyApply + +import Base: ==, +, -, *, /, convert, exp, show, adjoint, isless, hash + +export Op, OpSum, which_op, site, sites, params, Applied, expand + +##################################################################################### +# General functionality +# + +# Helper function to split a `Tuple` according to the function `f`. +# For example: +# +# julia> t = (1, "X", 1, 2, "Y", 2, "Z", 4) +# (1, "X", 1, 2, "Y", 2, "Z", 4) +# +# julia> split(x -> x isa AbstractString, t) +# [(1,), ("X", 1, 2), ("Y", 2), ("Z", 4)] +# +function split(f, t::Tuple) + n = findall(f, t) + nsplit = length(n) + 1 + s = Vector{Any}(undef, nsplit) + s[1] = t[1:(first(n) - 1)] + for i in 2:(nsplit - 1) + s[i] = t[n[i - 1]:(n[i] - 1)] + end + s[end] = t[last(n):end] + return s +end + +## XXX: Very long compile times: +## https://github.com/JuliaLang/julia/issues/45545 +## +## julia> using ITensors +## +## julia> @time ITensors.Ops.split(x -> x isa String, ("X", 1)) +## 7.588123 seconds (2.34 M allocations: 100.919 MiB, 1.71% gc time, 100.00% compilation time) +## ((), ("X", 1)) +## +## julia> @time ITensors.Ops.split(x -> x isa String, ("X", 1)) +## 0.042590 seconds (88.59 k allocations: 4.823 MiB, 19.13% gc time, 99.84% compilation time) +## ((), ("X", 1)) +## +## function split(f, t::Tuple) +## n = findall(f, t) +## ti = t[1:(first(n) - 1)] +## ts = ntuple(i -> t[n[i]:(n[i + 1] - 1)], length(n) - 1) +## tf = t[last(n):end] +## return ti, ts..., tf +## end + +struct Op + which_op + sites::Tuple + params::NamedTuple + function Op(which_op, site...; kwargs...) + return new(which_op, site, NamedTuple(kwargs)) + end +end + +which_op(o::Op) = o.which_op +name(o::Op) = which_op(o) +sites(o::Op) = o.sites +site(o::Op) = only(sites(o)) +params(o::Op) = o.params + +function (o1::Op == o2::Op) + return o1.which_op == o2.which_op && o1.sites == o2.sites && o1.params == o2.params +end + +function hash(o::Op, h::UInt) + return hash(which_op(o), hash(sites(o), hash(params(o), hash(:Op, h)))) +end + +# Version of `isless` defined for matrices +_isless(a, b) = isless(a, b) +_isless(a::AbstractMatrix, b::AbstractMatrix) = isless(hash(a), hash(b)) +_isless(a::AbstractString, b::AbstractMatrix) = true +_isless(a::AbstractMatrix, b::AbstractString) = !_isless(b, a) + +function isless(o1::Op, o2::Op) + if sites(o1) ≠ sites(o2) + return sites(o1) < sites(o2) + end + if which_op(o1) ≠ which_op(o2) + return _isless(which_op(o1), which_op(o2)) + end + return params(o1) < params(o2) +end + +function isless(o1::Prod{Op}, o2::Prod{Op}) + if length(o1) ≠ length(o2) + return length(o1) < length(o2) + end + for n in 1:length(o1) + if o1[n] ≠ o2[n] + return (o1[n] < o2[n]) + end + end + return false +end + +function isless(o1::Scaled{C1,Prod{Op}}, o2::Scaled{C2,Prod{Op}}) where {C1,C2} + if argument(o1) == argument(o2) + if coefficient(o1) ≈ coefficient(o2) + return false + else + c1 = coefficient(o1) + c2 = coefficient(o2) + #"lexicographic" ordering on complex numbers + return real(c1) < real(c2) || (real(c1) ≈ real(c2) && imag(c1) < imag(c2)) + end + end + return argument(o1) < argument(o2) +end + +## function Op(t::Tuple) +## which_op = first(t) +## site_params = Base.tail(t) +## if last(site_params) isa NamedTuple +## site = Base.front(site_params) +## params = last(site_params) +## else +## site = site_params +## params = (;) +## end +## return Op(which_op, site; params...) +## end + +## function Op(t::Tuple{WhichOp,NamedTuple,Vararg}) where {WhichOp} +## params = t[2] +## which_op = t[1] +## sites = t[3:end] +## return Op(which_op, sites...; params...) +## end + +function sites(a::Union{Sum,Prod}) + s = [] + for n in 1:length(a) + s = s ∪ sites(a[n]) + end + return map(identity, s) +end +sites(a::Scaled{C,<:Sum}) where {C} = sites(argument(a)) +sites(a::Scaled{C,<:Prod}) where {C} = sites(argument(a)) + +params(a::Scaled{C,<:Prod}) where {C} = params(only(argument(a))) + +which_op(a::Scaled{C,Op}) where {C} = which_op(argument(a)) +sites(a::Scaled{C,Op}) where {C} = sites(argument(a)) +params(a::Scaled{C,Op}) where {C} = params(argument(a)) + +# +# Op algebra +# + +function convert(::Type{Scaled{C1,Prod{Op}}}, o::Scaled{C2,Prod{Op}}) where {C1,C2} + c = convert(C1, coefficient(o)) + return c * argument(o) +end + +""" +An `OpSum` represents a sum of operator +terms. + +Often it is used to create matrix +product operator (`MPO`) approximation +of the sum of the terms in the `OpSum` oject. +Each term is a product of local operators +specified by names such as `"Sz"` or `"N"`, +times an optional coefficient which +can be real or complex. + +Which local operator names are available +is determined by the function `op` +associated with the `TagType` defined by +special Index tags, such as `"S=1/2"`, `"S=1"`, +`"Fermion"`, and `"Electron"`. +""" +const OpSum{C} = Sum{Scaled{C,Prod{Op}}} + +# This helps with in-place operations +OpSum() = OpSum{ComplexF64}() + +(o1::Op + o2::Op) = Applied(sum, ([o1, o2],)) +(o1::Op * o2::Op) = Applied(prod, ([o1, o2],)) +-(o::Op) = -one(Int) * o +(o1::Op - o2::Op) = o1 + (-o2) + +(c::Number * o::Op) = Applied(*, (c, o)) +(o::Op * c::Number) = Applied(*, (c, o)) +(o::Op / c::Number) = Applied(*, (inv(c), o)) + +(c::Number * o::Prod{Op}) = Applied(*, (c, o)) +(o::Prod{Op} * c::Number) = Applied(*, (c, o)) +(o::Prod{Op} / c::Number) = Applied(*, (inv(c), o)) + +# 1.3 * Op("X", 1) + Op("X", 2) +# 1.3 * Op("X", 1) * Op("X", 2) + Op("X", 3) +(co1::Scaled{C} + o2::Op) where {C} = co1 + one(C) * o2 + +# Op("X", 1) + 1.3 * Op("X", 2) +(o1::Op + co2::Scaled{C}) where {C} = one(C) * o1 + co2 + +(o1::Op * o2::Sum) = Applied(sum, (map(a -> o1 * a, only(o2.args)),)) +(o1::Sum * o2::Op) = Applied(sum, (map(a -> a * o2, only(o1.args)),)) + +# 1.3 * Op("X", 1) + Op("X", 2) * Op("X", 3) +# 1.3 * Op("X", 1) * Op("X", 2) + Op("X", 3) * Op("X", 4) +(co1::Scaled{C} + o2::Prod{Op}) where {C} = co1 + one(C) * o2 + +# 1.3 * Op("X", 1) * Op("X", 2) +(co1::Scaled{C} * o2::Op) where {C} = co1 * (one(C) * o2) + +exp(o::Op) = Applied(exp, (o,)) + +adjoint(o::Op) = Applied(adjoint, (o,)) +adjoint(o::LazyApply.Adjoint{Op}) = only(o.args) + +(o1::Exp{Op} * o2::Op) = Applied(prod, ([o1, o2],)) + +# +# Tuple interface +# + +const OpSumLike{C} = Union{ + Sum{Op}, + Sum{Scaled{C,Op}}, + Sum{Prod{Op}}, + Sum{Scaled{C,Prod{Op}}}, + Prod{Op}, + Scaled{C,Prod{Op}}, +} + +const WhichOp = Union{AbstractString,AbstractMatrix{<:Number}} + +# Make a `Scaled{C,Prod{Op}}` from a `Tuple` input, +# for example: +# +# (1.2, "X", 1, "Y", 2) -> 1.2 * Op("X", 1) * Op("Y", 2) +# +function op_term(a::Tuple{Number,Vararg}) + c = first(a) + return c * op_term(Base.tail(a)) +end + +function op_site(which_op, params::NamedTuple, sites...) + return Op(which_op, sites...; params...) +end + +function op_site(which_op, sites_params...) + if last(sites_params) isa NamedTuple + sites = Base.front(sites_params) + params = last(sites_params) + return Op(which_op, sites...; params...) + end + return Op(which_op, sites_params...) +end + +function op_term(a::Tuple{Vararg}) + a_split = split(x -> x isa WhichOp, a) + @assert isempty(first(a_split)) + popfirst!(a_split) + o = op_site(first(a_split)...) + popfirst!(a_split) + for aₙ in a_split + o *= op_site(aₙ...) + end + return o +end + +function (o1::OpSumLike + o2::Tuple) + return o1 + op_term(o2) +end + +function (o1::Tuple + o2::OpSumLike) + return op_term(o1) + o2 +end + +function (o1::OpSumLike - o2::Tuple) + return o1 - op_term(o2) +end + +function (o1::Tuple - o2::OpSumLike) + return op_term(o1) - o2 +end + +function (o1::OpSumLike * o2::Tuple) + return o1 * op_term(o2) +end + +function (o1::Tuple * o2::OpSumLike) + return op_term(o1) * o2 +end + +function show(io::IO, ::MIME"text/plain", o::Op) + print(io, which_op(o)) + print(io, sites(o)) + if !isempty(params(o)) + print(io, params(o)) + end + return nothing +end +show(io::IO, o::Op) = show(io, MIME("text/plain"), o) + +function show(io::IO, ::MIME"text/plain", o::Prod{Op}) + for n in 1:length(o) + print(io, o[n]) + if n < length(o) + print(io, " ") + end + end + return nothing +end +show(io::IO, o::Prod{Op}) = show(io, MIME("text/plain"), o) + +function show(io::IO, m::MIME"text/plain", o::Scaled{C,O}) where {C,O<:Union{Op,Prod{Op}}} + c = coefficient(o) + if isreal(c) + c = real(c) + end + print(io, c) + print(io, " ") + show(io, m, argument(o)) + return nothing +end +show(io::IO, o::Scaled{C,Prod{Op}}) where {C} = show(io, MIME("text/plain"), o) + +function show(io::IO, ::MIME"text/plain", o::LazyApply.Adjoint{Op}) + print(io, o') + print(io, "'") + return nothing +end +show(io::IO, o::LazyApply.Adjoint{Op}) = show(io, MIME("text/plain"), o) diff --git a/src/Ops/ops_itensor.jl b/src/Ops/ops_itensor.jl new file mode 100644 index 0000000000..e650376dcf --- /dev/null +++ b/src/Ops/ops_itensor.jl @@ -0,0 +1,73 @@ +using .SiteTypes: SiteTypes, op + +function SiteTypes.op(I::UniformScaling, s::Index...) + return I.λ * op("Id", s...) +end + +function ITensor(o::Op, s::Vector{<:Index}) + return op(o.which_op, map(n -> s[n], o.sites)...; o.params...) +end + +function ITensor(o::Scaled, s::Vector{<:Index}) + c = coefficient(o) + if isreal(c) + c = real(c) + end + return c * ITensor(argument(o), s) +end + +function ITensor(o::Prod, s::Vector{<:Index}) + T = ITensor(true) + for a in o.args[1] + Tₙ = ITensor(a, s) + # TODO: Implement this logic inside `apply` + if hascommoninds(T, Tₙ) + T = T(Tₙ) + else + T *= Tₙ + end + end + return T +end + +function ITensor(o::Sum, s::Vector{<:Index}) + T = ITensor() + for a in o.args[1] + T += ITensor(a, s) + end + return T +end + +function ITensor(o::Exp, s::Vector{<:Index}) + return exp(ITensor(argument(o), s)) +end + +function ITensor(o::LazyApply.Adjoint, s::Vector{<:Index}) + return swapprime(dag(ITensor(o', s)), 0 => 1) +end + +function Sum{ITensor}(o::Sum, s::Vector{<:Index}) + return Applied(sum, (map(oₙ -> ITensor(oₙ, s), only(o.args)),)) +end + +function Prod{ITensor}(o::Prod, s::Vector{<:Index}) + return Applied(prod, (map(oₙ -> ITensor(oₙ, s), only(o.args)),)) +end + +function Prod{ITensor}(o::Scaled{C,Prod{Op}}, s::Vector{<:Index}) where {C} + t = Prod{ITensor}(argument(o), s) + t1 = coefficient(o) * only(t.args)[1] + return Applied(prod, (vcat([t1], only(t.args)[2:end]),)) +end + +function apply(o::Prod{ITensor}, v::ITensor; kwargs...) + ov = v + for oₙ in reverse(only(o.args)) + ov = apply(oₙ, ov; kwargs...) + end + return ov +end + +function (o::Prod{ITensor})(v::ITensor; kwargs...) + return apply(o, v; kwargs...) +end diff --git a/src/Ops/trotter.jl b/src/Ops/trotter.jl new file mode 100644 index 0000000000..9c3b1a5641 --- /dev/null +++ b/src/Ops/trotter.jl @@ -0,0 +1,35 @@ +using ..LazyApply: Applied, Sum + +abstract type ExpAlgorithm end + +struct Exact <: ExpAlgorithm end + +struct Trotter{Order} <: ExpAlgorithm + nsteps::Int +end +Trotter{Order}() where {Order} = Trotter{Order}(1) +Base.one(::Trotter{Order}) where {Order} = Trotter{Order}(1) + +function Base.exp(o::Sum; alg::ExpAlgorithm=Exact()) + return exp(alg, o) +end + +function Base.exp(::Exact, o::Sum) + return Applied(prod, ([Applied(exp, (o,))],)) +end + +function exp_one_step(trotter::Trotter{1}, o::Sum) + exp_o = Applied(prod, (map(exp, reverse(only(o.args))),)) + return exp_o +end + +function exp_one_step(trotter::Trotter{2}, o::Sum) + exp_o_order_1 = exp_one_step(Trotter{1}(), o / 2) + exp_o = reverse(exp_o_order_1) * exp_o_order_1 + return exp_o +end + +function Base.exp(trotter::Trotter, o::Sum) + expδo = exp_one_step(one(trotter), o / trotter.nsteps) + return expδo^trotter.nsteps +end diff --git a/src/SiteTypes/SiteTypes.jl b/src/SiteTypes/SiteTypes.jl new file mode 100644 index 0000000000..f738e242ac --- /dev/null +++ b/src/SiteTypes/SiteTypes.jl @@ -0,0 +1,14 @@ +module SiteTypes +include("sitetype.jl") +include("SiteTypesChainRulesCoreExt.jl") +include("sitetypes/aliases.jl") +include("sitetypes/generic_sites.jl") +include("sitetypes/qubit.jl") +include("sitetypes/spinhalf.jl") +include("sitetypes/spinone.jl") +include("sitetypes/fermion.jl") +include("sitetypes/electron.jl") +include("sitetypes/tj.jl") +include("sitetypes/qudit.jl") +include("sitetypes/boson.jl") +end diff --git a/src/SiteTypes/SiteTypesChainRulesCoreExt.jl b/src/SiteTypes/SiteTypesChainRulesCoreExt.jl new file mode 100644 index 0000000000..625d64e6ab --- /dev/null +++ b/src/SiteTypes/SiteTypesChainRulesCoreExt.jl @@ -0,0 +1,7 @@ +module SiteTypesChainRulesCoreExt +using ChainRulesCore: @non_differentiable +using ..SiteTypes: SiteType, _sitetypes, has_fermion_string +@non_differentiable has_fermion_string(::AbstractString, ::Any) +@non_differentiable SiteType(::Any) +@non_differentiable _sitetypes(::Any) +end diff --git a/src/SiteTypes/sitetype.jl b/src/SiteTypes/sitetype.jl new file mode 100644 index 0000000000..712fa11df6 --- /dev/null +++ b/src/SiteTypes/sitetype.jl @@ -0,0 +1,827 @@ +using ChainRulesCore: @ignore_derivatives +using ..ITensors: + ITensors, Index, ITensor, itensor, dag, onehot, prime, product, swapprime, tags + +@eval struct SiteType{T} + (f::Type{<:SiteType})() = $(Expr(:new, :f)) +end + +# Note that the complicated definition of +# SiteType above is a workaround for performance +# issues when creating parameterized types +# in Julia 1.4 and 1.5-beta. Ideally we +# can just use the following in the future: +# struct SiteType{T} +# end + +""" +SiteType is a parameterized type which allows +making Index tags into Julia types. Use cases +include overloading functions such as `op`, +`siteinds`, and `state` which generate custom +operators, Index arrays, and IndexVals associated +with Index objects having a certain tag. + +To make a SiteType type, you can use the string +macro notation: `SiteType"MyTag"` + +To make a SiteType value or object, you can use +the notation: `SiteType("MyTag")` + +There are currently a few built-in site types +recognized by `jl`. The system is easily extensible +by users. To add new operators to an existing site type, +or to create new site types, you can follow the instructions +[here](https://itensor.github.io/jl/stable/examples/Physics.html). + +The current built-in site types are: + +- `SiteType"S=1/2"` (or `SiteType"S=½"`) +- `SiteType"S=1"` +- `SiteType"Qubit"` +- `SiteType"Qudit"` +- `SiteType"Boson"` +- `SiteType"Fermion"` +- `SiteType"tJ"` +- `SiteType"Electron"` + +# Examples + +Tags on indices get turned into SiteTypes internally, and then +we search for overloads of functions like `op` and `siteind`. +For example: + +```julia +julia> s = siteind("S=1/2") +(dim=2|id=862|"S=1/2,Site") + +julia> @show op("Sz", s); +op(s, "Sz") = ITensor ord=2 +Dim 1: (dim=2|id=862|"S=1/2,Site")' +Dim 2: (dim=2|id=862|"S=1/2,Site") +NDTensors.Dense{Float64,Array{Float64,1}} + 2×2 + 0.5 0.0 + 0.0 -0.5 + +julia> @show op("Sx", s); +op(s, "Sx") = ITensor ord=2 +Dim 1: (dim=2|id=862|"S=1/2,Site")' +Dim 2: (dim=2|id=862|"S=1/2,Site") +NDTensors.Dense{Float64,Array{Float64,1}} + 2×2 + 0.0 0.5 + 0.5 0.0 + +julia> @show op("Sy", s); +op(s, "Sy") = ITensor ord=2 +Dim 1: (dim=2|id=862|"S=1/2,Site")' +Dim 2: (dim=2|id=862|"S=1/2,Site") +NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}} + 2×2 + 0.0 + 0.0im -0.0 - 0.5im + 0.0 + 0.5im 0.0 + 0.0im + +julia> s = siteind("Electron") +(dim=4|id=734|"Electron,Site") + +julia> @show op("Nup", s); +op(s, "Nup") = ITensor ord=2 +Dim 1: (dim=4|id=734|"Electron,Site")' +Dim 2: (dim=4|id=734|"Electron,Site") +NDTensors.Dense{Float64,Array{Float64,1}} + 4×4 + 0.0 0.0 0.0 0.0 + 0.0 1.0 0.0 0.0 + 0.0 0.0 0.0 0.0 + 0.0 0.0 0.0 1.0 +``` + +Many operators are available, for example: + +- `SiteType"S=1/2"`: `"Sz"`, `"Sx"`, `"Sy"`, `"S+"`, `"S-"`, ... +- `SiteType"Electron"`: `"Nup"`, `"Ndn"`, `"Nupdn"`, `"Ntot"`, `"Cup"`, + `"Cdagup"`, `"Cdn"`, `"Cdagdn"`, `"Sz"`, `"Sx"`, `"Sy"`, `"S+"`, `"S-"`, ... +- ... + +You can view the source code for the internal SiteType definitions +and operators that are defined [here](https://github.com/ITensor/jl/tree/main/src/physics/site_types). +""" +SiteType(s::AbstractString) = SiteType{Symbol(s)}() + +SiteType(t::Integer) = SiteType{Symbol(t)}() + +tag(::SiteType{T}) where {T} = T + +macro SiteType_str(s) + return SiteType{Symbol(s)} +end + +# Keep TagType defined for backwards +# compatibility; will be deprecated later +const TagType = SiteType +macro TagType_str(s) + return TagType{Symbol(s)} +end + +#--------------------------------------- +# +# op system +# +#--------------------------------------- + +@eval struct OpName{Name} + (f::Type{<:OpName})() = $(Expr(:new, :f)) +end + +# Note that the complicated definition of +# OpName above is a workaround for performance +# issues when creating parameterized types +# in Julia 1.4 and 1.5-beta. Ideally we +# can just use the following in the future: +# struct OpName{Name} +# end + +""" +OpName is a parameterized type which allows +making strings into Julia types for the purpose +of representing operator names. +The main use of OpName is overloading the +`op!` method which generates operators +for indices with certain tags such as "S=1/2". + +To make a OpName type, you can use the string +macro notation: `OpName"MyTag"`. + +To make an OpName value or object, you can use +the notation: `OpName("myop")` +""" +OpName(s::AbstractString) = OpName{Symbol(s)}() +OpName(s::Symbol) = OpName{s}() +name(::OpName{N}) where {N} = N + +macro OpName_str(s) + return OpName{Symbol(s)} +end + +# Default implementations of op and op! +op(::OpName; kwargs...) = nothing +op(::OpName, ::SiteType; kwargs...) = nothing +op(::OpName, ::SiteType, ::Index...; kwargs...) = nothing +function op( + ::OpName, ::SiteType, ::SiteType, sitetypes_inds::Union{SiteType,Index}...; kwargs... +) + return nothing +end +op!(::ITensor, ::OpName, ::SiteType, ::Index...; kwargs...) = nothing +function op!( + ::ITensor, + ::OpName, + ::SiteType, + ::SiteType, + sitetypes_inds::Union{SiteType,Index}...; + kwargs..., +) + return nothing +end + +# Deprecated version, for backwards compatibility +op(::SiteType, ::Index, ::AbstractString; kwargs...) = nothing + +function _sitetypes(ts::Set) + return collect(SiteType.(ts)) +end + +_sitetypes(i::Index) = _sitetypes(tags(i)) + +""" + op(opname::String, s::Index; kwargs...) + +Return an ITensor corresponding to the operator +named `opname` for the Index `s`. The operator +is constructed by calling an overload of either +the `op` or `op!` methods which take a `SiteType` +argument that corresponds to one of the tags of +the Index `s` and an `OpName"opname"` argument +that corresponds to the input operator name. + +Operator names can be combined using the `"*"` +symbol, for example `"S+*S-"` or `"Sz*Sz*Sz"`. +The result is an ITensor made by forming each operator +then contracting them together in a way corresponding +to the usual operator product or matrix multiplication. + +The `op` system is used by the OpSum +system to convert operator names into ITensors, +and can be used directly such as for applying +operators to MPS. + +# Example + +```julia +s = Index(2, "Site,S=1/2") +Sz = op("Sz", s) +``` + +To see all of the operator names defined for the site types included with +ITensor, please view the [source code](https://github.com/ITensor/jl/tree/main/src/physics/site_types) +for each site type. Note that some site types such as "S=1/2" and "Qubit" +are aliases for each other and share operator definitions. +""" +function op(name::AbstractString, s::Index...; adjoint::Bool=false, kwargs...) + name = strip(name) + # TODO: filter out only commons tags + # if there are multiple indices + commontags_s = commontags(s...) + + # first we handle the + and - algebra, which requires a space between ops to avoid clashing + name_split = nothing + @ignore_derivatives name_split = String.(split(name, " ")) + oplocs = findall(x -> x ∈ ("+", "-"), name_split) + + if !isempty(oplocs) + @ignore_derivatives !isempty(kwargs) && + error("Lazy algebra on parametric gates not allowed") + + # the string representation of algebra ops: ex ["+", "-", "+"] + labels = name_split[oplocs] + # assign coefficients to each term: ex [+1, -1, +1] + coeffs = [1, [(-1)^Int(label == "-") for label in labels]...] + + # grad the name of each operator block separated by an algebra op, and do so by + # making sure blank spaces between opnames are kept when building the new block. + start, opnames = 0, String[] + for oploc in oplocs + finish = oploc + opnames = vcat( + opnames, [prod([name_split[k] * " " for k in (start + 1):(finish - 1)])] + ) + start = oploc + end + opnames = vcat( + opnames, [prod([name_split[k] * " " for k in (start + 1):length(name_split)])] + ) + + # build the vector of blocks and sum + op_list = [ + coeff * (op(opname, s...; kwargs...)) for (coeff, opname) in zip(coeffs, opnames) + ] + return sum(op_list) + end + + # the the multiplication come after + oploc = findfirst("*", name) + if !isnothing(oploc) + op1, op2 = nothing, nothing + @ignore_derivatives begin + op1 = name[1:prevind(name, oploc.start)] + op2 = name[nextind(name, oploc.start):end] + if !(op1[end] == ' ' && op2[1] == ' ') + @warn "($op1*$op2) composite op definition `A*B` deprecated: please use `A * B` instead (with spaces)" + end + end + return product(op(op1, s...; kwargs...), op(op2, s...; kwargs...)) + end + + common_stypes = _sitetypes(commontags_s) + @ignore_derivatives push!(common_stypes, SiteType("Generic")) + opn = OpName(name) + + # + # Try calling a function of the form: + # op(::OpName, ::SiteType, ::Index...; kwargs...) + # + for st in common_stypes + res = op(opn, st, s...; kwargs...) + if !isnothing(res) + adjoint && return swapprime(dag(res), 0 => 1) + return res + end + end + + # + # Try calling a function of the form: + # op(::OpName; kwargs...) + # for backward compatibility with previous + # gate system in PastaQ.jl + # + op_mat = op(opn; kwargs...) + if !isnothing(op_mat) + rs = reverse(s) + res = itensor(op_mat, prime.(rs)..., dag.(rs)...) + adjoint && return swapprime(dag(res), 0 => 1) + return res + end + # + # otherwise try calling a function of the form: + # op(::OpName, ::SiteType; kwargs...) + # which returns a Julia matrix + # + for st in common_stypes + op_mat = op(opn, st; kwargs...) + if !isnothing(op_mat) + rs = reverse(s) + #return itensor(op_mat, prime.(rs)..., dag.(rs)...) + res = itensor(op_mat, prime.(rs)..., dag.(rs)...) + adjoint && return swapprime(dag(res), 0 => 1) + return res + end + end + + # otherwise try calling a function of the form: + # op!(::ITensor, ::OpName, ::SiteType, ::Index...; kwargs...) + # + Op = ITensor(prime.(s)..., dag.(s)...) + for st in common_stypes + op!(Op, opn, st, s...; kwargs...) + if !isempty(Op) + adjoint && return swapprime(dag(Op), 0 => 1) + return Op + end + end + + if length(s) > 1 + # No overloads for common tags found. It might be a + # case of making an operator with mixed site types, + # searching for overloads like: + # op(::OpName, + # ::SiteType..., + # ::Index...; + # kwargs...) + # op!(::ITensor, ::OpName, + # ::SiteType..., + # ::Index...; + # kwargs...) + stypes = _sitetypes.(s) + + for st in Iterators.product(stypes...) + res = op(opn, st..., s...; kwargs...) + if !isnothing(res) + adjoint && return swapprime(dag(res), 0 => 1) + return res + end + end + + Op = ITensor(prime.(s)..., dag.(s)...) + for st in Iterators.product(stypes...) + op!(Op, opn, st..., s...; kwargs...) + if !isempty(Op) + adjoint && return swapprime(dag(Op), 0 => 1) + return Op + end + end + + throw( + ArgumentError( + "Overload of \"op\" or \"op!\" functions not found for operator name \"$name\" and Index tags: $(tags.(s)).", + ), + ) + end + + # + # otherwise try calling a function of the form: + # op(::SiteType, ::Index, ::AbstractString) + # + # (Note: this version is for backwards compatibility + # after version 0.1.10, and may be eventually + # deprecated) + # + for st in common_stypes + res = op(st, s[1], name; kwargs...) + if !isnothing(res) + adjoint && return dag(res) + return res + end + end + + return throw( + ArgumentError( + "Overload of \"op\" or \"op!\" functions not found for operator name \"$name\" and Index tags: $(tags.(s)).", + ), + ) +end + +op(name::AbstractString; kwargs...) = error("Must input indices when creating an `op`.") + +""" + op(X::AbstractArray, s::Index...) + op(M::Matrix, s::Index...) + +Given a matrix M and a set of indices s,t,... +return an operator ITensor with matrix elements given +by M and indices s, s', t, t' + +# Example + +```julia +julia> s = siteind("S=1/2") +(dim=2|id=575|"S=1/2,Site") + +julia> Sz = op([1/2 0; 0 -1/2],s) +ITensor ord=2 (dim=2|id=575|"S=1/2,Site")' (dim=2|id=575|"S=1/2,Site") +NDTensors.Dense{Float64, Vector{Float64}} + +julia> @show Sz +Sz = ITensor ord=2 +Dim 1: (dim=2|id=575|"S=1/2,Site")' +Dim 2: (dim=2|id=575|"S=1/2,Site") +NDTensors.Dense{Float64, Vector{Float64}} + 2×2 + 0.5 0.0 + 0.0 -0.5 +ITensor ord=2 (dim=2|id=575|"S=1/2,Site")' (dim=2|id=575|"S=1/2,Site") +NDTensors.Dense{Float64, Vector{Float64}} +``` +""" +op(X::AbstractArray, s::Index...) = itensor(X, prime.([s...]), dag.([s...])) + +op(opname, s::Vector{<:Index}; kwargs...) = op(opname, s...; kwargs...) + +op(s::Vector{<:Index}, opname; kwargs...) = op(opname, s...; kwargs...) + +# For backwards compatibility, version of `op` +# taking the arguments in the other order: +op(s::Index, opname; kwargs...) = op(opname, s; kwargs...) + +# To ease calling of other op overloads, +# allow passing a string as the op name +op(opname::AbstractString, t::SiteType; kwargs...) = op(OpName(opname), t; kwargs...) + +""" + op(opname::String,sites::Vector{<:Index},n::Int; kwargs...) + +Return an ITensor corresponding to the operator +named `opname` for the n'th Index in the array +`sites`. + +# Example + +```julia +s = siteinds("S=1/2", 4) +Sz2 = op("Sz", s, 2) +``` +""" +function op(opname, s::Vector{<:Index}, ns::NTuple{N,Integer}; kwargs...) where {N} + return op(opname, ntuple(n -> s[ns[n]], Val(N))...; kwargs...) +end + +function op(opname, s::Vector{<:Index}, ns::Vararg{Integer}; kwargs...) + return op(opname, s, ns; kwargs...) +end + +function op(s::Vector{<:Index}, opname, ns::Tuple{Vararg{Integer}}; kwargs...) + return op(opname, s, ns...; kwargs...) +end + +function op(s::Vector{<:Index}, opname, ns::Integer...; kwargs...) + return op(opname, s, ns; kwargs...) +end + +function op(s::Vector{<:Index}, opname, ns::Tuple{Vararg{Integer}}, kwargs::NamedTuple) + return op(opname, s, ns; kwargs...) +end + +function op(s::Vector{<:Index}, opname, ns::Integer, kwargs::NamedTuple) + return op(opname, s, (ns,); kwargs...) +end + +op(s::Vector{<:Index}, o::Tuple) = op(s, o...) + +op(o::Tuple, s::Vector{<:Index}) = op(s, o...) + +op(f::Function, args...; kwargs...) = f(op(args...; kwargs...)) + +function op( + s::Vector{<:Index}, + f::Function, + opname::AbstractString, + ns::Tuple{Vararg{Integer}}; + kwargs..., +) + return f(op(opname, s, ns...; kwargs...)) +end + +function op( + s::Vector{<:Index}, f::Function, opname::AbstractString, ns::Integer...; kwargs... +) + return f(op(opname, s, ns; kwargs...)) +end + +# Here, Ref is used to not broadcast over the vector of indices +# TODO: consider overloading broadcast for `op` with the example +# here: https://discourse.julialang.org/t/how-to-broadcast-over-only-certain-function-arguments/19274/5 +# so that `Ref` isn't needed. +ops(s::Vector{<:Index}, os::AbstractArray) = [op(oₙ, s) for oₙ in os] +ops(os::AbstractVector, s::Vector{<:Index}) = [op(oₙ, s) for oₙ in os] + +@doc """ + ops(s::Vector{<:Index}, os::Vector) + ops(os::Vector, s::Vector{<:Index}) + +Given a list of operators, create ITensors using the collection +of indices. + +# Examples + +```julia +s = siteinds("Qubit", 4) +os = [("H", 1), ("X", 2), ("CX", 2, 4)] +# gates = ops(s, os) +gates = ops(os, s) +``` +""" ops(::Vector{<:Index}, ::AbstractArray) + +#--------------------------------------- +# +# state system +# +#--------------------------------------- + +@eval struct StateName{Name} + (f::Type{<:StateName})() = $(Expr(:new, :f)) +end + +StateName(s::AbstractString) = StateName{Symbol(s)}() +StateName(s::Symbol) = StateName{s}() +name(::StateName{N}) where {N} = N + +macro StateName_str(s) + return StateName{Symbol(s)} +end + +state(::StateName, ::SiteType; kwargs...) = nothing +state(::StateName, ::SiteType, ::Index; kwargs...) = nothing +state!(::ITensor, ::StateName, ::SiteType, ::Index; kwargs...) = nothing + +# Syntax `state("Up", Index(2, "S=1/2"))` +state(sn::String, i::Index; kwargs...) = state(i, sn; kwargs...) + +""" + state(s::Index, name::String; kwargs...) + +Return an ITensor corresponding to the state +named `name` for the Index `s`. The returned +ITensor will have `s` as its only index. + +The terminology here is based on the idea of a +single-site state or wavefunction in physics. + +The `state` function is implemented for various +Index tags by overloading either the +`state` or `state!` methods which take a `SiteType` +argument corresponding to one of the tags of +the Index `s` and an `StateName"name"` argument +that corresponds to the input state name. + +The `state` system is used by the MPS type +to construct product-state MPS and for other purposes. + +# Example + +```julia +s = Index(2, "Site,S=1/2") +sup = state(s,"Up") +sdn = state(s,"Dn") +sxp = state(s,"X+") +sxm = state(s,"X-") +``` +""" +function state(s::Index, name::AbstractString; kwargs...)::ITensor + stypes = _sitetypes(s) + sname = StateName(name) + + # Try calling state(::StateName"Name",::SiteType"Tag",s::Index; kwargs...) + for st in stypes + v = state(sname, st, s; kwargs...) + if !isnothing(v) + if v isa ITensor + return v + else + # TODO: deprecate, only for backwards compatibility. + return itensor(v, s) + end + end + end + + # Try calling state!(::ITensor,::StateName"Name",::SiteType"Tag",s::Index;kwargs...) + T = ITensor(s) + for st in stypes + state!(T, sname, st, s; kwargs...) + !isempty(T) && return T + end + + # + # otherwise try calling a function of the form: + # state(::StateName"Name", ::SiteType"Tag"; kwargs...) + # which returns a Julia vector + # + for st in stypes + v = state(sname, st; kwargs...) + !isnothing(v) && return itensor(v, s) + end + + return throw( + ArgumentError( + "Overload of \"state\" or \"state!\" functions not found for state name \"$name\" and Index tags $(tags(s))", + ), + ) +end + +state(s::Index, n::Integer) = onehot(s => n) + +state(sset::Vector{<:Index}, j::Integer, st; kwargs...) = state(sset[j], st; kwargs...) + +#--------------------------------------- +# +# val system +# +#--------------------------------------- + +@eval struct ValName{Name} + (f::Type{<:ValName})() = $(Expr(:new, :f)) +end + +ValName(s::AbstractString) = ValName{Symbol(s)}() +ValName(s::Symbol) = ValName{s}() +name(::ValName{N}) where {N} = N + +macro ValName_str(s) + return ValName{Symbol(s)} +end + +val(::ValName, ::SiteType) = nothing +val(::AbstractString, ::SiteType) = nothing + +""" + val(s::Index, name::String) + +Return an integer corresponding to the `name` +of a certain value the Index `s` can take. +In other words, the `val` function maps strings +to specific integer values within the range `1:dim(s)`. + +The `val` function is implemented for various +Index tags by overloading methods named `val` +which take a `SiteType` argument corresponding to +one of the tags of the Index `s` and an `ValName"name"` +argument that corresponds to the input name. + +# Example + +```julia +s = Index(2, "Site,S=1/2") +val(s,"Up") == 1 +val(s,"Dn") == 2 + +s = Index(2, "Site,Fermion") +val(s,"Emp") == 1 +val(s,"Occ") == 2 +``` +""" +function val(s::Index, name::AbstractString)::Int + stypes = _sitetypes(s) + sname = ValName(name) + + # Try calling val(::StateName"Name",::SiteType"Tag",) + for st in stypes + res = val(sname, st) + !isnothing(res) && return res + end + + return throw( + ArgumentError("Overload of \"val\" function not found for Index tags $(tags(s))") + ) +end + +val(s::Index, n::Integer) = n + +val(sset::Vector{<:Index}, j::Integer, st) = val(sset[j], st) + +#--------------------------------------- +# +# siteind system +# +#--------------------------------------- + +space(st::SiteType; kwargs...) = nothing + +space(st::SiteType, n::Int; kwargs...) = space(st; kwargs...) + +function space_error_message(st::SiteType) + return "Overload of \"space\",\"siteind\", or \"siteinds\" functions not found for Index tag: $(tag(st))" +end + +function siteind(st::SiteType; addtags="", kwargs...) + sp = space(st; kwargs...) + isnothing(sp) && return nothing + return Index(sp, "Site, $(tag(st)), $addtags") +end + +function siteind(st::SiteType, n; kwargs...) + s = siteind(st; kwargs...) + !isnothing(s) && return addtags(s, "n=$n") + sp = space(st, n; kwargs...) + isnothing(sp) && error(space_error_message(st)) + return Index(sp, "Site, $(tag(st)), n=$n") +end + +siteind(tag::String; kwargs...) = siteind(SiteType(tag); kwargs...) + +siteind(tag::String, n; kwargs...) = siteind(SiteType(tag), n; kwargs...) + +# Special case of `siteind` where integer (dim) provided +# instead of a tag string +#siteind(d::Integer, n::Integer; kwargs...) = Index(d, "Site,n=$n") +function siteind(d::Integer, n::Integer; addtags="", kwargs...) + return Index(d, "Site,n=$n, $addtags") +end + +#--------------------------------------- +# +# siteinds system +# +#--------------------------------------- + +siteinds(::SiteType, N; kwargs...) = nothing + +""" + siteinds(tag::String, N::Integer; kwargs...) + +Create an array of `N` physical site indices of type `tag`. +Keyword arguments can be used to specify quantum number conservation, +see the `space` function corresponding to the site type `tag` for +supported keyword arguments. + +# Example + +```julia +N = 10 +s = siteinds("S=1/2", N; conserve_qns=true) +``` +""" +function siteinds(tag::String, N::Integer; kwargs...) + st = SiteType(tag) + + si = siteinds(st, N; kwargs...) + if !isnothing(si) + return si + end + + return [siteind(st, j; kwargs...) for j in 1:N] +end + +""" + siteinds(f::Function, N::Integer; kwargs...) + +Create an array of `N` physical site indices where the site type at site `n` is given +by `f(n)` (`f` should return a string). +""" +function siteinds(f::Function, N::Integer; kwargs...) + return [siteind(f(n), n; kwargs...) for n in 1:N] +end + +# Special case of `siteinds` where integer (dim) +# provided instead of a tag string +""" + siteinds(d::Integer, N::Integer; kwargs...) + +Create an array of `N` site indices, each of dimension `d`. + +# Keywords +- `addtags::String`: additional tags to be added to all indices +""" +function siteinds(d::Integer, N::Integer; kwargs...) + return [siteind(d, n; kwargs...) for n in 1:N] +end + +#--------------------------------------- +# +# has_fermion_string system +# +#--------------------------------------- + +has_fermion_string(operator::AbstractArray{<:Number}, s::Index; kwargs...)::Bool = false + +has_fermion_string(::OpName, ::SiteType) = nothing + +function has_fermion_string(opname::AbstractString, s::Index; kwargs...)::Bool + opname = strip(opname) + + # Interpret operator names joined by * + # as acting sequentially on the same site + starpos = findfirst(isequal('*'), opname) + if !isnothing(starpos) + op1 = opname[1:prevind(opname, starpos)] + op2 = opname[nextind(opname, starpos):end] + return xor(has_fermion_string(op1, s; kwargs...), has_fermion_string(op2, s; kwargs...)) + end + + Ntags = length(tags(s)) + stypes = _sitetypes(s) + opn = OpName(opname) + for st in stypes + res = has_fermion_string(opn, st) + !isnothing(res) && return res + end + return false +end diff --git a/src/SiteTypes/sitetypes/aliases.jl b/src/SiteTypes/sitetypes/aliases.jl new file mode 100644 index 0000000000..727ad7f990 --- /dev/null +++ b/src/SiteTypes/sitetypes/aliases.jl @@ -0,0 +1,40 @@ +alias(::OpName"c") = OpName"C"() +alias(::OpName"cdag") = OpName"Cdag"() +alias(::OpName"c†") = OpName"Cdag"() +alias(::OpName"n") = OpName"N"() +alias(::OpName"a") = OpName"A"() +alias(::OpName"adag") = OpName"Adag"() +alias(::OpName"a↑") = OpName"Aup"() +alias(::OpName"a↓") = OpName"Adn"() +alias(::OpName"a†↓") = OpName"Adagdn"() +alias(::OpName"a†↑") = OpName"Adagup"() +alias(::OpName"a†") = OpName"Adag"() +alias(::OpName"c↑") = OpName"Cup"() +alias(::OpName"c↓") = OpName"Cdn"() +alias(::OpName"c†↑") = OpName"Cdagup"() +alias(::OpName"c†↓") = OpName"Cdagdn"() +alias(::OpName"n↑") = OpName"Nup"() +alias(::OpName"n↓") = OpName"Ndn"() +alias(::OpName"n↑↓") = OpName"Nupdn"() +alias(::OpName"ntot") = OpName"Ntot"() +alias(::OpName"F↑") = OpName"Fup"() +alias(::OpName"F↓") = OpName"Fdn"() +alias(::OpName"I") = OpName"Id"() + +alias(::OpName"S²") = OpName"S2"() +alias(::OpName"Sᶻ") = OpName"Sz"() +alias(::OpName"Sʸ") = OpName"Sy"() +alias(::OpName"iSʸ") = OpName"iSy"() +alias(::OpName"Sˣ") = OpName"Sx"() +alias(::OpName"S⁻") = OpName"S-"() +alias(::OpName"Sminus") = OpName"S-"() +alias(::OpName"Sm") = OpName"S-"() +alias(::OpName"S⁺") = OpName"S+"() +alias(::OpName"Splus") = OpName"S+"() +alias(::OpName"Sp") = OpName"S+"() +alias(::OpName"projUp") = OpName"ProjUp"() +alias(::OpName"projDn") = OpName"ProjDn"() + +alias(::OpName"Proj0") = OpName"ProjUp"() +alias(::OpName"Proj1") = OpName"ProjDn"() +alias(::OpName"Rn̂") = OpName"Rn"() diff --git a/src/SiteTypes/sitetypes/boson.jl b/src/SiteTypes/sitetypes/boson.jl new file mode 100644 index 0000000000..387edad277 --- /dev/null +++ b/src/SiteTypes/sitetypes/boson.jl @@ -0,0 +1,32 @@ + +alias(::SiteType"Boson") = SiteType"Qudit"() + +""" + space(::SiteType"Boson"; + dim = 2, + conserve_qns = false, + conserve_number = false, + qnname_number = "Number") + +Create the Hilbert space for a site of type "Boson". + +Optionally specify the conserved symmetries and their quantum number labels. +""" +space(st::SiteType"Boson"; kwargs...) = space(alias(st); kwargs...) + +val(vn::ValName, st::SiteType"Boson") = val(vn, alias(st)) + +function state(sn::StateName, st::SiteType"Boson", s::Index; kwargs...) + return state(sn, alias(st), s; kwargs...) +end + +function op(on::OpName, st::SiteType"Boson", ds::Int...; kwargs...) + return op(on, alias(st), ds...; kwargs...) +end + +function op(on::OpName, st::SiteType"Boson", s1::Index, s_tail::Index...; kwargs...) + rs = reverse((s1, s_tail...)) + ds = dim.(rs) + opmat = op(on, st, ds...; kwargs...) + return itensor(opmat, prime.(rs)..., dag.(rs)...) +end diff --git a/src/SiteTypes/sitetypes/electron.jl b/src/SiteTypes/sitetypes/electron.jl new file mode 100644 index 0000000000..d8fd90848e --- /dev/null +++ b/src/SiteTypes/sitetypes/electron.jl @@ -0,0 +1,336 @@ +""" + space(::SiteType"Electron"; + conserve_qns = false, + conserve_sz = conserve_qns, + conserve_nf = conserve_qns, + conserve_nfparity = conserve_qns, + qnname_sz = "Sz", + qnname_nf = "Nf", + qnname_nfparity = "NfParity") + +Create the Hilbert space for a site of type "Electron". + +Optionally specify the conserved symmetries and their quantum number labels. +""" +function space( + ::SiteType"Electron"; + conserve_qns=false, + conserve_sz=conserve_qns, + conserve_nf=conserve_qns, + conserve_nfparity=conserve_qns, + qnname_sz="Sz", + qnname_nf="Nf", + qnname_nfparity="NfParity", + # Deprecated + conserve_parity=nothing, +) + if !isnothing(conserve_parity) + conserve_nfparity = conserve_parity + end + if conserve_sz && conserve_nf + return [ + QN((qnname_nf, 0, -1), (qnname_sz, 0)) => 1 + QN((qnname_nf, 1, -1), (qnname_sz, +1)) => 1 + QN((qnname_nf, 1, -1), (qnname_sz, -1)) => 1 + QN((qnname_nf, 2, -1), (qnname_sz, 0)) => 1 + ] + elseif conserve_nf + return [ + QN(qnname_nf, 0, -1) => 1 + QN(qnname_nf, 1, -1) => 2 + QN(qnname_nf, 2, -1) => 1 + ] + elseif conserve_sz + return [ + QN((qnname_sz, 0), (qnname_nfparity, 0, -2)) => 1 + QN((qnname_sz, +1), (qnname_nfparity, 1, -2)) => 1 + QN((qnname_sz, -1), (qnname_nfparity, 1, -2)) => 1 + QN((qnname_sz, 0), (qnname_nfparity, 0, -2)) => 1 + ] + elseif conserve_nfparity + return [ + QN(qnname_nfparity, 0, -2) => 1 + QN(qnname_nfparity, 1, -2) => 2 + QN(qnname_nfparity, 0, -2) => 1 + ] + end + return 4 +end + +val(::ValName"Emp", ::SiteType"Electron") = 1 +val(::ValName"Up", ::SiteType"Electron") = 2 +val(::ValName"Dn", ::SiteType"Electron") = 3 +val(::ValName"UpDn", ::SiteType"Electron") = 4 +val(::ValName"0", st::SiteType"Electron") = val(ValName("Emp"), st) +val(::ValName"↑", st::SiteType"Electron") = val(ValName("Up"), st) +val(::ValName"↓", st::SiteType"Electron") = val(ValName("Dn"), st) +val(::ValName"↑↓", st::SiteType"Electron") = val(ValName("UpDn"), st) + +state(::StateName"Emp", ::SiteType"Electron") = [1.0, 0, 0, 0] +state(::StateName"Up", ::SiteType"Electron") = [0.0, 1, 0, 0] +state(::StateName"Dn", ::SiteType"Electron") = [0.0, 0, 1, 0] +state(::StateName"UpDn", ::SiteType"Electron") = [0.0, 0, 0, 1] +state(::StateName"0", st::SiteType"Electron") = state(StateName("Emp"), st) +state(::StateName"↑", st::SiteType"Electron") = state(StateName("Up"), st) +state(::StateName"↓", st::SiteType"Electron") = state(StateName("Dn"), st) +state(::StateName"↑↓", st::SiteType"Electron") = state(StateName("UpDn"), st) + +function op(::OpName"Nup", ::SiteType"Electron") + return [ + 0.0 0.0 0.0 0.0 + 0.0 1.0 0.0 0.0 + 0.0 0.0 0.0 0.0 + 0.0 0.0 0.0 1.0 + ] +end +function op(on::OpName"n↑", st::SiteType"Electron") + return op(alias(on), st) +end + +function op(::OpName"Ndn", ::SiteType"Electron") + return [ + 0.0 0.0 0.0 0.0 + 0.0 0.0 0.0 0.0 + 0.0 0.0 1.0 0.0 + 0.0 0.0 0.0 1.0 + ] +end +function op(on::OpName"n↓", st::SiteType"Electron") + return op(alias(on), st) +end + +function op(::OpName"Nupdn", ::SiteType"Electron") + return [ + 0.0 0.0 0.0 0.0 + 0.0 0.0 0.0 0.0 + 0.0 0.0 0.0 0.0 + 0.0 0.0 0.0 1.0 + ] +end +function op(on::OpName"n↑↓", st::SiteType"Electron") + return op(alias(on), st) +end + +function op(::OpName"Ntot", ::SiteType"Electron") + return [ + 0.0 0.0 0.0 0.0 + 0.0 1.0 0.0 0.0 + 0.0 0.0 1.0 0.0 + 0.0 0.0 0.0 2.0 + ] +end +function op(on::OpName"ntot", st::SiteType"Electron") + return op(alias(on), st) +end + +function op(::OpName"Cup", ::SiteType"Electron") + return [ + 0.0 1.0 0.0 0.0 + 0.0 0.0 0.0 0.0 + 0.0 0.0 0.0 1.0 + 0.0 0.0 0.0 0.0 + ] +end +function op(on::OpName"c↑", st::SiteType"Electron") + return op(alias(on), st) +end + +function op(::OpName"Cdagup", ::SiteType"Electron") + return [ + 0.0 0.0 0.0 0.0 + 1.0 0.0 0.0 0.0 + 0.0 0.0 0.0 0.0 + 0.0 0.0 1.0 0.0 + ] +end +function op(on::OpName"c†↑", st::SiteType"Electron") + return op(alias(on), st) +end + +function op(::OpName"Cdn", ::SiteType"Electron") + return [ + 0.0 0.0 1.0 0.0 + 0.0 0.0 0.0 -1.0 + 0.0 0.0 0.0 0.0 + 0.0 0.0 0.0 0.0 + ] +end +function op(on::OpName"c↓", st::SiteType"Electron") + return op(alias(on), st) +end + +function op(::OpName"Cdagdn", ::SiteType"Electron") + return [ + 0.0 0.0 0.0 0.0 + 0.0 0.0 0.0 0.0 + 1.0 0.0 0.0 0.0 + 0.0 -1.0 0.0 0.0 + ] +end +function op(::OpName"c†↓", st::SiteType"Electron") + return op(OpName("Cdagdn"), st) +end + +function op(::OpName"Aup", ::SiteType"Electron") + return [ + 0.0 1.0 0.0 0.0 + 0.0 0.0 0.0 0.0 + 0.0 0.0 0.0 1.0 + 0.0 0.0 0.0 0.0 + ] +end +function op(::OpName"a↑", st::SiteType"Electron") + return op(OpName("Aup"), st) +end + +function op(::OpName"Adagup", ::SiteType"Electron") + return [ + 0.0 0.0 0.0 0.0 + 1.0 0.0 0.0 0.0 + 0.0 0.0 0.0 0.0 + 0.0 0.0 1.0 0.0 + ] +end +function op(::OpName"a†↑", st::SiteType"Electron") + return op(OpName("Adagup"), st) +end + +function op(::OpName"Adn", ::SiteType"Electron") + return [ + 0.0 0.0 1.0 0.0 + 0.0 0.0 0.0 1.0 + 0.0 0.0 0.0 0.0 + 0.0 0.0 0.0 0.0 + ] +end +function op(::OpName"a↓", st::SiteType"Electron") + return op(OpName("Adn"), st) +end + +function op(::OpName"Adagdn", ::SiteType"Electron") + return [ + 0.0 0.0 0.0 0.0 + 0.0 0.0 0.0 0.0 + 1.0 0.0 0.0 0.0 + 0.0 1.0 0.0 0.0 + ] +end +function op(::OpName"a†↓", st::SiteType"Electron") + return op(OpName("Adagdn"), st) +end + +function op(::OpName"F", ::SiteType"Electron") + return [ + 1.0 0.0 0.0 0.0 + 0.0 -1.0 0.0 0.0 + 0.0 0.0 -1.0 0.0 + 0.0 0.0 0.0 1.0 + ] +end + +function op(::OpName"Fup", ::SiteType"Electron") + return [ + 1.0 0.0 0.0 0.0 + 0.0 -1.0 0.0 0.0 + 0.0 0.0 1.0 0.0 + 0.0 0.0 0.0 -1.0 + ] +end +function op(::OpName"F↑", st::SiteType"Electron") + return op(OpName("Fup"), st) +end + +function op(::OpName"Fdn", ::SiteType"Electron") + return [ + 1.0 0.0 0.0 0.0 + 0.0 1.0 0.0 0.0 + 0.0 0.0 -1.0 0.0 + 0.0 0.0 0.0 -1.0 + ] +end +function op(::OpName"F↓", st::SiteType"Electron") + return op(OpName("Fdn"), st) +end + +function op(::OpName"Sz", ::SiteType"Electron") + #Op[s' => 2, s => 2] = +0.5 + #return Op[s' => 3, s => 3] = -0.5 + return [ + 0.0 0.0 0.0 0.0 + 0.0 0.5 0.0 0.0 + 0.0 0.0 -0.5 0.0 + 0.0 0.0 0.0 0.0 + ] +end + +function op(::OpName"Sᶻ", st::SiteType"Electron") + return op(OpName("Sz"), st) +end + +function op(::OpName"Sx", ::SiteType"Electron") + return [ + 0.0 0.0 0.0 0.0 + 0.0 0.0 0.5 0.0 + 0.0 0.5 0.0 0.0 + 0.0 0.0 0.0 0.0 + ] +end + +function op(::OpName"Sˣ", st::SiteType"Electron") + return op(OpName("Sx"), st) +end + +function op(::OpName"S+", ::SiteType"Electron") + return [ + 0.0 0.0 0.0 0.0 + 0.0 0.0 1.0 0.0 + 0.0 0.0 0.0 0.0 + 0.0 0.0 0.0 0.0 + ] +end + +function op(::OpName"S⁺", st::SiteType"Electron") + return op(OpName("S+"), st) +end +function op(::OpName"Sp", st::SiteType"Electron") + return op(OpName("S+"), st) +end +function op(::OpName"Splus", st::SiteType"Electron") + return op(OpName("S+"), st) +end + +function op(::OpName"S-", ::SiteType"Electron") + return [ + 0.0 0.0 0.0 0.0 + 0.0 0.0 0.0 0.0 + 0.0 1.0 0.0 0.0 + 0.0 0.0 0.0 0.0 + ] +end + +function op(::OpName"S⁻", st::SiteType"Electron") + return op(OpName("S-"), st) +end +function op(::OpName"Sm", st::SiteType"Electron") + return op(OpName("S-"), st) +end +function op(::OpName"Sminus", st::SiteType"Electron") + return op(OpName("S-"), st) +end + +has_fermion_string(::OpName"Cup", ::SiteType"Electron") = true +function has_fermion_string(on::OpName"c↑", st::SiteType"Electron") + return has_fermion_string(alias(on), st) +end +has_fermion_string(::OpName"Cdagup", ::SiteType"Electron") = true +function has_fermion_string(on::OpName"c†↑", st::SiteType"Electron") + return has_fermion_string(alias(on), st) +end +has_fermion_string(::OpName"Cdn", ::SiteType"Electron") = true +function has_fermion_string(on::OpName"c↓", st::SiteType"Electron") + return has_fermion_string(alias(on), st) +end +has_fermion_string(::OpName"Cdagdn", ::SiteType"Electron") = true +function has_fermion_string(on::OpName"c†↓", st::SiteType"Electron") + return has_fermion_string(alias(on), st) +end diff --git a/src/SiteTypes/sitetypes/fermion.jl b/src/SiteTypes/sitetypes/fermion.jl new file mode 100644 index 0000000000..008bdf0877 --- /dev/null +++ b/src/SiteTypes/sitetypes/fermion.jl @@ -0,0 +1,111 @@ + +""" + space(::SiteType"Fermion"; + conserve_qns=false, + conserve_nf=conserve_qns, + conserve_nfparity=conserve_qns, + qnname_nf = "Nf", + qnname_nfparity = "NfParity", + qnname_sz = "Sz", + conserve_sz = false) + +Create the Hilbert space for a site of type "Fermion". + +Optionally specify the conserved symmetries and their quantum number labels. +""" +function space( + ::SiteType"Fermion"; + conserve_qns=false, + conserve_nf=conserve_qns, + conserve_nfparity=conserve_qns, + qnname_nf="Nf", + qnname_nfparity="NfParity", + qnname_sz="Sz", + conserve_sz=false, + # Deprecated + conserve_parity=nothing, +) + if !isnothing(conserve_parity) + conserve_nfparity = conserve_parity + end + if conserve_sz == true + conserve_sz = "Up" + end + if conserve_nf && conserve_sz == "Up" + zer = QN((qnname_nf, 0, -1), (qnname_sz, 0)) => 1 + one = QN((qnname_nf, 1, -1), (qnname_sz, 1)) => 1 + return [zer, one] + elseif conserve_nf && conserve_sz == "Dn" + zer = QN((qnname_nf, 0, -1), (qnname_sz, 0)) => 1 + one = QN((qnname_nf, 1, -1), (qnname_sz, -1)) => 1 + return [zer, one] + elseif conserve_nfparity && conserve_sz == "Up" + zer = QN((qnname_nfparity, 0, -2), (qnname_sz, 0)) => 1 + one = QN((qnname_nfparity, 1, -2), (qnname_sz, 1)) => 1 + return [zer, one] + elseif conserve_nfparity && conserve_sz == "Dn" + zer = QN((qnname_nfparity, 0, -2), (qnname_sz, 0)) => 1 + one = QN((qnname_nfparity, 1, -2), (qnname_sz, -1)) => 1 + return [zer, one] + elseif conserve_nf + zer = QN(qnname_nf, 0, -1) => 1 + one = QN(qnname_nf, 1, -1) => 1 + return [zer, one] + elseif conserve_nfparity + zer = QN(qnname_nfparity, 0, -2) => 1 + one = QN(qnname_nfparity, 1, -2) => 1 + return [zer, one] + end + return 2 +end + +val(::ValName"Emp", ::SiteType"Fermion") = 1 +val(::ValName"Occ", ::SiteType"Fermion") = 2 +val(::ValName"0", st::SiteType"Fermion") = val(ValName("Emp"), st) +val(::ValName"1", st::SiteType"Fermion") = val(ValName("Occ"), st) + +state(::StateName"Emp", ::SiteType"Fermion") = [1.0 0.0] +state(::StateName"Occ", ::SiteType"Fermion") = [0.0 1.0] +state(::StateName"0", st::SiteType"Fermion") = state(StateName("Emp"), st) +state(::StateName"1", st::SiteType"Fermion") = state(StateName("Occ"), st) + +function op!(Op::ITensor, ::OpName"N", ::SiteType"Fermion", s::Index) + return Op[s' => 2, s => 2] = 1.0 +end +function op!(Op::ITensor, on::OpName"n", st::SiteType"Fermion", s::Index) + return op!(Op, alias(on), st, s) +end + +function op!(Op::ITensor, ::OpName"C", ::SiteType"Fermion", s::Index) + return Op[s' => 1, s => 2] = 1.0 +end +function op!(Op::ITensor, on::OpName"c", st::SiteType"Fermion", s::Index) + return op!(Op, alias(on), st, s) +end + +function op!(Op::ITensor, ::OpName"Cdag", ::SiteType"Fermion", s::Index) + return Op[s' => 2, s => 1] = 1.0 +end +function op!(Op::ITensor, on::OpName"c†", st::SiteType"Fermion", s::Index) + return op!(Op, alias(on), st, s) +end +function op!(Op::ITensor, on::OpName"cdag", st::SiteType"Fermion", s::Index) + return op!(Op, alias(on), st, s) +end + +function op!(Op::ITensor, ::OpName"F", ::SiteType"Fermion", s::Index) + Op[s' => 1, s => 1] = +1.0 + return Op[s' => 2, s => 2] = -1.0 +end + +has_fermion_string(::OpName"C", ::SiteType"Fermion") = true +function has_fermion_string(on::OpName"c", st::SiteType"Fermion") + return has_fermion_string(alias(on), st) +end +has_fermion_string(::OpName"Cdag", ::SiteType"Fermion") = true +function has_fermion_string(on::OpName"c†", st::SiteType"Fermion") + return has_fermion_string(alias(on), st) +end +function has_fermion_string(on::OpName"cdag", st::SiteType"Fermion") + return has_fermion_string(alias(on), st) +end diff --git a/src/SiteTypes/sitetypes/generic_sites.jl b/src/SiteTypes/sitetypes/generic_sites.jl new file mode 100644 index 0000000000..445d849848 --- /dev/null +++ b/src/SiteTypes/sitetypes/generic_sites.jl @@ -0,0 +1,48 @@ +using LinearAlgebra: I +using ..ITensors: ITensor, itensor, settensor! + +function op!( + o::ITensor, ::OpName"Id", ::SiteType"Generic", s1::Index, sn::Index...; eltype=Float64 +) + s = (s1, sn...) + n = prod(dim.(s)) + t = itensor(Matrix(one(eltype) * I, n, n), prime.(s)..., dag.(s)...) + return settensor!(o, tensor(t)) +end + +function op!(o::ITensor, on::OpName"I", st::SiteType"Generic", s::Index...; kwargs...) + return op!(o, alias(on), st, s...; kwargs...) +end + +function op!(o::ITensor, ::OpName"F", st::SiteType"Generic", s::Index; kwargs...) + return op!(o, OpName("Id"), st, s; kwargs...) +end + +function default_random_matrix(eltype::Type, s::Index...) + n = prod(dim.(s)) + return randn(eltype, n, n) +end + +# Haar-random unitary +# +# Reference: +# Section 4.6 +# http://math.mit.edu/~edelman/publications/random_matrix_theory.pdf +function op!( + o::ITensor, + ::OpName"RandomUnitary", + ::SiteType"Generic", + s1::Index, + sn::Index...; + eltype=ComplexF64, + random_matrix=default_random_matrix(eltype, s1, sn...), +) + s = (s1, sn...) + Q, _ = NDTensors.qr_positive(random_matrix) + t = itensor(Q, prime.(s)..., dag.(s)...) + return settensor!(o, tensor(t)) +end + +function op!(o::ITensor, ::OpName"randU", st::SiteType"Generic", s::Index...; kwargs...) + return op!(o, OpName("RandomUnitary"), st, s...; kwargs...) +end diff --git a/src/SiteTypes/sitetypes/qubit.jl b/src/SiteTypes/sitetypes/qubit.jl new file mode 100644 index 0000000000..e07bd547f0 --- /dev/null +++ b/src/SiteTypes/sitetypes/qubit.jl @@ -0,0 +1,502 @@ +using ..ITensors: ITensors + +# +# Qubit site type +# + +# Define Qubit space in terms of +# Qubit/2 space, but use different +# defaults for QN names + +""" + space(::SiteType"Qubit"; + conserve_qns = false, + conserve_parity = conserve_qns, + conserve_number = false, + qnname_parity = "Parity", + qnname_number = "Number") + +Create the Hilbert space for a site of type "Qubit". + +Optionally specify the conserved symmetries and their quantum number labels. +""" +function space( + ::SiteType"Qubit"; + conserve_qns=false, + conserve_parity=conserve_qns, + conserve_number=false, + qnname_parity="Parity", + qnname_number="Number", +) + if conserve_number && conserve_parity + return [ + QN((qnname_number, 0), (qnname_parity, 0, 2)) => 1, + QN((qnname_number, 1), (qnname_parity, 1, 2)) => 1, + ] + elseif conserve_number + return [QN(qnname_number, 0) => 1, QN(qnname_number, 1) => 1] + elseif conserve_parity + return [QN(qnname_parity, 0, 2) => 1, QN(qnname_parity, 1, 2) => 1] + end + return 2 +end + +val(::ValName"0", ::SiteType"Qubit") = 1 +val(::ValName"1", ::SiteType"Qubit") = 2 +val(::ValName"Up", ::SiteType"Qubit") = 1 +val(::ValName"Dn", ::SiteType"Qubit") = 2 +val(::ValName"↑", ::SiteType"Qubit") = 1 +val(::ValName"↓", ::SiteType"Qubit") = 2 + +state(::StateName"0", ::SiteType"Qubit") = [1.0, 0.0] +state(::StateName"1", ::SiteType"Qubit") = [0.0, 1.0] +state(::StateName"+", ::SiteType"Qubit") = [1.0, 1.0] / √2 +state(::StateName"-", ::SiteType"Qubit") = [1.0, -1.0] / √2 +state(::StateName"i", ::SiteType"Qubit") = [1.0, im] / √2 +state(::StateName"-i", ::SiteType"Qubit") = [1.0, -im] / √2 +state(::StateName"Up", t::SiteType"Qubit") = state(StateName("0"), t) +state(::StateName"Dn", t::SiteType"Qubit") = state(StateName("1"), t) +state(::StateName"↑", t::SiteType"Qubit") = state(StateName("0"), t) +state(::StateName"↓", t::SiteType"Qubit") = state(StateName("1"), t) + +# Pauli eingenstates +state(::StateName"X+", t::SiteType"Qubit") = state(StateName("+"), t) +state(::StateName"Xp", t::SiteType"Qubit") = state(StateName("+"), t) +state(::StateName"X-", t::SiteType"Qubit") = state(StateName("-"), t) +state(::StateName"Xm", t::SiteType"Qubit") = state(StateName("-"), t) + +state(::StateName"Y+", t::SiteType"Qubit") = state(StateName("i"), t) +state(::StateName"Yp", t::SiteType"Qubit") = state(StateName("i"), t) +state(::StateName"Y-", t::SiteType"Qubit") = state(StateName("-i"), t) +state(::StateName"Ym", t::SiteType"Qubit") = state(StateName("-i"), t) + +state(::StateName"Z+", t::SiteType"Qubit") = state(StateName("0"), t) +state(::StateName"Zp", t::SiteType"Qubit") = state(StateName("0"), t) +state(::StateName"Z-", t::SiteType"Qubit") = state(StateName("1"), t) +state(::StateName"Zm", t::SiteType"Qubit") = state(StateName("1"), t) + +# SIC-POVMs +state(::StateName"Tetra1", t::SiteType"Qubit") = state(StateName("Z+"), t) +state(::StateName"Tetra2", t::SiteType"Qubit") = [ + 1 / √3 + √2 / √3 +] +state(::StateName"Tetra3", t::SiteType"Qubit") = [ + 1 / √3 + √2 / √3 * exp(im * 2π / 3) +] +state(::StateName"Tetra4", t::SiteType"Qubit") = [ + 1 / √3 + √2 / √3 * exp(im * 4π / 3) +] + +# +# 1-Qubit gates +# +op(::OpName"X", ::SiteType"Qubit") = [ + 0 1 + 1 0 +] + +op(::OpName"σx", t::SiteType"Qubit") = op("X", t) + +op(::OpName"σ1", t::SiteType"Qubit") = op("X", t) + +op(::OpName"Y", ::SiteType"Qubit") = [ + 0.0 -1.0im + 1.0im 0.0 +] + +op(::OpName"σy", t::SiteType"Qubit") = op("Y", t) + +op(::OpName"σ2", t::SiteType"Qubit") = op("Y", t) + +op(::OpName"iY", ::SiteType"Qubit") = [ + 0 1 + -1 0 +] +op(::OpName"iσy", t::SiteType"Qubit") = op("iY", t) + +op(::OpName"iσ2", t::SiteType"Qubit") = op("iY", t) + +op(::OpName"Z", ::SiteType"Qubit") = [ + 1 0 + 0 -1 +] + +op(::OpName"σz", t::SiteType"Qubit") = op("Z", t) + +op(::OpName"σ3", t::SiteType"Qubit") = op("Z", t) + +function op(::OpName"√NOT", ::SiteType"Qubit") + return [ + (1 + im)/2 (1 - im)/2 + (1 - im)/2 (1 + im)/2 + ] +end + +op(::OpName"√X", t::SiteType"Qubit") = op("√NOT", t) + +op(::OpName"H", ::SiteType"Qubit") = [ + 1/sqrt(2) 1/sqrt(2) + 1/sqrt(2) -1/sqrt(2) +] + +# Rϕ with ϕ = π/2 +op(::OpName"Phase", ::SiteType"Qubit"; ϕ::Number=π / 2) = [ + 1 0 + 0 exp(im * ϕ) +] + +op(::OpName"P", t::SiteType"Qubit"; kwargs...) = op("Phase", t; kwargs...) + +op(::OpName"S", t::SiteType"Qubit") = op("Phase", t; ϕ=π / 2) + +## Rϕ with ϕ = π/4 +op(::OpName"π/8", ::SiteType"Qubit") = [ + 1 0 + 0 1 / sqrt(2)+im / sqrt(2) +] + +op(::OpName"T", t::SiteType"Qubit") = op("π/8", t) + +# Rotation around X-axis +function op(::OpName"Rx", ::SiteType"Qubit"; θ::Number) + return [ + cos(θ / 2) -im*sin(θ / 2) + -im*sin(θ / 2) cos(θ / 2) + ] +end + +# Rotation around Y-axis +function op(::OpName"Ry", ::SiteType"Qubit"; θ::Number) + return [ + cos(θ / 2) -sin(θ / 2) + sin(θ / 2) cos(θ / 2) + ] +end + +# Rotation around Z-axis +function op(::OpName"Rz", ::SiteType"Qubit"; θ=nothing, ϕ=nothing) + isone(count(isnothing, (θ, ϕ))) || error( + "Must specify the keyword argument `θ` (or the deprecated `ϕ`) when creating an Rz gate, but not both.", + ) + isnothing(θ) && (θ = ϕ) + return [ + exp(-im * θ / 2) 0 + 0 exp(im * θ / 2) + ] +end + +# Rotation around generic axis n̂ +function op(::OpName"Rn", ::SiteType"Qubit"; θ::Real, ϕ::Real, λ::Real) + return [ + cos(θ / 2) -exp(im * λ)*sin(θ / 2) + exp(im * ϕ)*sin(θ / 2) exp(im * (ϕ + λ))*cos(θ / 2) + ] +end + +function op(on::OpName"Rn̂", t::SiteType"Qubit"; kwargs...) + return op(alias(on), t; kwargs...) +end + +# +# 2-Qubit gates +# + +op(::OpName"CNOT", ::SiteType"Qubit") = [ + 1 0 0 0 + 0 1 0 0 + 0 0 0 1 + 0 0 1 0 +] + +op(::OpName"CX", t::SiteType"Qubit") = op("CNOT", t) + +op(::OpName"CY", ::SiteType"Qubit") = [ + 1 0 0 0 + 0 1 0 0 + 0 0 0 -im + 0 0 im 0 +] + +op(::OpName"CZ", ::SiteType"Qubit") = [ + 1 0 0 0 + 0 1 0 0 + 0 0 1 0 + 0 0 0 -1 +] + +function op(::OpName"CPHASE", ::SiteType"Qubit"; ϕ::Number) + return [ + 1 0 0 0 + 0 1 0 0 + 0 0 1 0 + 0 0 0 exp(im * ϕ) + ] +end +op(::OpName"Cphase", t::SiteType"Qubit"; kwargs...) = op("CPHASE", t; kwargs...) + +function op(::OpName"CRx", ::SiteType"Qubit"; θ::Number) + return [ + 1 0 0 0 + 0 1 0 0 + 0 0 cos(θ / 2) -im*sin(θ / 2) + 0 0 -im*sin(θ / 2) cos(θ / 2) + ] +end +op(::OpName"CRX", t::SiteType"Qubit"; kwargs...) = op("CRx", t; kwargs...) + +function op(::OpName"CRy", ::SiteType"Qubit"; θ::Number) + return [ + 1 0 0 0 + 0 1 0 0 + 0 0 cos(θ / 2) -sin(θ / 2) + 0 0 sin(θ / 2) cos(θ / 2) + ] +end +op(::OpName"CRY", t::SiteType"Qubit"; kwargs...) = op("CRy", t; kwargs...) + +function op(::OpName"CRz", ::SiteType"Qubit"; ϕ=nothing, θ=nothing) + isone(count(isnothing, (θ, ϕ))) || error( + "Must specify the keyword argument `θ` (or the deprecated `ϕ`) when creating a CRz gate, but not both.", + ) + isnothing(θ) && (θ = ϕ) + return [ + 1 0 0 0 + 0 1 0 0 + 0 0 exp(-im * θ / 2) 0 + 0 0 0 exp(im * θ / 2) + ] +end +op(::OpName"CRZ", t::SiteType"Qubit"; kwargs...) = op("CRz", t; kwargs...) + +function op(::OpName"CRn", ::SiteType"Qubit"; θ::Number, ϕ::Number, λ::Number) + return [ + 1 0 0 0 + 0 1 0 0 + 0 0 cos(θ / 2) -exp(im * λ)*sin(θ / 2) + 0 0 exp(im * ϕ)*sin(θ / 2) exp(im * (ϕ + λ))*cos(θ / 2) + ] +end +function op(::OpName"CRn̂", t::SiteType"Qubit"; kwargs...) + return op("CRn", t; kwargs...) +end + +op(::OpName"SWAP", ::SiteType"Qubit") = [ + 1 0 0 0 + 0 0 1 0 + 0 1 0 0 + 0 0 0 1 +] +op(::OpName"Swap", t::SiteType"Qubit") = op("SWAP", t) + +function op(::OpName"√SWAP", ::SiteType"Qubit") + return [ + 1 0 0 0 + 0 (1 + im)/2 (1 - im)/2 0 + 0 (1 - im)/2 (1 + im)/2 0 + 0 0 0 1 + ] +end +op(::OpName"√Swap", t::SiteType"Qubit") = op("√SWAP", t) + +op(::OpName"iSWAP", t::SiteType"Qubit") = [ + 1 0 0 0 + 0 0 im 0 + 0 im 0 0 + 0 0 0 1 +] +op(::OpName"iSwap", t::SiteType"Qubit") = op("iSWAP", t) + +function op(::OpName"√iSWAP", t::SiteType"Qubit") + return [ + 1 0 0 0 + 0 1/√2 im/√2 0 + 0 im/√2 1/√2 0 + 0 0 0 1 + ] +end +op(::OpName"√iSwap", t::SiteType"Qubit") = op("√iSWAP", t) + +# Ising (XX) coupling gate +function op(::OpName"Rxx", t::SiteType"Qubit"; ϕ::Number) + return [ + cos(ϕ) 0 0 -im*sin(ϕ) + 0 cos(ϕ) -im*sin(ϕ) 0 + 0 -im*sin(ϕ) cos(ϕ) 0 + -im*sin(ϕ) 0 0 cos(ϕ) + ] +end +op(::OpName"RXX", t::SiteType"Qubit"; kwargs...) = op("Rxx", t; kwargs...) + +# Ising (YY) coupling gate +function op(::OpName"Ryy", ::SiteType"Qubit"; ϕ::Number) + return [ + cos(ϕ) 0 0 im*sin(ϕ) + 0 cos(ϕ) -im*sin(ϕ) 0 + 0 -im*sin(ϕ) cos(ϕ) 0 + im*sin(ϕ) 0 0 cos(ϕ) + ] +end +op(::OpName"RYY", t::SiteType"Qubit"; kwargs...) = op("Ryy", t; kwargs...) + +# Ising (XY) coupling gate +function op(::OpName"Rxy", t::SiteType"Qubit"; ϕ::Number) + return [ + 1 0 0 0 + 0 cos(ϕ) -im*sin(ϕ) 0 + 0 -im*sin(ϕ) cos(ϕ) 0 + 0 0 0 1 + ] +end +op(::OpName"RXY", t::SiteType"Qubit"; kwargs...) = op("Rxy", t; kwargs...) + +# Ising (ZZ) coupling gate +function op(::OpName"Rzz", ::SiteType"Qubit"; ϕ::Number) + return [ + exp(-im * ϕ) 0 0 0 + 0 exp(im * ϕ) 0 0 + 0 0 exp(im * ϕ) 0 + 0 0 0 exp(-im * ϕ) + ] +end +op(::OpName"RZZ", t::SiteType"Qubit"; kwargs...) = op("Rzz", t; kwargs...) + +# +# 3-Qubit gates +# + +function op(::OpName"Toffoli", ::SiteType"Qubit") + return [ + 1 0 0 0 0 0 0 0 + 0 1 0 0 0 0 0 0 + 0 0 1 0 0 0 0 0 + 0 0 0 1 0 0 0 0 + 0 0 0 0 1 0 0 0 + 0 0 0 0 0 1 0 0 + 0 0 0 0 0 0 0 1 + 0 0 0 0 0 0 1 0 + ] +end + +op(::OpName"CCNOT", t::SiteType"Qubit") = op("Toffoli", t) + +op(::OpName"CCX", t::SiteType"Qubit") = op("Toffoli", t) + +op(::OpName"TOFF", t::SiteType"Qubit") = op("Toffoli", t) + +function op(::OpName"Fredkin", ::SiteType"Qubit") + return [ + 1 0 0 0 0 0 0 0 + 0 1 0 0 0 0 0 0 + 0 0 1 0 0 0 0 0 + 0 0 0 1 0 0 0 0 + 0 0 0 0 1 0 0 0 + 0 0 0 0 0 0 1 0 + 0 0 0 0 0 1 0 0 + 0 0 0 0 0 0 0 1 + ] +end + +op(::OpName"CSWAP", t::SiteType"Qubit") = op("Fredkin", t) +op(::OpName"CSwap", t::SiteType"Qubit") = op("Fredkin", t) + +op(::OpName"CS", t::SiteType"Qubit") = op("Fredkin", t) + +# +# 4-Qubit gates +# + +function op(::OpName"CCCNOT", ::SiteType"Qubit") + return [ + 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 + 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 + ] +end + +# spin-full operators +op(::OpName"Sz", ::SiteType"Qubit") = [ + 0.5 0.0 + 0.0 -0.5 +] + +op(on::OpName"Sᶻ", t::SiteType"Qubit") = op(alias(on), t) + +op(::OpName"S+", ::SiteType"Qubit") = [ + 0 1 + 0 0 +] + +op(on::OpName"S⁺", t::SiteType"Qubit") = op(alias(on), t) + +op(on::OpName"Splus", t::SiteType"Qubit") = op(alias(on), t) + +op(::OpName"S-", ::SiteType"Qubit") = [ + 0 0 + 1 0 +] + +op(on::OpName"S⁻", t::SiteType"Qubit") = op(alias(on), t) + +op(on::OpName"Sminus", t::SiteType"Qubit") = op(alias(on), t) + +op(::OpName"Sx", ::SiteType"Qubit") = [ + 0.0 0.5 + 0.5 0.0 +] + +op(on::OpName"Sˣ", t::SiteType"Qubit") = op(alias(on), t) + +op(::OpName"iSy", ::SiteType"Qubit") = [ + 0.0 0.5 + -0.5 0.0 +] + +op(on::OpName"iSʸ", t::SiteType"Qubit") = op(alias(on), t) + +op(::OpName"Sy", ::SiteType"Qubit") = [ + 0.0 -0.5im + 0.5im 0.0 +] + +op(on::OpName"Sʸ", t::SiteType"Qubit") = op(alias(on), t) + +op(::OpName"S2", ::SiteType"Qubit") = [ + 0.75 0.0 + 0.0 0.75 +] + +op(on::OpName"S²", t::SiteType"Qubit") = op(alias(on), t) + +op(::OpName"ProjUp", ::SiteType"Qubit") = [ + 1 0 + 0 0 +] + +op(on::OpName"projUp", t::SiteType"Qubit") = op(alias(on), t) + +op(on::OpName"Proj0", t::SiteType"Qubit") = op(alias(on), t) + +op(::OpName"ProjDn", ::SiteType"Qubit") = [ + 0 0 + 0 1 +] + +op(on::OpName"projDn", t::SiteType"Qubit") = op(alias(on), t) + +op(on::OpName"Proj1", t::SiteType"Qubit") = op(alias(on), t) diff --git a/src/SiteTypes/sitetypes/qudit.jl b/src/SiteTypes/sitetypes/qudit.jl new file mode 100644 index 0000000000..6644cf06e8 --- /dev/null +++ b/src/SiteTypes/sitetypes/qudit.jl @@ -0,0 +1,110 @@ +using ChainRulesCore: @non_differentiable + +""" + space(::SiteType"Qudit"; + dim = 2, + conserve_qns = false, + conserve_number = false, + qnname_number = "Number") + +Create the Hilbert space for a site of type "Qudit". + +Optionally specify the conserved symmetries and their quantum number labels. +""" +function space( + ::SiteType"Qudit"; + dim=2, + conserve_qns=false, + conserve_number=conserve_qns, + qnname_number="Number", +) + if conserve_number + return [QN(qnname_number, n - 1) => 1 for n in 1:dim] + end + return dim +end + +function val(::ValName{N}, ::SiteType"Qudit") where {N} + return parse(Int, String(N)) + 1 +end + +function state(::StateName{N}, ::SiteType"Qudit", s::Index) where {N} + n = parse(Int, String(N)) + st = zeros(dim(s)) + st[n + 1] = 1.0 + return itensor(st, s) +end + +# one-body operators +function op(::OpName"Id", ::SiteType"Qudit", ds::Int...) + d = prod(ds) + return Matrix(1.0I, d, d) +end +op(on::OpName"I", st::SiteType"Qudit", ds::Int...) = op(alias(on), st, ds...) +op(on::OpName"F", st::SiteType"Qudit", ds::Int...) = op(OpName"Id"(), st, ds...) + +function op(::OpName"Adag", ::SiteType"Qudit", d::Int) + mat = zeros(d, d) + for k in 1:(d - 1) + mat[k + 1, k] = √k + end + return mat +end +op(on::OpName"adag", st::SiteType"Qudit", d::Int) = op(alias(on), st, d) +op(on::OpName"a†", st::SiteType"Qudit", d::Int) = op(alias(on), st, d) + +function op(::OpName"A", ::SiteType"Qudit", d::Int) + mat = zeros(d, d) + for k in 1:(d - 1) + mat[k, k + 1] = √k + end + return mat +end +op(on::OpName"a", st::SiteType"Qudit", d::Int) = op(alias(on), st, d) + +function op(::OpName"N", ::SiteType"Qudit", d::Int) + mat = zeros(d, d) + for k in 1:d + mat[k, k] = k - 1 + end + return mat +end +op(on::OpName"n", st::SiteType"Qudit", d::Int) = op(alias(on), st, d) + +# two-body operators +function op(::OpName"ab", st::SiteType"Qudit", d1::Int, d2::Int) + return kron(op(OpName("a"), st, d1), op(OpName("a"), st, d2)) +end + +function op(::OpName"a†b", st::SiteType"Qudit", d1::Int, d2::Int) + return kron(op(OpName("a†"), st, d1), op(OpName("a"), st, d2)) +end + +function op(::OpName"ab†", st::SiteType"Qudit", d1::Int, d2::Int) + return kron(op(OpName("a"), st, d1), op(OpName("a†"), st, d2)) +end + +function op(::OpName"a†b†", st::SiteType"Qudit", d1::Int, d2::Int) + return kron(op(OpName("a†"), st, d1), op(OpName("a†"), st, d2)) +end + +# interface +function op(on::OpName, st::SiteType"Qudit", s1::Index, s_tail::Index...; kwargs...) + rs = reverse((s1, s_tail...)) + ds = dim.(rs) + opmat = op(on, st, ds...; kwargs...) + return itensor(opmat, prime.(rs)..., dag.(rs)...) +end + +function op(on::OpName, st::SiteType"Qudit"; kwargs...) + return error("`op` can't be called without indices or dimensions.") +end + +# Zygote +@non_differentiable op(::OpName"ab", ::SiteType"Qudit", ::Int, ::Int) +@non_differentiable op(::OpName"a†b", ::SiteType"Qudit", ::Int, ::Int) +@non_differentiable op(::OpName"ab†", ::SiteType"Qudit", ::Int, ::Int) +@non_differentiable op(::OpName"a†b†", ::SiteType"Qudit", ::Int, ::Int) +@non_differentiable op(::OpName"a", ::SiteType"Qudit", ::Int) +@non_differentiable op(::OpName"a†", ::SiteType"Qudit", ::Int) +@non_differentiable op(::OpName"N", ::SiteType"Qudit", ::Int) diff --git a/src/SiteTypes/sitetypes/spinhalf.jl b/src/SiteTypes/sitetypes/spinhalf.jl new file mode 100644 index 0000000000..c374d47c0d --- /dev/null +++ b/src/SiteTypes/sitetypes/spinhalf.jl @@ -0,0 +1,66 @@ + +""" + space(::SiteType"S=1/2"; + conserve_qns = false, + conserve_sz = conserve_qns, + conserve_szparity = false, + qnname_sz = "Sz", + qnname_szparity = "SzParity") + +Create the Hilbert space for a site of type "S=1/2". + +Optionally specify the conserved symmetries and their quantum number labels. +""" +function space( + ::SiteType"S=1/2"; + conserve_qns=false, + conserve_sz=conserve_qns, + conserve_szparity=false, + qnname_sz="Sz", + qnname_szparity="SzParity", +) + if conserve_sz && conserve_szparity + return [ + QN((qnname_sz, +1), (qnname_szparity, 1, 2)) => 1, + QN((qnname_sz, -1), (qnname_szparity, 0, 2)) => 1, + ] + elseif conserve_sz + return [QN(qnname_sz, +1) => 1, QN(qnname_sz, -1) => 1] + elseif conserve_szparity + return [QN(qnname_szparity, 1, 2) => 1, QN(qnname_szparity, 0, 2) => 1] + end + return 2 +end + +# Use Qubit definition of any operator/state +# called using S=1/2 SiteType +function val(vn::ValName, ::SiteType"S=1/2"; kwargs...) + return val(vn, SiteType("Qubit"); kwargs...) +end + +function state(sn::StateName, ::SiteType"S=1/2"; kwargs...) + return state(sn, SiteType("Qubit"); kwargs...) +end + +op(o::OpName, ::SiteType"S=1/2"; kwargs...) = op(o, SiteType("Qubit"); kwargs...) + +# Support the tag "SpinHalf" as equivalent to "S=1/2" +space(::SiteType"SpinHalf"; kwargs...) = space(SiteType("S=1/2"); kwargs...) + +val(name::ValName, ::SiteType"SpinHalf") = val(name, SiteType("S=1/2")) + +state(name::StateName, ::SiteType"SpinHalf") = state(name, SiteType("S=1/2")) + +function op(o::OpName, ::SiteType"SpinHalf"; kwargs...) + return op(o, SiteType("S=1/2"); kwargs...) +end + +# Support the tag "S=½" as equivalent to "S=1/2" + +space(::SiteType"S=½"; kwargs...) = space(SiteType("S=1/2"); kwargs...) + +val(name::ValName, ::SiteType"S=½") = val(name, SiteType("S=1/2")) + +state(name::StateName, ::SiteType"S=½") = state(name, SiteType("S=1/2")) + +op(o::OpName, ::SiteType"S=½"; kwargs...) = op(o, SiteType("S=1/2"); kwargs...) diff --git a/src/SiteTypes/sitetypes/spinone.jl b/src/SiteTypes/sitetypes/spinone.jl new file mode 100644 index 0000000000..1de6902243 --- /dev/null +++ b/src/SiteTypes/sitetypes/spinone.jl @@ -0,0 +1,143 @@ +using ..ITensors: complex!, QN + +alias(::SiteType"SpinOne") = SiteType"S=1"() + +""" + space(::SiteType"S=1"; + conserve_qns = false, + conserve_sz = conserve_qns, + qnname_sz = "Sz") + +Create the Hilbert space for a site of type "S=1". + +Optionally specify the conserved symmetries and their quantum number labels. +""" +function space( + ::SiteType"S=1"; conserve_qns=false, conserve_sz=conserve_qns, qnname_sz="Sz" +) + if conserve_sz + return [QN(qnname_sz, +2) => 1, QN(qnname_sz, 0) => 1, QN(qnname_sz, -2) => 1] + end + return 3 +end + +val(::ValName"Up", ::SiteType"S=1") = 1 +val(::ValName"Z0", ::SiteType"S=1") = 2 +val(::ValName"Dn", ::SiteType"S=1") = 3 + +val(::ValName"↑", st::SiteType"S=1") = 1 +val(::ValName"0", st::SiteType"S=1") = 2 +val(::ValName"↓", st::SiteType"S=1") = 3 + +val(::ValName"Z+", ::SiteType"S=1") = 1 +# -- Z0 is already defined above -- +val(::ValName"Z-", ::SiteType"S=1") = 3 + +state(::StateName"Up", ::SiteType"S=1") = [1.0, 0.0, 0.0] +state(::StateName"Z0", ::SiteType"S=1") = [0.0, 1.0, 0.0] +state(::StateName"Dn", ::SiteType"S=1") = [0.0, 0.0, 1.0] + +state(::StateName"↑", st::SiteType"S=1") = [1.0, 0.0, 0.0] +state(::StateName"0", st::SiteType"S=1") = [0.0, 1.0, 0.0] +state(::StateName"↓", st::SiteType"S=1") = [0.0, 0.0, 1.0] + +state(::StateName"Z+", st::SiteType"S=1") = [1.0, 0.0, 0.0] +# -- Z0 is already defined above -- +state(::StateName"Z-", st::SiteType"S=1") = [0.0, 0.0, 1.0] + +state(::StateName"X+", ::SiteType"S=1") = [1 / 2, 1 / sqrt(2), 1 / 2] +state(::StateName"X0", ::SiteType"S=1") = [-1 / sqrt(2), 0, 1 / sqrt(2)] +state(::StateName"X-", ::SiteType"S=1") = [1 / 2, -1 / sqrt(2), 1 / 2] + +state(::StateName"Y+", ::SiteType"S=1") = [-1 / 2, -im / sqrt(2), 1 / 2] +state(::StateName"Y0", ::SiteType"S=1") = [1 / sqrt(2), 0, 1 / sqrt(2)] +state(::StateName"Y-", ::SiteType"S=1") = [-1 / 2, im / sqrt(2), 1 / 2] + +op(::OpName"Sz", ::SiteType"S=1") = [ + 1.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 -1.0 +] + +op(on::OpName"Sᶻ", t::SiteType"S=1") = op(alias(on), t) + +op(::OpName"S+", ::SiteType"S=1") = [ + 0.0 √2 0.0 + 0.0 0.0 √2 + 0.0 0.0 0.0 +] + +op(on::OpName"S⁺", t::SiteType"S=1") = op(alias(on), t) +op(on::OpName"Splus", t::SiteType"S=1") = op(alias(on), t) +op(on::OpName"Sp", t::SiteType"S=1") = op(alias(on), t) + +op(::OpName"S-", ::SiteType"S=1") = [ + 0.0 0.0 0.0 + √2 0.0 0.0 + 0.0 √2 0.0 +] + +op(on::OpName"S⁻", t::SiteType"S=1") = op(alias(on), t) +op(on::OpName"Sminus", t::SiteType"S=1") = op(alias(on), t) +op(on::OpName"Sm", t::SiteType"S=1") = op(alias(on), t) + +op(::OpName"Sx", ::SiteType"S=1") = [ + 0.0 1/√2 0.0 + 1/√2 0.0 1/√2 + 0.0 1/√2 0.0 +] + +op(on::OpName"Sˣ", t::SiteType"S=1") = op(alias(on), t) + +op(::OpName"iSy", ::SiteType"S=1") = [ + 0.0 1/√2 0.0 + -1/√2 0.0 1/√2 + 0.0 -1/√2 0.0 +] + +op(on::OpName"iSʸ", t::SiteType"S=1") = op(alias(on), t) + +op(::OpName"Sy", ::SiteType"S=1") = [ + 0.0 -im/√2 0.0 + im/√2 0.0 -im/√2 + 0.0 im/√2 0.0 +] + +op(on::OpName"Sʸ", t::SiteType"S=1") = op(alias(on), t) + +op(::OpName"Sz2", ::SiteType"S=1") = [ + 1.0 0.0 0.0 + 0.0 0.0 0.0 + 0.0 0.0 1.0 +] + +op(::OpName"Sx2", ::SiteType"S=1") = [ + 0.5 0.0 0.5 + 0.0 1.0 0.0 + 0.5 0.0 0.5 +] + +op(::OpName"Sy2", ::SiteType"S=1") = [ + 0.5 0.0 -0.5 + 0.0 1.0 0.0 + -0.5 0.0 0.5 +] + +op(::OpName"S2", ::SiteType"S=1") = [ + 2.0 0.0 0.0 + 0.0 2.0 0.0 + 0.0 0.0 2.0 +] + +op(on::OpName"S²", t::SiteType"S=1") = op(alias(on), t) + +space(st::SiteType"SpinOne"; kwargs...) = space(alias(st); kwargs...) + +state(name::StateName, st::SiteType"SpinOne") = state(name, alias(st)) +val(name::ValName, st::SiteType"SpinOne") = val(name, alias(st)) + +function op!(Op::ITensor, o::OpName, st::SiteType"SpinOne", s::Index) + return op!(Op, o, alias(st), s) +end + +op(o::OpName, st::SiteType"SpinOne") = op(o, alias(st)) diff --git a/src/SiteTypes/sitetypes/tj.jl b/src/SiteTypes/sitetypes/tj.jl new file mode 100644 index 0000000000..063f46d721 --- /dev/null +++ b/src/SiteTypes/sitetypes/tj.jl @@ -0,0 +1,234 @@ + +""" + space(::SiteType"tJ"; + conserve_qns = false, + conserve_sz = conserve_qns, + conserve_nf = conserve_qns, + conserve_nfparity = conserve_qns, + qnname_sz = "Sz", + qnname_nf = "Nf", + qnname_nfparity = "NfParity") + +Create the Hilbert space for a site of type "tJ". + +Optionally specify the conserved symmetries and their quantum number labels. +""" +function space( + ::SiteType"tJ"; + conserve_qns=false, + conserve_sz=conserve_qns, + conserve_nf=conserve_qns, + conserve_nfparity=conserve_qns, + qnname_sz="Sz", + qnname_nf="Nf", + qnname_nfparity="NfParity", + # Deprecated + conserve_parity=nothing, +) + if !isnothing(conserve_parity) + conserve_nfparity = conserve_parity + end + if conserve_sz && conserve_nf + return [ + QN((qnname_nf, 0, -1), (qnname_sz, 0)) => 1 + QN((qnname_nf, 1, -1), (qnname_sz, +1)) => 1 + QN((qnname_nf, 1, -1), (qnname_sz, -1)) => 1 + ] + elseif conserve_nf + return [ + QN(qnname_nf, 0, -1) => 1 + QN(qnname_nf, 1, -1) => 2 + ] + elseif conserve_sz + return [ + QN((qnname_sz, 0), (qnname_nfparity, 0, -2)) => 1 + QN((qnname_sz, +1), (qnname_nfparity, 1, -2)) => 1 + QN((qnname_sz, -1), (qnname_nfparity, 1, -2)) => 1 + ] + elseif conserve_nfparity + return [ + QN(qnname_nfparity, 0, -2) => 1 + QN(qnname_nfparity, 1, -2) => 2 + ] + end + return 3 +end + +val(::ValName"Emp", ::SiteType"tJ") = 1 +val(::ValName"Up", ::SiteType"tJ") = 2 +val(::ValName"Dn", ::SiteType"tJ") = 3 +val(::ValName"0", st::SiteType"tJ") = val(ValName("Emp"), st) +val(::ValName"↑", st::SiteType"tJ") = val(ValName("Up"), st) +val(::ValName"↓", st::SiteType"tJ") = val(ValName("Dn"), st) + +state(::StateName"Emp", ::SiteType"tJ") = [1.0, 0, 0] +state(::StateName"Up", ::SiteType"tJ") = [0.0, 1, 0] +state(::StateName"Dn", ::SiteType"tJ") = [0.0, 0, 1] +state(::StateName"0", st::SiteType"tJ") = state(StateName("Emp"), st) +state(::StateName"↑", st::SiteType"tJ") = state(StateName("Up"), st) +state(::StateName"↓", st::SiteType"tJ") = state(StateName("Dn"), st) + +function op!(Op::ITensor, ::OpName"Nup", ::SiteType"tJ", s::Index) + return Op[s' => 2, s => 2] = 1.0 +end +function op!(Op::ITensor, on::OpName"n↑", st::SiteType"tJ", s::Index) + return op!(Op, alias(on), st, s) +end + +function op!(Op::ITensor, ::OpName"Ndn", ::SiteType"tJ", s::Index) + return Op[s' => 3, s => 3] = 1.0 +end +function op!(Op::ITensor, on::OpName"n↓", st::SiteType"tJ", s::Index) + return op!(Op, alias(on), st, s) +end + +function op!(Op::ITensor, ::OpName"Ntot", ::SiteType"tJ", s::Index) + Op[s' => 2, s => 2] = 1.0 + return Op[s' => 3, s => 3] = 1.0 +end +function op!(Op::ITensor, on::OpName"ntot", st::SiteType"tJ", s::Index) + return op!(Op, alias(on), st, s) +end + +function op!(Op::ITensor, ::OpName"Cup", ::SiteType"tJ", s::Index) + return Op[s' => 1, s => 2] = 1.0 +end +function op!(Op::ITensor, on::OpName"c↑", st::SiteType"tJ", s::Index) + return op!(Op, alias(on), st, s) +end + +function op!(Op::ITensor, ::OpName"Cdagup", ::SiteType"tJ", s::Index) + return Op[s' => 2, s => 1] = 1.0 +end +function op!(Op::ITensor, on::OpName"c†↑", st::SiteType"tJ", s::Index) + return op!(Op, alias(on), st, s) +end + +function op!(Op::ITensor, ::OpName"Cdn", ::SiteType"tJ", s::Index) + return Op[s' => 1, s => 3] = 1.0 +end +function op!(Op::ITensor, on::OpName"c↓", st::SiteType"tJ", s::Index) + return op!(Op, alias(on), st, s) +end + +function op!(Op::ITensor, ::OpName"Cdagdn", ::SiteType"tJ", s::Index) + return Op[s' => 3, s => 1] = 1.0 +end +function op!(Op::ITensor, on::OpName"c†↓", st::SiteType"tJ", s::Index) + return op!(Op, alias(on), st, s) +end + +function op!(Op::ITensor, ::OpName"Aup", ::SiteType"tJ", s::Index) + return Op[s' => 1, s => 2] = 1.0 +end +function op!(Op::ITensor, on::OpName"a↑", st::SiteType"tJ", s::Index) + return op!(Op, alias(on), st, s) +end + +function op!(Op::ITensor, ::OpName"Adagup", ::SiteType"tJ", s::Index) + return Op[s' => 2, s => 1] = 1.0 +end +function op!(Op::ITensor, on::OpName"a†↑", st::SiteType"tJ", s::Index) + return op!(Op, alias(on), st, s) +end + +function op!(Op::ITensor, ::OpName"Adn", ::SiteType"tJ", s::Index) + return Op[s' => 1, s => 3] = 1.0 +end +function op!(Op::ITensor, on::OpName"a↓", st::SiteType"tJ", s::Index) + return op!(Op, alias(on), st, s) +end + +function op!(Op::ITensor, ::OpName"Adagdn", ::SiteType"tJ", s::Index) + return Op[s' => 3, s => 1] = 1.0 +end +function op!(Op::ITensor, on::OpName"a†↓", st::SiteType"tJ", s::Index) + return op!(Op, alias(on), st, s) +end + +function op!(Op::ITensor, ::OpName"F", ::SiteType"tJ", s::Index) + Op[s' => 1, s => 1] = +1.0 + Op[s' => 2, s => 2] = -1.0 + return Op[s' => 3, s => 3] = -1.0 +end + +function op!(Op::ITensor, ::OpName"Fup", ::SiteType"tJ", s::Index) + Op[s' => 1, s => 1] = +1.0 + Op[s' => 2, s => 2] = -1.0 + return Op[s' => 3, s => 3] = +1.0 +end +function op!(Op::ITensor, on::OpName"F↑", st::SiteType"tJ", s::Index) + return op!(Op, alias(on), st, s) +end + +function op!(Op::ITensor, ::OpName"Fdn", ::SiteType"tJ", s::Index) + Op[s' => 1, s => 1] = +1.0 + Op[s' => 2, s => 2] = +1.0 + return Op[s' => 3, s => 3] = -1.0 +end +function op!(Op::ITensor, on::OpName"F↓", st::SiteType"tJ", s::Index) + return op!(Op, alias(on), st, s) +end + +function op!(Op::ITensor, ::OpName"Sz", ::SiteType"tJ", s::Index) + Op[s' => 2, s => 2] = +0.5 + return Op[s' => 3, s => 3] = -0.5 +end + +function op!(Op::ITensor, ::OpName"Sᶻ", st::SiteType"tJ", s::Index) + return op!(Op, OpName("Sz"), st, s) +end + +function op!(Op::ITensor, ::OpName"Sx", ::SiteType"tJ", s::Index) + Op[s' => 2, s => 3] = 0.5 + return Op[s' => 3, s => 2] = 0.5 +end + +function op!(Op::ITensor, ::OpName"Sˣ", st::SiteType"tJ", s::Index) + return op!(Op, OpName("Sx"), st, s) +end + +function op!(Op::ITensor, ::OpName"S+", ::SiteType"tJ", s::Index) + return Op[s' => 2, s => 3] = 1.0 +end + +function op!(Op::ITensor, ::OpName"S⁺", st::SiteType"tJ", s::Index) + return op!(Op, OpName("S+"), st, s) +end +function op!(Op::ITensor, ::OpName"Sp", st::SiteType"tJ", s::Index) + return op!(Op, OpName("S+"), st, s) +end +function op!(Op::ITensor, ::OpName"Splus", st::SiteType"tJ", s::Index) + return op!(Op, OpName("S+"), st, s) +end + +function op!(Op::ITensor, ::OpName"S-", ::SiteType"tJ", s::Index) + return Op[s' => 3, s => 2] = 1.0 +end + +function op!(Op::ITensor, ::OpName"S⁻", st::SiteType"tJ", s::Index) + return op!(Op, OpName("S-"), st, s) +end +function op!(Op::ITensor, ::OpName"Sm", st::SiteType"tJ", s::Index) + return op!(Op, OpName("S-"), st, s) +end +function op!(Op::ITensor, ::OpName"Sminus", st::SiteType"tJ", s::Index) + return op!(Op, OpName("S-"), st, s) +end + +has_fermion_string(::OpName"Cup", ::SiteType"tJ") = true +function has_fermion_string(on::OpName"c↑", st::SiteType"tJ") + return has_fermion_string(alias(on), st) +end +has_fermion_string(::OpName"Cdagup", ::SiteType"tJ") = true +function has_fermion_string(on::OpName"c†↑", st::SiteType"tJ") + return has_fermion_string(alias(on), st) +end +has_fermion_string(::OpName"Cdn", ::SiteType"tJ") = true +function has_fermion_string(on::OpName"c↓", st::SiteType"tJ") + return has_fermion_string(alias(on), st) +end +has_fermion_string(::OpName"Cdagdn", ::SiteType"tJ") = true +function has_fermion_string(on::OpName"c†↓", st::SiteType"tJ") + return has_fermion_string(alias(on), st) +end