diff --git a/.github/ISSUE_TEMPLATE/BUG_REPORT.md b/.github/ISSUE_TEMPLATE/BUG_REPORT.md index 45c1fb7..fd82368 100644 --- a/.github/ISSUE_TEMPLATE/BUG_REPORT.md +++ b/.github/ISSUE_TEMPLATE/BUG_REPORT.md @@ -1,6 +1,6 @@ --- -name: ITensorQuantumOperatorDefinitions.jl bug report -about: Create a bug report to help us improve ITensorQuantumOperatorDefinitions.jl +name: QuantumOperatorDefinitions.jl bug report +about: Create a bug report to help us improve QuantumOperatorDefinitions.jl title: "[BUG] YOUR SHORT DESCRIPTION OF THE BUG HERE" labels: ["bug"] assignees: '' @@ -55,8 +55,8 @@ If you provided a minimal code that demonstrates the bug or unexpected behavior, julia> versioninfo() [YOUR OUTPUT HERE] ``` - - Output from `using Pkg; Pkg.status("ITensorQuantumOperatorDefinitions")`: + - Output from `using Pkg; Pkg.status("QuantumOperatorDefinitions")`: ```julia -julia> using Pkg; Pkg.status("ITensorQuantumOperatorDefinitions") +julia> using Pkg; Pkg.status("QuantumOperatorDefinitions") [YOUR OUTPUT HERE] ``` diff --git a/.github/ISSUE_TEMPLATE/FEATURE_REQUEST.md b/.github/ISSUE_TEMPLATE/FEATURE_REQUEST.md index da7f0ca..515cf26 100644 --- a/.github/ISSUE_TEMPLATE/FEATURE_REQUEST.md +++ b/.github/ISSUE_TEMPLATE/FEATURE_REQUEST.md @@ -1,6 +1,6 @@ --- -name: ITensorQuantumOperatorDefinitions.jl feature request -about: Suggest an idea for ITensorQuantumOperatorDefinitions.jl +name: QuantumOperatorDefinitions.jl feature request +about: Suggest an idea for QuantumOperatorDefinitions.jl title: "[ENHANCEMENT] YOUR SHORT DESCRIPTION OF THE FEATURE REQUEST HERE" labels: ["enhancement"] assignees: '' diff --git a/.github/PULL_REQUEST_TEMPLATE.md b/.github/PULL_REQUEST_TEMPLATE.md index e8e9fb7..23edeee 100644 --- a/.github/PULL_REQUEST_TEMPLATE.md +++ b/.github/PULL_REQUEST_TEMPLATE.md @@ -33,7 +33,7 @@ Please give a summary of the tests that you added to verify your changes. # Checklist: -- [ ] My code follows the style guidelines of this project. Please run `using JuliaFormatter; format(".")` in the base directory of the repository (`~/.julia/dev/ITensorQuantumOperatorDefinitions`) to format your code according to our style guidelines. +- [ ] My code follows the style guidelines of this project. Please run `using JuliaFormatter; format(".")` in the base directory of the repository (`~/.julia/dev/QuantumOperatorDefinitions`) to format your code according to our style guidelines. - [ ] I have performed a self-review of my own code. - [ ] I have commented my code, particularly in hard-to-understand areas. - [ ] I have added tests that verify the behavior of the changes I made. diff --git a/Project.toml b/Project.toml index 8ff19e7..17f8383 100644 --- a/Project.toml +++ b/Project.toml @@ -1,15 +1,11 @@ -name = "ITensorQuantumOperatorDefinitions" -uuid = "fd9b415b-e710-4e2a-b407-cba023081494" +name = "QuantumOperatorDefinitions" +uuid = "826dd319-6fd5-459a-a990-3a4f214664bf" authors = ["ITensor developers and contributors"] -version = "0.1.2" +version = "0.1.0" [deps] -ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" -ITensorBase = "4795dd04-0d67-49bb-8f44-b89c448a1dc7" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" [compat] -ChainRulesCore = "1.25.1" -ITensorBase = "0.1.8" LinearAlgebra = "1.10" julia = "1.10" diff --git a/README.md b/README.md index 88b3f11..978d575 100644 --- a/README.md +++ b/README.md @@ -1,9 +1,9 @@ -# ITensorQuantumOperatorDefinitions.jl +# QuantumOperatorDefinitions.jl -[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://ITensor.github.io/ITensorQuantumOperatorDefinitions.jl/stable/) -[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://ITensor.github.io/ITensorQuantumOperatorDefinitions.jl/dev/) -[![Build Status](https://github.com/ITensor/ITensorQuantumOperatorDefinitions.jl/actions/workflows/Tests.yml/badge.svg?branch=main)](https://github.com/ITensor/ITensorQuantumOperatorDefinitions.jl/actions/workflows/Tests.yml?query=branch%3Amain) -[![Coverage](https://codecov.io/gh/ITensor/ITensorQuantumOperatorDefinitions.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/ITensor/ITensorQuantumOperatorDefinitions.jl) +[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://ITensor.github.io/QuantumOperatorDefinitions.jl/stable/) +[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://ITensor.github.io/QuantumOperatorDefinitions.jl/dev/) +[![Build Status](https://github.com/ITensor/QuantumOperatorDefinitions.jl/actions/workflows/Tests.yml/badge.svg?branch=main)](https://github.com/ITensor/QuantumOperatorDefinitions.jl/actions/workflows/Tests.yml?query=branch%3Amain) +[![Coverage](https://codecov.io/gh/ITensor/QuantumOperatorDefinitions.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/ITensor/QuantumOperatorDefinitions.jl) [![Code Style: Blue](https://img.shields.io/badge/code%20style-blue-4495d1.svg)](https://github.com/invenia/BlueStyle) [![Aqua](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl) @@ -26,50 +26,56 @@ if you want to use SSH credentials, which can make it so you don't have to enter Then, the package can be added as usual through the package manager: ```julia -julia> Pkg.add("ITensorQuantumOperatorDefinitions") +julia> Pkg.add("QuantumOperatorDefinitions") ``` ## Examples ````julia -using ITensorBase: ITensor, Index, tags -using ITensorQuantumOperatorDefinitions: - OpName, SiteType, StateName, ValName, op, siteind, siteinds, state, val +using QuantumOperatorDefinitions: OpName, SiteType, StateName, ⊗, controlled, op, state +using LinearAlgebra: Diagonal +using SparseArrays: SparseMatrixCSC, SparseVector using Test: @test -```` -States and operators as Julia arrays +@test state("0") == [1, 0] +@test state("1") == [0, 1] -````julia -@test val(ValName("Up"), SiteType("S=1/2")) == 1 -@test val(ValName("Dn"), SiteType("S=1/2")) == 2 -@test state(StateName("Up"), SiteType("S=1/2")) == [1, 0] -@test state(StateName("Dn"), SiteType("S=1/2")) == [0, 1] -@test op(OpName("X"), SiteType("S=1/2")) == [0 1; 1 0] -@test op(OpName("Z"), SiteType("S=1/2")) == [1 0; 0 -1] -```` +@test state(Float32, "0") == [1, 0] +@test eltype(state(Float32, "1")) === Float32 -Index constructors +@test Vector(StateName("0")) == [1, 0] +@test Vector(StateName("1")) == [0, 1] -````julia -@test siteind("S=1/2") isa Index -@test Int(length(siteind("S=1/2"))) == 2 -@test "S=1/2" in tags(siteind("S=1/2")) -@test siteinds("S=1/2", 4) isa Vector{<:Index} -@test length(siteinds("S=1/2", 4)) == 4 -@test all(s -> "S=1/2" in tags(s), siteinds("S=1/2", 4)) -```` +@test Vector{Float32}(StateName("0")) == [1, 0] +@test eltype(Vector{Float32}(StateName("0"))) === Float32 -States and operators as ITensors +@test state(SparseVector, "0") == [1, 0] +@test state(SparseVector, "0") isa SparseVector -````julia -s = Index(2, "S=1/2") -@test val(s, "Up") == 1 -@test val(s, "Dn") == 2 -@test state("Up", s) == ITensor([1, 0], (s,)) -@test state("Dn", s) == ITensor([0, 1], (s,)) -@test op("X", s) == ITensor([0 1; 1 0], (s', s)) -@test op("Z", s) == ITensor([1 0; 0 -1], (s', s)) +@test state("0", 3) == [1, 0, 0] +@test state("1", 3) == [0, 1, 0] +@test state("2", 3) == [0, 0, 1] + +@test Vector(StateName("0"), 3) == [1, 0, 0] +@test Vector(StateName("1"), 3) == [0, 1, 0] +@test Vector(StateName("2"), 3) == [0, 0, 1] + +@test op("X") == [0 1; 1 0] +@test op("Y") == [0 -im; im 0] +@test op("Z") == [1 0; 0 -1] + +@test op("Z") isa Diagonal + +@test op(Float32, "X") == [0 1; 1 0] +@test eltype(op(Float32, "X")) === Float32 +@test op(SparseMatrixCSC, "X") == [0 1; 1 0] +@test op(SparseMatrixCSC, "X") isa SparseMatrixCSC + +@test Matrix(OpName("X")) == [0 1; 1 0] +@test Matrix(OpName("Y")) == [0 -im; im 0] +@test Matrix(OpName("Z")) == [1 0; 0 -1] + +@test Matrix(OpName("Rx"; θ=π / 3)) ≈ [sin(π / 3) -cos(π / 3)*im; -cos(π / 3)*im sin(π / 3)] ```` --- diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl index edbe45d..a04a019 100644 --- a/benchmark/benchmarks.jl +++ b/benchmark/benchmarks.jl @@ -1,4 +1,4 @@ -using ITensorQuantumOperatorDefinitions +using QuantumOperatorDefinitions using BenchmarkTools SUITE = BenchmarkGroup() diff --git a/docs/Project.toml b/docs/Project.toml index cd64c85..41a0d99 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,6 +1,6 @@ [deps] Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -ITensorBase = "4795dd04-0d67-49bb-8f44-b89c448a1dc7" -ITensorQuantumOperatorDefinitions = "fd9b415b-e710-4e2a-b407-cba023081494" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" +QuantumOperatorDefinitions = "826dd319-6fd5-459a-a990-3a4f214664bf" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/docs/make.jl b/docs/make.jl index 26a28f4..de73136 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,21 +1,21 @@ -using ITensorQuantumOperatorDefinitions: ITensorQuantumOperatorDefinitions +using QuantumOperatorDefinitions: QuantumOperatorDefinitions using Documenter: Documenter, DocMeta, deploydocs, makedocs DocMeta.setdocmeta!( - ITensorQuantumOperatorDefinitions, + QuantumOperatorDefinitions, :DocTestSetup, - :(using ITensorQuantumOperatorDefinitions); + :(using QuantumOperatorDefinitions); recursive=true, ) include("make_index.jl") makedocs(; - modules=[ITensorQuantumOperatorDefinitions], + modules=[QuantumOperatorDefinitions], authors="ITensor developers and contributors", - sitename="ITensorQuantumOperatorDefinitions.jl", + sitename="QuantumOperatorDefinitions.jl", format=Documenter.HTML(; - canonical="https://ITensor.github.io/ITensorQuantumOperatorDefinitions.jl", + canonical="https://ITensor.github.io/QuantumOperatorDefinitions.jl", edit_link="main", assets=String[], ), @@ -23,7 +23,7 @@ makedocs(; ) deploydocs(; - repo="github.com/ITensor/ITensorQuantumOperatorDefinitions.jl", + repo="github.com/ITensor/QuantumOperatorDefinitions.jl", devbranch="main", push_preview=true, ) diff --git a/docs/make_index.jl b/docs/make_index.jl index 1f1eb5a..a1d332e 100644 --- a/docs/make_index.jl +++ b/docs/make_index.jl @@ -1,9 +1,9 @@ using Literate: Literate -using ITensorQuantumOperatorDefinitions: ITensorQuantumOperatorDefinitions +using QuantumOperatorDefinitions: QuantumOperatorDefinitions Literate.markdown( - joinpath(pkgdir(ITensorQuantumOperatorDefinitions), "examples", "README.jl"), - joinpath(pkgdir(ITensorQuantumOperatorDefinitions), "docs", "src"); + joinpath(pkgdir(QuantumOperatorDefinitions), "examples", "README.jl"), + joinpath(pkgdir(QuantumOperatorDefinitions), "docs", "src"); flavor=Literate.DocumenterFlavor(), name="index", ) diff --git a/docs/make_readme.jl b/docs/make_readme.jl index 82b10e9..6851a5c 100644 --- a/docs/make_readme.jl +++ b/docs/make_readme.jl @@ -1,9 +1,9 @@ using Literate: Literate -using ITensorQuantumOperatorDefinitions: ITensorQuantumOperatorDefinitions +using QuantumOperatorDefinitions: QuantumOperatorDefinitions Literate.markdown( - joinpath(pkgdir(ITensorQuantumOperatorDefinitions), "examples", "README.jl"), - joinpath(pkgdir(ITensorQuantumOperatorDefinitions)); + joinpath(pkgdir(QuantumOperatorDefinitions), "examples", "README.jl"), + joinpath(pkgdir(QuantumOperatorDefinitions)); flavor=Literate.CommonMarkFlavor(), name="README", ) diff --git a/docs/src/reference.md b/docs/src/reference.md index 36b8160..1c25104 100644 --- a/docs/src/reference.md +++ b/docs/src/reference.md @@ -1,5 +1,5 @@ # Reference ```@autodocs -Modules = [ITensorQuantumOperatorDefinitions] +Modules = [QuantumOperatorDefinitions] ``` diff --git a/examples/Project.toml b/examples/Project.toml index a685760..d78168c 100644 --- a/examples/Project.toml +++ b/examples/Project.toml @@ -1,3 +1,3 @@ [deps] -ITensorBase = "4795dd04-0d67-49bb-8f44-b89c448a1dc7" -ITensorQuantumOperatorDefinitions = "fd9b415b-e710-4e2a-b407-cba023081494" +QuantumOperatorDefinitions = "826dd319-6fd5-459a-a990-3a4f214664bf" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" diff --git a/examples/README.jl b/examples/README.jl index ee91a93..5b8b5b1 100644 --- a/examples/README.jl +++ b/examples/README.jl @@ -1,9 +1,9 @@ -# # ITensorQuantumOperatorDefinitions.jl +# # QuantumOperatorDefinitions.jl # -# [![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://ITensor.github.io/ITensorQuantumOperatorDefinitions.jl/stable/) -# [![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://ITensor.github.io/ITensorQuantumOperatorDefinitions.jl/dev/) -# [![Build Status](https://github.com/ITensor/ITensorQuantumOperatorDefinitions.jl/actions/workflows/Tests.yml/badge.svg?branch=main)](https://github.com/ITensor/ITensorQuantumOperatorDefinitions.jl/actions/workflows/Tests.yml?query=branch%3Amain) -# [![Coverage](https://codecov.io/gh/ITensor/ITensorQuantumOperatorDefinitions.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/ITensor/ITensorQuantumOperatorDefinitions.jl) +# [![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://ITensor.github.io/QuantumOperatorDefinitions.jl/stable/) +# [![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://ITensor.github.io/QuantumOperatorDefinitions.jl/dev/) +# [![Build Status](https://github.com/ITensor/QuantumOperatorDefinitions.jl/actions/workflows/Tests.yml/badge.svg?branch=main)](https://github.com/ITensor/QuantumOperatorDefinitions.jl/actions/workflows/Tests.yml?query=branch%3Amain) +# [![Coverage](https://codecov.io/gh/ITensor/QuantumOperatorDefinitions.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/ITensor/QuantumOperatorDefinitions.jl) # [![Code Style: Blue](https://img.shields.io/badge/code%20style-blue-4495d1.svg)](https://github.com/invenia/BlueStyle) # [![Aqua](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl) @@ -31,38 +31,53 @@ julia> Pkg.Registry.add(url="git@github.com:ITensor/ITensorRegistry.git") #= ```julia -julia> Pkg.add("ITensorQuantumOperatorDefinitions") +julia> Pkg.add("QuantumOperatorDefinitions") ``` =# # ## Examples -using ITensorBase: ITensor, Index, tags -using ITensorQuantumOperatorDefinitions: - OpName, SiteType, StateName, ValName, op, siteind, siteinds, state, val +using QuantumOperatorDefinitions: OpName, SiteType, StateName, ⊗, controlled, op, state +using LinearAlgebra: Diagonal +using SparseArrays: SparseMatrixCSC, SparseVector using Test: @test -# States and operators as Julia arrays -@test val(ValName("Up"), SiteType("S=1/2")) == 1 -@test val(ValName("Dn"), SiteType("S=1/2")) == 2 -@test state(StateName("Up"), SiteType("S=1/2")) == [1, 0] -@test state(StateName("Dn"), SiteType("S=1/2")) == [0, 1] -@test op(OpName("X"), SiteType("S=1/2")) == [0 1; 1 0] -@test op(OpName("Z"), SiteType("S=1/2")) == [1 0; 0 -1] - -# Index constructors -@test siteind("S=1/2") isa Index -@test Int(length(siteind("S=1/2"))) == 2 -@test "S=1/2" in tags(siteind("S=1/2")) -@test siteinds("S=1/2", 4) isa Vector{<:Index} -@test length(siteinds("S=1/2", 4)) == 4 -@test all(s -> "S=1/2" in tags(s), siteinds("S=1/2", 4)) - -# States and operators as ITensors -s = Index(2, "S=1/2") -@test val(s, "Up") == 1 -@test val(s, "Dn") == 2 -@test state("Up", s) == ITensor([1, 0], (s,)) -@test state("Dn", s) == ITensor([0, 1], (s,)) -@test op("X", s) == ITensor([0 1; 1 0], (s', s)) -@test op("Z", s) == ITensor([1 0; 0 -1], (s', s)) +@test state("0") == [1, 0] +@test state("1") == [0, 1] + +@test state(Float32, "0") == [1, 0] +@test eltype(state(Float32, "1")) === Float32 + +@test Vector(StateName("0")) == [1, 0] +@test Vector(StateName("1")) == [0, 1] + +@test Vector{Float32}(StateName("0")) == [1, 0] +@test eltype(Vector{Float32}(StateName("0"))) === Float32 + +@test state(SparseVector, "0") == [1, 0] +@test state(SparseVector, "0") isa SparseVector + +@test state("0", 3) == [1, 0, 0] +@test state("1", 3) == [0, 1, 0] +@test state("2", 3) == [0, 0, 1] + +@test Vector(StateName("0"), 3) == [1, 0, 0] +@test Vector(StateName("1"), 3) == [0, 1, 0] +@test Vector(StateName("2"), 3) == [0, 0, 1] + +@test op("X") == [0 1; 1 0] +@test op("Y") == [0 -im; im 0] +@test op("Z") == [1 0; 0 -1] + +@test op("Z") isa Diagonal + +@test op(Float32, "X") == [0 1; 1 0] +@test eltype(op(Float32, "X")) === Float32 +@test op(SparseMatrixCSC, "X") == [0 1; 1 0] +@test op(SparseMatrixCSC, "X") isa SparseMatrixCSC + +@test Matrix(OpName("X")) == [0 1; 1 0] +@test Matrix(OpName("Y")) == [0 -im; im 0] +@test Matrix(OpName("Z")) == [1 0; 0 -1] + +@test Matrix(OpName("Rx"; θ=π / 3)) ≈ [sin(π / 3) -cos(π / 3)*im; -cos(π / 3)*im sin(π / 3)] diff --git a/src/ITensorQuantumOperatorDefinitions.jl b/src/ITensorQuantumOperatorDefinitions.jl deleted file mode 100644 index 11c65a8..0000000 --- a/src/ITensorQuantumOperatorDefinitions.jl +++ /dev/null @@ -1,29 +0,0 @@ -module ITensorQuantumOperatorDefinitions - -include("sitetype.jl") -include("space.jl") -include("val.jl") -include("state.jl") -include("op.jl") -include("has_fermion_string.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") - -include("ITensorQuantumOperatorDefinitionsChainRulesCoreExt.jl") - -include("itensor/siteinds.jl") -include("itensor/val.jl") -include("itensor/state.jl") -include("itensor/op.jl") -include("itensor/has_fermion_string.jl") - -end diff --git a/src/ITensorQuantumOperatorDefinitionsChainRulesCoreExt.jl b/src/ITensorQuantumOperatorDefinitionsChainRulesCoreExt.jl deleted file mode 100644 index c420db9..0000000 --- a/src/ITensorQuantumOperatorDefinitionsChainRulesCoreExt.jl +++ /dev/null @@ -1,10 +0,0 @@ -module ITensorQuantumOperatorDefinitionsChainRulesCoreExt - -using ChainRulesCore: @non_differentiable -using ..ITensorQuantumOperatorDefinitions: 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/QuantumOperatorDefinitions.jl b/src/QuantumOperatorDefinitions.jl new file mode 100644 index 0000000..f653b6a --- /dev/null +++ b/src/QuantumOperatorDefinitions.jl @@ -0,0 +1,14 @@ +module QuantumOperatorDefinitions + +include("sitetype.jl") +include("state.jl") +include("op.jl") +include("has_fermion_string.jl") +include("sitetypes/qubit.jl") +include("sitetypes/spinone.jl") +include("sitetypes/fermion.jl") +include("sitetypes/electron.jl") +include("sitetypes/tj.jl") +include("sitetypes/qudit.jl") + +end diff --git a/src/itensor/has_fermion_string.jl b/src/itensor/has_fermion_string.jl deleted file mode 100644 index 22330b2..0000000 --- a/src/itensor/has_fermion_string.jl +++ /dev/null @@ -1,25 +0,0 @@ -using ITensorBase: Index - -has_fermion_string(operator::AbstractArray{<:Number}, s::Index; kwargs...)::Bool = false - -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/itensor/op.jl b/src/itensor/op.jl deleted file mode 100644 index 40a8074..0000000 --- a/src/itensor/op.jl +++ /dev/null @@ -1,138 +0,0 @@ -using ITensorBase: Index, ITensor, dag, prime, tags - -op(::OpName, ::SiteType, ::Index...; kwargs...) = nothing -function op( - ::OpName, ::SiteType, ::SiteType, sitetypes_inds::Union{SiteType,Index}...; kwargs... -) - return nothing -end - -_sitetypes(i::Index) = _sitetypes(tags(i)) - -function commontags(i::Index...) - return union(tags.(i)...) -end - -op(X::AbstractArray, s::Index...) = ITensor(X, (prime.(s)..., dag.(s)...)) - -# TODO: Delete these. -op(opname, s::Vector{<:Index}; kwargs...) = op(opname, s...; kwargs...) -op(s::Vector{<:Index}, opname; kwargs...) = op(opname, s...; kwargs...) - -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) - for st in common_stypes - op_mat = op(opn, st; 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 - end - return throw( - ArgumentError( - "Overload of \"op\" or \"op!\" functions not found for operator name \"$name\" and Index tags: $(tags.(s)).", - ), - ) -end - -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...) - -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] diff --git a/src/itensor/siteinds.jl b/src/itensor/siteinds.jl deleted file mode 100644 index be81f46..0000000 --- a/src/itensor/siteinds.jl +++ /dev/null @@ -1,78 +0,0 @@ -using ITensorBase: Index, addtags - -function siteind(st::SiteType; addtags="", kwargs...) - sp = space(st; kwargs...) - isnothing(sp) && return nothing - return Index(sp, "Site, $(value(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, $(value(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(::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 diff --git a/src/itensor/state.jl b/src/itensor/state.jl deleted file mode 100644 index 52b355f..0000000 --- a/src/itensor/state.jl +++ /dev/null @@ -1,24 +0,0 @@ -using ITensorBase: ITensor, Index, onehot - -state(::StateName, ::SiteType, ::Index; kwargs...) = nothing - -# Syntax `state("Up", Index(2, "S=1/2"))` -state(sn::String, i::Index; kwargs...) = state(i, sn; kwargs...) - -function state(s::Index, name::AbstractString; kwargs...)::ITensor - stypes = _sitetypes(s) - sname = StateName(name) - for st in stypes - v = state(sname, st; kwargs...) - !isnothing(v) && return ITensor(v, (s,)) - end - return throw( - ArgumentError( - "Overload of \"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...) diff --git a/src/itensor/val.jl b/src/itensor/val.jl deleted file mode 100644 index 2bc938c..0000000 --- a/src/itensor/val.jl +++ /dev/null @@ -1,20 +0,0 @@ -using ITensorBase: Index - -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) diff --git a/src/op.jl b/src/op.jl index 529be1b..b909026 100644 --- a/src/op.jl +++ b/src/op.jl @@ -1,24 +1,480 @@ -# TODO: Add `params<:NamedTuple` field. -struct OpName{Name} end -OpName(s::AbstractString) = OpName{Symbol(s)}() -OpName(s::Symbol) = OpName{s}() +struct OpName{Name,Params} + params::Params +end +params(n::OpName) = getfield(n, :params) + +Base.getproperty(n::OpName, name::Symbol) = getfield(params(n), name) + +OpName{N}(params) where {N} = OpName{N,typeof(params)}(params) +OpName{N}(; kwargs...) where {N} = OpName{N}((; kwargs...)) + +# This compiles operator expressions, such as: +# ```julia +# opexpr("X + Y") == OpName("X") + OpName("Y") +# opexpr("Ry{θ=π/2}") == OpName("Ry"; θ=π/2) +# ``` +function opexpr(n::String) + return opexpr(Meta.parse(n)) +end +opexpr(n::Number) = n +function opexpr(n::Symbol) + n === :im && return im + n === :π && return π + return OpName{n}() +end +function opexpr(ex::Expr) + if Meta.isexpr(ex, :call) + return eval(ex.args[1])(opexpr.(ex.args[2:end])...) + end + if Meta.isexpr(ex, :curly) + # Syntax for parametrized gates, i.e. + # `opexpr("Ry{θ=π/2}")`. + params = ex.args[2:end] + kwargs = Dict( + map(params) do param + @assert Meta.isexpr(param, :(=)) + return param.args[1] => eval(param.args[2]) + end, + ) + return OpName{ex.args[1]}(; kwargs...) + end + return error("Can't parse expression $ex.") +end + +OpName(s::AbstractString; kwargs...) = OpName{Symbol(s)}(; kwargs...) +OpName(s::Symbol; kwargs...) = OpName{s}(; kwargs...) name(::OpName{N}) where {N} = N macro OpName_str(s) return OpName{Symbol(s)} end -# Default implementations of op -op(::OpName; kwargs...) = nothing -op(::OpName, ::SiteType; kwargs...) = nothing +function op_alias_expr(name1, name2, pars...) + return :(function alias(n::OpName{Symbol($name1)}) + return OpName{Symbol($name2)}(; params(n)..., $(esc.(pars)...)) + end) +end +macro op_alias(name1, name2, pars...) + return op_alias_expr(name1, name2, pars...) +end + +alias(n::OpName) = n +function op_convert( + arrtype::Type{<:AbstractArray{<:Any,N}}, + domain_size::Tuple{Vararg{Integer}}, + a::AbstractArray{<:Any,N}, +) where {N} + # TODO: Check the dimensions. + return convert(arrtype, a) +end +function op_convert( + arrtype::Type{<:AbstractArray}, domain_size::Tuple{Vararg{Integer}}, a::AbstractArray +) + # TODO: Check the dimensions. + return convert(arrtype, a) +end +function op_convert( + arrtype::Type{<:AbstractArray{<:Any,N}}, + domain_size::Tuple{Vararg{Integer}}, + a::AbstractArray, +) where {N} + size = (domain_size..., domain_size...) + @assert length(size) == N + return convert(arrtype, reshape(a, size)) +end +function (arrtype::Type{<:AbstractArray})(n::OpName, ts::Tuple{Vararg{SiteType}}) + return op_convert(arrtype, length.(ts), AbstractArray(n, ts)) +end +function (arrtype::Type{<:AbstractArray})(n::OpName, domain_size::Tuple{Vararg{Integer}}) + return op_convert(arrtype, domain_size, AbstractArray(n, Int.(domain_size))) +end +function (arrtype::Type{<:AbstractArray})(n::OpName, domain_size::Integer...) + return arrtype(n, domain_size) +end +(arrtype::Type{<:AbstractArray})(n::OpName, ts::SiteType...) = arrtype(n, ts) +Base.AbstractArray(n::OpName, ts::SiteType...) = AbstractArray(n, ts) +function Base.AbstractArray(n::OpName, ts::Tuple{Vararg{SiteType}}) + n′ = alias(n) + ts′ = alias.(ts) + if n′ == n && ts′ == ts + return AbstractArray(n′, length.(ts′)) + end + return AbstractArray(n′, ts′) +end +function Base.AbstractArray(n::OpName, domain_size::Tuple{Vararg{Int}}) + n′ = alias(n) + if n′ == n + error("Not implemented.") + end + return AbstractArray(n′, domain_size) +end + +# TODO: Decide on this. +function Base.AbstractArray(n::OpName) + return AbstractArray(n, ntuple(Returns(default_sitetype()), nsites(n))) +end +function (arrtype::Type{<:AbstractArray})(n::OpName) + return arrtype(n, ntuple(Returns(default_sitetype()), nsites(n))) +end + +function op(arrtype::Type{<:AbstractArray}, n::String, domain...; kwargs...) + return arrtype(OpName(n; kwargs...), domain...) +end +function op(elt::Type{<:Number}, n::String, domain...; kwargs...) + return op(AbstractArray{elt}, n, domain...; kwargs...) +end +function op(n::String, domain...; kwargs...) + return op(AbstractArray, n, domain...; kwargs...) +end + +# `Int(ndims // 2)`, i.e. `ndims_domain`/`ndims_codomain`. +function nsites(n::Union{StateName,OpName}) + n′ = alias(n) + if n′ == n + # Default value, assume 1-site operator/state. + return 1 + end + return nsites(n′) +end + +## TODO: Delete. +## # Default implementations of op +## op(::OpName; kwargs...) = nothing +## op(::OpName, ::SiteType; kwargs...) = nothing function _sitetypes(ts::Set) return collect(SiteType, SiteType.(ts)) end -op(name::AbstractString; kwargs...) = error("Must input indices when creating an `op`.") +## TODO: Delete. +## op(name::AbstractString; kwargs...) = error("Must input indices when creating an `op`.") # 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...) +## TODO: Bring this back? +## op(opname::AbstractString, t::SiteType; kwargs...) = op(OpName(opname), t; kwargs...) + +# TODO: Bring this back? +# op(f::Function, args...; kwargs...) = f(op(args...; kwargs...)) + +using LinearAlgebra: Diagonal +function Base.AbstractArray(::OpName"Id", domain_size::Tuple{Int}) + return Diagonal(trues(only(domain_size))) +end +function Base.AbstractArray(n::OpName"Id", domain_size::Tuple{Int,Vararg{Int}}) + return Base.AbstractArray(kron(ntuple(Returns(n), length(domain_size))...), domain_size) +end +@op_alias "I" "Id" +@op_alias "σ0" "Id" +@op_alias "σ⁰" "Id" +@op_alias "σ₀" "Id" +# TODO: Is this a good definition? +@op_alias "F" "Id" + +# TODO: Define as `::Tuple{Int}`. +function Base.AbstractArray(n::OpName"a†", domain_size::Tuple{Int}) + d = only(domain_size) + a = zeros(d, d) + for k in 1:(d - 1) + a[k + 1, k] = √k + end + return a +end +@op_alias "Adag" "a†" +@op_alias "adag" "a†" +alias(::OpName"a") = OpName"a†"()' +@op_alias "A" "a" + +alias(::OpName"n") = OpName"a†"() * OpName"a"() +@op_alias "N" "n" + +# `cis(x) = exp(im * x)` +alias(n::OpName"Phase") = cis(n.θ * OpName"n"()) +@op_alias "PHASE" "Phase" +@op_alias "P" "Phase" +@op_alias "π/8" "Phase" θ = π / 4 +@op_alias "T" "π/8" +@op_alias "S" "Phase" θ = π / 2 + +alias(::OpName"aa") = OpName("a") ⊗ OpName("a") +alias(::OpName"a†a") = OpName("a†") ⊗ OpName("a") +alias(::OpName"aa†") = OpName("a") ⊗ OpName("a†") +alias(::OpName"a†a†") = OpName("a†") ⊗ OpName("a†") + +δ(x, y) = (x == y) ? 1 : 0 + +# See: +# https://en.wikipedia.org/wiki/Spin_(physics)#Higher_spins +# https://en.wikipedia.org/wiki/Pauli_matrices +# https://en.wikipedia.org/wiki/Generalizations_of_Pauli_matrices +# https://en.wikipedia.org/wiki/Generalized_Clifford_algebra +# https://github.com/QuantumKitHub/MPSKitModels.jl/blob/v0.4.0/src/operators/spinoperators.jl +function Base.AbstractArray(n::OpName"σ⁺", domain_size::Tuple{Int}) + d = only(domain_size) + s = (d - 1) / 2 + return [2 * δ(i + 1, j) * √((s + 1) * (i + j - 1) - i * j) for i in 1:d, j in 1:d] +end +alias(::OpName"S⁺") = OpName("σ⁺") / 2 +@op_alias "S+" "S⁺" +@op_alias "Splus" "S+" +@op_alias "Sp" "S+" + +alias(::OpName"σ⁻") = OpName"σ⁺"()' +alias(::OpName"S⁻") = OpName("σ⁻") / 2 +@op_alias "S-" "S⁻" +@op_alias "Sminus" "S-" +@op_alias "Sm" "S-" + +alias(::OpName"X") = (OpName"σ⁺"() + OpName"σ⁻"()) / 2 +@op_alias "σx" "X" +@op_alias "σˣ" "X" +@op_alias "σₓ" "X" +@op_alias "σ1" "X" +@op_alias "σ¹" "X" +@op_alias "σ₁" "X" +@op_alias "σy" "Y" +alias(::OpName"iX") = im * OpName"X"() +alias(::OpName"√X") = √OpName"X"() +@op_alias "√NOT" "√X" +alias(n::OpName"Sx") = OpName"X"() / 2 +@op_alias "Sˣ" "Sx" +@op_alias "Sₓ" "Sx" +alias(::OpName"Sx2") = OpName"Sx"()^2 -op(f::Function, args...; kwargs...) = f(op(args...; kwargs...)) +# Rotation around X-axis +# exp(-im * n.θ / 2 * X) +alias(n::OpName"Rx") = cis(-(n.θ / 2) * OpName"X"()) + +alias(::OpName"Y") = -im * (OpName"σ⁺"() - OpName"σ⁻"()) / 2 +# TODO: No subsript `\_y` available +# in unicode. +@op_alias "σʸ" "Y" +@op_alias "σ2" "Y" +@op_alias "σ²" "Y" +@op_alias "σ₂" "Y" +alias(::OpName"iY") = (OpName"σ⁺"() - OpName"σ⁻"()) / 2 +@op_alias "iσy" "iY" +# TODO: No subsript `\_y` available +# in unicode. +@op_alias "iσʸ" "iY" +@op_alias "iσ2" "iY" +@op_alias "iσ²" "iY" +@op_alias "iσ₂" "iY" +alias(n::OpName"Sy") = OpName"Y"() / 2 +@op_alias "Sʸ" "Sy" +alias(n::OpName"iSy") = OpName"iY"() / 2 +@op_alias "iSʸ" "iSy" +alias(::OpName"Sy2") = -OpName"iSy"()^2 + +# Rotation around Y-axis +# exp(-im * n.θ / 2 * Y) +alias(n::OpName"Ry") = exp(-(n.θ / 2) * OpName"iY"()) + +# Ising (XX) coupling gate +# exp(-im * θ/2 * X ⊗ X) +alias(n::OpName"Rxx") = exp(-im * (n.θ / 2) * OpName"X"() ⊗ OpName"X"()) +@op_alias "RXX" "Rxx" + +# Ising (YY) coupling gate +# exp(-im * θ/2 * Y ⊗ Y) +alias(n::OpName"Ryy") = exp(-im * (n.θ / 2) * OpName"Y"() ⊗ OpName"Y"()) +@op_alias "RYY" "Ryy" + +# Ising (ZZ) coupling gate +# exp(-im * θ/2 * Z ⊗ Z) +alias(n::OpName"Rzz") = exp(-im * (n.θ / 2) * OpName"Z"() ⊗ OpName"Z"()) +@op_alias "RZZ" "Rzz" + +## TODO: Check this definition and see if it is worth defining this. +## # Ising (XY) coupling gate +## # exp(-im * θ/2 * X ⊗ Y) +## alias(n::OpName"Rxy") = exp(-im * (n.θ / 2) * OpName"X"() ⊗ OpName"Y"()) +## @op_alias "RXY" "Rxy" + +function Base.AbstractArray(n::OpName"σᶻ", domain_size::Tuple{Int}) + d = only(domain_size) + s = (d - 1) / 2 + return Diagonal([2 * (s + 1 - i) for i in 1:d]) +end +# TODO: No subsript `\_z` available +# in unicode. +@op_alias "Z" "σᶻ" +@op_alias "σ3" "Z" +@op_alias "σ³" "Z" +@op_alias "σ₃" "Z" +@op_alias "σz" "Z" +alias(::OpName"iZ") = im * OpName"Z"() +alias(n::OpName"Sz") = OpName"Z"() / 2 +@op_alias "Sᶻ" "Sz" +# TODO: Make sure it ends up real, using `S⁺` and `S⁻`, +# define `Sx²`, `Sy²`, and `Sz²`, etc. +# Should be equal to: +# ```julia +# α * I = s * (s + 1) * I +# = (d - 1) * (d + 1) / 4 * I +# ``` +# where `s = (d - 1) / 2`. See https://en.wikipedia.org/wiki/Spin_(physics)#Higher_spins. +alias(n::OpName"S2") = OpName"Sx2"() + OpName"Sy2"() + OpName"Sz2"() +@op_alias "S²" "S2" +alias(::OpName"Sz2") = OpName"Sz"()^2 + +# Rotation around Z-axis +# exp(-im * n.θ / 2 * Z) +alias(n::OpName"Rz") = exp(-im * (n.θ / 2) * OpName"Z"()) + +using LinearAlgebra: eigen +function Base.AbstractArray(n::OpName"H", domain_size::Tuple{Int}) + Λ, H = eigen(AbstractArray(OpName("X"), domain_size)) + p = sortperm(Λ; rev=true) + return H[:, p] +end + +using LinearAlgebra: Diagonal +nsites(::OpName"SWAP") = 2 +function Base.AbstractArray(::OpName"SWAP", domain_size::Tuple{Int,Int}) + I_matrix = Diagonal(trues(prod(domain_size))) + I_array = reshape(I_matrix, (domain_size..., domain_size...)) + SWAP_array = permutedims(I_array, (2, 1, 3, 4)) + SWAP_matrix = reshape(SWAP_array, (prod(domain_size), prod(domain_size))) + return SWAP_matrix +end +@op_alias "Swap" "SWAP" +alias(::OpName"√SWAP") = √(OpName"SWAP"()) +@op_alias "√Swap" "√SWAP" + +using LinearAlgebra: diagind +nsites(::OpName"iSWAP") = 2 +function Base.AbstractArray(::OpName"iSWAP", domain_size::Tuple{Int,Int}) + swap = AbstractArray(OpName"SWAP"(), domain_size) + iswap = im * swap + iswap[diagind(iswap)] .*= -im + return iswap +end +@op_alias "iSwap" "iSWAP" +alias(::OpName"√iSWAP") = √(OpName"iSWAP"()) +@op_alias "√iSwap" "√iSWAP" + +## TODO: Bring back these definitions. +## 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 + +# Unary operations +nsites(n::OpName"f") = nsites(n.op) +function Base.AbstractArray(n::OpName"f", domain_size::Tuple{Vararg{Int}}) + return n.f(AbstractArray(n.op, domain_size)) +end + +for f in ( + :(Base.sqrt), + :(Base.real), + :(Base.imag), + :(Base.complex), + :(Base.exp), + :(Base.cis), + :(Base.cos), + :(Base.sin), + :(Base.adjoint), + :(Base.:+), + :(Base.:-), +) + @eval begin + $f(n::OpName) = OpName"f"(; f=$f, op=n) + end +end + +nsites(n::OpName"^") = nsites(n.op) +function Base.AbstractArray(n::OpName"^", domain_size::Tuple{Vararg{Int}}) + return AbstractArray(n.op, domain_size)^n.exponent +end +Base.:^(n::OpName, exponent) = OpName"^"(; op=n, exponent) + +nsites(n::OpName"kron") = nsites(n.op1) + nsites(n.op2) +function Base.AbstractArray(n::OpName"kron", domain_size::Tuple{Vararg{Int}}) + domain_size1 = domain_size[1:nsites(n.op1)] + domain_size2 = domain_size[(nsites(n.op1) + 1):end] + @assert length(domain_size2) == nsites(n.op2) + return kron(AbstractArray(n.op1, domain_size1), AbstractArray(n.op2, domain_size2)) +end +Base.kron(n1::OpName, n2::OpName) = OpName"kron"(; op1=n1, op2=n2) +⊗(n1::OpName, n2::OpName) = kron(n1, n2) + +function nsites(n::OpName"+") + @assert nsites(n.op1) == nsites(n.op2) + return nsites(n.op1) +end +function Base.AbstractArray(n::OpName"+", domain_size::Tuple{Vararg{Int}}) + return AbstractArray(n.op1, domain_size) + AbstractArray(n.op2, domain_size) +end +Base.:+(n1::OpName, n2::OpName) = OpName"+"(; op1=n1, op2=n2) +Base.:-(n1::OpName, n2::OpName) = n1 + (-n2) + +function nsites(n::OpName"*") + @assert nsites(n.op1) == nsites(n.op2) + return nsites(n.op1) +end +function Base.AbstractArray(n::OpName"*", domain_size::Tuple{Vararg{Int}}) + return AbstractArray(n.op1, domain_size) * AbstractArray(n.op2, domain_size) +end +Base.:*(n1::OpName, n2::OpName) = OpName"*"(; op1=n1, op2=n2) + +nsites(n::OpName"scaled") = nsites(n.op) +function Base.AbstractArray(n::OpName"scaled", domain_size::Tuple{Vararg{Int}}) + return AbstractArray(n.op, domain_size) * n.c +end +function Base.:*(c::Number, n::OpName) + return OpName"scaled"(; op=n, c) +end +function Base.:*(n::OpName, c::Number) + return OpName"scaled"(; op=n, c) +end +function Base.:/(n::OpName, c::Number) + return OpName"scaled"(; op=n, c=inv(c)) +end + +function Base.:*(c::Number, n::OpName"scaled") + return OpName"scaled"(; op=n.op, c=(c * n.c)) +end +function Base.:*(n::OpName"scaled", c::Number) + return OpName"scaled"(; op=n.op, c=(n.c * c)) +end +function Base.:/(n::OpName"scaled", c::Number) + return OpName"scaled"(; op=n.op, c=(n.c / c)) +end + +controlled(n::OpName; ncontrol=1) = OpName"Controlled"(; ncontrol, op=n) + +# Expand the operator in a new basis. +using LinearAlgebra: ⋅ +function expand(v::AbstractArray, basis) + gramian = [basisᵢ ⋅ basisⱼ for basisᵢ in basis, basisⱼ in basis] + vbasis = [basisᵢ ⋅ v for basisᵢ in basis] + return gramian \ vbasis +end +function expand(n::OpName, basis, ts...) + return expand(AbstractArray(n, ts...), map(basisᵢ -> AbstractArray(basisᵢ, ts...), basis)) +end diff --git a/src/sitetype.jl b/src/sitetype.jl index d6459cd..2d06f9b 100644 --- a/src/sitetype.jl +++ b/src/sitetype.jl @@ -1,13 +1,23 @@ -using ChainRulesCore: @ignore_derivatives +struct SiteType{T,Params} + params::Params +end +params(t::SiteType) = getfield(t, :params) + +Base.getproperty(t::SiteType, name::Symbol) = getfield(params(t), name) -# TODO: Need to define or replace. -# using ITensorBase: product, swapprime +SiteType{N}(params) where {N} = SiteType{N,typeof(params)}(params) +SiteType{N}(; kwargs...) where {N} = SiteType{N}((; kwargs...)) -# TODO: Add `params<:NamedTuple` field. -struct SiteType{T} end -SiteType(s::AbstractString) = SiteType{Symbol(s)}() -SiteType(i::Integer) = SiteType{Symbol(i)}() +SiteType(s::AbstractString; kwargs...) = SiteType{Symbol(s)}(; kwargs...) +SiteType(i::Integer; kwargs...) = SiteType{Symbol(i)}(; kwargs...) value(::SiteType{T}) where {T} = T macro SiteType_str(s) return SiteType{Symbol(s)} end + +alias(t::SiteType) = t +Base.AbstractUnitRange(t::SiteType) = Base.OneTo(length(t)) +Base.size(t::SiteType) = (length(t),) +Base.size(t::SiteType, dim::Integer) = size(t)[dim] +Base.axes(t::SiteType) = (AbstractUnitRange(t),) +Base.axes(t::SiteType, dim::Integer) = axes(t)[dim] diff --git a/src/sitetypes/aliases.jl b/src/sitetypes/aliases.jl deleted file mode 100644 index 727ad7f..0000000 --- a/src/sitetypes/aliases.jl +++ /dev/null @@ -1,40 +0,0 @@ -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/boson.jl b/src/sitetypes/boson.jl deleted file mode 100644 index 6a011a4..0000000 --- a/src/sitetypes/boson.jl +++ /dev/null @@ -1,31 +0,0 @@ -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/electron.jl b/src/sitetypes/electron.jl index d8fd908..782a85b 100644 --- a/src/sitetypes/electron.jl +++ b/src/sitetypes/electron.jl @@ -1,81 +1,26 @@ -""" - 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") +Base.length(::SiteType"Electron") = 4 -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 +Base.AbstractArray(::StateName"Emp", ::Tuple{SiteType"Electron"}) = [1.0, 0, 0, 0] +Base.AbstractArray(::StateName"Up", ::Tuple{SiteType"Electron"}) = [0.0, 1, 0, 0] +Base.AbstractArray(::StateName"Dn", ::Tuple{SiteType"Electron"}) = [0.0, 0, 1, 0] +Base.AbstractArray(::StateName"UpDn", ::Tuple{SiteType"Electron"}) = [0.0, 0, 0, 1] +# TODO: Use aliasing. +function Base.AbstractArray(::StateName"0", st::Tuple{SiteType"Electron"}) + return AbstractArray(StateName("Emp"), st) +end +function Base.AbstractArray(::StateName"↑", st::Tuple{SiteType"Electron"}) + return AbstractArray(StateName("Up"), st) +end +function Base.AbstractArray(::StateName"↓", st::Tuple{SiteType"Electron"}) + return AbstractArray(StateName("Dn"), st) +end +function Base.AbstractArray(::StateName"↑↓", st::Tuple{SiteType"Electron"}) + return AbstractArray(StateName("UpDn"), st) 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") +# I ⊗ n +# TODO: Define as `AbstractArray(OpName"I"() ⊗ OpName"n"(), (SiteType("Fermion"), SiteType("Fermion")))`? +function Base.AbstractArray(::OpName"Nup", ::Tuple{SiteType"Electron"}) return [ 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 @@ -83,11 +28,12 @@ function op(::OpName"Nup", ::SiteType"Electron") 0.0 0.0 0.0 1.0 ] end -function op(on::OpName"n↑", st::SiteType"Electron") - return op(alias(on), st) +function Base.AbstractArray(on::OpName"n↑", st::Tuple{SiteType"Electron"}) + return AbstractArray(alias(on), st) end -function op(::OpName"Ndn", ::SiteType"Electron") +# n ⊗ I +function Base.AbstractArray(::OpName"Ndn", ::Tuple{SiteType"Electron"}) return [ 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 @@ -95,11 +41,12 @@ function op(::OpName"Ndn", ::SiteType"Electron") 0.0 0.0 0.0 1.0 ] end -function op(on::OpName"n↓", st::SiteType"Electron") - return op(alias(on), st) +function Base.AbstractArray(on::OpName"n↓", st::Tuple{SiteType"Electron"}) + return AbstractArray(alias(on), st) end -function op(::OpName"Nupdn", ::SiteType"Electron") +# n ⊗ n +function Base.AbstractArray(::OpName"Nupdn", ::Tuple{SiteType"Electron"}) return [ 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 @@ -107,11 +54,12 @@ function op(::OpName"Nupdn", ::SiteType"Electron") 0.0 0.0 0.0 1.0 ] end -function op(on::OpName"n↑↓", st::SiteType"Electron") - return op(alias(on), st) +function Base.AbstractArray(on::OpName"n↑↓", st::Tuple{SiteType"Electron"}) + return AbstractArray(alias(on), st) end -function op(::OpName"Ntot", ::SiteType"Electron") +# I ⊗ n + n ⊗ I +function Base.AbstractArray(::OpName"Ntot", ::Tuple{SiteType"Electron"}) return [ 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 @@ -119,11 +67,12 @@ function op(::OpName"Ntot", ::SiteType"Electron") 0.0 0.0 0.0 2.0 ] end -function op(on::OpName"ntot", st::SiteType"Electron") - return op(alias(on), st) +function Base.AbstractArray(on::OpName"ntot", st::Tuple{SiteType"Electron"}) + return AbstractArray(alias(on), st) end -function op(::OpName"Cup", ::SiteType"Electron") +# I ⊗ c +function Base.AbstractArray(::OpName"Cup", ::Tuple{SiteType"Electron"}) return [ 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 @@ -131,11 +80,12 @@ function op(::OpName"Cup", ::SiteType"Electron") 0.0 0.0 0.0 0.0 ] end -function op(on::OpName"c↑", st::SiteType"Electron") - return op(alias(on), st) +function Base.AbstractArray(on::OpName"c↑", st::Tuple{SiteType"Electron"}) + return AbstractArray(alias(on), st) end -function op(::OpName"Cdagup", ::SiteType"Electron") +# I ⊗ c† +function Base.AbstractArray(::OpName"Cdagup", ::Tuple{SiteType"Electron"}) return [ 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 @@ -143,11 +93,12 @@ function op(::OpName"Cdagup", ::SiteType"Electron") 0.0 0.0 1.0 0.0 ] end -function op(on::OpName"c†↑", st::SiteType"Electron") - return op(alias(on), st) +function Base.AbstractArray(on::OpName"c†↑", st::Tuple{SiteType"Electron"}) + return AbstractArray(alias(on), st) end -function op(::OpName"Cdn", ::SiteType"Electron") +# c ⊗ F +function Base.AbstractArray(::OpName"Cdn", ::Tuple{SiteType"Electron"}) return [ 0.0 0.0 1.0 0.0 0.0 0.0 0.0 -1.0 @@ -155,11 +106,12 @@ function op(::OpName"Cdn", ::SiteType"Electron") 0.0 0.0 0.0 0.0 ] end -function op(on::OpName"c↓", st::SiteType"Electron") - return op(alias(on), st) +function Base.AbstractArray(on::OpName"c↓", st::Tuple{SiteType"Electron"}) + return AbstractArray(alias(on), st) end -function op(::OpName"Cdagdn", ::SiteType"Electron") +# c† ⊗ F +function Base.AbstractArray(::OpName"Cdagdn", ::Tuple{SiteType"Electron"}) return [ 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 @@ -167,11 +119,12 @@ function op(::OpName"Cdagdn", ::SiteType"Electron") 0.0 -1.0 0.0 0.0 ] end -function op(::OpName"c†↓", st::SiteType"Electron") - return op(OpName("Cdagdn"), st) +function Base.AbstractArray(::OpName"c†↓", st::Tuple{SiteType"Electron"}) + return AbstractArray(OpName("Cdagdn"), st) end -function op(::OpName"Aup", ::SiteType"Electron") +# I ⊗ a +function Base.AbstractArray(::OpName"Aup", ::Tuple{SiteType"Electron"}) return [ 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 @@ -179,11 +132,12 @@ function op(::OpName"Aup", ::SiteType"Electron") 0.0 0.0 0.0 0.0 ] end -function op(::OpName"a↑", st::SiteType"Electron") - return op(OpName("Aup"), st) +function Base.AbstractArray(::OpName"a↑", st::Tuple{SiteType"Electron"}) + return AbstractArray(OpName("Aup"), st) end -function op(::OpName"Adagup", ::SiteType"Electron") +# I ⊗ a† +function Base.AbstractArray(::OpName"Adagup", ::Tuple{SiteType"Electron"}) return [ 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 @@ -191,11 +145,12 @@ function op(::OpName"Adagup", ::SiteType"Electron") 0.0 0.0 1.0 0.0 ] end -function op(::OpName"a†↑", st::SiteType"Electron") - return op(OpName("Adagup"), st) +function Base.AbstractArray(::OpName"a†↑", st::Tuple{SiteType"Electron"}) + return AbstractArray(OpName("Adagup"), st) end -function op(::OpName"Adn", ::SiteType"Electron") +# a ⊗ I +function Base.AbstractArray(::OpName"Adn", ::Tuple{SiteType"Electron"}) return [ 0.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0 @@ -203,11 +158,12 @@ function op(::OpName"Adn", ::SiteType"Electron") 0.0 0.0 0.0 0.0 ] end -function op(::OpName"a↓", st::SiteType"Electron") - return op(OpName("Adn"), st) +function Base.AbstractArray(::OpName"a↓", st::Tuple{SiteType"Electron"}) + return AbstractArray(OpName("Adn"), st) end -function op(::OpName"Adagdn", ::SiteType"Electron") +# a† ⊗ I +function Base.AbstractArray(::OpName"Adagdn", ::Tuple{SiteType"Electron"}) return [ 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 @@ -215,11 +171,12 @@ function op(::OpName"Adagdn", ::SiteType"Electron") 0.0 1.0 0.0 0.0 ] end -function op(::OpName"a†↓", st::SiteType"Electron") - return op(OpName("Adagdn"), st) +function Base.AbstractArray(::OpName"a†↓", st::Tuple{SiteType"Electron"}) + return AbstractArray(OpName("Adagdn"), st) end -function op(::OpName"F", ::SiteType"Electron") +# F ⊗ F +function Base.AbstractArray(::OpName"F", ::Tuple{SiteType"Electron"}) return [ 1.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 @@ -228,7 +185,8 @@ function op(::OpName"F", ::SiteType"Electron") ] end -function op(::OpName"Fup", ::SiteType"Electron") +# I ⊗ F +function Base.AbstractArray(::OpName"Fup", ::Tuple{SiteType"Electron"}) return [ 1.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 @@ -236,11 +194,12 @@ function op(::OpName"Fup", ::SiteType"Electron") 0.0 0.0 0.0 -1.0 ] end -function op(::OpName"F↑", st::SiteType"Electron") - return op(OpName("Fup"), st) +function Base.AbstractArray(::OpName"F↑", st::Tuple{SiteType"Electron"}) + return AbstractArray(OpName("Fup"), st) end -function op(::OpName"Fdn", ::SiteType"Electron") +# F ⊗ I +function Base.AbstractArray(::OpName"Fdn", ::Tuple{SiteType"Electron"}) return [ 1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 @@ -248,13 +207,12 @@ function op(::OpName"Fdn", ::SiteType"Electron") 0.0 0.0 0.0 -1.0 ] end -function op(::OpName"F↓", st::SiteType"Electron") - return op(OpName("Fdn"), st) +function Base.AbstractArray(::OpName"F↓", st::Tuple{SiteType"Electron"}) + return AbstractArray(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 +function Base.AbstractArray(::OpName"Sz", ::Tuple{SiteType"Electron"}) + # cat(falses(1, 1), Matrix(OpName("Sz")), falses(1, 1); dims=(1, 2)) return [ 0.0 0.0 0.0 0.0 0.0 0.5 0.0 0.0 @@ -262,12 +220,12 @@ function op(::OpName"Sz", ::SiteType"Electron") 0.0 0.0 0.0 0.0 ] end - -function op(::OpName"Sᶻ", st::SiteType"Electron") - return op(OpName("Sz"), st) +function Base.AbstractArray(::OpName"Sᶻ", st::Tuple{SiteType"Electron"}) + return AbstractArray(OpName("Sz"), st) end -function op(::OpName"Sx", ::SiteType"Electron") +function Base.AbstractArray(::OpName"Sx", ::Tuple{SiteType"Electron"}) + # cat(falses(1, 1), Matrix(OpName("Sx")), falses(1, 1); dims=(1, 2)) return [ 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.0 @@ -276,11 +234,12 @@ function op(::OpName"Sx", ::SiteType"Electron") ] end -function op(::OpName"Sˣ", st::SiteType"Electron") - return op(OpName("Sx"), st) +function Base.AbstractArray(::OpName"Sˣ", st::Tuple{SiteType"Electron"}) + return AbstractArray(OpName("Sx"), st) end -function op(::OpName"S+", ::SiteType"Electron") +function Base.AbstractArray(::OpName"S+", ::Tuple{SiteType"Electron"}) + # cat(falses(1, 1), Matrix(OpName("S+")), falses(1, 1); dims=(1, 2)) return [ 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 @@ -289,17 +248,18 @@ function op(::OpName"S+", ::SiteType"Electron") ] end -function op(::OpName"S⁺", st::SiteType"Electron") - return op(OpName("S+"), st) +function Base.AbstractArray(::OpName"S⁺", st::Tuple{SiteType"Electron"}) + return AbstractArray(OpName("S+"), st) end -function op(::OpName"Sp", st::SiteType"Electron") - return op(OpName("S+"), st) +function Base.AbstractArray(::OpName"Sp", st::Tuple{SiteType"Electron"}) + return AbstractArray(OpName("S+"), st) end -function op(::OpName"Splus", st::SiteType"Electron") - return op(OpName("S+"), st) +function Base.AbstractArray(::OpName"Splus", st::Tuple{SiteType"Electron"}) + return AbstractArray(OpName("S+"), st) end -function op(::OpName"S-", ::SiteType"Electron") +function Base.AbstractArray(::OpName"S-", ::Tuple{SiteType"Electron"}) + # cat(falses(1, 1), Matrix(OpName("S-")), falses(1, 1); dims=(1, 2)) return [ 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 @@ -308,29 +268,44 @@ function op(::OpName"S-", ::SiteType"Electron") ] end -function op(::OpName"S⁻", st::SiteType"Electron") - return op(OpName("S-"), st) +function Base.AbstractArray(::OpName"S⁻", st::Tuple{SiteType"Electron"}) + return AbstractArray(OpName("S-"), st) end -function op(::OpName"Sm", st::SiteType"Electron") - return op(OpName("S-"), st) +function Base.AbstractArray(::OpName"Sm", st::Tuple{SiteType"Electron"}) + return AbstractArray(OpName("S-"), st) end -function op(::OpName"Sminus", st::SiteType"Electron") - return op(OpName("S-"), st) +function Base.AbstractArray(::OpName"Sminus", st::Tuple{SiteType"Electron"}) + return AbstractArray(OpName("S-"), st) end -has_fermion_string(::OpName"Cup", ::SiteType"Electron") = true -function has_fermion_string(on::OpName"c↑", st::SiteType"Electron") +@op_alias "a↑" "Aup" +@op_alias "a↓" "Adn" +@op_alias "a†↓" "Adagdn" +@op_alias "a†↑" "Adagup" +@op_alias "c↑" "Cup" +@op_alias "c↓" "Cdn" +@op_alias "c†↑" "Cdagup" +@op_alias "c†↓" "Cdagdn" +@op_alias "n↑" "Nup" +@op_alias "n↓" "Ndn" +@op_alias "n↑↓" "Nupdn" +@op_alias "ntot" "Ntot" +@op_alias "F↑" "Fup" +@op_alias "F↓" "Fdn" + +has_fermion_string(::OpName"Cup", ::Tuple{SiteType"Electron"}) = true +function has_fermion_string(on::OpName"c↑", st::Tuple{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") +has_fermion_string(::OpName"Cdagup", ::Tuple{SiteType"Electron"}) = true +function has_fermion_string(on::OpName"c†↑", st::Tuple{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") +has_fermion_string(::OpName"Cdn", ::Tuple{SiteType"Electron"}) = true +function has_fermion_string(on::OpName"c↓", st::Tuple{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") +has_fermion_string(::OpName"Cdagdn", ::Tuple{SiteType"Electron"}) = true +function has_fermion_string(on::OpName"c†↓", st::Tuple{SiteType"Electron"}) return has_fermion_string(alias(on), st) end diff --git a/src/sitetypes/fermion.jl b/src/sitetypes/fermion.jl index f9735c2..797b584 100644 --- a/src/sitetypes/fermion.jl +++ b/src/sitetypes/fermion.jl @@ -1,101 +1,27 @@ -""" - 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) +Base.length(::SiteType"Fermion") = 2 -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 +Base.AbstractArray(::StateName"Emp", ::SiteType"Fermion") = [1.0 0.0] +Base.AbstractArray(::StateName"Occ", ::SiteType"Fermion") = [0.0 1.0] +function Base.AbstractArray(::StateName"0", st::SiteType"Fermion") + return AbstractArray(StateName("Emp"), st) +end +function Base.AbstractArray(::StateName"1", st::SiteType"Fermion") + return AbstractArray(StateName("Occ"), st) 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) +function Base.AbstractArray(::OpName"F", ::Tuple{SiteType"Fermion"}) + return [ + 1 0 + 0 -1 + ] +end -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) +@op_alias "c" "a" +@op_alias "C" "c" -## 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 +@op_alias "c†" "a†" +@op_alias "Cdag" "c†" +@op_alias "cdag" "c†" has_fermion_string(::OpName"C", ::SiteType"Fermion") = true function has_fermion_string(on::OpName"c", st::SiteType"Fermion") diff --git a/src/sitetypes/generic_sites.jl b/src/sitetypes/generic_sites.jl deleted file mode 100644 index 9f06047..0000000 --- a/src/sitetypes/generic_sites.jl +++ /dev/null @@ -1,51 +0,0 @@ -## using ITensorBase: ITensor, itensor -## using LinearAlgebra: I -## -## # TODO: Need to define or replace. -## # using ITensorBase: 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/qubit.jl b/src/sitetypes/qubit.jl index b0e50af..fef6d8f 100644 --- a/src/sitetypes/qubit.jl +++ b/src/sitetypes/qubit.jl @@ -1,481 +1,164 @@ -# -# Qubit site type -# - -# Define Qubit space in terms of -# Qubit/2 space, but use different -# defaults for QN names - -""" - space(::SiteType"Qubit") - -Create the Hilbert space for a site of type "Qubit". -""" -function space(::SiteType"Qubit") - 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) +using LinearAlgebra: I + +alias(::SiteType"S=1/2") = SiteType"Qubit"() +alias(::SiteType"S=½") = SiteType"Qubit"() +alias(::SiteType"SpinHalf=1/2") = SiteType"Qubit"() + +Base.length(::SiteType"Qubit") = 2 + +# `eigvecs(Z)` +Base.AbstractArray(::StateName"0", ::Tuple{SiteType"Qubit"}) = [1, 0] +Base.AbstractArray(::StateName"1", ::Tuple{SiteType"Qubit"}) = [0, 1] +@state_alias "Up" "0" +@state_alias "↑" "0" +@state_alias "Z+" "0" +@state_alias "Zp" "0" +@state_alias "↓" "1" +@state_alias "Dn" "1" +@state_alias "Z-" "1" +@state_alias "Zm" "1" + +# `eigvecs(X)` +Base.AbstractArray(::StateName"+", ::Tuple{SiteType"Qubit"}) = [1, 1] / √2 +Base.AbstractArray(::StateName"-", ::Tuple{SiteType"Qubit"}) = [1, -1] / √2 +@state_alias "X+" "+" +@state_alias "Xp" "+" +@state_alias "X-" "-" +@state_alias "Xm" "-" + +# `eigvecs(Y)` +Base.AbstractArray(::StateName"i", ::Tuple{SiteType"Qubit"}) = [1, im] / √2 +Base.AbstractArray(::StateName"-i", ::Tuple{SiteType"Qubit"}) = [1, -im] / √2 +@state_alias "Y+" "i" +@state_alias "Yp" "i" +@state_alias "Y-" "-i" +@state_alias "Ym" "-i" # SIC-POVMs -state(::StateName"Tetra1", t::SiteType"Qubit") = state(StateName("Z+"), t) -state(::StateName"Tetra2", t::SiteType"Qubit") = [ - 1 / √3 - √2 / √3 +Base.AbstractArray(::StateName"Tetra0", ::Tuple{SiteType"Qubit"}) = [ + 1 + 0 ] -state(::StateName"Tetra3", t::SiteType"Qubit") = [ +Base.AbstractArray(::StateName"Tetra2", ::Tuple{SiteType"Qubit"}) = [ 1 / √3 - √2 / √3 * exp(im * 2π / 3) -] -state(::StateName"Tetra4", t::SiteType"Qubit") = [ - 1 / √3 - √2 / √3 * exp(im * 4π / 3) -] - -op(::OpName"Id", ::SiteType"Qubit") = [ - 1 0 - 0 1 -] - -# -# 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) + √2 / √3 ] - -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(θ) && (θ = ϕ) +function Base.AbstractArray(::StateName"Tetra3", ::Tuple{SiteType"Qubit"}) return [ - exp(-im * θ / 2) 0 - 0 exp(im * θ / 2) + 1 / √3 + √2 / √3 * exp(im * 2π / 3) ] end - -# Rotation around generic axis n̂ -function op(::OpName"Rn", ::SiteType"Qubit"; θ::Real, ϕ::Real, λ::Real) +function Base.AbstractArray(::StateName"Tetra4", ::Tuple{SiteType"Qubit"}) return [ - cos(θ / 2) -exp(im * λ)*sin(θ / 2) - exp(im * ϕ)*sin(θ / 2) exp(im * (ϕ + λ))*cos(θ / 2) + 1 / √3 + √2 / √3 * exp(im * 4π / 3) ] 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 +# TODO: Write as `(I + σᶻ) / 2`? +Base.AbstractArray(::OpName"ProjUp", ::Tuple{SiteType"Qubit"}) = [ + 1 0 + 0 0 ] +@op_alias "projUp" "ProjUp" +@op_alias "Proj↑" "ProjUp" +@op_alias "proj↑" "ProjUp" +@op_alias "Proj0" "ProjUp" +@op_alias "proj0" "ProjUp" -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 +# TODO: Define as `σ⁺ * σ−`? +# TODO: Write as `(I - σᶻ) / 2`? +Base.AbstractArray(::OpName"ProjDn", ::Tuple{SiteType"Qubit"}) = [ + 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_alias "projDn" "ProjDn" +@op_alias "Proj↓" "ProjDn" +@op_alias "proj↓" "ProjDn" +@op_alias "Proj1" "ProjDn" +@op_alias "proj1" "ProjDn" -op(::OpName"iSWAP", t::SiteType"Qubit") = [ - 1 0 0 0 - 0 0 im 0 - 0 im 0 0 - 0 0 0 1 +# TODO: Determine a general spin definition, such as +# `eigvecs(X)`. +Base.AbstractArray(::OpName"H", ::Tuple{SiteType"Qubit"}) = [ + 1/√2 1/√2 + 1/√2 -1/√2 ] -op(::OpName"iSwap", t::SiteType"Qubit") = op("iSWAP", t) -function op(::OpName"√iSWAP", t::SiteType"Qubit") +# Rotation around generic axis n̂ +# exp(-im * n.θ / 2 * n̂ ⋅ σ⃗) +#= +TODO: Define R-gate when `λ == -ϕ`, i.e.: +```julia +function Base.AbstractArray(n::OpName"R", ::Tuple{SiteType"Qubit"}) return [ - 1 0 0 0 - 0 1/√2 im/√2 0 - 0 im/√2 1/√2 0 - 0 0 0 1 + cos(n.θ / 2) -exp(-im * n.ϕ)*sin(n.θ / 2) + exp(im * n.ϕ)*sin(n.θ / 2) cos(n.θ / 2) ] end -op(::OpName"√iSwap", t::SiteType"Qubit") = op("√iSWAP", t) - -# Ising (XX) coupling gate -function op(::OpName"Rxx", t::SiteType"Qubit"; ϕ::Number) +``` +or: +```julia +alias(n::OpName"R") = OpName"Rn"(; θ=n.θ, ϕ=n.ϕ, λ=-n.ϕ) +=# +# https://docs.quantum.ibm.com/api/qiskit/qiskit.circuit.library.UGate +function Base.AbstractArray(n::OpName"Rn", ::Tuple{SiteType"Qubit"}) return [ - cos(ϕ) 0 0 -im*sin(ϕ) - 0 cos(ϕ) -im*sin(ϕ) 0 - 0 -im*sin(ϕ) cos(ϕ) 0 - -im*sin(ϕ) 0 0 cos(ϕ) + cos(n.θ / 2) -exp(im * n.λ)*sin(n.θ / 2) + exp(im * n.ϕ)*sin(n.θ / 2) exp(im * (n.ϕ + n.λ))*cos(n.θ / 2) ] end -op(::OpName"RXX", t::SiteType"Qubit"; kwargs...) = op("Rxx", t; kwargs...) +@op_alias "Rn̂" "Rn" -# Ising (YY) coupling gate -function op(::OpName"Ryy", ::SiteType"Qubit"; ϕ::Number) +# TODO: Generalize to `"Qudit"` and other SiteTypes. +# https://docs.quantum.ibm.com/api/qiskit/qiskit.circuit.library.UCGate +nsites(n::OpName"Controlled") = get(params(n), :ncontrol, 1) + nsites(n.op) +function Base.AbstractArray(n::OpName"Controlled", ts::Tuple{Vararg{SiteType"Qubit"}}) + # Number of target qubits. + nt = nsites(n.op) + # Number of control qubits. + nc = get(params(n), :ncontrol, length(ts) - nt) + @assert length(ts) == nc + nt return [ - cos(ϕ) 0 0 im*sin(ϕ) - 0 cos(ϕ) -im*sin(ϕ) 0 - 0 -im*sin(ϕ) cos(ϕ) 0 - im*sin(ϕ) 0 0 cos(ϕ) + I(2^nc) falses(2^nc, 2^nt) + falses(2^nt, 2^nc) AbstractArray(n.op, ts[(nc + 1):end]) ] 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 - ] +@op_alias "CNOT" "Controlled" op = OpName"X"() +@op_alias "CX" "Controlled" op = OpName"X"() +@op_alias "CY" "Controlled" op = OpName"Y"() +@op_alias "CZ" "Controlled" op = OpName"Z"() +function alias(n::OpName"CPhase") + return controlled(OpName"Phase"(; params(n)...)) 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 * ϕ) - ] +@op_alias "CPHASE" "CPhase" +@op_alias "Cphase" "CPhase" +function alias(n::OpName"CRx") + return controlled(OpName"Rx"(; params(n)...)) 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 - ] +@op_alias "CRX" "CRx" +function Base.AbstractArray(::OpName"CRy") + return controlled(OpName"Ry"(; params(n)...)) 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 - ] +@op_alias "CRY" "CRy" +function Base.AbstractArray(::OpName"CRz") + return controlled(OpName"Rz"(; params(n)...)) 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 - ] +@op_alias "CRZ" "CRz" +function Base.AbstractArray(::OpName"CRn") + return controlled(; op=OpName"Rn"(; params(n)...)) end +@op_alias "CRn̂" "CRn" -# 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_alias "CCNOT" "Controlled" ncontrol = 2 op = OpName"X"() +@op_alias "Toffoli" "CCNOT" +@op_alias "CCX" "CCNOT" +@op_alias "TOFF" "CCNOT" -op(on::OpName"projDn", t::SiteType"Qubit") = op(alias(on), t) +@op_alias "CSWAP" "Controlled" ncontrol = 2 op = OpName"SWAP"() +@op_alias "Fredkin" "CSWAP" +@op_alias "CSwap" "CSWAP" +@op_alias "CS" "CSWAP" -op(on::OpName"Proj1", t::SiteType"Qubit") = op(alias(on), t) +@op_alias "CCCNOT" "Controlled" ncontrol = 3 op = OpName"X"() diff --git a/src/sitetypes/qudit.jl b/src/sitetypes/qudit.jl index 2ae2489..ad8463d 100644 --- a/src/sitetypes/qudit.jl +++ b/src/sitetypes/qudit.jl @@ -1,97 +1,2 @@ -using ChainRulesCore: @non_differentiable - -""" - space(::SiteType"Qudit") - -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) - return dim -end - -function val(::ValName{N}, ::SiteType"Qudit") where {N} - return parse(Int, String(N)) + 1 -end - -function state(::StateName{N}, ::SiteType"Qudit") where {N} - n = parse(Int, String(N)) - st = zeros(dim(s)) - st[n + 1] = 1.0 - return st -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) +Base.length(t::SiteType"Qudit") = t.length +alias(::SiteType"Boson") = SiteType"Qudit"() diff --git a/src/sitetypes/spinhalf.jl b/src/sitetypes/spinhalf.jl deleted file mode 100644 index c374d47..0000000 --- a/src/sitetypes/spinhalf.jl +++ /dev/null @@ -1,66 +0,0 @@ - -""" - 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/spinone.jl b/src/sitetypes/spinone.jl index 91766f1..87d612c 100644 --- a/src/sitetypes/spinone.jl +++ b/src/sitetypes/spinone.jl @@ -1,150 +1,28 @@ -# TODO: Need to define or replace. -# using ITensorBase: complex!, QN +Base.length(::SiteType"S=1") = 3 alias(::SiteType"SpinOne") = SiteType"S=1"() -""" - space(::SiteType"S=1"; - conserve_qns = false, - conserve_sz = conserve_qns, - qnname_sz = "Sz") +# TODO: Decide on these names, use `alias`. +# TODO: Make a more general `SiteType"Spin`/`SiteType"S"` +# with a `spin` field that can be set to a rational number +# `1//2, `2//2`, `3//2`, etc. that maps to `Qudit` +# of length `2 * spin + 1`. +Base.AbstractArray(::StateName"Up", ::SiteType"S=1") = [1, 0, 0] +Base.AbstractArray(::StateName"Z0", ::SiteType"S=1") = [0, 1, 0] +Base.AbstractArray(::StateName"Dn", ::SiteType"S=1") = [0, 0, 1] -Create the Hilbert space for a site of type "S=1". +Base.AbstractArray(::StateName"↑", st::SiteType"S=1") = [1, 0, 0] +Base.AbstractArray(::StateName"0", st::SiteType"S=1") = [0, 1, 0] +Base.AbstractArray(::StateName"↓", st::SiteType"S=1") = [0, 0, 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] +Base.AbstractArray(::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"Id", ::SiteType"S=1") = [ - 1.0 0.0 0.0 - 0.0 1.0 0.0 - 0.0 0.0 1.0 -] - -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)) +Base.AbstractArray(::StateName"Z-", st::SiteType"S=1") = [0.0, 0.0, 1.0] -## function op!(Op::ITensor, o::OpName, st::SiteType"SpinOne", s::Index) -## return op!(Op, o, alias(st), s) -## end +Base.AbstractArray(::StateName"X+", ::SiteType"S=1") = [1 / 2, 1 / sqrt(2), 1 / 2] +Base.AbstractArray(::StateName"X0", ::SiteType"S=1") = [-1 / sqrt(2), 0, 1 / sqrt(2)] +Base.AbstractArray(::StateName"X-", ::SiteType"S=1") = [1 / 2, -1 / sqrt(2), 1 / 2] -op(o::OpName, st::SiteType"SpinOne") = op(o, alias(st)) +Base.AbstractArray(::StateName"Y+", ::SiteType"S=1") = [-1 / 2, -im / sqrt(2), 1 / 2] +Base.AbstractArray(::StateName"Y0", ::SiteType"S=1") = [1 / sqrt(2), 0, 1 / sqrt(2)] +Base.AbstractArray(::StateName"Y-", ::SiteType"S=1") = [-1 / 2, im / sqrt(2), 1 / 2] diff --git a/src/sitetypes/tj.jl b/src/sitetypes/tj.jl index e603f19..92bebce 100644 --- a/src/sitetypes/tj.jl +++ b/src/sitetypes/tj.jl @@ -1,72 +1,13 @@ -""" - 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") +Base.length(::SiteType"tJ") = 3 -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) +Base.AbstractArray(::StateName"Emp", ::SiteType"tJ") = [1.0, 0, 0] +Base.AbstractArray(::StateName"Up", ::SiteType"tJ") = [0.0, 1, 0] +Base.AbstractArray(::StateName"Dn", ::SiteType"tJ") = [0.0, 0, 1] +Base.AbstractArray(::StateName"0", st::SiteType"tJ") = AbstractArray(StateName("Emp"), st) +Base.AbstractArray(::StateName"↑", st::SiteType"tJ") = AbstractArray(StateName("Up"), st) +Base.AbstractArray(::StateName"↓", st::SiteType"tJ") = AbstractArray(StateName("Dn"), st) +# TODO: Update these to the latest syntax. ## function op!(Op::ITensor, ::OpName"Nup", ::SiteType"tJ", s::Index) ## return Op[s' => 2, s => 2] = 1.0 ## end diff --git a/src/space.jl b/src/space.jl deleted file mode 100644 index 78f8fad..0000000 --- a/src/space.jl +++ /dev/null @@ -1,7 +0,0 @@ -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: $(value(st))" -end diff --git a/src/state.jl b/src/state.jl index 1e2ccb8..d0cb08b 100644 --- a/src/state.jl +++ b/src/state.jl @@ -1,13 +1,91 @@ -@eval struct StateName{Name} - (f::Type{<:StateName})() = $(Expr(:new, :f)) +struct StateName{Name,Params} + params::Params end +params(n::StateName) = getfield(n, :params) -StateName(s::AbstractString) = StateName{Symbol(s)}() -StateName(s::Symbol) = StateName{s}() -name(::StateName{N}) where {N} = N +Base.getproperty(n::StateName, name::Symbol) = getfield(params(n), name) + +StateName{N}(params) where {N} = StateName{N,typeof(params)}(params) +StateName{N}(; kwargs...) where {N} = StateName{N}((; kwargs...)) +StateName(s::AbstractString; kwargs...) = StateName{Symbol(s)}(; kwargs...) +StateName(s::Symbol; kwargs...) = StateName{s}(; kwargs...) +name(::StateName{N}) where {N} = N macro StateName_str(s) return StateName{Symbol(s)} end -state(::StateName, ::SiteType; kwargs...) = nothing +function state_alias_expr(name1, name2, pars...) + return :(function alias(n::StateName{Symbol($name1)}) + return StateName{Symbol($name2)}(; params(n)..., $(esc.(pars)...)) + end) +end +macro state_alias(name1, name2, params...) + return state_alias_expr(name1, name2) +end + +# TODO: Decide on this. +default_sitetype() = SiteType"Qubit"() + +alias(n::StateName) = n +function (arrtype::Type{<:AbstractArray})(n::StateName, ts::Tuple{Vararg{SiteType}}) + # TODO: Define `state_convert` to handle reshaping multisite states + # to higher order arrays. + return convert(arrtype, AbstractArray(n, ts)) +end +function (arrtype::Type{<:AbstractArray})(n::StateName, domain_size::Tuple{Vararg{Integer}}) + # TODO: Define `state_convert` to handle reshaping multisite states + # to higher order arrays. + return convert(arrtype, AbstractArray(n, Int.(domain_size))) +end +function (arrtype::Type{<:AbstractArray})(n::StateName, domain_size::Integer...) + return arrtype(n, domain_size) +end +(arrtype::Type{<:AbstractArray})(n::StateName, ts::SiteType...) = arrtype(n, ts) +function Base.AbstractArray(n::StateName, ts::Tuple{Vararg{SiteType}}) + n′ = alias(n) + ts′ = alias.(ts) + if n′ == n && ts′ == ts + return AbstractArray(n′, length.(ts′)) + end + return AbstractArray(n′, ts′) +end +function Base.AbstractArray(n::StateName, domain_size::Tuple{Vararg{Int}}) + n′ = alias(n) + if n′ == n + error("Not implemented.") + end + return AbstractArray(n′, domain_size) +end + +# TODO: Decide on this. +function Base.AbstractArray(n::StateName) + return AbstractArray(n, ntuple(Returns(default_sitetype()), nsites(n))) +end +function (arrtype::Type{<:AbstractArray})(n::StateName) + return arrtype(n, ntuple(Returns(default_sitetype()), nsites(n))) +end + +function state(arrtype::Type{<:AbstractArray}, n::String, domain...; kwargs...) + return arrtype(StateName(n; kwargs...), domain...) +end +function state(elt::Type{<:Number}, n::String, domain...; kwargs...) + return state(AbstractArray{elt}, n, domain...; kwargs...) +end +function state(n::String, domain...; kwargs...) + return state(AbstractArray, n, domain...; kwargs...) +end + +# TODO: Add this. +## function Base.Integer(::StateName{N}) where {N} +## return parse(Int, String(N)) +## end + +# TODO: Define as `::Tuple{Int}`. +function Base.AbstractArray(n::StateName{N}, domain_size::Tuple{Int}) where {N} + # TODO: Use `Integer(n)`. + n = parse(Int, String(N)) + a = falses(only(domain_size)) + a[n + 1] = one(Bool) + return a +end diff --git a/src/val.jl b/src/val.jl deleted file mode 100644 index f45634e..0000000 --- a/src/val.jl +++ /dev/null @@ -1,14 +0,0 @@ -@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 diff --git a/test/Project.toml b/test/Project.toml index 236b1e3..32b3272 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,8 +1,10 @@ [deps] Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" ITensorBase = "4795dd04-0d67-49bb-8f44-b89c448a1dc7" -ITensorQuantumOperatorDefinitions = "fd9b415b-e710-4e2a-b407-cba023081494" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +QuantumOperatorDefinitions = "826dd319-6fd5-459a-a990-3a4f214664bf" SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/test/basics/test_basics.jl b/test/basics/test_basics.jl deleted file mode 100644 index 26fe749..0000000 --- a/test/basics/test_basics.jl +++ /dev/null @@ -1,6 +0,0 @@ -using ITensorQuantumOperatorDefinitions: ITensorQuantumOperatorDefinitions -using Test: @test, @testset - -@testset "ITensorQuantumOperatorDefinitions" begin - # Tests go here. -end diff --git a/test/test_aqua.jl b/test/test_aqua.jl index 9b7d862..eedff41 100644 --- a/test/test_aqua.jl +++ b/test/test_aqua.jl @@ -1,7 +1,7 @@ -using ITensorQuantumOperatorDefinitions: ITensorQuantumOperatorDefinitions +using QuantumOperatorDefinitions: QuantumOperatorDefinitions using Aqua: Aqua using Test: @testset @testset "Code quality (Aqua.jl)" begin - Aqua.test_all(ITensorQuantumOperatorDefinitions; ambiguities=false) + Aqua.test_all(QuantumOperatorDefinitions; ambiguities=false) end diff --git a/test/test_basics.jl b/test/test_basics.jl new file mode 100644 index 0000000..9ccf717 --- /dev/null +++ b/test/test_basics.jl @@ -0,0 +1,132 @@ +using QuantumOperatorDefinitions: OpName, SiteType, ⊗, expand, op, opexpr, state, nsites +using LinearAlgebra: Diagonal +using Test: @test, @testset + +const real_elts = (Float32, Float64) +const complex_elts = map(elt -> Complex{elt}, real_elts) +const elts = (real_elts..., complex_elts...) + +@testset "QuantumOperatorDefinitions" begin + @testset "Qubit/Qudit" begin + # https://en.wikipedia.org/wiki/Pauli_matrices + # https://en.wikipedia.org/wiki/Spin_(physics)#Higher_spins + for (t, len, Xmat, Ymat, Zmat, Nmat, SWAPmat, iSWAPmat, RXXmat, RYYmat, RZZmat) in ( + ( + SiteType("Qubit"), + 2, + [0 1; 1 0], + [0 -im; im 0], + [1 0; 0 -1], + [0 0; 0 1], + [1 0 0 0; 0 0 1 0; 0 1 0 0; 0 0 0 1], + [1 0 0 0; 0 0 im 0; 0 im 0 0; 0 0 0 1], + (_, θ) -> [ + cos(θ / 2) 0 0 -im*sin(θ / 2) + 0 cos(θ / 2) -im*sin(θ / 2) 0 + 0 -im*sin(θ / 2) cos(θ / 2) 0 + -im*sin(θ / 2) 0 0 cos(θ / 2) + ], + (_, θ) -> [ + cos(θ / 2) 0 0 im*sin(θ / 2) + 0 cos(θ / 2) -im*sin(θ / 2) 0 + 0 -im*sin(θ / 2) cos(θ / 2) 0 + im*sin(θ / 2) 0 0 cos(θ / 2) + ], + (_, θ) -> + Diagonal([exp(-im * θ / 2), exp(im * θ / 2), exp(im * θ / 2), exp(-im * θ / 2)]), + ), + ( + SiteType("Qudit"; length=3), + 3, + √2 * [0 1 0; 1 0 1; 0 1 0], + √2 * [0 -im 0; im 0 -im; 0 im 0], + 2 * [1 0 0; 0 0 0; 0 0 -1], + [0 0 0; 0 1 0; 0 0 2], + [ + 1 0 0 0 0 0 0 0 0 + 0 0 0 1 0 0 0 0 0 + 0 0 0 0 0 0 1 0 0 + 0 1 0 0 0 0 0 0 0 + 0 0 0 0 1 0 0 0 0 + 0 0 0 0 0 0 0 1 0 + 0 0 1 0 0 0 0 0 0 + 0 0 0 0 0 1 0 0 0 + 0 0 0 0 0 0 0 0 1 + ], + [ + 1 0 0 0 0 0 0 0 0 + 0 0 0 im 0 0 0 0 0 + 0 0 0 0 0 0 im 0 0 + 0 im 0 0 0 0 0 0 0 + 0 0 0 0 1 0 0 0 0 + 0 0 0 0 0 0 0 im 0 + 0 0 im 0 0 0 0 0 0 + 0 0 0 0 0 im 0 0 0 + 0 0 0 0 0 0 0 0 1 + ], + (O, θ) -> exp(-im * (θ / 2) * kron(O, O)), + (O, θ) -> exp(-im * (θ / 2) * kron(O, O)), + (O, θ) -> exp(-im * (θ / 2) * kron(O, O)), + ), + ) + @test length(t) == len + @test size(t) == (len,) + @test size(t, 1) == len + @test axes(t) == (Base.OneTo(len),) + @test axes(t, 1) == Base.OneTo(len) + for (o, nbits, elts, ref) in ( + (OpName("X"), 1, elts, Xmat), + (OpName("σˣ"), 1, elts, Xmat), + (OpName("√X"), 1, complex_elts, √Xmat), + (OpName("iX"), 1, complex_elts, im * Xmat), + (OpName("Y"), 1, complex_elts, Ymat), + (OpName("σʸ"), 1, complex_elts, Ymat), + (OpName("iY"), 1, elts, im * Ymat), + (OpName("Z"), 1, elts, Zmat), + (OpName("σᶻ"), 1, elts, Zmat), + (OpName("iZ"), 1, complex_elts, im * Zmat), + (OpName("N"), 1, elts, Nmat), + (OpName("n"), 1, elts, Nmat), + (OpName("Phase"; θ=π / 3), 1, complex_elts, exp(im * π / 3 * Nmat)), + (OpName("π/8"), 1, complex_elts, exp(im * π / 4 * Nmat)), + (OpName("Rx"; θ=π / 3), 1, complex_elts, exp(-im * π / 6 * Xmat)), + (OpName("Ry"; θ=π / 3), 1, complex_elts, exp(-im * π / 6 * Ymat)), + (OpName("Rz"; θ=π / 3), 1, complex_elts, exp(-im * π / 6 * Zmat)), + (OpName("SWAP"), 2, elts, SWAPmat), + (OpName("√SWAP"), 2, complex_elts, √SWAPmat), + (OpName("iSWAP"), 2, complex_elts, iSWAPmat), + (OpName("√iSWAP"), 2, complex_elts, √iSWAPmat), + (OpName("Rxx"; θ=π / 3), 2, complex_elts, RXXmat(Xmat, π / 3)), + (OpName("RXX"; θ=π / 3), 2, complex_elts, RXXmat(Xmat, π / 3)), + (OpName("Ryy"; θ=π / 3), 2, complex_elts, RYYmat(Ymat, π / 3)), + (OpName("RYY"; θ=π / 3), 2, complex_elts, RYYmat(Ymat, π / 3)), + (OpName("Rzz"; θ=π / 3), 2, complex_elts, RZZmat(Zmat, π / 3)), + (OpName("RZZ"; θ=π / 3), 2, complex_elts, RZZmat(Zmat, π / 3)), + ) + @test nsites(o) == nbits + for arraytype in (AbstractArray, AbstractMatrix, Array, Matrix) + for elt in elts + ts = ntuple(Returns(t), nbits) + lens = ntuple(Returns(len), nbits) + for domain in (ts, (ts,), lens, (lens,)) + @test arraytype(o, domain...) ≈ ref + @test arraytype{elt}(o, domain...) ≈ ref + @test eltype(arraytype{elt}(o, domain...)) === elt + end + end + end + end + end + end + @testset "Parsing" begin + # TODO: Should this be called in the `OpName`/`op` constructor? + @test Matrix(opexpr("X * Y")) == op("X") * op("Y") + @test Matrix(opexpr("X * Y + Z")) == op("X") * op("Y") + op("Z") + @test Matrix(opexpr("X * Y + 2 * Z")) == op("X") * op("Y") + 2 * op("Z") + @test Matrix(opexpr("exp(im * (X * Y + 2 * Z))")) == + exp(im * (op("X") * op("Y") + 2 * op("Z"))) + @test Matrix(opexpr("exp(im * (X ⊗ Y + Z ⊗ Z))")) == + exp(im * (kron(op("X"), op("Y")) + kron(op("Z"), op("Z")))) + @test Matrix(opexpr("Ry{θ=π/2}")) == op("Ry"; θ=π / 2) + end +end