diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml deleted file mode 100644 index 9613e05..0000000 --- a/.JuliaFormatter.toml +++ /dev/null @@ -1 +0,0 @@ -style = "yas" \ No newline at end of file diff --git a/.github/workflows/FormatCheck.yml b/.github/workflows/FormatCheck.yml index 439646d..383c502 100644 --- a/.github/workflows/FormatCheck.yml +++ b/.github/workflows/FormatCheck.yml @@ -1,63 +1,15 @@ -name: FormatCheck +name: 'Format' on: - push: - branches: - - 'main' - - 'master' - - 'release-' - tags: '*' - pull_request: + pull_request_target: + paths: ['**/*.jl'] + types: [opened, synchronize, reopened, ready_for_review] -jobs: - build: - runs-on: ${{ matrix.os }} - strategy: - matrix: - version: - - '1' # automatically expands to the latest stable 1.x release of Julia - os: - - ubuntu-latest - arch: - - x64 - steps: - - uses: julia-actions/setup-julia@latest - with: - version: ${{ matrix.version }} - arch: ${{ matrix.arch }} - - - uses: actions/checkout@v5 - - name: Install JuliaFormatter and format - run: | - julia -e 'using Pkg; Pkg.add(PackageSpec(name="JuliaFormatter", version="v1"))' - julia -e 'using JuliaFormatter; format(".", verbose=true)' - - name: Format check - id: format - run: | - julia -e ' - out = Cmd(`git diff --name-only`) |> read |> String - if out == "" - exit(0) - else - @error "Some files have not been formatted !!!" - write(stdout, out) - exit(1) - end' - - name: Create pull request - if: ${{ failure() && steps.format.conclusion == 'failure' }} - id: cpr - uses: peter-evans/create-pull-request@v7 - with: - token: ${{ secrets.GITHUB_TOKEN }} - commit-message: Format .jl files - title: 'Automatic JuliaFormatter.jl run' - base: ${{ github.head_ref }} - branch: auto-juliaformatter-pr - delete-branch: true - labels: formatting, automated pr, no changelog +permissions: + contents: read + actions: write + pull-requests: write - - name: Check outputs - if: ${{ success() && steps.cpr.conclusion == 'success' }} - run: | - echo "Pull Request Number - ${{ steps.cpr.outputs.pull-request-number }}" - echo "Pull Request URL - ${{ steps.cpr.outputs.pull-request-url }}" +jobs: + formatcheck: + uses: "QuantumKitHub/QuantumKitHubActions/.github/workflows/FormatCheck.yml@main" diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml new file mode 100644 index 0000000..3e2823c --- /dev/null +++ b/.pre-commit-config.yaml @@ -0,0 +1,5 @@ +repos: + - repo: https://github.com/fredrikekre/runic-pre-commit + rev: v2.0.1 + hooks: + - id: runic diff --git a/docs/make.jl b/docs/make.jl index 3ff4160..937ed67 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -2,7 +2,7 @@ if Base.active_project() != joinpath(@__DIR__, "Project.toml") using Pkg Pkg.activate(@__DIR__) - Pkg.develop(PackageSpec(; path=(@__DIR__) * "/../")) + Pkg.develop(PackageSpec(; path = (@__DIR__) * "/../")) Pkg.resolve() Pkg.instantiate() end @@ -11,18 +11,24 @@ using Documenter using MPSKitModels makedocs(; - modules=[MPSKitModels], - sitename="MPSKitModels.jl", - authors="Maarten Vandamme", - format=Documenter.HTML(; - prettyurls=get(ENV, "CI", nothing) == "true", - mathengine=MathJax()), - pages=["Home" => "index.md", - "Manual" => ["man/operators.md", - "man/mpoham.md", - "man/lattices.md", - "man/models.md"], - "Index" => "package_index.md"], - checkdocs=:public) + modules = [MPSKitModels], + sitename = "MPSKitModels.jl", + authors = "Maarten Vandamme", + format = Documenter.HTML(; + prettyurls = get(ENV, "CI", nothing) == "true", + mathengine = MathJax() + ), + pages = [ + "Home" => "index.md", + "Manual" => [ + "man/operators.md", + "man/mpoham.md", + "man/lattices.md", + "man/models.md", + ], + "Index" => "package_index.md", + ], + checkdocs = :public +) -deploydocs(; repo="github.com/QuantumKitHub/MPSKitModels.jl.git") +deploydocs(; repo = "github.com/QuantumKitHub/MPSKitModels.jl.git") diff --git a/src/lattices/chains.jl b/src/lattices/chains.jl index c92573c..a3228dd 100644 --- a/src/lattices/chains.jl +++ b/src/lattices/chains.jl @@ -5,7 +5,7 @@ A one dimensional infinite lattice with a unit cell containing `L` sites. """ struct InfiniteChain <: AbstractLattice{1} L::Int - function InfiniteChain(L::Integer=1) + function InfiniteChain(L::Integer = 1) return L > 0 ? new(L) : error("period should be positive ($L)") end end @@ -19,7 +19,7 @@ A one-dimensional lattice of length `L` """ struct FiniteChain <: AbstractLattice{1} L::Int - function FiniteChain(L::Integer=1) + function FiniteChain(L::Integer = 1) return L > 0 ? new(L) : error("length should be positive ($L)") end end @@ -27,7 +27,7 @@ Base.axes(chain::FiniteChain) = (1:(chain.L),) Base.isfinite(::Type{FiniteChain}) = true Base.isfinite(::FiniteChain) = true -const Chain = Union{InfiniteChain,FiniteChain} +const Chain = Union{InfiniteChain, FiniteChain} vertices(chain::Chain) = LatticePoint.(1:(chain.L), Ref(chain)) nearest_neighbours(chain::InfiniteChain) = map(v -> v => v + 1, vertices(chain)) @@ -44,12 +44,12 @@ function bipartition(chain::Chain) end linearize_index(chain::InfiniteChain, i::Int) = i -linearize_index(chain::Chain, i::NTuple{1,Int}) = linearize_index(chain, i...) +linearize_index(chain::Chain, i::NTuple{1, Int}) = linearize_index(chain, i...) function linearize_index(chain::FiniteChain, i::Int) 0 < i <= chain.L || throw(BoundsError("lattice point out of bounds")) return i end -function LinearAlgebra.norm(p::LatticePoint{1,<:Union{FiniteChain,InfiniteChain}}) +function LinearAlgebra.norm(p::LatticePoint{1, <:Union{FiniteChain, InfiniteChain}}) return abs(p.coordinates[1]) end diff --git a/src/lattices/latticepoints.jl b/src/lattices/latticepoints.jl index 39917de..2d0b1c2 100644 --- a/src/lattices/latticepoints.jl +++ b/src/lattices/latticepoints.jl @@ -3,16 +3,16 @@ represents an `N`-dimensional point on a `G` lattice. """ -struct LatticePoint{N,G<:AbstractLattice{N}} - coordinates::NTuple{N,Int} +struct LatticePoint{N, G <: AbstractLattice{N}} + coordinates::NTuple{N, Int} lattice::G - function LatticePoint(coordinates::NTuple{N,Int}, lattice::AbstractLattice{N}) where {N} + function LatticePoint(coordinates::NTuple{N, Int}, lattice::AbstractLattice{N}) where {N} checkbounds(lattice, coordinates...) - return new{N,typeof(lattice)}(coordinates, lattice) + return new{N, typeof(lattice)}(coordinates, lattice) end end -function LatticePoint(ind::Int, lattice::G) where {G<:AbstractLattice{1}} +function LatticePoint(ind::Int, lattice::G) where {G <: AbstractLattice{1}} return LatticePoint((ind,), lattice) end @@ -21,47 +21,47 @@ function Base.show(io::IO, p::LatticePoint) end function Base.show(io::IO, ::MIME"text/plain", p::LatticePoint) - if get(io, :typeinfo, Any) === typeof(p) + return if get(io, :typeinfo, Any) === typeof(p) print(io, [p.coordinates...]) else print(io, p.lattice, [p.coordinates...]) end end -function Base.getindex(lattice::AbstractLattice{N}, inds::Vararg{Int,N}) where {N} +function Base.getindex(lattice::AbstractLattice{N}, inds::Vararg{Int, N}) where {N} return LatticePoint(inds, lattice) end Base.to_index(p::LatticePoint) = linearize_index(p) linearize_index(p::LatticePoint) = linearize_index(p.lattice, p.coordinates...) -function Base.:+(i::LatticePoint{N,G}, j::LatticePoint{N,G}) where {N,G} +function Base.:+(i::LatticePoint{N, G}, j::LatticePoint{N, G}) where {N, G} i.lattice == j.lattice || throw(ArgumentError("lattices should be equal")) return LatticePoint(i.coordinates .+ j.coordinates, i.lattice) end -function Base.:-(i::LatticePoint{N,G}, j::LatticePoint{N,G}) where {N,G} +function Base.:-(i::LatticePoint{N, G}, j::LatticePoint{N, G}) where {N, G} i.lattice == j.lattice || throw(ArgumentError("lattices should be equal")) return LatticePoint(i.coordinates .- j.coordinates, i.lattice) end -function Base.:+(i::LatticePoint{N}, j::NTuple{N,Int}) where {N} +function Base.:+(i::LatticePoint{N}, j::NTuple{N, Int}) where {N} return LatticePoint(i.coordinates .+ j, i.lattice) end Base.:+(i::LatticePoint{1}, j::Int) = LatticePoint(i.coordinates .+ j, i.lattice) -Base.:+(i::NTuple{N,Int}, j::LatticePoint{N}) where {N} = j + i +Base.:+(i::NTuple{N, Int}, j::LatticePoint{N}) where {N} = j + i Base.:+(i::Int, j::LatticePoint{1}) = j + i -function Base.:-(i::LatticePoint{N}, j::NTuple{N,Int}) where {N} +function Base.:-(i::LatticePoint{N}, j::NTuple{N, Int}) where {N} return LatticePoint(i.coordinates .- j, i.lattice) end Base.:-(i::LatticePoint{1}, j::Int) = LatticePoint(i.coordinates .- j, i.lattice) -function Base.:-(i::NTuple{N,Int}, j::LatticePoint{N}) where {N} +function Base.:-(i::NTuple{N, Int}, j::LatticePoint{N}) where {N} return LatticePoint(i .- j.coordinates, j.lattice) end Base.:-(i::Int, j::LatticePoint{1}) = LatticePoint(i .- j, j.lattice) -Base.isless(i::L, j::L) where {L<:LatticePoint} = linearize_index(i) < linearize_index(j) -function Base.isfinite(::Union{LatticePoint{N,G},Type{<:LatticePoint{N,G}}}) where {N,G} +Base.isless(i::L, j::L) where {L <: LatticePoint} = linearize_index(i) < linearize_index(j) +function Base.isfinite(::Union{LatticePoint{N, G}, Type{<:LatticePoint{N, G}}}) where {N, G} return isfinite(G) end -latticetype(::Union{LatticePoint{N,G},Type{<:LatticePoint{N,G}}}) where {N,G} = G +latticetype(::Union{LatticePoint{N, G}, Type{<:LatticePoint{N, G}}}) where {N, G} = G diff --git a/src/lattices/lattices.jl b/src/lattices/lattices.jl index fa69c8b..f0b8eb1 100644 --- a/src/lattices/lattices.jl +++ b/src/lattices/lattices.jl @@ -44,11 +44,13 @@ Base.length(L::AbstractLattice) = length(vertices(L)) Base.isfinite(L::AbstractLattice) = isfinite(typeof(L)) Base.iterate(L::AbstractLattice) = iterate(vertices(L)) -function Base.checkbounds(L::AbstractLattice{N}, inds::Vararg{Int,N}) where {N} +function Base.checkbounds(L::AbstractLattice{N}, inds::Vararg{Int, N}) where {N} return checkbounds(Bool, L, inds...) || throw(BoundsError(L, inds)) end -function Base.checkbounds(::Type{Bool}, L::AbstractLattice{N}, - inds::Vararg{Int,N}) where {N} +function Base.checkbounds( + ::Type{Bool}, L::AbstractLattice{N}, + inds::Vararg{Int, N} + ) where {N} return Base.checkbounds_indices(Bool, axes(L), inds) end diff --git a/src/lattices/snakepattern.jl b/src/lattices/snakepattern.jl index 2866799..b2fde49 100644 --- a/src/lattices/snakepattern.jl +++ b/src/lattices/snakepattern.jl @@ -3,7 +3,7 @@ Represents a given lattice with a linear order that is provided by `pattern`. """ -struct SnakePattern{N,G<:AbstractLattice{N},F} <: AbstractLattice{N} +struct SnakePattern{N, G <: AbstractLattice{N}, F} <: AbstractLattice{N} lattice::G pattern::F end @@ -11,7 +11,7 @@ end SnakePattern(lattice) = SnakePattern(lattice, identity) Base.axes(lattice::SnakePattern) = axes(lattice.lattice) -Base.isfinite(::Type{SnakePattern{N,G}}) where {N,G} = isfinite(G) +Base.isfinite(::Type{SnakePattern{N, G}}) where {N, G} = isfinite(G) function linearize_index(snake::SnakePattern, i...) return snake.pattern(linearize_index(snake.lattice, i...)) diff --git a/src/lattices/squarelattice.jl b/src/lattices/squarelattice.jl index 1bdb7b5..51a771e 100644 --- a/src/lattices/squarelattice.jl +++ b/src/lattices/squarelattice.jl @@ -8,7 +8,7 @@ This representes an `L` by `N÷L` rectangular patch. struct FiniteStrip <: AbstractLattice{2} L::Int N::Int - function FiniteStrip(L::Integer, N::Integer=L) + function FiniteStrip(L::Integer, N::Integer = L) N > 0 || throw(ArgumentError("period should be positive")) mod(N, L) == 0 || throw(ArgumentError("period should be a multiple of circumference")) @@ -27,7 +27,7 @@ An infinite strip with `L` sites per rung and `N` sites per unit cell. struct InfiniteStrip <: AbstractLattice{2} L::Int N::Int - function InfiniteStrip(L::Integer, N::Integer=L) + function InfiniteStrip(L::Integer, N::Integer = L) N > 0 || throw(ArgumentError("period should be positive")) mod(N, L) == 0 || throw(ArgumentError("period should be a multiple of circumference")) @@ -46,7 +46,7 @@ A cylinder with circumference `L` and `N` sites in total. struct FiniteCylinder <: AbstractLattice{2} L::Int N::Int - function FiniteCylinder(L::Integer, N::Integer=L) + function FiniteCylinder(L::Integer, N::Integer = L) N > 0 || throw(ArgumentError("period should be positive")) mod(N, L) == 0 || throw(ArgumentError("period should be a multiple of circumference")) @@ -67,7 +67,7 @@ An infinite cylinder with `L` sites per rung and `N` sites per unit cell. struct InfiniteCylinder <: AbstractLattice{2} L::Int N::Int - function InfiniteCylinder(L::Integer, N::Integer=L) + function InfiniteCylinder(L::Integer, N::Integer = L) N > 0 || throw(ArgumentError("period should be positive")) mod(N, L) == 0 || throw(ArgumentError("period should be a multiple of circumference")) @@ -86,7 +86,7 @@ A finite helix with `L` sites per rung and `N` sites in total. struct FiniteHelix <: AbstractLattice{2} L::Int N::Int - function FiniteHelix(L::Integer, N::Integer=1) + function FiniteHelix(L::Integer, N::Integer = 1) N > 0 || error("period should be positive") return new(L, N) end @@ -103,7 +103,7 @@ An infinite helix with `L` sites per rung and `N` sites per unit cell. struct InfiniteHelix <: AbstractLattice{2} L::Int N::Int - function InfiniteHelix(L::Integer, N::Integer=1) + function InfiniteHelix(L::Integer, N::Integer = 1) N > 0 || error("period should be positive") return new(L, N) end @@ -140,49 +140,63 @@ function linearize_index(helix::InfiniteHelix, i::Int, j::Int) return mod1(i, helix.L) + helix.L * (j + (i - 1) ÷ helix.L - 1) end -function vertices(lattice::Union{FiniteStrip,InfiniteStrip,FiniteCylinder,InfiniteCylinder}) - return (LatticePoint((i, j), lattice) for i in 1:(lattice.L), - j in 1:(lattice.N ÷ lattice.L)) +function vertices(lattice::Union{FiniteStrip, InfiniteStrip, FiniteCylinder, InfiniteCylinder}) + return ( + LatticePoint((i, j), lattice) for i in 1:(lattice.L), + j in 1:(lattice.N ÷ lattice.L) + ) end -function vertices(lattice::Union{FiniteHelix,InfiniteHelix}) +function vertices(lattice::Union{FiniteHelix, InfiniteHelix}) return (LatticePoint((i, 1), lattice) for i in 1:(lattice.N)) end function nearest_neighbours(lattice::FiniteStrip) rows = lattice.L cols = lattice.N ÷ lattice.L - horizontal = (LatticePoint((i, j), lattice) => LatticePoint((i, j + 1), lattice) - for i in 1:rows, j in 1:(cols - 1)) - vertical = (LatticePoint((i, j), lattice) => LatticePoint((i + 1, j), lattice) - for i in 1:(rows - 1), j in 1:cols) + horizontal = ( + LatticePoint((i, j), lattice) => LatticePoint((i, j + 1), lattice) + for i in 1:rows, j in 1:(cols - 1) + ) + vertical = ( + LatticePoint((i, j), lattice) => LatticePoint((i + 1, j), lattice) + for i in 1:(rows - 1), j in 1:cols + ) return [horizontal..., vertical...] end function nearest_neighbours(lattice::FiniteCylinder) rows = lattice.L cols = lattice.N ÷ lattice.L - horizontal = (LatticePoint((i, j), lattice) => LatticePoint((i, j + 1), lattice) - for i in 1:rows, j in 1:(cols - 1)) - vertical = (LatticePoint((i, j), lattice) => LatticePoint((i + 1, j), lattice) - for i in 1:rows, j in 1:cols) + horizontal = ( + LatticePoint((i, j), lattice) => LatticePoint((i, j + 1), lattice) + for i in 1:rows, j in 1:(cols - 1) + ) + vertical = ( + LatticePoint((i, j), lattice) => LatticePoint((i + 1, j), lattice) + for i in 1:rows, j in 1:cols + ) return [horizontal..., vertical...] end function nearest_neighbours(lattice::FiniteHelix) rows = lattice.L cols = lattice.N ÷ lattice.L - horizontal = (LatticePoint((i, j), lattice) => LatticePoint((i, j + 1), lattice) - for i in 1:rows, j in 1:(cols - 1)) - vertical = (LatticePoint((i, j), lattice) => LatticePoint((i + 1, j), lattice) - for i in 1:rows, j in 1:cols if (i != rows && j != cols)) + horizontal = ( + LatticePoint((i, j), lattice) => LatticePoint((i, j + 1), lattice) + for i in 1:rows, j in 1:(cols - 1) + ) + vertical = ( + LatticePoint((i, j), lattice) => LatticePoint((i + 1, j), lattice) + for i in 1:rows, j in 1:cols if (i != rows && j != cols) + ) return [horizontal..., vertical...] end -function nearest_neighbours(lattice::Union{InfiniteStrip,InfiniteCylinder,InfiniteHelix}) +function nearest_neighbours(lattice::Union{InfiniteStrip, InfiniteCylinder, InfiniteHelix}) V = vertices(lattice) - neighbours = Pair{eltype(V),eltype(V)}[] + neighbours = Pair{eltype(V), eltype(V)}[] for v in V push!(neighbours, v => v + (0, 1)) if v.coordinates[1] < lattice.L || - lattice isa InfiniteCylinder || - lattice isa InfiniteHelix + lattice isa InfiniteCylinder || + lattice isa InfiniteHelix push!(neighbours, v => v + (1, 0)) end end @@ -190,22 +204,27 @@ function nearest_neighbours(lattice::Union{InfiniteStrip,InfiniteCylinder,Infini end function next_nearest_neighbours(lattice::AbstractLattice{2}) - diag1 = (i => i + (1, 1) for i in vertices(lattice) if checkbounds(Bool, lattice, - (i.coordinates .+ - (1, 1))...)) - diag2 = (i => i + (1, -1) for i in vertices(lattice) if checkbounds(Bool, lattice, - (i.coordinates .+ - (1, -1))...)) + diag1 = ( + i => i + (1, 1) for i in vertices(lattice) if checkbounds( + Bool, lattice, (i.coordinates .+ (1, 1))... + ) + ) + diag2 = ( + i => i + (1, -1) for i in vertices(lattice) if checkbounds( + Bool, lattice, (i.coordinates .+ (1, -1))... + ) + ) return [diag1..., diag2...] end -function next_nearest_neighbours(lattice::Union{InfiniteStrip,InfiniteCylinder, - InfiniteHelix}) +function next_nearest_neighbours( + lattice::Union{InfiniteStrip, InfiniteCylinder, InfiniteHelix} + ) V = vertices(lattice) - neighbours = Pair{eltype(V),eltype(V)}[] + neighbours = Pair{eltype(V), eltype(V)}[] for v in V if v.coordinates[1] < lattice.L || - lattice isa InfiniteCylinder || - lattice isa InfiniteHelix + lattice isa InfiniteCylinder || + lattice isa InfiniteHelix push!(neighbours, v => v + (1, 1)) end if v.coordinates[1] > 1 || lattice isa InfiniteCylinder || lattice isa InfiniteHelix @@ -215,12 +234,14 @@ function next_nearest_neighbours(lattice::Union{InfiniteStrip,InfiniteCylinder, return neighbours end -LinearAlgebra.norm(p::LatticePoint{2,InfiniteStrip}) = LinearAlgebra.norm(p.coordinates) -function LinearAlgebra.norm(p::LatticePoint{2,InfiniteCylinder}) - return min(sqrt(mod(p.coordinates[1], p.lattice.L)^2 + p.coordinates[2]^2), - sqrt(mod(-p.coordinates[1], p.lattice.L)^2 + p.coordinates[2]^2)) +LinearAlgebra.norm(p::LatticePoint{2, InfiniteStrip}) = LinearAlgebra.norm(p.coordinates) +function LinearAlgebra.norm(p::LatticePoint{2, InfiniteCylinder}) + return min( + sqrt(mod(p.coordinates[1], p.lattice.L)^2 + p.coordinates[2]^2), + sqrt(mod(-p.coordinates[1], p.lattice.L)^2 + p.coordinates[2]^2) + ) end -function LinearAlgebra.norm(p::LatticePoint{2,InfiniteHelix}) +function LinearAlgebra.norm(p::LatticePoint{2, InfiniteHelix}) x₁ = mod(p.coordinates[1], p.lattice.L) y₁ = p.coordinates[2] + (p.coordinates[1] ÷ p.lattice.L) x₂ = mod(-p.coordinates[1], p.lattice.L) diff --git a/src/lattices/triangularlattice.jl b/src/lattices/triangularlattice.jl index c942b2a..1175ec6 100644 --- a/src/lattices/triangularlattice.jl +++ b/src/lattices/triangularlattice.jl @@ -7,7 +7,7 @@ The y-axis is aligned along an edge of the hexagons, and the circumference is `` struct HoneycombYC <: AbstractLattice{2} L::Int N::Int - function HoneycombYC(L::Integer, N::Integer=L) + function HoneycombYC(L::Integer, N::Integer = L) (L > 0 && N > 0) || throw(ArgumentError("period and length should be strictly positive")) mod(L, 4) == 0 || throw(ArgumentError("period must be a multiple of 4")) @@ -20,11 +20,11 @@ end Base.isfinite(::Type{HoneycombYC}) = false # TODO: do proper boundscheck -function Base.checkbounds(::Type{Bool}, lattice::HoneycombYC, inds::Vararg{Int,2}) +function Base.checkbounds(::Type{Bool}, lattice::HoneycombYC, inds::Vararg{Int, 2}) return true end -function LinearAlgebra.norm(p::LatticePoint{2,HoneycombYC}) +function LinearAlgebra.norm(p::LatticePoint{2, HoneycombYC}) x = p.coordinates[1] + p.coordinates[2] * cos(2π / 6) y = p.coordinates[2] * sin(2π / 6) x₁ = mod(x, 3p.lattice.L ÷ 4) @@ -65,7 +65,7 @@ end function nearest_neighbours(lattice::HoneycombYC) V = vertices(lattice) - neighbours = Pair{eltype(V),eltype(V)}[] + neighbours = Pair{eltype(V), eltype(V)}[] for v in V I = linearize_index(v) if mod(I, 4) == 1 diff --git a/src/models/hamiltonians.jl b/src/models/hamiltonians.jl index 087490d..7d36c28 100644 --- a/src/models/hamiltonians.jl +++ b/src/models/hamiltonians.jl @@ -22,37 +22,41 @@ function transverse_field_ising end function transverse_field_ising(lattice::AbstractLattice; kwargs...) return transverse_field_ising(ComplexF64, Trivial, lattice; kwargs...) end -function transverse_field_ising(S::Type{<:Sector}, - lattice::AbstractLattice=InfiniteChain(1); - kwargs...) +function transverse_field_ising( + S::Type{<:Sector}, lattice::AbstractLattice = InfiniteChain(1); + kwargs... + ) return transverse_field_ising(ComplexF64, S, lattice; kwargs...) end function transverse_field_ising(T::Type{<:Number}, lattice::AbstractLattice; kwargs...) return transverse_field_ising(T, Trivial, lattice; kwargs...) end -function transverse_field_ising(T::Type{<:Number}, - S::Type{<:Sector}, - lattice::AbstractLattice; - kwargs...) +function transverse_field_ising( + T::Type{<:Number}, S::Type{<:Sector}, lattice::AbstractLattice; + kwargs... + ) throw(ArgumentError("`symmetry` must be either `Trivial`, `Z2Irrep` or `FermionParity`")) end -function transverse_field_ising(T::Type{<:Number}=ComplexF64, - S::Union{Type{Trivial},Type{Z2Irrep}}=Trivial, - lattice::AbstractLattice=InfiniteChain(1); - J=1.0, g=1.0) +function transverse_field_ising( + T::Type{<:Number} = ComplexF64, S::Union{Type{Trivial}, Type{Z2Irrep}} = Trivial, + lattice::AbstractLattice = InfiniteChain(1); + J = 1.0, g = 1.0 + ) ZZ = scale!(σᶻᶻ(T, S), -J) X = scale!(σˣ(T, S), g * -J) return @mpoham begin sum(nearest_neighbours(lattice)) do (i, j) - return ZZ{i,j} + return ZZ{i, j} end + sum(vertices(lattice)) do i return X{i} end end end -function transverse_field_ising(T::Type{<:Number}, ::Type{fℤ₂}, - lattice::AbstractLattice=InfiniteChain(1); - J=1.0, g=1.0) +function transverse_field_ising( + T::Type{<:Number}, ::Type{fℤ₂}, + lattice::AbstractLattice = InfiniteChain(1); + J = 1.0, g = 1.0 + ) hop = add!(c_plusmin(T), c_minplus(T)) sc = add!(c_plusplus(T), c_minmin(T)) twosite = add!(hop, sc, J, -J) @@ -63,7 +67,7 @@ function transverse_field_ising(T::Type{<:Number}, ::Type{fℤ₂}, return @mpoham begin sum(nearest_neighbours(lattice)) do (i, j) - return twosite{i,j} + return twosite{i, j} end + sum(vertices(lattice)) do i return onesite{i} end @@ -89,16 +93,17 @@ function kitaev_model end function kitaev_model(lattice::AbstractLattice; kwargs...) return kitaev_model(ComplexF64, lattice; kwargs...) end -function kitaev_model(elt::Type{<:Number}=ComplexF64, - lattice::AbstractLattice=InfiniteChain(1); - t=1.0, mu=1.0, Delta=1.0) +function kitaev_model( + elt::Type{<:Number} = ComplexF64, lattice::AbstractLattice = InfiniteChain(1); + t = 1.0, mu = 1.0, Delta = 1.0 + ) TB = rmul!(c_plusmin(elt) + c_minplus(elt), -t / 2) # tight-binding term SC = rmul!(c_plusplus(elt) + c_minmin(elt), Delta / 2) # superconducting term CP = rmul!(c_number(elt), -mu) # chemical potential term return @mpoham begin sum(nearest_neighbours(lattice)) do (i, j) - return (TB + SC){i,j} + return (TB + SC){i, j} end + sum(vertices(lattice)) do i return CP{i} end @@ -127,20 +132,23 @@ function heisenberg_XXX end function heisenberg_XXX(lattice::AbstractLattice; kwargs...) return heisenberg_XXX(ComplexF64, Trivial, lattice; kwargs...) end -function heisenberg_XXX(symmetry::Type{<:Sector}, lattice::AbstractLattice=InfiniteChain(1); - kwargs...) +function heisenberg_XXX( + symmetry::Type{<:Sector}, lattice::AbstractLattice = InfiniteChain(1); + kwargs... + ) return heisenberg_XXX(ComplexF64, symmetry, lattice; kwargs...) end function heisenberg_XXX(elt::Type{<:Number}, lattice::AbstractLattice; kwargs...) return heisenberg_XXX(elt, Trivial, lattice; kwargs...) end -function heisenberg_XXX(T::Type{<:Number}=ComplexF64, - symmetry::Type{<:Sector}=Trivial, - lattice::AbstractLattice=InfiniteChain(1); - J::Real=1.0, spin::Real=1) - term = rmul!(S_exchange(T, symmetry; spin=spin), J) +function heisenberg_XXX( + T::Type{<:Number} = ComplexF64, symmetry::Type{<:Sector} = Trivial, + lattice::AbstractLattice = InfiniteChain(1); + J::Real = 1.0, spin::Real = 1 + ) + term = rmul!(S_exchange(T, symmetry; spin = spin), J) return @mpoham sum(nearest_neighbours(lattice)) do (i, j) - return term{i,j} + return term{i, j} end end @@ -159,22 +167,25 @@ function heisenberg_XXZ end function heisenberg_XXZ(lattice::AbstractLattice; kwargs...) return heisenberg_XXZ(ComplexF64, Trivial, lattice; kwargs...) end -function heisenberg_XXZ(symmetry::Type{<:Sector}, lattice::AbstractLattice=InfiniteChain(1); - kwargs...) +function heisenberg_XXZ( + symmetry::Type{<:Sector}, lattice::AbstractLattice = InfiniteChain(1); + kwargs... + ) return heisenberg_XXZ(ComplexF64, symmetry, lattice; kwargs...) end function heisenberg_XXZ(elt::Type{<:Number}, lattice::AbstractLattice; kwargs...) return heisenberg_XXZ(elt, Trivial, lattice; kwargs...) end -function heisenberg_XXZ(elt::Type{<:Number}=ComplexF64, - symmetry::Type{<:Sector}=Trivial, - lattice::AbstractLattice=InfiniteChain(1); - J=1.0, Delta=1.0, spin=1) - term = rmul!(S_xx(elt, symmetry; spin=spin), J) + - rmul!(S_yy(elt, symmetry; spin=spin), J) + - rmul!(S_zz(elt, symmetry; spin=spin), Delta * J) +function heisenberg_XXZ( + elt::Type{<:Number} = ComplexF64, symmetry::Type{<:Sector} = Trivial, + lattice::AbstractLattice = InfiniteChain(1); + J = 1.0, Delta = 1.0, spin = 1 + ) + term = rmul!(S_xx(elt, symmetry; spin = spin), J) + + rmul!(S_yy(elt, symmetry; spin = spin), J) + + rmul!(S_zz(elt, symmetry; spin = spin), Delta * J) return @mpoham sum(nearest_neighbours(lattice)) do (i, j) - return term{i,j} + return term{i, j} end end @@ -193,14 +204,15 @@ function heisenberg_XYZ end function heisenberg_XYZ(lattice::AbstractLattice; kwargs...) return heisenberg_XYZ(ComplexF64, lattice; kwargs...) end -function heisenberg_XYZ(T::Type{<:Number}=ComplexF64, - lattice::AbstractLattice=InfiniteChain(1); - Jx=1.0, Jy=1.0, Jz=1.0, spin=1) - term = rmul!(S_xx(T, Trivial; spin=spin), Jx) + - rmul!(S_yy(T, Trivial; spin=spin), Jy) + - rmul!(S_zz(T, Trivial; spin=spin), Jz) +function heisenberg_XYZ( + T::Type{<:Number} = ComplexF64, lattice::AbstractLattice = InfiniteChain(1); + Jx = 1.0, Jy = 1.0, Jz = 1.0, spin = 1 + ) + term = rmul!(S_xx(T, Trivial; spin = spin), Jx) + + rmul!(S_yy(T, Trivial; spin = spin), Jy) + + rmul!(S_zz(T, Trivial; spin = spin), Jz) return @mpoham sum(nearest_neighbours(lattice)) do (i, j) - return term{i,j} + return term{i, j} end end @@ -220,21 +232,26 @@ function bilinear_biquadratic_model end function bilinear_biquadratic_model(lattice::AbstractLattice; kwargs...) return bilinear_biquadratic_model(ComplexF64, Trivial, lattice; kwargs...) end -function bilinear_biquadratic_model(symmetry::Type{<:Sector}, - lattice::AbstractLattice=InfiniteChain(1); kwargs...) +function bilinear_biquadratic_model( + symmetry::Type{<:Sector}, lattice::AbstractLattice = InfiniteChain(1); + kwargs... + ) return bilinear_biquadratic_model(ComplexF64, symmetry, lattice; kwargs...) end -function bilinear_biquadratic_model(elt::Type{<:Number}, lattice::AbstractLattice; - kwargs...) +function bilinear_biquadratic_model( + elt::Type{<:Number}, lattice::AbstractLattice; + kwargs... + ) return bilinear_biquadratic_model(elt, Trivial, lattice; kwargs...) end -function bilinear_biquadratic_model(elt::Type{<:Number}=ComplexF64, - symmetry::Type{<:Sector}=Trivial, - lattice::AbstractLattice=InfiniteChain(1); - spin=1, J=1.0, θ=0.0) +function bilinear_biquadratic_model( + elt::Type{<:Number} = ComplexF64, symmetry::Type{<:Sector} = Trivial, + lattice::AbstractLattice = InfiniteChain(1); + spin = 1, J = 1.0, θ = 0.0 + ) return @mpoham sum(nearest_neighbours(lattice)) do (i, j) - return J * cos(θ) * S_exchange(elt, symmetry; spin=spin){i,j} + - J * sin(θ) * (S_exchange(elt, symmetry; spin=spin)^2){i,j} + return J * cos(θ) * S_exchange(elt, symmetry; spin = spin){i, j} + + J * sin(θ) * (S_exchange(elt, symmetry; spin = spin)^2){i, j} end end @@ -257,23 +274,30 @@ function quantum_potts end function quantum_potts(lattice::AbstractLattice; kwargs...) return quantum_potts(ComplexF64, Trivial, lattice; kwargs...) end -function quantum_potts(symmetry::Type{<:Sector}, - lattice::AbstractLattice=InfiniteChain(1); kwargs...) +function quantum_potts( + symmetry::Type{<:Sector}, lattice::AbstractLattice = InfiniteChain(1); + kwargs... + ) return quantum_potts(ComplexF64, symmetry, lattice; kwargs...) end -function quantum_potts(elt::Type{<:Number}, lattice::AbstractLattice; - kwargs...) +function quantum_potts( + elt::Type{<:Number}, lattice::AbstractLattice; + kwargs... + ) return quantum_potts(elt, Trivial, lattice; kwargs...) end -function quantum_potts(elt::Type{<:Number}=ComplexF64, - symmetry::Type{<:Sector}=Trivial, - lattice::AbstractLattice=InfiniteChain(1); - q=3, J=1.0, g=1.0) - return @mpoham sum(sum(nearest_neighbours(lattice)) do (i, j) - return -J * (potts_ZZ(elt, symmetry; q)^k){i,j} - end - sum(vertices(lattice)) do i - return g * (potts_field(elt, symmetry; q)^k){i} - end for k in 1:(q - 1)) +function quantum_potts( + elt::Type{<:Number} = ComplexF64, symmetry::Type{<:Sector} = Trivial, + lattice::AbstractLattice = InfiniteChain(1); + q = 3, J = 1.0, g = 1.0 + ) + return @mpoham sum( + sum(nearest_neighbours(lattice)) do (i, j) + return -J * (potts_ZZ(elt, symmetry; q)^k){i, j} + end - sum(vertices(lattice)) do i + return g * (potts_field(elt, symmetry; q)^k){i} + end for k in 1:(q - 1) + ) end #=========================================================================================== @@ -298,27 +322,29 @@ function hubbard_model end function hubbard_model(lattice::AbstractLattice; kwargs...) return hubbard_model(ComplexF64, Trivial, Trivial, lattice; kwargs...) end -function hubbard_model(particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; - kwargs...) +function hubbard_model( + particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; + kwargs... + ) return hubbard_model(ComplexF64, particle_symmetry, spin_symmetry; kwargs...) end function hubbard_model(elt::Type{<:Number}, lattice::AbstractLattice; kwargs...) return hubbard_model(elt, Trivial, Trivial, lattice; kwargs...) end -function hubbard_model(T::Type{<:Number}=ComplexF64, - particle_symmetry::Type{<:Sector}=Trivial, - spin_symmetry::Type{<:Sector}=Trivial, - lattice::AbstractLattice=InfiniteChain(1); - t=1.0, U=1.0, mu=0.0, n::Integer=0) +function hubbard_model( + T::Type{<:Number} = ComplexF64, particle_symmetry::Type{<:Sector} = Trivial, + spin_symmetry::Type{<:Sector} = Trivial, lattice::AbstractLattice = InfiniteChain(1); + t = 1.0, U = 1.0, mu = 0.0, n::Integer = 0 + ) hopping = e⁺e⁻(T, particle_symmetry, spin_symmetry) + - e⁻e⁺(T, particle_symmetry, spin_symmetry) + e⁻e⁺(T, particle_symmetry, spin_symmetry) interaction_term = nꜛnꜜ(T, particle_symmetry, spin_symmetry) N = e_number(T, particle_symmetry, spin_symmetry) return @mpoham begin sum(nearest_neighbours(lattice)) do (i, j) - return -t * hopping{i,j} + return -t * hopping{i, j} end + - sum(vertices(lattice)) do i + sum(vertices(lattice)) do i return U * interaction_term{i} - mu * N{i} end end @@ -344,24 +370,27 @@ function bose_hubbard_model end function bose_hubbard_model(lattice::AbstractLattice; kwargs...) return bose_hubbard_model(ComplexF64, Trivial, lattice; kwargs...) end -function bose_hubbard_model(symmetry::Type{<:Sector}, - lattice::AbstractLattice=InfiniteChain(1); kwargs...) +function bose_hubbard_model( + symmetry::Type{<:Sector}, lattice::AbstractLattice = InfiniteChain(1); + kwargs... + ) return bose_hubbard_model(ComplexF64, symmetry, lattice; kwargs...) end -function bose_hubbard_model(elt::Type{<:Number}=ComplexF64, - symmetry::Type{<:Sector}=Trivial, - lattice::AbstractLattice=InfiniteChain(1); - cutoff::Integer=5, t=1.0, U=1.0, mu=0.0, n=0) - hopping_term = a_plusmin(elt, symmetry; cutoff=cutoff) + - a_minplus(elt, symmetry; cutoff=cutoff) - N = a_number(elt, symmetry; cutoff=cutoff) +function bose_hubbard_model( + elt::Type{<:Number} = ComplexF64, symmetry::Type{<:Sector} = Trivial, + lattice::AbstractLattice = InfiniteChain(1); + cutoff::Integer = 5, t = 1.0, U = 1.0, mu = 0.0, n = 0 + ) + hopping_term = a_plusmin(elt, symmetry; cutoff = cutoff) + + a_minplus(elt, symmetry; cutoff = cutoff) + N = a_number(elt, symmetry; cutoff = cutoff) interaction_term = contract_onesite(N, N - id(domain(N))) H = @mpoham begin sum(nearest_neighbours(lattice)) do (i, j) - return -t * hopping_term{i,j} + return -t * hopping_term{i, j} end + - sum(vertices(lattice)) do i + sum(vertices(lattice)) do i return U / 2 * interaction_term{i} - mu * N{i} end end @@ -401,27 +430,31 @@ function tj_model end function tj_model(lattice::AbstractLattice; kwargs...) return tj_model(ComplexF64, Trivial, Trivial, lattice; kwargs...) end -function tj_model(particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; - kwargs...) +function tj_model( + particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; + kwargs... + ) return tj_model(ComplexF64, particle_symmetry, spin_symmetry; kwargs...) end function tj_model(elt::Type{<:Number}, lattice::AbstractLattice; kwargs...) return tj_model(elt, Trivial, Trivial, lattice; kwargs...) end -function tj_model(T::Type{<:Number}=ComplexF64, - particle_symmetry::Type{<:Sector}=Trivial, - spin_symmetry::Type{<:Sector}=Trivial, - lattice::AbstractLattice=InfiniteChain(1); - t=2.5, J=1.0, mu=0.0, slave_fermion::Bool=false) +function tj_model( + T::Type{<:Number} = ComplexF64, particle_symmetry::Type{<:Sector} = Trivial, + spin_symmetry::Type{<:Sector} = Trivial, lattice::AbstractLattice = InfiniteChain(1); + t = 2.5, J = 1.0, mu = 0.0, slave_fermion::Bool = false + ) hopping = TJOperators.e_plusmin(T, particle_symmetry, spin_symmetry; slave_fermion) + - TJOperators.e_minplus(T, particle_symmetry, spin_symmetry; slave_fermion) + TJOperators.e_minplus(T, particle_symmetry, spin_symmetry; slave_fermion) num = TJOperators.e_number(T, particle_symmetry, spin_symmetry; slave_fermion) - heisenberg = TJOperators.S_exchange(T, particle_symmetry, spin_symmetry; - slave_fermion) - - (1 / 4) * (num ⊗ num) + heisenberg = TJOperators.S_exchange( + T, particle_symmetry, spin_symmetry; + slave_fermion + ) - + (1 / 4) * (num ⊗ num) return @mpoham begin sum(nearest_neighbours(lattice)) do (i, j) - return (-t) * hopping{i,j} + J * heisenberg{i,j} + return (-t) * hopping{i, j} + J * heisenberg{i, j} end + sum(vertices(lattice)) do i return (-mu) * num{i} end diff --git a/src/models/quantum_chemistry.jl b/src/models/quantum_chemistry.jl index 4291b02..de0dd0f 100644 --- a/src/models/quantum_chemistry.jl +++ b/src/models/quantum_chemistry.jl @@ -29,29 +29,33 @@ where ``s`` and ``t`` are spin indices, which can be ``\\uparrow`` or ``\\downar - MPSKit does not contain many required algorithms in quantum chemistry (orbital ordering/optimization) - MPOHamiltonian is not well suited for quantum chemistry """ -function quantum_chemistry_hamiltonian(E0, K, V, Elt=ComplexF64) +function quantum_chemistry_hamiltonian(E0, K, V, Elt = ComplexF64) return mapped_quantum_chemistry_hamiltonian(E0, K, V, Elt)[1] end -function mapped_quantum_chemistry_hamiltonian(E0, K, V, Elt=ComplexF64) +function mapped_quantum_chemistry_hamiltonian(E0, K, V, Elt = ComplexF64) basis_size = size(K, 1) half_basis_size = Int(ceil(basis_size / 2)) # the phsyical space - psp = Vect[(Irrep[U₁] ⊠ Irrep[SU₂] ⊠ FermionParity)]((0, 0, 0) => 1, - (1, 1 // 2, 1) => 1, - (2, 0, 0) => 1) - - ap = TensorMap(ones, Elt, - psp * - Vect[(Irrep[U₁] ⊠ Irrep[SU₂] ⊠ FermionParity)]((-1, 1 // 2, 1) => 1), - psp) + psp = Vect[(Irrep[U₁] ⊠ Irrep[SU₂] ⊠ FermionParity)]( + (0, 0, 0) => 1, + (1, 1 // 2, 1) => 1, + (2, 0, 0) => 1 + ) + + ap = TensorMap( + ones, Elt, + psp * Vect[(Irrep[U₁] ⊠ Irrep[SU₂] ⊠ FermionParity)]((-1, 1 // 2, 1) => 1), + psp + ) blocks(ap)[(U₁(0) ⊠ SU₂(0) ⊠ FermionParity(0))] .*= -sqrt(2) blocks(ap)[(U₁(1) ⊠ SU₂(1 // 2) ⊠ FermionParity(1))] .*= 1 - bm = TensorMap(ones, Elt, psp, - Vect[(Irrep[U₁] ⊠ Irrep[SU₂] ⊠ FermionParity)]((-1, 1 // 2, 1) => 1) * - psp) + bm = TensorMap( + ones, Elt, psp, + Vect[(Irrep[U₁] ⊠ Irrep[SU₂] ⊠ FermionParity)]((-1, 1 // 2, 1) => 1) * psp + ) blocks(bm)[(U₁(0) ⊠ SU₂(0) ⊠ FermionParity(0))] .*= sqrt(2) blocks(bm)[(U₁(1) ⊠ SU₂(1 // 2) ⊠ FermionParity(1))] .*= -1 @@ -63,17 +67,17 @@ function mapped_quantum_chemistry_hamiltonian(E0, K, V, Elt=ComplexF64) @plansor b_derp[-1 -2; -3] := bp[1; 2 -2] * τ[-3 -1; 2 1] Lmap_ap_to_bp = inv(ap' * ap) * ap' * b_derp - @assert norm(ap * Lmap_ap_to_bp - b_derp) < 1e-12 + @assert norm(ap * Lmap_ap_to_bp - b_derp) < 1.0e-12 @plansor b_derp[-1 -2; -3] := bm[1; 2 -2] * τ[-3 -1; 2 1] Lmap_am_to_bm = inv(am' * am) * am' * b_derp - @assert norm(am * Lmap_am_to_bm - b_derp) < 1e-12 + @assert norm(am * Lmap_am_to_bm - b_derp) < 1.0e-12 Rmap_bp_to_ap = transpose(Lmap_am_to_bm', (2,), (1,)) Rmap_bm_to_am = transpose(Lmap_ap_to_bp', (2,), (1,)) @plansor a_derp[-1 -2; -3] := bp[1; -1 2] * Rmap_bp_to_ap[1; 3] * τ[2 3; -3 -2] - @assert norm(a_derp - ap) < 1e-12 + @assert norm(a_derp - ap) < 1.0e-12 @plansor a_derp[-1 -2; -3] := bm[1; -1 2] * Rmap_bm_to_am[1; 3] * τ[2 3; -3 -2] - @assert norm(a_derp - am) < 1e-12 + @assert norm(a_derp - am) < 1.0e-12 h_pm = TensorMap(ones, Elt, psp, psp) blocks(h_pm)[(U₁(0) ⊠ SU₂(0) ⊠ FermionParity(0))] .= 0 @@ -83,12 +87,12 @@ function mapped_quantum_chemistry_hamiltonian(E0, K, V, Elt=ComplexF64) @plansor o_derp[-1 -2; -3 -4] := am[-1 1; -3] * ap[1 -2; -4] h_pm_derp = transpose(h_pm, (2, 1), ()) Lmap_apam_to_pm = inv(o_derp' * o_derp) * o_derp' * h_pm_derp - @assert norm(o_derp * Lmap_apam_to_pm - h_pm_derp) < 1e-12 + @assert norm(o_derp * Lmap_apam_to_pm - h_pm_derp) < 1.0e-12 @plansor o_derp[-1 -2; -3 -4] := bm[-1; -3 1] * bp[-2; 1 -4] h_pm_derp2 = transpose(h_pm, (), (2, 1)) Rmap_bpbm_to_pm = h_pm_derp2 * o_derp' * inv(o_derp * o_derp') - @assert norm(transpose(h_pm, (), (2, 1)) - Rmap_bpbm_to_pm * o_derp) < 1e-12 + @assert norm(transpose(h_pm, (), (2, 1)) - Rmap_bpbm_to_pm * o_derp) < 1.0e-12 h_ppmm = h_pm * h_pm - h_pm @@ -118,7 +122,7 @@ function mapped_quantum_chemistry_hamiltonian(E0, K, V, Elt=ComplexF64) indmap_2R[pm1, end - j + 1, pm2, end - i + 1] = cnt end - hamdat = convert(Array{Any,3}, fill(EmptyVal(), basis_size, cnt + 1, cnt + 1))#Array{Any,3}(missing,basis_size,cnt+2,cnt+2); + hamdat = convert(Array{Any, 3}, fill(EmptyVal(), basis_size, cnt + 1, cnt + 1)) #Array{Any,3}(missing,basis_size,cnt+2,cnt+2); hamdat[:, 1, 1] .+= Elt(1) hamdat[:, end, end] .+= Elt(1) @@ -157,32 +161,28 @@ function mapped_quantum_chemistry_hamiltonian(E0, K, V, Elt=ComplexF64) # indmap_2 onsite part # we need pp, mm, pm - pp_f = isometry(fuse(_lastspace(ap)' * _lastspace(ap)'), - _lastspace(ap)' * _lastspace(ap)') - mm_f = isometry(fuse(_lastspace(am)' * _lastspace(am)'), - _lastspace(am)' * _lastspace(am)') - mp_f = isometry(fuse(_lastspace(am)' * _lastspace(ap)'), - _lastspace(am)' * _lastspace(ap)') - pm_f = isometry(fuse(_lastspace(ap)' * _lastspace(am)'), - _lastspace(ap)' * _lastspace(am)') + pp_f = isometry(fuse(_lastspace(ap)' * _lastspace(ap)'), _lastspace(ap)' * _lastspace(ap)') + mm_f = isometry(fuse(_lastspace(am)' * _lastspace(am)'), _lastspace(am)' * _lastspace(am)') + mp_f = isometry(fuse(_lastspace(am)' * _lastspace(ap)'), _lastspace(am)' * _lastspace(ap)') + pm_f = isometry(fuse(_lastspace(ap)' * _lastspace(am)'), _lastspace(ap)' * _lastspace(am)') @plansor ut_apap[-1 -2; -3 -4] := ut[-1] * ap[-3 1; 3] * ap[1 -2; 4] * - conj(pp_f[-4; 3 4]) + conj(pp_f[-4; 3 4]) @plansor ut_amam[-1 -2; -3 -4] := ut[-1] * am[-3 1; 3] * am[1 -2; 4] * - conj(mm_f[-4; 3 4]) + conj(mm_f[-4; 3 4]) @plansor ut_amap[-1 -2; -3 -4] := ut[-1] * am[-3 1; 3] * ap[1 -2; 4] * - conj(mp_f[-4; 3 4]) + conj(mp_f[-4; 3 4]) @plansor ut_apam[-1 -2; -3 -4] := ut[-1] * ap[-3 1; 3] * am[1 -2; 4] * - conj(pm_f[-4; 3 4]) + conj(pm_f[-4; 3 4]) @plansor bpbp_ut[-1 -2; -3 -4] := mm_f[-1; 1 2] * bp[1; -3 3] * bp[2; 3 -2] * - conj(ut[-4]) + conj(ut[-4]) @plansor bmbm_ut[-1 -2; -3 -4] := pp_f[-1; 1 2] * bm[1; -3 3] * bm[2; 3 -2] * - conj(ut[-4]) + conj(ut[-4]) @plansor bmbp_ut[-1 -2; -3 -4] := pm_f[-1; 1 2] * bm[1; -3 3] * bp[2; 3 -2] * - conj(ut[-4]) + conj(ut[-4]) @plansor bpbm_ut[-1 -2; -3 -4] := mp_f[-1; 1 2] * bp[1; -3 3] * bm[2; 3 -2] * - conj(ut[-4]) + conj(ut[-4]) for i in 1:basis_size if i < half_basis_size @@ -214,23 +214,15 @@ function mapped_quantum_chemistry_hamiltonian(E0, K, V, Elt=ComplexF64) iso_pp = isomorphism(_lastspace(ap)', _lastspace(ap)') iso_mm = isomorphism(_lastspace(am)', _lastspace(am)') - @plansor p_ap[-1 -2; -3 -4] := iso_pp[-1; 1] * τ[1 2; -3 3] * ap[2 -2; 4] * - conj(pp_f[-4; 3 4]) - @plansor m_ap[-1 -2; -3 -4] := iso_mm[-1; 1] * τ[1 2; -3 3] * ap[2 -2; 4] * - conj(mp_f[-4; 3 4]) - @plansor p_am[-1 -2; -3 -4] := iso_pp[-1; 1] * τ[1 2; -3 3] * am[2 -2; 4] * - conj(pm_f[-4; 3 4]) - @plansor m_am[-1 -2; -3 -4] := iso_mm[-1; 1] * τ[1 2; -3 3] * am[2 -2; 4] * - conj(mm_f[-4; 3 4]) - - @plansor bp_p[-1 -2; -3 -4] := bp[2; -3 3] * iso_mm[1; -4] * τ[4 -2; 3 1] * - mm_f[-1; 2 4] - @plansor bm_p[-1 -2; -3 -4] := bm[2; -3 3] * iso_mm[1; -4] * τ[4 -2; 3 1] * - pm_f[-1; 2 4] - @plansor bm_m[-1 -2; -3 -4] := bm[2; -3 3] * iso_pp[1; -4] * τ[4 -2; 3 1] * - pp_f[-1; 2 4] - @plansor bp_m[-1 -2; -3 -4] := bp[2; -3 3] * iso_pp[1; -4] * τ[4 -2; 3 1] * - mp_f[-1; 2 4] + @plansor p_ap[-1 -2; -3 -4] := iso_pp[-1; 1] * τ[1 2; -3 3] * ap[2 -2; 4] * conj(pp_f[-4; 3 4]) + @plansor m_ap[-1 -2; -3 -4] := iso_mm[-1; 1] * τ[1 2; -3 3] * ap[2 -2; 4] * conj(mp_f[-4; 3 4]) + @plansor p_am[-1 -2; -3 -4] := iso_pp[-1; 1] * τ[1 2; -3 3] * am[2 -2; 4] * conj(pm_f[-4; 3 4]) + @plansor m_am[-1 -2; -3 -4] := iso_mm[-1; 1] * τ[1 2; -3 3] * am[2 -2; 4] * conj(mm_f[-4; 3 4]) + + @plansor bp_p[-1 -2; -3 -4] := bp[2; -3 3] * iso_mm[1; -4] * τ[4 -2; 3 1] * mm_f[-1; 2 4] + @plansor bm_p[-1 -2; -3 -4] := bm[2; -3 3] * iso_mm[1; -4] * τ[4 -2; 3 1] * pm_f[-1; 2 4] + @plansor bm_m[-1 -2; -3 -4] := bm[2; -3 3] * iso_pp[1; -4] * τ[4 -2; 3 1] * pp_f[-1; 2 4] + @plansor bp_m[-1 -2; -3 -4] := bp[2; -3 3] * iso_pp[1; -4] * τ[4 -2; 3 1] * mp_f[-1; 2 4] for i in 1:basis_size, j in (i + 1):basis_size if j < half_basis_size @@ -326,18 +318,15 @@ function mapped_quantum_chemistry_hamiltonian(E0, K, V, Elt=ComplexF64) # 2|2 for i in 1:basis_size, j in (i + 1):basis_size # p p | . . m m - @plansor __mm[-1 -2; -3 -4] := pp_f[-1; 1 2] * bm[1; -3 3] * bm[2; 3 -2] * - conj(ut[-4]) + @plansor __mm[-1 -2; -3 -4] := pp_f[-1; 1 2] * bm[1; -3 3] * bm[2; 3 -2] * conj(ut[-4]) hamdat[j, map_2[1, i, 1, i], end] += V[i, i, j, j] * __mm # m m | p p . . - @plansor __pp[-1 -2; -3 -4] := mm_f[-1; 1 2] * bp[1; -3 3] * bp[2; 3 -2] * - conj(ut[-4]) + @plansor __pp[-1 -2; -3 -4] := mm_f[-1; 1 2] * bp[1; -3 3] * bp[2; 3 -2] * conj(ut[-4]) hamdat[j, map_2[2, i, 2, i], end] += V[j, j, i, i] * __pp # p m | . p . m - @plansor _p_m[-1 -2; -3 -4] := mp_f[-1; 1 2] * bp[1; -3 3] * bm[2; 3 -2] * - conj(ut[-4]) + @plansor _p_m[-1 -2; -3 -4] := mp_f[-1; 1 2] * bp[1; -3 3] * bm[2; 3 -2] * conj(ut[-4]) hamdat[j, map_2[2, i, 1, i], end] += V[i, j, i, j] * _p_m hamdat[i, 1, end] -= V[i, j, i, j] * h_pm @@ -464,26 +453,23 @@ function mapped_quantum_chemistry_hamiltonian(E0, K, V, Elt=ComplexF64) for i in 1:basis_size, j in (i + 1):basis_size, k in (j + 1):basis_size # j i m m @plansor jimm[-1 -2; -3 -4] := bm[1; -3 2] * bm[3; 2 -2] * τ[4 5; 1 3] * - pp_f[-1; 4 5] * conj(ut[-4]) + pp_f[-1; 4 5] * conj(ut[-4]) hamdat[k, map_2[1, i, 1, j], end] += V[j, i, k, k] * jimm # i j m m @plansor ijmm[-1 -2; -3 -4] := bm[1; -3 2] * - bm[3; 2 -2] * - τ[4 5; 1 3] * - permute(pp_f, (1,), (3, 2))[-1; 4 5] * - conj(ut[-4]) + bm[3; 2 -2] * + τ[4 5; 1 3] * + permute(pp_f, (1,), (3, 2))[-1; 4 5] * + conj(ut[-4]) hamdat[k, map_2[1, i, 1, j], end] += V[i, j, k, k] * ijmm # j p i m @plansor jpim[-1 -2; -3 -4] := bm[1; -3 2] * bp[3; 2 -2] * τ[4 5; 1 3] * - mp_f[-1; 4 5] * conj(ut[-4]) + mp_f[-1; 4 5] * conj(ut[-4]) hamdat[k, map_2[2, i, 1, j], end] += V[j, k, i, k] * jpim # i p j m - @plansor ipjm[-1 -2; -3 -4] := bm[1; -3 2] * - bp[3; 2 -2] * - τ[4 5; 1 3] * - permute(pm_f, (1,), (3, 2))[-1; 4 5] * - conj(ut[-4]) + @plansor ipjm[-1 -2; -3 -4] := bm[1; -3 2] * bp[3; 2 -2] * + τ[4 5; 1 3] * permute(pm_f, (1,), (3, 2))[-1; 4 5] * conj(ut[-4]) hamdat[k, map_2[1, i, 2, j], end] += V[i, k, j, k] * ipjm # j p m i @@ -494,12 +480,11 @@ function mapped_quantum_chemistry_hamiltonian(E0, K, V, Elt=ComplexF64) hamdat[k, map_3[1, i, 2, j], end] += V[i, k, k, j] * ipmj # p p j i - @plansor ppji[-1 -2; -3 -4] := mm_f[-1; 1 2] * bp[1; -3 3] * bp[2; 3 -2] * - conj(ut[-4]) + @plansor ppji[-1 -2; -3 -4] := mm_f[-1; 1 2] * bp[1; -3 3] * bp[2; 3 -2] * conj(ut[-4]) hamdat[k, map_2[2, i, 2, j], end] += V[k, k, j, i] * ppji # p p i j @plansor ppij[-1 -2; -3 -4] := permute(mm_f, (1,), (3, 2))[-1; 1 2] * bp[1; -3 3] * - bp[2; 3 -2] * conj(ut[-4]) + bp[2; 3 -2] * conj(ut[-4]) hamdat[k, map_2[2, i, 2, j], end] += V[k, k, i, j] * ppij # p . . m @@ -508,7 +493,7 @@ function mapped_quantum_chemistry_hamiltonian(E0, K, V, Elt=ComplexF64) pijm = ipmj hamdat[k, map_3[1, i, 2, j], end] += V[k, i, j, k] * pijm - # p . m . + # p . m . pjmi = jpim hamdat[k, map_2[2, i, 1, j], end] += V[k, j, k, i] * pjmi pimj = ipjm @@ -518,10 +503,8 @@ function mapped_quantum_chemistry_hamiltonian(E0, K, V, Elt=ComplexF64) # 1|1|1|1 # (i,j) in map_2, 1 in map_4, 1 onsite - for i in 1:basis_size, - j in (i + 1):basis_size, - k in (j + 1):basis_size, - l in (k + 1):basis_size + for i in 1:basis_size, j in (i + 1):basis_size, + k in (j + 1):basis_size, l in (k + 1):basis_size # p j i R @plansor pjiR[-1 -2; -3 -4] := ut[-1] * ap[-3 -2; -4] hamdat[k, map_3[2, i, 1, j], map_4[2, l]] += V[k, j, i, l] * pjiR @@ -560,7 +543,7 @@ function mapped_quantum_chemistry_hamiltonian(E0, K, V, Elt=ComplexF64) # i j m R @plansor ijmR[-1 -2; -3 -4] := permute(pp_f, (1,), (3, 2))[-1; 1 2] * - τ[3 2; -4 -2] * bm[1; -3 3] + τ[3 2; -4 -2] * bm[1; -3 3] hamdat[k, map_2[1, i, 1, j], map_4[2, l]] += V[i, j, k, l] * ijmR # j i R m @@ -569,7 +552,7 @@ function mapped_quantum_chemistry_hamiltonian(E0, K, V, Elt=ComplexF64) # i j R m @plansor ijRm[-1 -2; -3 -4] := permute(pp_f, (1,), (3, 2))[-1; 1 2] * - τ[1 3; -3 -4] * bm[2; 3 -2] + τ[1 3; -3 -4] * bm[2; 3 -2] hamdat[k, map_2[1, i, 1, j], map_4[2, l]] += V[i, j, l, k] * ijRm # j p i R @@ -578,7 +561,7 @@ function mapped_quantum_chemistry_hamiltonian(E0, K, V, Elt=ComplexF64) # i p j R @plansor ipjR[-1 -2; -3 -4] := permute(pm_f, (1,), (3, 2))[-1; 1 2] * - τ[3 2; -4 -2] * bp[1; -3 3] + τ[3 2; -4 -2] * bp[1; -3 3] hamdat[k, map_2[1, i, 2, j], map_4[2, l]] += V[i, k, j, l] * ipjR # j R i m @@ -587,7 +570,7 @@ function mapped_quantum_chemistry_hamiltonian(E0, K, V, Elt=ComplexF64) # i R j m @plansor iRjm[-1 -2; -3 -4] := permute(pm_f, (1,), (3, 2))[-1; 1 2] * - τ[1 3; -3 -4] * bm[2; 3 -2] + τ[1 3; -3 -4] * bm[2; 3 -2] hamdat[k, map_2[1, i, 2, j], map_4[1, l]] += V[i, l, j, k] * iRjm # p j R i @@ -596,7 +579,7 @@ function mapped_quantum_chemistry_hamiltonian(E0, K, V, Elt=ComplexF64) # p i R j @plansor piRj[-1 -2; -3 -4] := permute(pm_f, (1,), (3, 2))[-1; 1 2] * - τ[3 2; -4 -2] * bp[1; -3 3] + τ[3 2; -4 -2] * bp[1; -3 3] hamdat[k, map_2[1, i, 2, j], map_4[2, l]] += V[k, i, l, j] * piRj # R j m i @@ -605,7 +588,7 @@ function mapped_quantum_chemistry_hamiltonian(E0, K, V, Elt=ComplexF64) # R i m j @plansor Rimj[-1 -2; -3 -4] := permute(pm_f, (1,), (3, 2))[-1; 1 2] * - τ[1 3; -3 -4] * bm[2; 3 -2] + τ[1 3; -3 -4] * bm[2; 3 -2] hamdat[k, map_2[1, i, 2, j], map_4[1, l]] += V[l, i, k, j] * Rimj # R p j i @@ -614,7 +597,7 @@ function mapped_quantum_chemistry_hamiltonian(E0, K, V, Elt=ComplexF64) # R p i j @plansor Rpij[-1 -2; -3 -4] := permute(mm_f, (1,), (3, 2))[-1; 1 2] * - τ[1 3; -3 -4] * bp[2; 3 -2] + τ[1 3; -3 -4] * bp[2; 3 -2] hamdat[k, map_2[2, i, 2, j], map_4[1, l]] += V[l, k, i, j] * Rpij # p R j i @@ -623,7 +606,7 @@ function mapped_quantum_chemistry_hamiltonian(E0, K, V, Elt=ComplexF64) # p R i j @plansor pRij[-1 -2; -3 -4] := permute(mm_f, (1,), (3, 2))[-1; 1 2] * - τ[3 2; -4 -2] * bp[1; -3 3] + τ[3 2; -4 -2] * bp[1; -3 3] hamdat[k, map_2[2, i, 2, j], map_4[1, l]] += V[k, l, i, j] * pRij end diff --git a/src/models/transfermatrices.jl b/src/models/transfermatrices.jl index c9ff600..0f7fc54 100644 --- a/src/models/transfermatrices.jl +++ b/src/models/transfermatrices.jl @@ -17,8 +17,10 @@ function classical_ising end function classical_ising(symmetry::Type{<:Sector}; kwargs...) return classical_ising(ComplexF64, symmetry; kwargs...) end -function classical_ising(elt::Type{<:Number}=ComplexF64, ::Type{Trivial}=Trivial; - beta=log(1 + sqrt(2)) / 2) +function classical_ising( + elt::Type{<:Number} = ComplexF64, ::Type{Trivial} = Trivial; + beta = log(1 + sqrt(2)) / 2 + ) t = elt[exp(beta) exp(-beta); exp(-beta) exp(beta)] r = eigen(t) @@ -33,14 +35,14 @@ function classical_ising(elt::Type{<:Number}=ComplexF64, ::Type{Trivial}=Trivial return InfiniteMPO([TensorMap(o, ℂ^2 * ℂ^2, ℂ^2 * ℂ^2)]) end -function classical_ising(elt::Type{<:Number}, ::Type{Z2Irrep}; beta=log(1 + sqrt(2)) / 2) +function classical_ising(elt::Type{<:Number}, ::Type{Z2Irrep}; beta = log(1 + sqrt(2)) / 2) x = cosh(beta) y = sinh(beta) sec = ℤ₂Space(0 => 1, 1 => 1) mpo = zeros(elt, sec * sec, sec * sec) - block(mpo, Irrep[ℤ₂](0)) .= [2x^2 2x*y; 2x*y 2y^2] - block(mpo, Irrep[ℤ₂](1)) .= [2x*y 2x*y; 2x*y 2x*y] + block(mpo, Irrep[ℤ₂](0)) .= [2x^2 2x * y; 2x * y 2y^2] + block(mpo, Irrep[ℤ₂](1)) .= [2x * y 2x * y; 2x * y 2x * y] return InfiniteMPO([mpo]) end @@ -57,15 +59,19 @@ MPO for the partition function of the two-dimensional six vertex model. """ function sixvertex end sixvertex(symmetry::Type{<:Sector}; kwargs...) = sixvertex(ComplexF64, symmetry; kwargs...) -function sixvertex(elt::Type{<:Number}=ComplexF64, ::Type{Trivial}=Trivial; a=1.0, b=1.0, - c=1.0) - d = elt[a 0 0 0 - 0 c b 0 - 0 b c 0 - 0 0 0 a] +function sixvertex( + elt::Type{<:Number} = ComplexF64, ::Type{Trivial} = Trivial; + a = 1.0, b = 1.0, c = 1.0 + ) + d = elt[ + a 0 0 0 + 0 c b 0 + 0 b c 0 + 0 0 0 a + ] return InfiniteMPO([permute(TensorMap(d, ℂ^2 ⊗ ℂ^2, ℂ^2 ⊗ ℂ^2), ((1, 2), (4, 3)))]) end -function sixvertex(elt::Type{<:Number}, ::Type{U1Irrep}; a=1.0, b=1.0, c=1.0) +function sixvertex(elt::Type{<:Number}, ::Type{U1Irrep}; a = 1.0, b = 1.0, c = 1.0) pspace = U1Space(-1 // 2 => 1, 1 // 2 => 1) mpo = zeros(elt, pspace ⊗ pspace, pspace ⊗ pspace) block(mpo, Irrep[U₁](0)) .= [b c; c b] @@ -73,7 +79,7 @@ function sixvertex(elt::Type{<:Number}, ::Type{U1Irrep}; a=1.0, b=1.0, c=1.0) block(mpo, Irrep[U₁](-1)) .= reshape([a], (1, 1)) return InfiniteMPO([permute(mpo, ((1, 2), (4, 3)))]) end -function sixvertex(elt::Type{<:Number}, ::Type{CU1Irrep}; a=1.0, b=1.0, c=1.0) +function sixvertex(elt::Type{<:Number}, ::Type{CU1Irrep}; a = 1.0, b = 1.0, c = 1.0) pspace = CU1Space(1 // 2 => 1) mpo = zeros(elt, pspace ⊗ pspace, pspace ⊗ pspace) block(mpo, Irrep[CU₁](0, 0)) .= reshape([b + c], (1, 1)) @@ -91,7 +97,7 @@ end MPO for the partition function of the two-dimensional hard hexagon model. """ -function hard_hexagon(elt::Type{<:Number}=ComplexF64) +function hard_hexagon(elt::Type{<:Number} = ComplexF64) P = Vect[FibonacciAnyon](:τ => 1) O = ones(elt, P ⊗ P ← P ⊗ P) block(O, FibonacciAnyon(:I)) .*= 0 @@ -107,13 +113,14 @@ end MPO for the partition function of the two-dimensional discrete clock model with ``q`` states. """ -function qstate_clock(elt::Type{<:Number}=ComplexF64, ::Type{Trivial}=Trivial; - beta::Number=1.0, q::Integer=3) +function qstate_clock( + elt::Type{<:Number} = ComplexF64, ::Type{Trivial} = Trivial; + beta::Number = 1.0, q::Integer = 3 + ) comega(d) = cos(2 * pi * d / q) O = zeros(elt, q, q, q, q) for i in 1:q, j in 1:q, k in 1:q, l in 1:q - O[i, j, k, l] = exp(beta * - (comega(i - j) + comega(j - k) + comega(k - l) + comega(l - i))) + O[i, j, k, l] = exp(beta * (comega(i - j) + comega(j - k) + comega(k - l) + comega(l - i))) end return InfiniteMPO([TensorMap(O, ℂ^q * ℂ^q, ℂ^q * ℂ^q)]) diff --git a/src/operators/bosonoperators.jl b/src/operators/bosonoperators.jl index 5ffbf4b..0f6c4ac 100644 --- a/src/operators/bosonoperators.jl +++ b/src/operators/bosonoperators.jl @@ -9,7 +9,7 @@ a_plus(; kwargs...) = a_plus(ComplexF64, Trivial; kwargs...) a_plus(elt::Type{<:Number}; kwargs...) = a_plus(elt, Trivial; kwargs...) a_plus(symm::Type{<:Sector}; kwargs...) = a_plus(ComplexF64, symm; kwargs...) -function a_plus(elt::Type{<:Number}, ::Type{Trivial}; cutoff::Integer=5) +function a_plus(elt::Type{<:Number}, ::Type{Trivial}; cutoff::Integer = 5) a⁺ = zeros(elt, ComplexSpace(cutoff + 1), ComplexSpace(cutoff + 1)) for n in 1:cutoff a⁺[n + 1, n] = sqrt(n) @@ -17,7 +17,7 @@ function a_plus(elt::Type{<:Number}, ::Type{Trivial}; cutoff::Integer=5) return a⁺ end -function a_plus(elt::Type{<:Number}, ::Type{U1Irrep}; cutoff::Integer=5, side=:L) +function a_plus(elt::Type{<:Number}, ::Type{U1Irrep}; cutoff::Integer = 5, side = :L) pspace = U1Space(n => 1 for n in 0:cutoff) if side === :L vspace = U1Space(1 => 1) @@ -56,7 +56,7 @@ a_min(; kwargs...) = a_min(ComplexF64, Trivial; kwargs...) a_min(elt::Type{<:Number}; kwargs...) = a_min(elt, Trivial; kwargs...) a_min(symm::Type{<:Sector}; kwargs...) = a_min(ComplexF64, symm; kwargs...) -function a_min(elt::Type{<:Number}, ::Type{Trivial}; cutoff::Integer=5) +function a_min(elt::Type{<:Number}, ::Type{Trivial}; cutoff::Integer = 5) a⁻ = zeros(elt, ComplexSpace(cutoff + 1), ComplexSpace(cutoff + 1)) for n in 1:cutoff a⁻[n, n + 1] = sqrt(n) @@ -64,7 +64,7 @@ function a_min(elt::Type{<:Number}, ::Type{Trivial}; cutoff::Integer=5) return a⁻ end -function a_min(elt::Type{<:Number}, ::Type{U1Irrep}; cutoff::Integer=5, side=:L) +function a_min(elt::Type{<:Number}, ::Type{U1Irrep}; cutoff::Integer = 5, side = :L) pspace = U1Space(n => 1 for n in 0:cutoff) if side === :L vspace = U1Space(-1 => 1) @@ -93,23 +93,31 @@ end const a⁻ = a_min a_plusmin(symmetry::Type{<:Sector}; kwargs...) = a_plusmin(ComplexF64, symmetry; kwargs...) -function a_plusmin(elt::Type{<:Number}=ComplexF64, ::Type{Trivial}=Trivial; - cutoff::Integer=5) - return contract_twosite(a⁺(elt; cutoff=cutoff), a⁻(elt; cutoff=cutoff)) +function a_plusmin( + elt::Type{<:Number} = ComplexF64, ::Type{Trivial} = Trivial; + cutoff::Integer = 5 + ) + return contract_twosite(a⁺(elt; cutoff = cutoff), a⁻(elt; cutoff = cutoff)) end -function a_plusmin(elt::Type{<:Number}, ::Type{U1Irrep}; cutoff::Integer=5) - return contract_twosite(a⁺(elt, U1Irrep; cutoff=cutoff, side=:L), - a⁻(elt, U1Irrep; cutoff=cutoff, side=:R)) +function a_plusmin(elt::Type{<:Number}, ::Type{U1Irrep}; cutoff::Integer = 5) + return contract_twosite( + a⁺(elt, U1Irrep; cutoff = cutoff, side = :L), + a⁻(elt, U1Irrep; cutoff = cutoff, side = :R) + ) end a_minplus(symmetry::Type{<:Sector}; kwargs...) = a_minplus(ComplexF64, symmetry; kwargs...) -function a_minplus(elt::Type{<:Number}=ComplexF64, ::Type{Trivial}=Trivial; - cutoff::Integer=5) - return contract_twosite(a⁻(elt; cutoff=cutoff), a⁺(elt; cutoff=cutoff)) +function a_minplus( + elt::Type{<:Number} = ComplexF64, ::Type{Trivial} = Trivial; + cutoff::Integer = 5 + ) + return contract_twosite(a⁻(elt; cutoff = cutoff), a⁺(elt; cutoff = cutoff)) end -function a_minplus(elt::Type{<:Number}, ::Type{U1Irrep}; cutoff::Integer=5) - return contract_twosite(a⁻(elt, U1Irrep; cutoff=cutoff, side=:L), - a⁺(elt, U1Irrep; cutoff=cutoff, side=:R)) +function a_minplus(elt::Type{<:Number}, ::Type{U1Irrep}; cutoff::Integer = 5) + return contract_twosite( + a⁻(elt, U1Irrep; cutoff = cutoff, side = :L), + a⁺(elt, U1Irrep; cutoff = cutoff, side = :R) + ) end """ @@ -122,7 +130,7 @@ a_number(; kwargs...) = a_number(ComplexF64, Trivial; kwargs...) a_number(elt::Type{<:Number}; kwargs...) = a_number(elt, Trivial; kwargs...) a_number(symm::Type{<:Sector}; kwargs...) = a_number(ComplexF64, symm; kwargs...) -function a_number(elt::Type{<:Number}, ::Type{Trivial}; cutoff::Integer=5) +function a_number(elt::Type{<:Number}, ::Type{Trivial}; cutoff::Integer = 5) N = zeros(elt, ComplexSpace(cutoff + 1), ComplexSpace(cutoff + 1)) for n in 0:cutoff N[n + 1, n + 1] = n @@ -130,7 +138,7 @@ function a_number(elt::Type{<:Number}, ::Type{Trivial}; cutoff::Integer=5) return N end -function a_number(elt::Type{<:Number}, ::Type{U1Irrep}; cutoff::Integer=5) +function a_number(elt::Type{<:Number}, ::Type{U1Irrep}; cutoff::Integer = 5) pspace = U1Space(n => 1 for n in 0:cutoff) N = zeros(elt, pspace, pspace) for (c, b) in blocks(N) diff --git a/src/operators/fermionoperators.jl b/src/operators/fermionoperators.jl index d3ef213..8f148f9 100644 --- a/src/operators/fermionoperators.jl +++ b/src/operators/fermionoperators.jl @@ -8,14 +8,14 @@ Fermionic creation operator. """ -function c_plus(elt::Type{<:Number}=ComplexF64; side=:L) +function c_plus(elt::Type{<:Number} = ComplexF64; side = :L) vspace = Vect[fℤ₂](1 => 1) if side === :L pspace = Vect[fℤ₂](0 => 1, 1 => 1) c⁺ = zeros(elt, pspace ← pspace ⊗ vspace) block(c⁺, fℤ₂(1)) .= one(elt) elseif side === :R - C = c_plus(elt; side=:L) + C = c_plus(elt; side = :L) F = isomorphism(storagetype(C), vspace, flip(vspace)) @planar c⁺[-1 -2; -3] := C[-2; 1 2] * τ[1 2; 3 -3] * F[3; -1] else @@ -31,13 +31,13 @@ const c⁺ = c_plus Fermionic annihilation operator. """ -function c_min(elt::Type{<:Number}=ComplexF64; side=:L) +function c_min(elt::Type{<:Number} = ComplexF64; side = :L) if side === :L - C = c_plus(elt; side=:L)' + C = c_plus(elt; side = :L)' F = isomorphism(flip(space(C, 2)), space(C, 2)) @planar c⁻[-1; -2 -3] := C[-1 1; -2] * F[-3; 1] elseif side === :R - c⁻ = permute(c_plus(elt; side=:L)', ((2, 1), (3,))) + c⁻ = permute(c_plus(elt; side = :L)', ((2, 1), (3,))) else throw(ArgumentError("invalid side `:$side`, expected `:L` or `:R`")) end @@ -46,13 +46,13 @@ end const c⁻ = c_min -c_plusmin(elt=ComplexF64) = contract_twosite(c⁺(elt; side=:L), c⁻(elt; side=:R)) +c_plusmin(elt = ComplexF64) = contract_twosite(c⁺(elt; side = :L), c⁻(elt; side = :R)) const c⁺c⁻ = c_plusmin -c_minplus(elt=ComplexF64) = contract_twosite(c⁻(elt; side=:L), c⁺(elt; side=:R)) +c_minplus(elt = ComplexF64) = contract_twosite(c⁻(elt; side = :L), c⁺(elt; side = :R)) const c⁻c⁺ = c_minplus -c_plusplus(elt=ComplexF64) = contract_twosite(c⁺(elt; side=:L), c⁺(elt; side=:R)) +c_plusplus(elt = ComplexF64) = contract_twosite(c⁺(elt; side = :L), c⁺(elt; side = :R)) const c⁺c⁺ = c_plusplus -c_minmin(elt=ComplexF64) = contract_twosite(c⁻(elt; side=:L), c⁻(elt; side=:R)) +c_minmin(elt = ComplexF64) = contract_twosite(c⁻(elt; side = :L), c⁻(elt; side = :R)) const c⁻c⁻ = c_minmin """ @@ -60,7 +60,7 @@ const c⁻c⁻ = c_minmin Fermionic number operator. """ -function c_number(elt::Type{<:Number}=ComplexF64) +function c_number(elt::Type{<:Number} = ComplexF64) pspace = Vect[fℤ₂](0 => 1, 1 => 1) n = zeros(elt, pspace ← pspace) block(n, fℤ₂(1)) .= one(elt) diff --git a/src/operators/hubbardoperators.jl b/src/operators/hubbardoperators.jl index 2bf4204..1255850 100644 --- a/src/operators/hubbardoperators.jl +++ b/src/operators/hubbardoperators.jl @@ -19,7 +19,7 @@ export nꜛ, nꜜ, nꜛnꜜ Return the local hilbert space for a Hubbard-type model with the given particle and spin symmetries. The possible symmetries are `Trivial`, `U1Irrep`, and `SU2Irrep`, for both particle number and spin. """ -function hubbard_space(::Type{Trivial}=Trivial, ::Type{Trivial}=Trivial) +function hubbard_space(::Type{Trivial} = Trivial, ::Type{Trivial} = Trivial) return Vect[FermionParity](0 => 2, 1 => 2) end function hubbard_space(::Type{Trivial}, ::Type{U1Irrep}) @@ -32,12 +32,15 @@ function hubbard_space(::Type{U1Irrep}, ::Type{Trivial}) return Vect[FermionParity ⊠ U1Irrep]((0, 0) => 1, (1, 1) => 2, (0, 2) => 1) end function hubbard_space(::Type{U1Irrep}, ::Type{U1Irrep}) - return Vect[FermionParity ⊠ U1Irrep ⊠ U1Irrep]((0, 0, 0) => 1, (1, 1, 1 // 2) => 1, - (1, 1, -1 // 2) => 1, (0, 2, 0) => 1) + return Vect[FermionParity ⊠ U1Irrep ⊠ U1Irrep]( + (0, 0, 0) => 1, (1, 1, 1 // 2) => 1, + (1, 1, -1 // 2) => 1, (0, 2, 0) => 1 + ) end function hubbard_space(::Type{U1Irrep}, ::Type{SU2Irrep}) - return Vect[FermionParity ⊠ U1Irrep ⊠ SU2Irrep]((0, 0, 0) => 1, (1, 1, 1 // 2) => 1, - (0, 2, 0) => 1) + return Vect[FermionParity ⊠ U1Irrep ⊠ SU2Irrep]( + (0, 0, 0) => 1, (1, 1, 1 // 2) => 1, (0, 2, 0) => 1 + ) end function hubbard_space(::Type{SU2Irrep}, ::Type{Trivial}) return Vect[FermionParity ⊠ SU2Irrep]((0, 0) => 2, (1, 1 // 2) => 1) @@ -49,14 +52,16 @@ function hubbard_space(::Type{SU2Irrep}, ::Type{SU2Irrep}) return Vect[FermionParity ⊠ SU2Irrep ⊠ SU2Irrep]((1, 1 // 2, 1 // 2) => 1) end -function single_site_operator(T, particle_symmetry::Type{<:Sector}, - spin_symmetry::Type{<:Sector}) +function single_site_operator( + T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector} + ) V = hubbard_space(particle_symmetry, spin_symmetry) return zeros(T, V ← V) end -function two_site_operator(T, particle_symmetry::Type{<:Sector}, - spin_symmetry::Type{<:Sector}) +function two_site_operator( + T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector} + ) V = hubbard_space(particle_symmetry, spin_symmetry) return zeros(T, V ⊗ V ← V ⊗ V) end @@ -214,7 +219,7 @@ This is the sum of `e_plusmin_up` and `e_plusmin_down`. e_plusmin(P::Type{<:Sector}, S::Type{<:Sector}) = e_plusmin(ComplexF64, P, S) function e_plusmin(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}) return e_plusmin_up(T, particle_symmetry, spin_symmetry) + - e_plusmin_down(T, particle_symmetry, spin_symmetry) + e_plusmin_down(T, particle_symmetry, spin_symmetry) end function e_plusmin(T, ::Type{Trivial}, ::Type{SU2Irrep}) t = two_site_operator(T, Trivial, SU2Irrep) @@ -270,7 +275,7 @@ const e⁻e⁺ = e_minplus Return the one-body operator that counts the number of spin-up electrons. """ e_number_up(P::Type{<:Sector}, S::Type{<:Sector}) = e_number_up(ComplexF64, P, S) -function e_number_up(T::Type{<:Number}, ::Type{Trivial}=Trivial, ::Type{Trivial}=Trivial) +function e_number_up(T::Type{<:Number}, ::Type{Trivial} = Trivial, ::Type{Trivial} = Trivial) t = single_site_operator(T, Trivial, Trivial) I = sectortype(t) t[(I(1), I(1))][1, 1] = 1 @@ -321,7 +326,7 @@ const nꜛ = e_number_up Return the one-body operator that counts the number of spin-down electrons. """ e_number_down(P::Type{<:Sector}, S::Type{<:Sector}) = e_number_down(ComplexF64, P, S) -function e_number_down(T::Type{<:Number}, ::Type{Trivial}=Trivial, ::Type{Trivial}=Trivial) +function e_number_down(T::Type{<:Number}, ::Type{Trivial} = Trivial, ::Type{Trivial} = Trivial) t = single_site_operator(T, Trivial, Trivial) I = sectortype(t) t[(I(1), I(1))][2, 2] = 1 @@ -374,7 +379,7 @@ Return the one-body operator that counts the number of particles. e_number(P::Type{<:Sector}, S::Type{<:Sector}) = e_number(ComplexF64, P, S) function e_number(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}) return e_number_up(T, particle_symmetry, spin_symmetry) + - e_number_down(T, particle_symmetry, spin_symmetry) + e_number_down(T, particle_symmetry, spin_symmetry) end function e_number(T, ::Type{Trivial}, ::Type{SU2Irrep}) t = single_site_operator(T, Trivial, SU2Irrep) @@ -398,10 +403,12 @@ const n = e_number Return the one-body operator that counts the number of doubly occupied sites. """ e_number_updown(P::Type{<:Sector}, S::Type{<:Sector}) = e_number_updown(ComplexF64, P, S) -function e_number_updown(T, particle_symmetry::Type{<:Sector}, - spin_symmetry::Type{<:Sector}) +function e_number_updown( + T, particle_symmetry::Type{<:Sector}, + spin_symmetry::Type{<:Sector} + ) return e_number_up(T, particle_symmetry, spin_symmetry) * - e_number_down(T, particle_symmetry, spin_symmetry) + e_number_down(T, particle_symmetry, spin_symmetry) end function e_number_updown(T, ::Type{Trivial}, ::Type{SU2Irrep}) t = single_site_operator(T, Trivial, SU2Irrep) diff --git a/src/operators/localoperators.jl b/src/operators/localoperators.jl index 88f8bb1..a52e701 100644 --- a/src/operators/localoperators.jl +++ b/src/operators/localoperators.jl @@ -8,30 +8,30 @@ operator is represented as a vector of `MPOTens, instantiate_operatoror`s, each - `opp::Vector{T}`: `N`-body operator represented by an MPO. - `inds::Vector{G}`: `N` site indices. """ -struct LocalOperator{T<:AbstractTensorMap{<:Number,<:Any,2,2},G<:LatticePoint} +struct LocalOperator{T <: AbstractTensorMap{<:Number, <:Any, 2, 2}, G <: LatticePoint} opp::Vector{T} inds::Vector{G} - function LocalOperator{T,G}(O::Vector{T}, - inds::Vector{G}) where {T<:AbstractTensorMap{<:Number,<:Any, - 2,2}, - G<:LatticePoint} + function LocalOperator{T, G}( + O::Vector{T}, inds::Vector{G} + ) where {T <: AbstractTensorMap{<:Number, <:Any, 2, 2}, G <: LatticePoint} length(O) == length(inds) || throw(ArgumentError("number of operators and indices should be the same")) issorted(inds) && allunique(inds) || throw(ArgumentError("indices should be ascending and unique")) allequal(getfield.(inds, :lattice)) || throw(ArgumentError("points should be defined on the same lattice")) - return new{T,G}(O, inds) + return new{T, G}(O, inds) end end -function LocalOperator(t::AbstractTensorMap{<:Number,<:Any,N,N}, - inds::Vararg{G,N}) where {N,G<:LatticePoint} +function LocalOperator( + t::AbstractTensorMap{<:Number, <:Any, N, N}, inds::Vararg{G, N} + ) where {N, G <: LatticePoint} p = TupleTools.sortperm(linearize_index.(inds)) t = permute(t, (p, p .+ N)) t_mpo = collect(MPSKit.decompose_localmpo(MPSKit.add_util_leg(t))) - return LocalOperator{eltype(t_mpo),G}(t_mpo, collect(getindex.(Ref(inds), p))) + return LocalOperator{eltype(t_mpo), G}(t_mpo, collect(getindex.(Ref(inds), p))) end function MPSKit.instantiate_operator(lattice, O::LocalOperator) @@ -54,9 +54,9 @@ end # O.inds[p]) # end -Base.copy(O::LocalOperator{T,G}) where {T,G} = LocalOperator{T,G}(copy(O.opp), copy(O.inds)) -function Base.deepcopy(O::LocalOperator{T,G}) where {T,G} - return LocalOperator{T,G}(deepcopy(O.opp), deepcopy(O.inds)) +Base.copy(O::LocalOperator{T, G}) where {T, G} = LocalOperator{T, G}(copy(O.opp), copy(O.inds)) +function Base.deepcopy(O::LocalOperator{T, G}) where {T, G} + return LocalOperator{T, G}(deepcopy(O.opp), deepcopy(O.inds)) end # Linear Algebra @@ -76,7 +76,7 @@ Base.:*(a::Number, b::LocalOperator) = lmul!(a, deepcopy(b)) Base.:/(a::LocalOperator, b::Number) = a * inv(b) Base.:\(a::Number, b::LocalOperator) = inv(a) * b -function Base.:*(a::LocalOperator{T₁,G}, b::LocalOperator{T₂,G}) where {T₁,T₂,G} +function Base.:*(a::LocalOperator{T₁, G}, b::LocalOperator{T₂, G}) where {T₁, T₂, G} inds = sort!(union(a.inds, b.inds)) T = promote_type(T₁, T₂) operators = Vector{T}(undef, length(inds)) @@ -92,26 +92,24 @@ function Base.:*(a::LocalOperator{T₁,G}, b::LocalOperator{T₂,G}) where {T₁ right_vspace_A = isnothing(i_A) ? left_vspace_A : space(a.opp[i_A], 4)' right_vspace_B = isnothing(i_B) ? left_vspace_B : space(b.opp[i_B], 4)' - left_fuse = unitary(M, fuse(left_vspace_B, left_vspace_A), - left_vspace_B ⊗ left_vspace_A) - right_fuse = unitary(M, fuse(right_vspace_B, right_vspace_A), - right_vspace_B ⊗ right_vspace_A) + left_fuse = unitary( + M, fuse(left_vspace_B, left_vspace_A), + left_vspace_B ⊗ left_vspace_A + ) + right_fuse = unitary( + M, fuse(right_vspace_B, right_vspace_A), + right_vspace_B ⊗ right_vspace_A + ) if !isnothing(i_A) && !isnothing(i_B) @plansor operators[i][-1 -2; -3 -4] := b.opp[i_B][1 2; -3 4] * - a.opp[i_A][3 -2; 2 5] * - left_fuse[-1; 1 3] * - conj(right_fuse[-4; 4 5]) + a.opp[i_A][3 -2; 2 5] * left_fuse[-1; 1 3] * conj(right_fuse[-4; 4 5]) elseif !isnothing(i_A) @plansor operators[i][-1 -2; -3 -4] := τ[1 2; -3 4] * - a.opp[i_A][3 -2; 2 5] * - left_fuse[-1; 1 3] * - conj(right_fuse[-4; 4 5]) + a.opp[i_A][3 -2; 2 5] * left_fuse[-1; 1 3] * conj(right_fuse[-4; 4 5]) elseif !isnothing(i_B) @plansor operators[i][-1 -2; -3 -4] := b.opp[i_B][1 2; -3 4] * - τ[3 -2; 2 5] * - left_fuse[-1; 1 3] * - conj(right_fuse[-4; 4 5]) + τ[3 -2; 2 5] * left_fuse[-1; 1 3] * conj(right_fuse[-4; 4 5]) else error("this should not happen") end @@ -120,14 +118,14 @@ function Base.:*(a::LocalOperator{T₁,G}, b::LocalOperator{T₂,G}) where {T₁ left_vspace_B = right_vspace_B end - return LocalOperator{T,G}(operators, inds) + return LocalOperator{T, G}(operators, inds) end lattice(O::LocalOperator) = first(O.inds).lattice latticetype(O::LocalOperator) = latticetype(typeof(O)) -latticetype(::Type{<:LocalOperator{T,G}}) where {T,G} = G -tensortype(::Union{O,Type{O}}) where {T,O<:LocalOperator{T}} = T -function TensorKit.spacetype(O::Union{LocalOperator,Type{<:LocalOperator}}) +latticetype(::Type{<:LocalOperator{T, G}}) where {T, G} = G +tensortype(::Union{O, Type{O}}) where {T, O <: LocalOperator{T}} = T +function TensorKit.spacetype(O::Union{LocalOperator, Type{<:LocalOperator}}) return spacetype(tensortype(O)) end @@ -138,7 +136,7 @@ end Lazy sum of local operators. """ -struct SumOfLocalOperators{L<:LocalOperator} +struct SumOfLocalOperators{L <: LocalOperator} opps::Vector{L} function SumOfLocalOperators(opps::Vector{<:LocalOperator}) if length(opps) > 1 @@ -163,7 +161,7 @@ Base.:*(a::SumOfLocalOperators, b::Number) = SumOfLocalOperators(a.opps .* b) Base.:\(a::Number, b::SumOfLocalOperators) = SumOfLocalOperators(a .\ b.opps) Base.:/(a::SumOfLocalOperators, b::Number) = SumOfLocalOperators(a.opps ./ b) -const LocalOrSumOfLocal = Union{LocalOperator,SumOfLocalOperators} +const LocalOrSumOfLocal = Union{LocalOperator, SumOfLocalOperators} Base.:-(a::LocalOrSumOfLocal) = -1 * a Base.:-(a::LocalOrSumOfLocal, b::LocalOrSumOfLocal) = a + (-b) @@ -174,4 +172,4 @@ lattice(Os::SumOfLocalOperators) = lattice(first(Os.opps)) latticetype(Os::SumOfLocalOperators) = latticetype(typeof(Os)) latticetype(::Type{<:SumOfLocalOperators{L}}) where {L} = latticetype(L) -tensortype(::Union{Os,Type{<:Os}}) where {L,Os<:SumOfLocalOperators{L}} = tensortype(L) +tensortype(::Union{Os, Type{<:Os}}) where {L, Os <: SumOfLocalOperators{L}} = tensortype(L) diff --git a/src/operators/mpoham.jl b/src/operators/mpoham.jl index cae1b63..61b1f62 100644 --- a/src/operators/mpoham.jl +++ b/src/operators/mpoham.jl @@ -28,8 +28,7 @@ end function process_geometries_sugar(ex) @capture(ex, (((-Inf):Inf) | ((-∞):∞))) && return :(vertices(InfiniteChain())) - @capture(ex, (((-Inf):step_:(Inf | -∞)):step_:∞)) && - return :(vertices(InfiniteChain($step))) + @capture(ex, (((-Inf):step_:(Inf | -∞)):step_:∞)) && return :(vertices(InfiniteChain($step))) return ex end @@ -67,9 +66,8 @@ Attempt to automatically deduce the physical spaces for all sites of the lattice """ function deduce_pspaces(opps::SumOfLocalOperators) S = spacetype(opps) - MissingS = Union{S,Missing} - pspaces = MPSKit.PeriodicVector{MissingS}(missing, - length(lattice(opps))) + MissingS = Union{S, Missing} + pspaces = MPSKit.PeriodicVector{MissingS}(missing, length(lattice(opps))) for opp in opps.opps for (ind, O) in zip(opp.inds, opp.opp) lind = linearize_index(ind) @@ -100,9 +98,9 @@ end # define a partial order on local operators, sorting them by starting site # and then by decreasing length. -function _isless(a::L, b::L) where {L<:LocalOperator} +function _isless(a::L, b::L) where {L <: LocalOperator} return first(a.inds) == first(b.inds) ? length(a.inds) < length(b.inds) : - first(a.inds) < first(b.inds) + first(a.inds) < first(b.inds) end MPSKit.MPOHamiltonian(o::LocalOperator) = MPOHamiltonian(SumOfLocalOperators([o])) diff --git a/src/operators/spinoperators.jl b/src/operators/spinoperators.jl index 07858ef..e2e1fe5 100644 --- a/src/operators/spinoperators.jl +++ b/src/operators/spinoperators.jl @@ -6,8 +6,10 @@ end function _pauliterm(spin, i::U1Irrep, j::U1Irrep) -spin <= i.charge <= spin || return 0.0 -spin <= j.charge <= spin || return 0.0 - return sqrt((spin + 1) * (i.charge + j.charge + 2 * spin + 1) - - (i.charge + spin + 1) * (j.charge + spin + 1)) / 2.0 + return sqrt( + (spin + 1) * (i.charge + j.charge + 2 * spin + 1) - + (i.charge + spin + 1) * (j.charge + spin + 1) + ) / 2.0 end """ @@ -15,7 +17,7 @@ end the spinmatrices according to [Wikipedia](https://en.wikipedia.org/wiki/Spin_(physics)#Higher_spins). """ -function spinmatrices(s::Union{Rational{Int},Int}, elt=ComplexF64) +function spinmatrices(s::Union{Rational{Int}, Int}, elt = ComplexF64) N = Int(2 * s) Sx = zeros(elt, N + 1, N + 1) @@ -57,13 +59,13 @@ S_x(; kwargs...) = S_x(ComplexF64, Trivial; kwargs...) S_x(elt::Type{<:Number}; kwargs...) = S_x(elt, Trivial; kwargs...) S_x(symm::Type{<:Sector}; kwargs...) = S_x(ComplexF64, symm; kwargs...) -function S_x(elt::Type{<:Number}, ::Type{Trivial}; spin=1 // 2) +function S_x(elt::Type{<:Number}, ::Type{Trivial}; spin = 1 // 2) S_x_mat, _, _ = spinmatrices(spin, elt) pspace = ComplexSpace(size(S_x_mat, 1)) return TensorMap(S_x_mat, pspace ← pspace) end -function S_x(elt::Type{<:Number}, ::Type{Z2Irrep}; spin=1 // 2) +function S_x(elt::Type{<:Number}, ::Type{Z2Irrep}; spin = 1 // 2) spin == 1 // 2 || error("not implemented") pspace = Z2Space(0 => 1, 1 => 1) X = zeros(elt, pspace, pspace) @@ -72,7 +74,7 @@ function S_x(elt::Type{<:Number}, ::Type{Z2Irrep}; spin=1 // 2) return X end -function S_x(elt::Type{<:Number}, ::Type{U1Irrep}; spin=1 // 2, side=:L) +function S_x(elt::Type{<:Number}, ::Type{U1Irrep}; spin = 1 // 2, side = :L) pspace = U1Space(i => 1 for i in (-spin):spin) vspace = U1Space(1 => 1, -1 => 1) if side == :L @@ -119,13 +121,13 @@ S_y(; kwargs...) = S_y(ComplexF64, Trivial; kwargs...) S_y(elt::Type{<:Complex}; kwargs...) = S_y(elt, Trivial; kwargs...) S_y(symm::Type{<:Sector}; kwargs...) = S_y(ComplexF64, symm; kwargs...) -function S_y(elt::Type{<:Complex}, ::Type{Trivial}; spin=1 // 2) +function S_y(elt::Type{<:Complex}, ::Type{Trivial}; spin = 1 // 2) _, Y, _, _ = spinmatrices(spin, elt) pspace = ComplexSpace(size(Y, 1)) return TensorMap(Y, pspace ← pspace) end -function S_y(elt::Type{<:Complex}, ::Type{Z2Irrep}; spin=1 // 2, side=:L) +function S_y(elt::Type{<:Complex}, ::Type{Z2Irrep}; spin = 1 // 2, side = :L) spin == 1 // 2 || error("not implemented") pspace = Z2Space(0 => 1, 1 => 1) vspace = Z2Space(1 => 1) @@ -143,7 +145,7 @@ function S_y(elt::Type{<:Complex}, ::Type{Z2Irrep}; spin=1 // 2, side=:L) return Y end -function S_y(elt::Type{<:Complex}, ::Type{U1Irrep}; spin=1 // 2, side=:L) +function S_y(elt::Type{<:Complex}, ::Type{U1Irrep}; spin = 1 // 2, side = :L) pspace = U1Space(i => 1 for i in (-spin):spin) vspace = U1Space(1 => 1, -1 => 1) if side == :L @@ -195,13 +197,13 @@ S_z(; kwargs...) = S_z(ComplexF64, Trivial; kwargs...) S_z(elt::Type{<:Number}; kwargs...) = S_z(elt, Trivial; kwargs...) S_z(symm::Type{<:Sector}; kwargs...) = S_z(ComplexF64, symm; kwargs...) -function S_z(elt::Type{<:Number}, ::Type{Trivial}; spin=1 // 2) +function S_z(elt::Type{<:Number}, ::Type{Trivial}; spin = 1 // 2) _, _, S_z_mat = spinmatrices(spin, elt) pspace = ComplexSpace(size(S_z_mat, 1)) return TensorMap(S_z_mat, pspace ← pspace) end -function S_z(elt::Type{<:Number}, ::Type{Z2Irrep}; spin=1 // 2, side=:L) +function S_z(elt::Type{<:Number}, ::Type{Z2Irrep}; spin = 1 // 2, side = :L) spin == 1 // 2 || error("Z2 symmetry only implemented for spin 1 // 2") pspace = Z2Space(0 => 1, 1 => 1) vspace = Z2Space(1 => 1) @@ -219,7 +221,7 @@ function S_z(elt::Type{<:Number}, ::Type{Z2Irrep}; spin=1 // 2, side=:L) return Z end -function S_z(elt::Type{<:Number}, ::Type{U1Irrep}; spin=1 // 2) +function S_z(elt::Type{<:Number}, ::Type{U1Irrep}; spin = 1 // 2) charges = U1Irrep.((-spin):spin) pspace = U1Space((v => 1 for v in charges)) Z = zeros(elt, pspace ← pspace) @@ -251,12 +253,12 @@ S_plus(; kwargs...) = S_plus(ComplexF64, Trivial; kwargs...) S_plus(elt::Type{<:Number}; kwargs...) = S_plus(elt, Trivial; kwargs...) S_plus(symm::Type{<:Sector}; kwargs...) = S_plus(ComplexF64, symm; kwargs...) -function S_plus(elt::Type{<:Number}, ::Type{Trivial}; spin=1 // 2) - S⁺ = S_x(elt, Trivial; spin=spin) + 1im * S_y(complex(elt), Trivial; spin=spin) +function S_plus(elt::Type{<:Number}, ::Type{Trivial}; spin = 1 // 2) + S⁺ = S_x(elt, Trivial; spin = spin) + 1im * S_y(complex(elt), Trivial; spin = spin) return elt <: Real ? real(S⁺) : S⁺ end -function S_plus(elt::Type{<:Number}, ::Type{Z2Irrep}; spin=1 // 2, side=:L) +function S_plus(elt::Type{<:Number}, ::Type{Z2Irrep}; spin = 1 // 2, side = :L) spin == 1 // 2 || error("Z2 symmetry only implemented for spin 1 // 2") pspace = Z2Space(0 => 1, 1 => 1) vspace = Z2Space(0 => 1, 1 => 1) @@ -274,7 +276,7 @@ function S_plus(elt::Type{<:Number}, ::Type{Z2Irrep}; spin=1 // 2, side=:L) return S⁺ end -function S_plus(elt::Type{<:Number}, ::Type{U1Irrep}; spin=1 // 2, side=:L) +function S_plus(elt::Type{<:Number}, ::Type{U1Irrep}; spin = 1 // 2, side = :L) pspace = U1Space(i => 1 for i in (-spin):spin) if side == :L vspace = U1Space(1 => 1) @@ -316,12 +318,12 @@ S_min(; kwargs...) = S_min(ComplexF64, Trivial; kwargs...) S_min(elt::Type{<:Number}; kwargs...) = S_min(elt, Trivial; kwargs...) S_min(symm::Type{<:Sector}; kwargs...) = S_min(ComplexF64, symm; kwargs...) -function S_min(elt::Type{<:Number}, ::Type{Trivial}; spin=1 // 2) - S⁻ = S_x(elt, Trivial; spin=spin) - 1im * S_y(complex(elt), Trivial; spin=spin) +function S_min(elt::Type{<:Number}, ::Type{Trivial}; spin = 1 // 2) + S⁻ = S_x(elt, Trivial; spin = spin) - 1im * S_y(complex(elt), Trivial; spin = spin) return elt <: Real ? real(S⁻) : S⁻ end -function S_min(elt::Type{<:Number}, ::Type{Z2Irrep}; spin=1 // 2, side=:L) +function S_min(elt::Type{<:Number}, ::Type{Z2Irrep}; spin = 1 // 2, side = :L) spin == 1 // 2 || error("Z2 symmetry only implemented for spin 1 // 2") pspace = Z2Space(0 => 1, 1 => 1) vspace = Z2Space(0 => 1, 1 => 1) @@ -339,7 +341,7 @@ function S_min(elt::Type{<:Number}, ::Type{Z2Irrep}; spin=1 // 2, side=:L) return S⁻ end -function S_min(elt::Type{<:Number}, ::Type{U1Irrep}; spin=1 // 2, side=:L) +function S_min(elt::Type{<:Number}, ::Type{U1Irrep}; spin = 1 // 2, side = :L) pspace = U1Space(i => 1 for i in (-spin):spin) if side == :L vspace = U1Space(-1 => 1) @@ -402,14 +404,18 @@ for (L, R) in ((:x, :x), (:y, :y), (:z, :z), (:plus, :min), (:min, :plus)) ($f)(elt::Type{<:Number}; kwargs...) = ($f)(elt, Trivial; kwargs...) ($f)(symmetry::Type{<:Sector}; kwargs...) = ($f)(ComplexF64, symmetry; kwargs...) - function ($f)(elt::Type{<:Number}, ::Type{Trivial}; spin=1 // 2) - return contract_twosite($(fₗ)(elt, Trivial; spin=spin), - $(fᵣ)(elt, Trivial; spin=spin)) + function ($f)(elt::Type{<:Number}, ::Type{Trivial}; spin = 1 // 2) + return contract_twosite( + $(fₗ)(elt, Trivial; spin = spin), + $(fᵣ)(elt, Trivial; spin = spin) + ) end - function ($f)(elt::Type{<:Number}, symmetry::Type{<:Sector}; spin=1 // 2) - return contract_twosite($(fₗ)(elt, symmetry; spin=spin, side=:L), - $(fᵣ)(elt, symmetry; spin=spin, side=:R)) + function ($f)(elt::Type{<:Number}, symmetry::Type{<:Sector}; spin = 1 // 2) + return contract_twosite( + $(fₗ)(elt, symmetry; spin = spin, side = :L), + $(fᵣ)(elt, symmetry; spin = spin, side = :R) + ) end const $f_unicode = $f @@ -420,11 +426,11 @@ for (L, R) in ((:x, :x), (:y, :y), (:z, :z), (:plus, :min), (:min, :plus)) end end -function S_xx(elt::Type{<:Number}, ::Type{Z2Irrep}; spin=1 // 2) - return contract_twosite(S_x(elt, Z2Irrep; spin=spin), S_x(elt, Z2Irrep; spin=spin)) +function S_xx(elt::Type{<:Number}, ::Type{Z2Irrep}; spin = 1 // 2) + return contract_twosite(S_x(elt, Z2Irrep; spin = spin), S_x(elt, Z2Irrep; spin = spin)) end -function S_zz(elt::Type{<:Number}, ::Type{U1Irrep}; spin=1 // 2) - return contract_twosite(S_z(elt, U1Irrep; spin=spin), S_z(elt, U1Irrep; spin=spin)) +function S_zz(elt::Type{<:Number}, ::Type{U1Irrep}; spin = 1 // 2) + return contract_twosite(S_z(elt, U1Irrep; spin = spin), S_z(elt, U1Irrep; spin = spin)) end """ @@ -441,14 +447,16 @@ S_exchange(elt::Type{<:Number}; kwargs...) = S_exchange(elt, Trivial; kwargs...) function S_exchange(symmetry::Type{<:Sector}; kwargs...) return S_exchange(ComplexF64, symmetry; kwargs...) end -function S_exchange(elt::Type{<:Number}, symmetry::Type{<:Sector}; spin=1 // 2) +function S_exchange(elt::Type{<:Number}, symmetry::Type{<:Sector}; spin = 1 // 2) elt_complex = complex(elt) - SS = (S_plusmin(elt_complex, symmetry; spin=spin) + - S_minplus(elt_complex, symmetry; spin=spin)) / 2 + - S_zz(elt_complex, symmetry; spin=spin) + SS = ( + S_plusmin(elt_complex, symmetry; spin = spin) + + S_minplus(elt_complex, symmetry; spin = spin) + ) / 2 + + S_zz(elt_complex, symmetry; spin = spin) return elt <: Real ? real(SS) : SS end -function S_exchange(elt::Type{<:Number}, ::Type{SU2Irrep}; spin=1 // 2) +function S_exchange(elt::Type{<:Number}, ::Type{SU2Irrep}; spin = 1 // 2) pspace = SU2Space(spin => 1) aspace = SU2Space(1 => 1) @@ -480,12 +488,12 @@ function potts_ZZ(symmetry::Type{<:Sector}; kwargs...) return potts_ZZ(ComplexF64, symmetry; kwargs...) end -function potts_ZZ(elt::Type{<:Number}, ::Type{Trivial}; q=3) - Z = potts_Z(elt, Trivial; q=q) +function potts_ZZ(elt::Type{<:Number}, ::Type{Trivial}; q = 3) + Z = potts_Z(elt, Trivial; q = q) return Z ⊗ Z' end -function potts_ZZ(elt::Type{<:Number}, ::Type{ZNIrrep{Q}}; q=Q) where {Q} +function potts_ZZ(elt::Type{<:Number}, ::Type{ZNIrrep{Q}}; q = Q) where {Q} @assert q == Q "q must match the irrep charge" pspace = Vect[ZNIrrep{Q}](i => 1 for i in 0:(Q - 1)) ZZ = zeros(elt, pspace ⊗ pspace ← pspace ⊗ pspace) @@ -502,7 +510,7 @@ end the Weyl-Heisenberg matrices according to [Wikipedia](https://en.wikipedia.org/wiki/Generalizations_of_Pauli_matrices#Sylvester's_generalized_Pauli_matrices_(non-Hermitian)). """ -function weyl_heisenberg_matrices(Q::Int, elt=ComplexF64) +function weyl_heisenberg_matrices(Q::Int, elt = ComplexF64) U = zeros(elt, Q, Q) # clock matrix V = zeros(elt, Q, Q) # shift matrix W = zeros(elt, Q, Q) # DFT @@ -527,7 +535,7 @@ function potts_Z end potts_Z(; kwargs...) = potts_Z(ComplexF64, Trivial; kwargs...) potts_Z(elt::Type{<:Complex}; kwargs...) = potts_Z(elt, Trivial; kwargs...) potts_Z(symm::Type{<:Sector}; kwargs...) = potts_Z(ComplexF64, symm; kwargs...) -function potts_Z(elt::Type{<:Number}, ::Type{Trivial}; q=3) +function potts_Z(elt::Type{<:Number}, ::Type{Trivial}; q = 3) U, _, _ = weyl_heisenberg_matrices(q, elt) Z = TensorMap(U, ComplexSpace(q) ← ComplexSpace(q)) return Z @@ -543,13 +551,13 @@ function potts_X end potts_X(; kwargs...) = potts_X(ComplexF64, Trivial; kwargs...) potts_X(elt::Type{<:Complex}; kwargs...) = potts_X(elt, Trivial; kwargs...) potts_X(symm::Type{<:Sector}; kwargs...) = potts_X(ComplexF64, symm; kwargs...) -function potts_X(elt::Type{<:Number}, ::Type{Trivial}; q=3) +function potts_X(elt::Type{<:Number}, ::Type{Trivial}; q = 3) _, V, _ = weyl_heisenberg_matrices(q, elt) X = TensorMap(V, ComplexSpace(q) ← ComplexSpace(q)) return X end -function potts_X(elt::Type{<:Number}, ::Type{ZNIrrep{Q}}; q=Q) where {Q} +function potts_X(elt::Type{<:Number}, ::Type{ZNIrrep{Q}}; q = Q) where {Q} @assert q == Q "q must match the irrep charge" pspace = Vect[ZNIrrep{Q}](i => 1 for i in 0:(Q - 1)) X = zeros(elt, pspace ← pspace) diff --git a/src/operators/tjoperators.jl b/src/operators/tjoperators.jl index 96b796d..298c732 100644 --- a/src/operators/tjoperators.jl +++ b/src/operators/tjoperators.jl @@ -33,49 +33,59 @@ The possible symmetries are Setting `slave_fermion = true` switches to the slave-fermion basis. """ -function tj_space(::Type{Trivial}=Trivial, ::Type{Trivial}=Trivial; - slave_fermion::Bool=false) +function tj_space( + ::Type{Trivial} = Trivial, ::Type{Trivial} = Trivial; + slave_fermion::Bool = false + ) return slave_fermion ? Vect[FermionParity](0 => 2, 1 => 1) : - Vect[FermionParity](0 => 1, 1 => 2) + Vect[FermionParity](0 => 1, 1 => 2) end -function tj_space(::Type{Trivial}, ::Type{U1Irrep}; slave_fermion::Bool=false) +function tj_space(::Type{Trivial}, ::Type{U1Irrep}; slave_fermion::Bool = false) return if slave_fermion Vect[FermionParity ⊠ U1Irrep]((1, 0) => 1, (0, 1 // 2) => 1, (0, -1 // 2) => 1) else Vect[FermionParity ⊠ U1Irrep]((0, 0) => 1, (1, 1 // 2) => 1, (1, -1 // 2) => 1) end end -function tj_space(::Type{Trivial}, ::Type{SU2Irrep}; slave_fermion::Bool=false) +function tj_space(::Type{Trivial}, ::Type{SU2Irrep}; slave_fermion::Bool = false) return error("Not implemented") end -function tj_space(::Type{U1Irrep}, ::Type{Trivial}; slave_fermion::Bool=false) +function tj_space(::Type{U1Irrep}, ::Type{Trivial}; slave_fermion::Bool = false) return if slave_fermion Vect[FermionParity ⊠ U1Irrep]((1, 0) => 1, (0, 1) => 2) else Vect[FermionParity ⊠ U1Irrep]((0, 0) => 1, (1, 1) => 2) end end -function tj_space(::Type{U1Irrep}, ::Type{U1Irrep}; slave_fermion::Bool=false) +function tj_space(::Type{U1Irrep}, ::Type{U1Irrep}; slave_fermion::Bool = false) return if slave_fermion - Vect[FermionParity ⊠ U1Irrep ⊠ U1Irrep]((1, 0, 0) => 1, (0, 1, 1 // 2) => 1, - (0, 1, -1 // 2) => 1) + Vect[FermionParity ⊠ U1Irrep ⊠ U1Irrep]( + (1, 0, 0) => 1, (0, 1, 1 // 2) => 1, + (0, 1, -1 // 2) => 1 + ) else - Vect[FermionParity ⊠ U1Irrep ⊠ U1Irrep]((0, 0, 0) => 1, (1, 1, 1 // 2) => 1, - (1, 1, -1 // 2) => 1) + Vect[FermionParity ⊠ U1Irrep ⊠ U1Irrep]( + (0, 0, 0) => 1, (1, 1, 1 // 2) => 1, + (1, 1, -1 // 2) => 1 + ) end end -function tj_space(::Type{U1Irrep}, ::Type{SU2Irrep}; slave_fermion::Bool=false) +function tj_space(::Type{U1Irrep}, ::Type{SU2Irrep}; slave_fermion::Bool = false) return error("Not implemented") end -function single_site_operator(T, particle_symmetry::Type{<:Sector}, - spin_symmetry::Type{<:Sector}; slave_fermion::Bool=false) +function single_site_operator( + T, particle_symmetry::Type{<:Sector}, + spin_symmetry::Type{<:Sector}; slave_fermion::Bool = false + ) V = tj_space(particle_symmetry, spin_symmetry; slave_fermion) return zeros(T, V ← V) end -function two_site_operator(T, particle_symmetry::Type{<:Sector}, - spin_symmetry::Type{<:Sector}; slave_fermion::Bool=false) +function two_site_operator( + T, particle_symmetry::Type{<:Sector}, + spin_symmetry::Type{<:Sector}; slave_fermion::Bool = false + ) V = tj_space(particle_symmetry, spin_symmetry; slave_fermion) return zeros(T, V ⊗ V ← V ⊗ V) end @@ -86,10 +96,10 @@ end Return the two-body operator ``e†_{1,↑}, e_{2,↑}`` that creates a spin-up electron at the first site and annihilates a spin-up electron at the second. The only nonzero matrix element corresponds to `|↑0⟩ <-- |0↑⟩`. """ -function e_plusmin_up(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool=false) +function e_plusmin_up(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool = false) return e_plusmin_up(ComplexF64, P, S; slave_fermion) end -function e_plusmin_up(T, ::Type{Trivial}, ::Type{Trivial}; slave_fermion::Bool=false) +function e_plusmin_up(T, ::Type{Trivial}, ::Type{Trivial}; slave_fermion::Bool = false) t = two_site_operator(T, Trivial, Trivial; slave_fermion) I = sectortype(t) (h, b, sgn) = slave_fermion ? (1, 0, -1) : (0, 1, 1) @@ -103,31 +113,31 @@ function e_plusmin_up(T, ::Type{Trivial}, ::Type{Trivial}; slave_fermion::Bool=f t[(I(b), I(h), dual(I(h)), dual(I(b)))][1, 1, 1, 1] = sgn * 1 return t end -function e_plusmin_up(T, ::Type{Trivial}, ::Type{U1Irrep}; slave_fermion::Bool=false) +function e_plusmin_up(T, ::Type{Trivial}, ::Type{U1Irrep}; slave_fermion::Bool = false) t = two_site_operator(T, Trivial, U1Irrep; slave_fermion) I = sectortype(t) (h, b, sgn) = slave_fermion ? (1, 0, -1) : (0, 1, 1) t[(I(b, 1 // 2), I(h, 0), dual(I(h, 0)), dual(I(b, 1 // 2)))] .= sgn * 1 return t end -function e_plusmin_up(T, ::Type{Trivial}, ::Type{SU2Irrep}; slave_fermion::Bool=false) +function e_plusmin_up(T, ::Type{Trivial}, ::Type{SU2Irrep}; slave_fermion::Bool = false) return error("Not implemented") end -function e_plusmin_up(T, ::Type{U1Irrep}, ::Type{Trivial}; slave_fermion::Bool=false) +function e_plusmin_up(T, ::Type{U1Irrep}, ::Type{Trivial}; slave_fermion::Bool = false) t = two_site_operator(T, U1Irrep, Trivial; slave_fermion) I = sectortype(t) (h, b, sgn) = slave_fermion ? (1, 0, -1) : (0, 1, 1) t[(I(b, 1), I(h, 0), dual(I(h, 0)), dual(I(b, 1)))][1, 1, 1, 1] = sgn * 1 return t end -function e_plusmin_up(T, ::Type{U1Irrep}, ::Type{U1Irrep}; slave_fermion::Bool=false) +function e_plusmin_up(T, ::Type{U1Irrep}, ::Type{U1Irrep}; slave_fermion::Bool = false) t = two_site_operator(T, U1Irrep, U1Irrep; slave_fermion) I = sectortype(t) (h, b, sgn) = slave_fermion ? (1, 0, -1) : (0, 1, 1) t[(I(b, 1, 1 // 2), I(h, 0, 0), dual(I(h, 0, 0)), dual(I(b, 1, 1 // 2)))] .= sgn * 1 return t end -function e_plusmin_up(T, ::Type{U1Irrep}, ::Type{SU2Irrep}; slave_fermion::Bool=false) +function e_plusmin_up(T, ::Type{U1Irrep}, ::Type{SU2Irrep}; slave_fermion::Bool = false) return error("Not implemented") end const e⁺e⁻ꜛ = e_plusmin_up @@ -138,41 +148,41 @@ const e⁺e⁻ꜛ = e_plusmin_up Return the two-body operator ``e†_{1,↓}, e_{2,↓}`` that creates a spin-down electron at the first site and annihilates a spin-down electron at the second. The only nonzero matrix element corresponds to `|↓0⟩ <-- |0↓⟩`. """ -function e_plusmin_down(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool=false) +function e_plusmin_down(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool = false) return e_plusmin_down(ComplexF64, P, S; slave_fermion) end -function e_plusmin_down(T, ::Type{Trivial}, ::Type{Trivial}; slave_fermion::Bool=false) +function e_plusmin_down(T, ::Type{Trivial}, ::Type{Trivial}; slave_fermion::Bool = false) t = two_site_operator(T, Trivial, Trivial; slave_fermion) I = sectortype(t) (h, b, sgn) = slave_fermion ? (1, 0, -1) : (0, 1, 1) t[(I(b), I(h), dual(I(h)), dual(I(b)))][2, 1, 1, 2] = sgn * 1 return t end -function e_plusmin_down(T, ::Type{Trivial}, ::Type{U1Irrep}; slave_fermion::Bool=false) +function e_plusmin_down(T, ::Type{Trivial}, ::Type{U1Irrep}; slave_fermion::Bool = false) t = two_site_operator(T, Trivial, U1Irrep; slave_fermion) I = sectortype(t) (h, b, sgn) = slave_fermion ? (1, 0, -1) : (0, 1, 1) t[(I(b, -1 // 2), I(h, 0), dual(I(h, 0)), dual(I(b, -1 // 2)))] .= sgn * 1 return t end -function e_plusmin_down(T, ::Type{Trivial}, ::Type{SU2Irrep}; slave_fermion::Bool=false) +function e_plusmin_down(T, ::Type{Trivial}, ::Type{SU2Irrep}; slave_fermion::Bool = false) return error("Not implemented") end -function e_plusmin_down(T, ::Type{U1Irrep}, ::Type{Trivial}; slave_fermion::Bool=false) +function e_plusmin_down(T, ::Type{U1Irrep}, ::Type{Trivial}; slave_fermion::Bool = false) t = two_site_operator(T, U1Irrep, Trivial; slave_fermion) I = sectortype(t) (h, b, sgn) = slave_fermion ? (1, 0, -1) : (0, 1, 1) t[(I(b, 1), I(h, 0), dual(I(h, 0)), dual(I(b, 1)))][2, 1, 1, 2] = sgn * 1 return t end -function e_plusmin_down(T, ::Type{U1Irrep}, ::Type{U1Irrep}; slave_fermion::Bool=false) +function e_plusmin_down(T, ::Type{U1Irrep}, ::Type{U1Irrep}; slave_fermion::Bool = false) t = two_site_operator(T, U1Irrep, U1Irrep; slave_fermion) I = sectortype(t) (h, b, sgn) = slave_fermion ? (1, 0, -1) : (0, 1, 1) t[(I(b, 1, -1 // 2), I(h, 0, 0), dual(I(h, 0, 0)), dual(I(b, 1, -1 // 2)))] .= sgn * 1 return t end -function e_plusmin_down(T, ::Type{U1Irrep}, ::Type{SU2Irrep}; slave_fermion::Bool=false) +function e_plusmin_down(T, ::Type{U1Irrep}, ::Type{SU2Irrep}; slave_fermion::Bool = false) return error("Not implemented") end const e⁺e⁻ꜜ = e_plusmin_down @@ -185,11 +195,13 @@ Return the Hermitian conjugate of `e_plusmin_up`, i.e. It annihilates a spin-up electron at the first site and creates a spin-up electron at the second. The only nonzero matrix element corresponds to `|0↑⟩ <-- |↑0⟩`. """ -function e_minplus_up(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool=false) +function e_minplus_up(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool = false) return e_minplus_up(ComplexF64, P, S; slave_fermion) end -function e_minplus_up(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; - slave_fermion::Bool=false) +function e_minplus_up( + T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; + slave_fermion::Bool = false + ) return copy(adjoint(e_plusmin_up(T, particle_symmetry, spin_symmetry; slave_fermion))) end const e⁻⁺ꜛ = e_minplus_up @@ -202,11 +214,13 @@ Return the Hermitian conjugate of `e_plusmin_down`, i.e. It annihilates a spin-down electron at the first site and creates a spin-down electron at the second. The only nonzero matrix element corresponds to `|0↓⟩ <-- |↓0⟩`. """ -function e_minplus_down(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool=false) +function e_minplus_down(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool = false) return e_minplus_down(ComplexF64, P, S; slave_fermion) end -function e_minplus_down(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; - slave_fermion::Bool=false) +function e_minplus_down( + T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; + slave_fermion::Bool = false + ) return copy(adjoint(e_plusmin_down(T, particle_symmetry, spin_symmetry; slave_fermion))) end const e⁻e⁺ꜜ = e_minplus_down @@ -217,13 +231,15 @@ const e⁻e⁺ꜜ = e_minplus_down Return the two-body operator that creates a particle at the first site and annihilates a particle at the second. This is the sum of `e_plusmin_up` and `e_plusmin_down`. """ -function e_plusmin(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool=false) +function e_plusmin(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool = false) return e_plusmin(ComplexF64, P, S; slave_fermion) end -function e_plusmin(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; - slave_fermion::Bool=false) +function e_plusmin( + T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; + slave_fermion::Bool = false + ) return e_plusmin_up(T, particle_symmetry, spin_symmetry; slave_fermion) + - e_plusmin_down(T, particle_symmetry, spin_symmetry; slave_fermion) + e_plusmin_down(T, particle_symmetry, spin_symmetry; slave_fermion) end const e⁺e⁻ = e_plusmin @@ -233,11 +249,13 @@ const e⁺e⁻ = e_plusmin Return the two-body operator that annihilates a particle at the first site and creates a particle at the second. This is the sum of `e_minplus_up` and `e_minplus_down`. """ -function e_minplus(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool=false) +function e_minplus(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool = false) return e_minplus(ComplexF64, P, S; slave_fermion) end -function e_minplus(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; - slave_fermion::Bool=false) +function e_minplus( + T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; + slave_fermion::Bool = false + ) return copy(adjoint(e_plusmin(T, particle_symmetry, spin_symmetry; slave_fermion))) end const e⁻e⁺ = e_minplus @@ -248,24 +266,24 @@ const e⁻e⁺ = e_minplus Return the two-body operator ``e_{1,↑} e_{2,↓}`` that annihilates a spin-up particle at the first site and a spin-down particle at the second site. The only nonzero matrix element corresponds to `|00⟩ <-- |↑↓⟩`. """ -function e_minmin_updown(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool=false) +function e_minmin_updown(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool = false) return e_minmin_updown(ComplexF64, P, S; slave_fermion) end -function e_minmin_updown(T, ::Type{Trivial}, ::Type{Trivial}; slave_fermion::Bool=false) +function e_minmin_updown(T, ::Type{Trivial}, ::Type{Trivial}; slave_fermion::Bool = false) t = two_site_operator(T, Trivial, Trivial; slave_fermion) I = sectortype(t) (h, b, sgn) = slave_fermion ? (1, 0, -1) : (0, 1, 1) t[(I(h), I(h), dual(I(b)), dual(I(b)))][1, 1, 1, 2] = -sgn * 1 return t end -function e_minmin_updown(T, ::Type{Trivial}, ::Type{U1Irrep}; slave_fermion::Bool=false) +function e_minmin_updown(T, ::Type{Trivial}, ::Type{U1Irrep}; slave_fermion::Bool = false) t = two_site_operator(T, Trivial, U1Irrep; slave_fermion) I = sectortype(t) (h, b, sgn) = slave_fermion ? (1, 0, -1) : (0, 1, 1) t[(I(h, 0), I(h, 0), dual(I(b, 1 // 2)), dual(I(b, -1 // 2)))] .= -sgn * 1 return t end -function e_minmin_updown(T, ::Type{U1Irrep}, ::Type{<:Sector}; slave_fermion::Bool=false) +function e_minmin_updown(T, ::Type{U1Irrep}, ::Type{<:Sector}; slave_fermion::Bool = false) throw(ArgumentError("`e_minmin_updown` is not symmetric under `U1Irrep` particle symmetry")) end @@ -275,24 +293,24 @@ end Return the two-body operator ``e_{1,↓} e_{2,↑}`` that annihilates a spin-down particle at the first site and a spin-up particle at the second site. The only nonzero matrix element corresponds to `|00⟩ <-- |↓↑⟩`. """ -function e_minmin_downup(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool=false) +function e_minmin_downup(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool = false) return e_minmin_downup(ComplexF64, P, S; slave_fermion) end -function e_minmin_downup(T, ::Type{Trivial}, ::Type{Trivial}; slave_fermion::Bool=false) +function e_minmin_downup(T, ::Type{Trivial}, ::Type{Trivial}; slave_fermion::Bool = false) t = two_site_operator(T, Trivial, Trivial; slave_fermion) I = sectortype(t) (h, b, sgn) = slave_fermion ? (1, 0, -1) : (0, 1, 1) t[(I(h), I(h), dual(I(b)), dual(I(b)))][1, 1, 2, 1] = -sgn * 1 return t end -function e_minmin_downup(T, ::Type{Trivial}, ::Type{U1Irrep}; slave_fermion::Bool=false) +function e_minmin_downup(T, ::Type{Trivial}, ::Type{U1Irrep}; slave_fermion::Bool = false) t = two_site_operator(T, Trivial, U1Irrep; slave_fermion) I = sectortype(t) (h, b, sgn) = slave_fermion ? (1, 0, -1) : (0, 1, 1) t[(I(h, 0), I(h, 0), dual(I(b, -1 // 2)), dual(I(b, 1 // 2)))] .= -sgn * 1 return t end -function e_minmin_downup(T, ::Type{U1Irrep}, ::Type{<:Sector}; slave_fermion::Bool=false) +function e_minmin_downup(T, ::Type{U1Irrep}, ::Type{<:Sector}; slave_fermion::Bool = false) throw(ArgumentError("`e_minmin_downup` is not symmetric under `U1Irrep` particle symmetry")) end @@ -301,13 +319,17 @@ end Return the two-body singlet operator ``(e_{1,↓} e_{2,↑} - e_{1,↓} e_{2,↑}) / sqrt(2)``. """ -function e_singlet(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool=false) +function e_singlet(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool = false) return e_singlet(ComplexF64, P, S; slave_fermion) end -function e_singlet(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; - slave_fermion::Bool=false) - return (e_minmin_updown(T, particle_symmetry, spin_symmetry; slave_fermion) - - e_minmin_downup(T, particle_symmetry, spin_symmetry; slave_fermion)) / sqrt(2) +function e_singlet( + T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; + slave_fermion::Bool = false + ) + return ( + e_minmin_updown(T, particle_symmetry, spin_symmetry; slave_fermion) - + e_minmin_downup(T, particle_symmetry, spin_symmetry; slave_fermion) + ) / sqrt(2) end """ @@ -315,42 +337,44 @@ end Return the one-body operator that counts the number of spin-up electrons. """ -function e_number_up(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool=false) +function e_number_up(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool = false) return e_number_up(ComplexF64, P, S; slave_fermion) end -function e_number_up(T::Type{<:Number}, ::Type{Trivial}, ::Type{Trivial}; - slave_fermion::Bool=false) +function e_number_up( + T::Type{<:Number}, ::Type{Trivial}, ::Type{Trivial}; + slave_fermion::Bool = false + ) t = single_site_operator(T, Trivial, Trivial; slave_fermion) I = sectortype(t) b = slave_fermion ? 0 : 1 t[(I(b), dual(I(b)))][1, 1] = 1 return t end -function e_number_up(T, ::Type{Trivial}, ::Type{U1Irrep}; slave_fermion::Bool=false) +function e_number_up(T, ::Type{Trivial}, ::Type{U1Irrep}; slave_fermion::Bool = false) t = single_site_operator(T, Trivial, U1Irrep; slave_fermion) I = sectortype(t) b = slave_fermion ? 0 : 1 t[(I(b, 1 // 2), dual(I(b, 1 // 2)))][1, 1] = 1 return t end -function e_number_up(T, ::Type{Trivial}, ::Type{SU2Irrep}; slave_fermion::Bool=false) +function e_number_up(T, ::Type{Trivial}, ::Type{SU2Irrep}; slave_fermion::Bool = false) throw(ArgumentError("`e_number_up` is not symmetric under `SU2Irrep` spin symmetry")) end -function e_number_up(T, ::Type{U1Irrep}, ::Type{Trivial}; slave_fermion::Bool=false) +function e_number_up(T, ::Type{U1Irrep}, ::Type{Trivial}; slave_fermion::Bool = false) t = single_site_operator(T, U1Irrep, Trivial; slave_fermion) I = sectortype(t) b = slave_fermion ? 0 : 1 t[(I(b, 1), dual(I(b, 1)))][1, 1] = 1 return t end -function e_number_up(T, ::Type{U1Irrep}, ::Type{U1Irrep}; slave_fermion::Bool=false) +function e_number_up(T, ::Type{U1Irrep}, ::Type{U1Irrep}; slave_fermion::Bool = false) t = single_site_operator(T, U1Irrep, U1Irrep; slave_fermion) I = sectortype(t) b = slave_fermion ? 0 : 1 t[(I(b, 1, 1 // 2), dual(I(b, 1, 1 // 2)))] .= 1 return t end -function e_number_up(T, ::Type{U1Irrep}, ::Type{SU2Irrep}; slave_fermion::Bool=false) +function e_number_up(T, ::Type{U1Irrep}, ::Type{SU2Irrep}; slave_fermion::Bool = false) throw(ArgumentError("`e_number_up` is not symmetric under `SU2Irrep` spin symmetry")) end const nꜛ = e_number_up @@ -360,42 +384,44 @@ const nꜛ = e_number_up Return the one-body operator that counts the number of spin-down electrons. """ -function e_number_down(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool=false) +function e_number_down(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool = false) return e_number_down(ComplexF64, P, S; slave_fermion) end -function e_number_down(T::Type{<:Number}, ::Type{Trivial}, ::Type{Trivial}; - slave_fermion::Bool=false) +function e_number_down( + T::Type{<:Number}, ::Type{Trivial}, ::Type{Trivial}; + slave_fermion::Bool = false + ) t = single_site_operator(T, Trivial, Trivial; slave_fermion) I = sectortype(t) b = slave_fermion ? 0 : 1 t[(I(b), dual(I(b)))][2, 2] = 1 return t end -function e_number_down(T, ::Type{Trivial}, ::Type{U1Irrep}; slave_fermion::Bool=false) +function e_number_down(T, ::Type{Trivial}, ::Type{U1Irrep}; slave_fermion::Bool = false) t = single_site_operator(T, Trivial, U1Irrep; slave_fermion) I = sectortype(t) b = slave_fermion ? 0 : 1 t[(I(b, -1 // 2), dual(I(b, -1 // 2)))][1, 1] = 1 return t end -function e_number_down(T, ::Type{Trivial}, ::Type{SU2Irrep}; slave_fermion::Bool=false) +function e_number_down(T, ::Type{Trivial}, ::Type{SU2Irrep}; slave_fermion::Bool = false) throw(ArgumentError("`e_number_down` is not symmetric under `SU2Irrep` spin symmetry")) end -function e_number_down(T, ::Type{U1Irrep}, ::Type{Trivial}; slave_fermion::Bool=false) +function e_number_down(T, ::Type{U1Irrep}, ::Type{Trivial}; slave_fermion::Bool = false) t = single_site_operator(T, U1Irrep, Trivial; slave_fermion) I = sectortype(t) b = slave_fermion ? 0 : 1 t[(I(b, 1), dual(I(b, 1)))][2, 2] = 1 return t end -function e_number_down(T, ::Type{U1Irrep}, ::Type{U1Irrep}; slave_fermion::Bool=false) +function e_number_down(T, ::Type{U1Irrep}, ::Type{U1Irrep}; slave_fermion::Bool = false) t = single_site_operator(T, U1Irrep, U1Irrep; slave_fermion) I = sectortype(t) b = slave_fermion ? 0 : 1 t[(I(b, 1, -1 // 2), dual(I(b, 1, -1 // 2)))] .= 1 return t end -function e_number_down(T, ::Type{U1Irrep}, ::Type{SU2Irrep}; slave_fermion::Bool=false) +function e_number_down(T, ::Type{U1Irrep}, ::Type{SU2Irrep}; slave_fermion::Bool = false) throw(ArgumentError("`e_number_down` is not symmetric under `SU2Irrep` spin symmetry")) end const nꜜ = e_number_down @@ -405,13 +431,15 @@ const nꜜ = e_number_down Return the one-body operator that counts the number of particles. """ -function e_number(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool=false) +function e_number(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool = false) return e_number(ComplexF64, P, S; slave_fermion) end -function e_number(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; - slave_fermion::Bool=false) +function e_number( + T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; + slave_fermion::Bool = false + ) return e_number_up(T, particle_symmetry, spin_symmetry; slave_fermion) + - e_number_down(T, particle_symmetry, spin_symmetry; slave_fermion) + e_number_down(T, particle_symmetry, spin_symmetry; slave_fermion) end const n = e_number @@ -420,11 +448,13 @@ const n = e_number Return the one-body operator that counts the number of holes. """ -function e_number_hole(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool=false) +function e_number_hole(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool = false) return e_number_hole(ComplexF64, P, S; slave_fermion) end -function e_number_hole(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; - slave_fermion::Bool=false) +function e_number_hole( + T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; + slave_fermion::Bool = false + ) iden = TensorKit.id(tj_space(particle_symmetry, spin_symmetry; slave_fermion)) return iden - e_number(T, particle_symmetry, spin_symmetry; slave_fermion) end @@ -435,11 +465,13 @@ const nʰ = e_number_hole Return the one-body spin-1/2 x-operator on the electrons. """ -function S_x(P::Type{<:Sector}=Trivial, S::Type{<:Sector}=Trivial; - slave_fermion::Bool=false) +function S_x( + P::Type{<:Sector} = Trivial, S::Type{<:Sector} = Trivial; + slave_fermion::Bool = false + ) return S_x(ComplexF64, P, S; slave_fermion) end -function S_x(T::Type{<:Number}, ::Type{Trivial}, ::Type{Trivial}; slave_fermion::Bool=false) +function S_x(T::Type{<:Number}, ::Type{Trivial}, ::Type{Trivial}; slave_fermion::Bool = false) t = single_site_operator(T, Trivial, Trivial; slave_fermion) I = sectortype(t) b = slave_fermion ? 0 : 1 @@ -447,7 +479,7 @@ function S_x(T::Type{<:Number}, ::Type{Trivial}, ::Type{Trivial}; slave_fermion: t[(I(b), dual(I(b)))][2, 1] = 0.5 return t end -function S_x(T::Type{<:Number}, ::Type{U1Irrep}, ::Type{Trivial}; slave_fermion::Bool=false) +function S_x(T::Type{<:Number}, ::Type{U1Irrep}, ::Type{Trivial}; slave_fermion::Bool = false) t = single_site_operator(T, U1Irrep, Trivial; slave_fermion) I = sectortype(t) b = slave_fermion ? 0 : 1 @@ -462,11 +494,13 @@ const Sˣ = S_x Return the one-body spin-1/2 x-operator on the electrons (only defined for `Trivial` symmetry). """ -function S_y(P::Type{<:Sector}=Trivial, S::Type{<:Sector}=Trivial; - slave_fermion::Bool=false) +function S_y( + P::Type{<:Sector} = Trivial, S::Type{<:Sector} = Trivial; + slave_fermion::Bool = false + ) return S_y(ComplexF64, P, S; slave_fermion) end -function S_y(T::Type{<:Number}, ::Type{Trivial}, ::Type{Trivial}; slave_fermion::Bool=false) +function S_y(T::Type{<:Number}, ::Type{Trivial}, ::Type{Trivial}; slave_fermion::Bool = false) t = single_site_operator(T, Trivial, Trivial; slave_fermion) I = sectortype(t) b = slave_fermion ? 0 : 1 @@ -474,7 +508,7 @@ function S_y(T::Type{<:Number}, ::Type{Trivial}, ::Type{Trivial}; slave_fermion: t[(I(b), dual(I(b)))][2, 1] = 0.5im return t end -function S_y(T::Type{<:Number}, ::Type{U1Irrep}, ::Type{Trivial}; slave_fermion::Bool=false) +function S_y(T::Type{<:Number}, ::Type{U1Irrep}, ::Type{Trivial}; slave_fermion::Bool = false) t = single_site_operator(T, U1Irrep, Trivial; slave_fermion) I = sectortype(t) b = slave_fermion ? 0 : 1 @@ -489,11 +523,13 @@ const Sʸ = S_y Return the one-body spin-1/2 z-operator on the electrons. """ -function S_z(P::Type{<:Sector}=Trivial, S::Type{<:Sector}=Trivial; - slave_fermion::Bool=false) +function S_z( + P::Type{<:Sector} = Trivial, S::Type{<:Sector} = Trivial; + slave_fermion::Bool = false + ) return S_z(ComplexF64, P, S; slave_fermion) end -function S_z(T::Type{<:Number}, ::Type{Trivial}, ::Type{Trivial}; slave_fermion::Bool=false) +function S_z(T::Type{<:Number}, ::Type{Trivial}, ::Type{Trivial}; slave_fermion::Bool = false) t = single_site_operator(T, Trivial, Trivial; slave_fermion) I = sectortype(t) b = slave_fermion ? 0 : 1 @@ -501,7 +537,7 @@ function S_z(T::Type{<:Number}, ::Type{Trivial}, ::Type{Trivial}; slave_fermion: t[(I(b), dual(I(b)))][2, 2] = -0.5 return t end -function S_z(T::Type{<:Number}, ::Type{Trivial}, ::Type{U1Irrep}; slave_fermion::Bool=false) +function S_z(T::Type{<:Number}, ::Type{Trivial}, ::Type{U1Irrep}; slave_fermion::Bool = false) t = single_site_operator(T, Trivial, U1Irrep; slave_fermion) I = sectortype(t) b = slave_fermion ? 0 : 1 @@ -509,7 +545,7 @@ function S_z(T::Type{<:Number}, ::Type{Trivial}, ::Type{U1Irrep}; slave_fermion: t[(I(b, -1 // 2), dual(I(b, -1 // 2)))] .= -0.5 return t end -function S_z(T::Type{<:Number}, ::Type{U1Irrep}, ::Type{Trivial}; slave_fermion::Bool=false) +function S_z(T::Type{<:Number}, ::Type{U1Irrep}, ::Type{Trivial}; slave_fermion::Bool = false) t = single_site_operator(T, U1Irrep, Trivial; slave_fermion) I = sectortype(t) b = slave_fermion ? 0 : 1 @@ -517,7 +553,7 @@ function S_z(T::Type{<:Number}, ::Type{U1Irrep}, ::Type{Trivial}; slave_fermion: t[(I(b, 1), dual(I(b, 1)))][2, 2] = -0.5 return t end -function S_z(T::Type{<:Number}, ::Type{U1Irrep}, ::Type{U1Irrep}; slave_fermion::Bool=false) +function S_z(T::Type{<:Number}, ::Type{U1Irrep}, ::Type{U1Irrep}; slave_fermion::Bool = false) t = single_site_operator(T, U1Irrep, U1Irrep; slave_fermion) I = sectortype(t) b = slave_fermion ? 0 : 1 @@ -532,13 +568,15 @@ const Sᶻ = S_z Return the spin-plus operator. """ -function S_plus(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool=false) +function S_plus(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool = false) return S_plus(ComplexF64, P, S; slave_fermion) end -function S_plus(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; - slave_fermion::Bool=false) +function S_plus( + T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; + slave_fermion::Bool = false + ) return S_x(T, particle_symmetry, spin_symmetry; slave_fermion) + - 1im * S_y(T, particle_symmetry, spin_symmetry; slave_fermion) + 1im * S_y(T, particle_symmetry, spin_symmetry; slave_fermion) end const S⁺ = S_plus @@ -547,13 +585,15 @@ const S⁺ = S_plus Return the spin-minus operator. """ -function S_min(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool=false) +function S_min(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool = false) return S_min(ComplexF64, P, S; slave_fermion) end -function S_min(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; - slave_fermion::Bool=false) +function S_min( + T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; + slave_fermion::Bool = false + ) return S_x(T, particle_symmetry, spin_symmetry; slave_fermion) - - 1im * S_y(T, particle_symmetry, spin_symmetry; slave_fermion) + 1im * S_y(T, particle_symmetry, spin_symmetry; slave_fermion) end const S⁻ = S_min @@ -563,31 +603,31 @@ const S⁻ = S_min Return the two-body operator S⁺S⁻. The only nonzero matrix element corresponds to `|↑↓⟩ <-- |↓↑⟩`. """ -function S_plusmin(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool=false) +function S_plusmin(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool = false) return S_plusmin(ComplexF64, P, S; slave_fermion) end -function S_plusmin(T, ::Type{Trivial}, ::Type{Trivial}; slave_fermion::Bool=false) +function S_plusmin(T, ::Type{Trivial}, ::Type{Trivial}; slave_fermion::Bool = false) t = two_site_operator(T, Trivial, Trivial; slave_fermion) I = sectortype(t) b = slave_fermion ? 0 : 1 t[(I(b), I(b), dual(I(b)), dual(I(b)))][1, 2, 2, 1] = 1 return t end -function S_plusmin(T, ::Type{Trivial}, ::Type{U1Irrep}; slave_fermion::Bool=false) +function S_plusmin(T, ::Type{Trivial}, ::Type{U1Irrep}; slave_fermion::Bool = false) t = two_site_operator(T, Trivial, U1Irrep; slave_fermion) I = sectortype(t) b = slave_fermion ? 0 : 1 t[(I(b, 1 // 2), I(b, -1 // 2), dual(I(b, -1 // 2)), dual(I(b, 1 // 2)))] .= 1 return t end -function S_plusmin(T, ::Type{U1Irrep}, ::Type{Trivial}; slave_fermion::Bool=false) +function S_plusmin(T, ::Type{U1Irrep}, ::Type{Trivial}; slave_fermion::Bool = false) t = two_site_operator(T, U1Irrep, Trivial; slave_fermion) I = sectortype(t) b = slave_fermion ? 0 : 1 t[(I(b, 1), I(b, 1), dual(I(b, 1)), dual(I(b, 1)))][1, 2, 2, 1] = 1 return t end -function S_plusmin(T, ::Type{U1Irrep}, ::Type{U1Irrep}; slave_fermion::Bool=false) +function S_plusmin(T, ::Type{U1Irrep}, ::Type{U1Irrep}; slave_fermion::Bool = false) t = two_site_operator(T, U1Irrep, U1Irrep; slave_fermion) I = sectortype(t) b = slave_fermion ? 0 : 1 @@ -601,11 +641,13 @@ end Return the two-body operator S⁻S⁺. The only nonzero matrix element corresponds to `|↓↑⟩ <-- |↑↓⟩`. """ -function S_minplus(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool=false) +function S_minplus(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool = false) return S_minplus(ComplexF64, P, S; slave_fermion) end -function S_minplus(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; - slave_fermion::Bool=false) +function S_minplus( + T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; + slave_fermion::Bool = false + ) return copy(adjoint(S_plusmin(T, particle_symmetry, spin_symmetry; slave_fermion))) end @@ -614,16 +656,18 @@ end Return the spin exchange operator S⋅S. """ -function S_exchange(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool=false) +function S_exchange(P::Type{<:Sector}, S::Type{<:Sector}; slave_fermion::Bool = false) return S_exchange(ComplexF64, P, S; slave_fermion) end -function S_exchange(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; - slave_fermion::Bool=false) +function S_exchange( + T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; + slave_fermion::Bool = false + ) Sz = S_z(T, particle_symmetry, spin_symmetry; slave_fermion) - return (1 / 2) * (S_plusmin(T, particle_symmetry, spin_symmetry; slave_fermion) - + - S_minplus(T, particle_symmetry, spin_symmetry; slave_fermion)) + - Sz ⊗ Sz + return (1 / 2) * ( + S_plusmin(T, particle_symmetry, spin_symmetry; slave_fermion) + + S_minplus(T, particle_symmetry, spin_symmetry; slave_fermion) + ) + Sz ⊗ Sz end end diff --git a/src/utility.jl b/src/utility.jl index baaef27..a80fdfa 100644 --- a/src/utility.jl +++ b/src/utility.jl @@ -12,13 +12,17 @@ end contract two single-site operators inta a single-site operator. """ -function contract_onesite(L::AbstractTensorMap{<:Number,<:Any,1,2}, - R::AbstractTensorMap{<:Number,<:Any,2,1}) +function contract_onesite( + L::AbstractTensorMap{<:Number, <:Any, 1, 2}, + R::AbstractTensorMap{<:Number, <:Any, 2, 1} + ) @plansor H[-1; -2] := L[-1; 1 2] * τ[1 2; 3 4] * R[3 4; -2] return H end -function contract_onesite(L::AbstractTensorMap{<:Number,<:Any,1,1}, - R::AbstractTensorMap{<:Number,<:Any,1,1}) +function contract_onesite( + L::AbstractTensorMap{<:Number, <:Any, 1, 1}, + R::AbstractTensorMap{<:Number, <:Any, 1, 1} + ) return L * R end @@ -27,13 +31,17 @@ end contract two single-site operators into a two-site operator. """ -function contract_twosite(L::AbstractTensorMap{<:Number,<:Any,1,2}, - R::AbstractTensorMap{<:Number,<:Any,2,1}) +function contract_twosite( + L::AbstractTensorMap{<:Number, <:Any, 1, 2}, + R::AbstractTensorMap{<:Number, <:Any, 2, 1} + ) @plansor H[-1 -2; -3 -4] := L[-1; -3 1] * R[1 -2; -4] return H end -function contract_twosite(L::AbstractTensorMap{<:Any,<:Any,1,1}, - R::AbstractTensorMap{<:Any,<:Any,1,1}) +function contract_twosite( + L::AbstractTensorMap{<:Any, <:Any, 1, 1}, + R::AbstractTensorMap{<:Any, <:Any, 1, 1} + ) return L ⊗ R end @@ -42,8 +50,8 @@ end Split a two-site operator into two single-site operators with a connecting auxiliary leg. """ -function split_twosite(O::AbstractTensorMap{<:Any,<:Any,2,2}) - U, S, V, = tsvd(O, ((3, 1), (4, 2)); trunc=truncbelow(eps(real(scalartype(O))))) +function split_twosite(O::AbstractTensorMap{<:Any, <:Any, 2, 2}) + U, S, V, = tsvd(O, ((3, 1), (4, 2)); trunc = truncbelow(eps(real(scalartype(O))))) sqrtS = sqrt(S) @plansor L[p'; p a] := U[p p'; 1] * sqrtS[1; a] @plansor R[a p'; p] := sqrtS[a; 1] * V[1; p p'] diff --git a/test/bose_hubbard.jl b/test/bose_hubbard.jl index 6616129..4def832 100644 --- a/test/bose_hubbard.jl +++ b/test/bose_hubbard.jl @@ -15,16 +15,18 @@ lattice = InfiniteChain() Vspace = U1Space(0 => 8, 1 => 6, -1 => 6, 2 => 4, -2 => 4, 3 => 2, -3 => 2) -alg = VUMPS(; maxiter=25, verbosity=0) +alg = VUMPS(; maxiter = 25, verbosity = 0) # compare against higher-order analytic expansion from https://arxiv.org/pdf/1507.06426 -function exact_bose_hubbard_energy(; t=1.0, U=1.0) +function exact_bose_hubbard_energy(; t = 1.0, U = 1.0) J = t / U E = 4 * U * - (-J^2 + J^4 + 68 / 9 * J^6 - 1267 / 81 * J^8 + 44171 / 1458 * J^10 - - 4902596 / 6561 * J^12 - - 8020902135607 / 2645395200 * J^14 - 32507578587517774813 / 466647713280000 * J^16) + ( + -J^2 + J^4 + 68 / 9 * J^6 - 1267 / 81 * J^8 + 44171 / 1458 * J^10 - + 4902596 / 6561 * J^12 - + 8020902135607 / 2645395200 * J^14 - 32507578587517774813 / 466647713280000 * J^16 + ) return E end @@ -34,7 +36,7 @@ end @testset "Bose-Hubbard ground state" begin H = bose_hubbard_model(symmetry, lattice; cutoff, t, U, mu, n) ψ = InfiniteMPS([physicalspace(H, 1)], [Vspace]) - @test imag(expectation_value(ψ, H)) ≈ 0 atol = 1e-10 + @test imag(expectation_value(ψ, H)) ≈ 0 atol = 1.0e-10 ψ, envs, δ = find_groundstate(ψ, H, alg) - @test expectation_value(ψ, H, envs) ≈ exact_bose_hubbard_energy(; t, U) atol = 1e-2 + @test expectation_value(ψ, H, envs) ≈ exact_bose_hubbard_energy(; t, U) atol = 1.0e-2 end diff --git a/test/bosonoperators.jl b/test/bosonoperators.jl index 0b6ae7c..0955afa 100644 --- a/test/bosonoperators.jl +++ b/test/bosonoperators.jl @@ -7,17 +7,17 @@ cutoff = 3 elt = ComplexF64 @testset "non-symmetric bosonic operators" begin - raising = a_plus(; cutoff=cutoff) - lowering = a_min(; cutoff=cutoff) + raising = a_plus(; cutoff = cutoff) + lowering = a_min(; cutoff = cutoff) @test raising' ≈ lowering - @test a_number(; cutoff=cutoff) ≈ raising * lowering atol = 1e-4 + @test a_number(; cutoff = cutoff) ≈ raising * lowering atol = 1.0e-4 end @testset "U1-symmetric bosonic operators" begin - @test convert(Array, a_number(U1Irrep; cutoff=cutoff)) ≈ - convert(Array, a_number(; cutoff=cutoff)) - @test permute(a_plus(U1Irrep; cutoff=cutoff, side=:L)', ((2, 1), (3,))) ≈ - a_min(U1Irrep; cutoff=cutoff, side=:R) - @test permute(a_min(U1Irrep; cutoff=cutoff, side=:L)', ((2, 1), (3,))) ≈ - a_plus(U1Irrep; cutoff=cutoff, side=:R) + @test convert(Array, a_number(U1Irrep; cutoff = cutoff)) ≈ + convert(Array, a_number(; cutoff = cutoff)) + @test permute(a_plus(U1Irrep; cutoff = cutoff, side = :L)', ((2, 1), (3,))) ≈ + a_min(U1Irrep; cutoff = cutoff, side = :R) + @test permute(a_min(U1Irrep; cutoff = cutoff, side = :L)', ((2, 1), (3,))) ≈ + a_plus(U1Irrep; cutoff = cutoff, side = :R) end diff --git a/test/classical_ising.jl b/test/classical_ising.jl index 9daaaa6..4fe4ed2 100644 --- a/test/classical_ising.jl +++ b/test/classical_ising.jl @@ -8,7 +8,7 @@ using MPSKitModels beta = 0.6 Vspaces = [ComplexSpace(12), Z2Space(0 => 6, 1 => 6)] -alg = VUMPS(; tol=1e-8, maxiter=25, verbosity=1) +alg = VUMPS(; tol = 1.0e-8, maxiter = 25, verbosity = 1) """ classical_ising_free_energy(; beta) @@ -20,11 +20,12 @@ for the free energy of the 2D classical Ising Model with partition function \\mathcal{Z}(\\beta) = \\sum_{\\{s\\}} \\exp(-\\beta H(s)) \\text{ with } H(s) = - \\sum_{\\langle i, j \\rangle} s_i s_j ``` """ -function classical_ising_free_energy(; beta=log(1 + sqrt(2)) / 2) +function classical_ising_free_energy(; beta = log(1 + sqrt(2)) / 2) k = 1 / sinh(2 * beta)^2 - F = quadgk(theta -> log(cosh(2 * beta)^2 + - 1 / k * sqrt(1 + k^2 - 2 * k * cos(2 * theta))), - 0, pi)[1] + F = quadgk( + θ -> log(cosh(2 * beta)^2 + 1 / k * sqrt(1 + k^2 - 2 * k * cos(2 * θ))), + 0, pi + )[1] return -1 / beta * (log(2) / 2 + 1 / (2 * pi) * F) end @@ -37,5 +38,5 @@ end λ = expectation_value(psi, O, envs) f = -log(λ) / beta f_exact = classical_ising_free_energy(; beta) - @test abs(f - f_exact) < 1e-4 + @test abs(f - f_exact) < 1.0e-4 end diff --git a/test/fermionoperators.jl b/test/fermionoperators.jl index 212c48e..7a49268 100644 --- a/test/fermionoperators.jl +++ b/test/fermionoperators.jl @@ -10,10 +10,10 @@ using MPSKitModels: contract_twosite, contract_onesite # {cᵢ, cⱼ†} = δᵢⱼ @testset "simple fermions" begin - cc = contract_twosite(c⁻(; side=:L), c⁻(; side=:R)) - cc⁺ = contract_twosite(c⁻(; side=:L), c⁺(; side=:R)) - c⁺c = contract_twosite(c⁺(; side=:L), c⁻(; side=:R)) - c⁺c⁺ = contract_twosite(c⁺(; side=:L), c⁺(; side=:R)) + cc = contract_twosite(c⁻(; side = :L), c⁻(; side = :R)) + cc⁺ = contract_twosite(c⁻(; side = :L), c⁺(; side = :R)) + c⁺c = contract_twosite(c⁺(; side = :L), c⁻(; side = :R)) + c⁺c⁺ = contract_twosite(c⁺(; side = :L), c⁺(; side = :R)) @test cc ≈ -permute(cc, ((2, 1), (4, 3))) @test c⁺c⁺ ≈ -permute(c⁺c⁺, ((2, 1), (4, 3))) @@ -27,5 +27,5 @@ using MPSKitModels: contract_twosite, contract_onesite @test (c⁺c + cc⁺)' ≈ cc⁺ + c⁺c @test (c⁺c - cc⁺)' ≈ cc⁺ - c⁺c - @test c_number() ≈ contract_onesite(c⁺(; side=:L), c⁻(; side=:R)) + @test c_number() ≈ contract_onesite(c⁺(; side = :L), c⁻(; side = :R)) end diff --git a/test/heisenberg.jl b/test/heisenberg.jl index 7cc6c94..d20b642 100644 --- a/test/heisenberg.jl +++ b/test/heisenberg.jl @@ -1,43 +1,53 @@ using MPSKit using TensorKit -alg = VUMPS(; maxiter=25, verbosity=0) +alg = VUMPS(; maxiter = 25, verbosity = 0) E₀ = -1.401484014561 E₁ = 0.41047925 @testset "xxx" begin H = heisenberg_XXX() ψ = InfiniteMPS([ComplexSpace(3)], [ComplexSpace(48)]) - @test imag(expectation_value(ψ, H)) ≈ 0 atol = 1e-10 + @test imag(expectation_value(ψ, H)) ≈ 0 atol = 1.0e-10 ψ, envs, δ = find_groundstate(ψ, H, alg) - @test E₀ ≈ expectation_value(ψ, H, envs) atol = 1e-2 + @test E₀ ≈ expectation_value(ψ, H, envs) atol = 1.0e-2 ΔEs, qps = excitations(H, QuasiparticleAnsatz(), Float64(pi), ψ, envs) - @test E₁ ≈ first(ΔEs) atol = 1e-2 + @test E₁ ≈ first(ΔEs) atol = 1.0e-2 end @testset "xxx SU2" begin H = heisenberg_XXX(SU2Irrep) ψ = InfiniteMPS([Rep[SU₂](1 => 1)], [Rep[SU₂](1 // 2 => 5, 3 // 2 => 5, 5 // 2 => 1)]) - @test imag(expectation_value(ψ, H)) ≈ 0 atol = 1e-10 + @test imag(expectation_value(ψ, H)) ≈ 0 atol = 1.0e-10 ψ, envs, δ = find_groundstate(ψ, H, alg) - @test E₀ ≈ expectation_value(ψ, H, envs) atol = 1e-2 + @test E₀ ≈ expectation_value(ψ, H, envs) atol = 1.0e-2 - ΔEs, qps = excitations(H, QuasiparticleAnsatz(), Float64(pi), ψ, envs; - sector=SU2Irrep(1)) - @test E₁ ≈ first(ΔEs) atol = 1e-2 + ΔEs, qps = excitations( + H, QuasiparticleAnsatz(), Float64(pi), ψ, envs; + sector = SU2Irrep(1) + ) + @test E₁ ≈ first(ΔEs) atol = 1.0e-2 end @testset "xxx U1" begin H = heisenberg_XXX(U1Irrep) - ψ = InfiniteMPS([Rep[U₁](0 => 1, 1 => 1, -1 => 1)], - [Rep[U₁](1 // 2 => 10, -1 // 2 => 10, 3 // 2 => 5, -3 // 2 => 5, - 5 // 2 => 3, -5 // 2 => 3)]) - @test imag(expectation_value(ψ, H)) ≈ 0 atol = 1e-10 + ψ = InfiniteMPS( + [Rep[U₁](0 => 1, 1 => 1, -1 => 1)], + [ + Rep[U₁]( + 1 // 2 => 10, -1 // 2 => 10, 3 // 2 => 5, -3 // 2 => 5, + 5 // 2 => 3, -5 // 2 => 3 + ), + ] + ) + @test imag(expectation_value(ψ, H)) ≈ 0 atol = 1.0e-10 ψ, envs, δ = find_groundstate(ψ, H, alg) - @test E₀ ≈ expectation_value(ψ, H, envs) atol = 1e-2 + @test E₀ ≈ expectation_value(ψ, H, envs) atol = 1.0e-2 - ΔEs, qps = excitations(H, QuasiparticleAnsatz(), Float64(pi), ψ, envs; - sector=U1Irrep(1)) - @test E₁ ≈ first(ΔEs) atol = 1e-2 + ΔEs, qps = excitations( + H, QuasiparticleAnsatz(), Float64(pi), ψ, envs; + sector = U1Irrep(1) + ) + @test E₁ ≈ first(ΔEs) atol = 1.0e-2 end diff --git a/test/hubbardoperators.jl b/test/hubbardoperators.jl index 55163fa..ea67ebb 100644 --- a/test/hubbardoperators.jl +++ b/test/hubbardoperators.jl @@ -3,38 +3,40 @@ using TensorKit using MPSKitModels.HubbardOperators using LinearAlgebra: eigvals -implemented_symmetries = [(Trivial, Trivial), (Trivial, U1Irrep), (Trivial, SU2Irrep), - (U1Irrep, Trivial), (U1Irrep, U1Irrep), (U1Irrep, SU2Irrep)] +implemented_symmetries = [ + (Trivial, Trivial), (Trivial, U1Irrep), (Trivial, SU2Irrep), + (U1Irrep, Trivial), (U1Irrep, U1Irrep), (U1Irrep, SU2Irrep), +] @testset "basic properties" begin for particle_symmetry in (Trivial, U1Irrep, SU2Irrep), - spin_symmetry in (Trivial, U1Irrep, SU2Irrep) + spin_symmetry in (Trivial, U1Irrep, SU2Irrep) if (particle_symmetry, spin_symmetry) in implemented_symmetries # test hermiticity @test e_plusmin(particle_symmetry, spin_symmetry)' ≈ - e_minplus(particle_symmetry, spin_symmetry) + e_minplus(particle_symmetry, spin_symmetry) if spin_symmetry !== SU2Irrep @test e_plusmin_down(particle_symmetry, spin_symmetry)' ≈ - e_minplus_down(particle_symmetry, spin_symmetry) + e_minplus_down(particle_symmetry, spin_symmetry) @test e_plusmin_up(particle_symmetry, spin_symmetry)' ≈ - e_minplus_up(particle_symmetry, spin_symmetry) + e_minplus_up(particle_symmetry, spin_symmetry) @test e_plusmin_down(particle_symmetry, spin_symmetry)' ≈ - e_minplus_down(particle_symmetry, spin_symmetry) + e_minplus_down(particle_symmetry, spin_symmetry) @test e_plusmin_up(particle_symmetry, spin_symmetry)' ≈ - e_minplus_up(particle_symmetry, spin_symmetry) + e_minplus_up(particle_symmetry, spin_symmetry) end # test number operator if spin_symmetry !== SU2Irrep @test e_number(particle_symmetry, spin_symmetry) ≈ - e_number_up(particle_symmetry, spin_symmetry) + - e_number_down(particle_symmetry, spin_symmetry) + e_number_up(particle_symmetry, spin_symmetry) + + e_number_down(particle_symmetry, spin_symmetry) @test e_number_updown(particle_symmetry, spin_symmetry) ≈ - e_number_up(particle_symmetry, spin_symmetry) * - e_number_down(particle_symmetry, spin_symmetry) ≈ - e_number_down(particle_symmetry, spin_symmetry) * - e_number_up(particle_symmetry, spin_symmetry) + e_number_up(particle_symmetry, spin_symmetry) * + e_number_down(particle_symmetry, spin_symmetry) ≈ + e_number_down(particle_symmetry, spin_symmetry) * + e_number_up(particle_symmetry, spin_symmetry) end else @test_broken e_plusmin(particle_symmetry, spin_symmetry) @@ -44,20 +46,22 @@ implemented_symmetries = [(Trivial, Trivial), (Trivial, U1Irrep), (Trivial, SU2I end function hubbard_hamiltonian(particle_symmetry, spin_symmetry; t, U, mu, L) - hopping = t * (e_plusmin(particle_symmetry, spin_symmetry) + - e_minplus(particle_symmetry, spin_symmetry)) + hopping = t * ( + e_plusmin(particle_symmetry, spin_symmetry) + + e_minplus(particle_symmetry, spin_symmetry) + ) interaction = U * e_number_updown(particle_symmetry, spin_symmetry) chemical_potential = mu * e_number(particle_symmetry, spin_symmetry) I = id(hubbard_space(particle_symmetry, spin_symmetry)) H = sum(1:(L - 1)) do i - return reduce(⊗, insert!(collect(Any, fill(I, L - 2)), i, hopping)) - end + + return reduce(⊗, insert!(collect(Any, fill(I, L - 2)), i, hopping)) + end + sum(1:L) do i - return reduce(⊗, insert!(collect(Any, fill(I, L - 1)), i, interaction)) - end + + return reduce(⊗, insert!(collect(Any, fill(I, L - 1)), i, interaction)) + end + sum(1:L) do i - return reduce(⊗, insert!(collect(Any, fill(I, L - 1)), i, chemical_potential)) - end + return reduce(⊗, insert!(collect(Any, fill(I, L - 1)), i, chemical_potential)) + end return H end diff --git a/test/lattices.jl b/test/lattices.jl index 036111f..be54920 100644 --- a/test/lattices.jl +++ b/test/lattices.jl @@ -220,7 +220,7 @@ end NN = nearest_neighbours(lattice) @test all(zip(NN, base_NN)) do (x, y) return linearize_index.(Tuple(x)) == - linearize_index.(Tuple(y)) + linearize_index.(Tuple(y)) end pattern(i) = length(lattice) - i + 1 @@ -238,6 +238,6 @@ end NN = nearest_neighbours(lattice) @test all(zip(NN, base_NN)) do (x, y) return linearize_index.(Tuple(x)) == - pattern.(linearize_index.(Tuple(y))) + pattern.(linearize_index.(Tuple(y))) end end diff --git a/test/mpoham.jl b/test/mpoham.jl index bb4884e..48b52bd 100644 --- a/test/mpoham.jl +++ b/test/mpoham.jl @@ -3,7 +3,7 @@ using MPSKitModels, MPSKit, TensorKit lattice = InfiniteChain(1) H1 = @mpoham begin sum(nearest_neighbours(lattice)) do (i, j) - return (σˣˣ() + σʸʸ()){i,j} + return (σˣˣ() + σʸʸ()){i, j} end end @@ -17,16 +17,16 @@ end ZZ = S_zz() lattice = InfiniteCylinder(5) - H = @mpoham sum(ZZ{i,j} for (i, j) in nearest_neighbours(lattice)) + H = @mpoham sum(ZZ{i, j} for (i, j) in nearest_neighbours(lattice)) @test length(H) == length(lattice) end @testset "deduce_pspaces" begin # Not fully defining the pspaces should still work lattice = FiniteChain(5) - H = @mpoham S_zz(){lattice[2],lattice[3]} + H = @mpoham S_zz(){lattice[2], lattice[3]} @test unique(MPSKit.physicalspace(H))[1] == ComplexSpace(2) - @test_throws Exception @mpoham σˣ(){lattice[1]} + σˣ(; spin=3 // 2){lattice[2]} + @test_throws Exception @mpoham σˣ(){lattice[1]} + σˣ(; spin = 3 // 2){lattice[2]} end diff --git a/test/quantum_potts.jl b/test/quantum_potts.jl index 07bef37..cd36d6a 100644 --- a/test/quantum_potts.jl +++ b/test/quantum_potts.jl @@ -9,14 +9,16 @@ using MPSKitModels qs = [3, 4, 5] symmetries = [Trivial, ZNIrrep] Vspaces = [ComplexSpace(12), Z2Space(0 => 6, 1 => 6)] -alg = VUMPS(; maxiter=25, verbosity=0) +alg = VUMPS(; maxiter = 25, verbosity = 0) # https://iopscience.iop.org/article/10.1088/0305-4470/14/11/020/meta -function quantum_potts_energy(Q::Int; maxiter::Int=1000) +function quantum_potts_energy(Q::Int; maxiter::Int = 1000) Q == 3 && return -(4 / 3 + 2sqrt(3) / π) Q == 4 && return 2 - 8 * log(2) - summation = sum((-1)^n / (sqrt(Q) / 2 - cosh((2 * n + 1) * acosh(sqrt(Q) / 2))) - for n in 1:maxiter) + summation = sum( + (-1)^n / (sqrt(Q) / 2 - cosh((2 * n + 1) * acosh(sqrt(Q) / 2))) + for n in 1:maxiter + ) limit = 2 - Q - sqrt(Q) * (Q - 4) * summation return limit end @@ -28,13 +30,10 @@ _virtualspace(::Type{ZNIrrep}, q::Int, χ::Int) = Rep[ℤ{q}](i => χ for i in 0 ## Test -@testset "$q-state Potts with $(_sectortype(symmetry, q))) symmetry" for q in qs, - symmetry in - symmetries - +@testset "$q-state Potts with $(_sectortype(symmetry, q))) symmetry" for q in qs, symmetry in symmetries H = quantum_potts(_sectortype(symmetry, q); q) ψ = InfiniteMPS(physicalspace(H, 1), _virtualspace(symmetry, q, χ)) - @test imag(expectation_value(ψ, H)) ≈ 0 atol = 1e-10 + @test imag(expectation_value(ψ, H)) ≈ 0 atol = 1.0e-10 ψ, envs, δ = find_groundstate(ψ, H, alg) - @test quantum_potts_energy(q) ≈ expectation_value(ψ, H, envs) atol = 1e-2 + @test quantum_potts_energy(q) ≈ expectation_value(ψ, H, envs) atol = 1.0e-2 end diff --git a/test/sixvertex.jl b/test/sixvertex.jl index 4385ee6..8bd88e0 100644 --- a/test/sixvertex.jl +++ b/test/sixvertex.jl @@ -3,35 +3,39 @@ using TensorKit using Test F₀ = (4 / 3)^(3 / 2) -alg = VUMPS(; maxiter=25, verbosity=0) +alg = VUMPS(; maxiter = 25, verbosity = 0) @testset "ℤ₁" begin mpo = @inferred sixvertex() ψ = InfiniteMPS(ℂ^2, ℂ^20) ψ, _ = leading_boundary(ψ, mpo, alg) F = prod(expectation_value(ψ, mpo)) - @test imag(F) < 1e-3 - @test real(F) ≈ F₀ atol = 1e-2 + @test imag(F) < 1.0e-3 + @test real(F) ≈ F₀ atol = 1.0e-2 end @testset "U₁" begin mpo = @inferred sixvertex(ComplexF64, U1Irrep) mpo2 = repeat(mpo, 2, 2) - vspaces = [U1Space(0 => 20, 1 => 10, -1 => 10, 2 => 5, -2 => 5), - U1Space(1 // 2 => 15, -1 // 2 => 15, 3 // 2 => 5, -3 // 2 => 5)] + vspaces = [ + U1Space(0 => 20, 1 => 10, -1 => 10, 2 => 5, -2 => 5), + U1Space(1 // 2 => 15, -1 // 2 => 15, 3 // 2 => 5, -3 // 2 => 5), + ] ψ = MultilineMPS(repeat(physicalspace(mpo), 2, 2), [vspaces circshift(vspaces, 1)]) ψ, _ = leading_boundary(ψ, mpo2, alg) F = prod(expectation_value(ψ, mpo2)) - @test abs(F)^(1 / 4) ≈ F₀ atol = 1e-2 + @test abs(F)^(1 / 4) ≈ F₀ atol = 1.0e-2 end @testset "CU₁" begin mpo = @inferred sixvertex(ComplexF64, CU1Irrep) mpo2 = repeat(mpo, 2, 2) - vspaces = [CU1Space((0, 0) => 10, (0, 1) => 10, (1, 2) => 5, (2, 2) => 5), - CU1Space((1 // 2, 2) => 15, (3 // 2, 2) => 5)] + vspaces = [ + CU1Space((0, 0) => 10, (0, 1) => 10, (1, 2) => 5, (2, 2) => 5), + CU1Space((1 // 2, 2) => 15, (3 // 2, 2) => 5), + ] ψ = MultilineMPS(repeat(physicalspace(mpo), 2, 2), [vspaces circshift(vspaces, 1)]) ψ, _ = leading_boundary(ψ, mpo2, alg) F = prod(expectation_value(ψ, mpo2)) - @test abs(F)^(1 / 4) ≈ F₀ atol = 1e-2 + @test abs(F)^(1 / 4) ≈ F₀ atol = 1.0e-2 end diff --git a/test/spinoperators.jl b/test/spinoperators.jl index 7011c66..2a94608 100644 --- a/test/spinoperators.jl +++ b/test/spinoperators.jl @@ -13,17 +13,17 @@ end @testset "non-symmetric spin $(Int(2S))/2 operators" for S in (1 // 2):(1 // 2):4 # inferrability - X = @inferred S_x(; spin=S) - Y = @inferred S_y(; spin=S) - Z = @inferred S_z(; spin=S) - S⁺ = @inferred S_plus(; spin=S) - S⁻ = @inferred S_min(; spin=S) - S⁺⁻ = @inferred S_plusmin(; spin=S) - S⁻⁺ = @inferred S_minplus(; spin=S) - XX = @inferred S_xx(; spin=S) - YY = @inferred S_yy(; spin=S) - ZZ = @inferred S_zz(; spin=S) - SS = @inferred S_exchange(; spin=S) + X = @inferred S_x(; spin = S) + Y = @inferred S_y(; spin = S) + Z = @inferred S_z(; spin = S) + S⁺ = @inferred S_plus(; spin = S) + S⁻ = @inferred S_min(; spin = S) + S⁺⁻ = @inferred S_plusmin(; spin = S) + S⁻⁺ = @inferred S_minplus(; spin = S) + XX = @inferred S_xx(; spin = S) + YY = @inferred S_yy(; spin = S) + ZZ = @inferred S_zz(; spin = S) + SS = @inferred S_exchange(; spin = S) Svec = [X Y Z] # operators should be hermitian @@ -37,7 +37,7 @@ end # commutation relations for i in 1:3, j in 1:3 @test Svec[i] * Svec[j] - Svec[j] * Svec[i] ≈ - sum(im * ε[i, j, k] * Svec[k] for k in 1:3) + sum(im * ε[i, j, k] * Svec[k] for k in 1:3) end # definition of +- @@ -62,22 +62,22 @@ end @test H * convert(Array, S_x()) * H' ≈ convert(Array, S_x(Z2Irrep)) for S in (S_y, S_z, S_plus, S_min) array1 = H * convert(Array, S()) * H' - arrayL = reshape(sum(convert(Array, S(Z2Irrep; side=:L)); dims=3), 2, 2) - arrayR = reshape(sum(convert(Array, S(Z2Irrep; side=:R)); dims=1), 2, 2) + arrayL = reshape(sum(convert(Array, S(Z2Irrep; side = :L)); dims = 3), 2, 2) + arrayR = reshape(sum(convert(Array, S(Z2Irrep; side = :R)); dims = 1), 2, 2) @test array1 ≈ arrayL @test array1 ≈ arrayR end # inferrability X = @inferred S_x(Z2Irrep) - YL = @constinferred S_y(Z2Irrep; side=:L) - YR = @constinferred S_y(Z2Irrep; side=:R) - ZL = @constinferred S_z(Z2Irrep; side=:L) - ZR = @constinferred S_z(Z2Irrep; side=:R) - S⁺L = @constinferred S_plus(Z2Irrep; side=:L) - S⁺R = @constinferred S_plus(Z2Irrep; side=:R) - S⁻L = @constinferred S_min(Z2Irrep; side=:L) - S⁻R = @constinferred S_min(Z2Irrep; side=:R) + YL = @constinferred S_y(Z2Irrep; side = :L) + YR = @constinferred S_y(Z2Irrep; side = :R) + ZL = @constinferred S_z(Z2Irrep; side = :L) + ZR = @constinferred S_z(Z2Irrep; side = :R) + S⁺L = @constinferred S_plus(Z2Irrep; side = :L) + S⁺R = @constinferred S_plus(Z2Irrep; side = :R) + S⁻L = @constinferred S_min(Z2Irrep; side = :L) + S⁻R = @constinferred S_min(Z2Irrep; side = :R) S⁺⁻ = @inferred S_plusmin(Z2Irrep) S⁻⁺ = @inferred S_minplus(Z2Irrep) XX = @inferred S_xx(Z2Irrep) @@ -91,51 +91,45 @@ end @test permute(S⁻L', ((2, 1), (3,))) ≈ S⁺R # composite operators - @test (S⁺⁻ + S⁻⁺) / 2 ≈ XX + YY rtol = 1e-3 + @test (S⁺⁻ + S⁻⁺) / 2 ≈ XX + YY rtol = 1.0e-3 end @testset "U1-symmetric spin $(Int(2spin))/2 operators" for spin in (1 // 2):(1 // 2):4 # array conversion N = Int(2spin + 1) - p = sortperm(reverse((-spin):spin); by=x -> abs(x - 0.1)) # sort as 0, 1, -1, 2, -2, ... + p = sortperm(reverse((-spin):spin); by = x -> abs(x - 0.1)) # sort as 0, 1, -1, 2, -2, ... H = one(zeros(N, N))[p, :] - @test H * convert(Array, S_z(; spin=spin)) * H' ≈ - convert(Array, S_z(U1Irrep; spin=spin)) + @test H * convert(Array, S_z(; spin = spin)) * H' ≈ + convert(Array, S_z(U1Irrep; spin = spin)) for S in (S_x, S_y, S_plus, S_min) - array1 = convert(Array, S(; spin=spin)) - arrayL = H' * - reshape(sum(convert(Array, S(U1Irrep; side=:L, spin=spin)); dims=3), N, - N) * - H - arrayR = H' * - reshape(sum(convert(Array, S(U1Irrep; side=:R, spin=spin)); dims=1), N, - N) * - H + array1 = convert(Array, S(; spin = spin)) + arrayL = H' * reshape(sum(convert(Array, S(U1Irrep; side = :L, spin = spin)); dims = 3), N, N) * H + arrayR = H' * reshape(sum(convert(Array, S(U1Irrep; side = :R, spin = spin)); dims = 1), N, N) * H @test array1 ≈ arrayL @test array1 ≈ arrayR end # # hermiticity - @test S_z(U1Irrep; spin=spin)' ≈ S_z(U1Irrep; spin=spin) - @test permute(S_x(U1Irrep; spin=spin, side=:L)', ((2, 1), (3,))) ≈ - S_x(U1Irrep; spin=spin, side=:R) - @test permute(S_y(U1Irrep; spin=spin, side=:L)', ((2, 1), (3,))) ≈ - S_y(U1Irrep; spin=spin, side=:R) - @test permute(S_plus(U1Irrep; spin=spin, side=:L)', ((2, 1), (3,))) ≈ - S_min(U1Irrep; spin=spin, side=:R) - @test permute(S_min(U1Irrep; spin=spin, side=:L)', ((2, 1), (3,))) ≈ - S_plus(U1Irrep; spin=spin, side=:R) + @test S_z(U1Irrep; spin = spin)' ≈ S_z(U1Irrep; spin = spin) + @test permute(S_x(U1Irrep; spin = spin, side = :L)', ((2, 1), (3,))) ≈ + S_x(U1Irrep; spin = spin, side = :R) + @test permute(S_y(U1Irrep; spin = spin, side = :L)', ((2, 1), (3,))) ≈ + S_y(U1Irrep; spin = spin, side = :R) + @test permute(S_plus(U1Irrep; spin = spin, side = :L)', ((2, 1), (3,))) ≈ + S_min(U1Irrep; spin = spin, side = :R) + @test permute(S_min(U1Irrep; spin = spin, side = :L)', ((2, 1), (3,))) ≈ + S_plus(U1Irrep; spin = spin, side = :R) # # composite operators - @test (S_plusmin(U1Irrep; spin=spin) + S_minplus(U1Irrep; spin=spin)) / 2 ≈ - S_xx(U1Irrep; spin=spin) + S_yy(U1Irrep; spin=spin) rtol = 1e-3 - @test S_exchange(U1Irrep; spin=spin) ≈ - S_xx(U1Irrep; spin=spin) + S_yy(U1Irrep; spin=spin) + S_zz(U1Irrep; spin=spin) rtol = 1e-3 + @test (S_plusmin(U1Irrep; spin = spin) + S_minplus(U1Irrep; spin = spin)) / 2 ≈ + S_xx(U1Irrep; spin = spin) + S_yy(U1Irrep; spin = spin) rtol = 1.0e-3 + @test S_exchange(U1Irrep; spin = spin) ≈ + S_xx(U1Irrep; spin = spin) + S_yy(U1Irrep; spin = spin) + S_zz(U1Irrep; spin = spin) rtol = 1.0e-3 end @testset "Real eltype" begin for operator in (S_x, S_z, S_xx, S_zz, S_plusmin, S_minplus, S_exchange), - symmetry in (Trivial, Z2Irrep, U1Irrep) + symmetry in (Trivial, Z2Irrep, U1Irrep) @test real(operator(ComplexF64, symmetry)) ≈ operator(Float64, symmetry) end @@ -145,9 +139,9 @@ end # potts_ZZ test? @testset "non-symmetric Q-state potts operators" for Q in 3:5 # inferrability - X = @inferred potts_X(; q=Q) - Z = @inferred potts_Z(; q=Q) - ZZ = @inferred potts_ZZ(; q=Q) + X = @inferred potts_X(; q = Q) + Z = @inferred potts_Z(; q = Q) + ZZ = @inferred potts_ZZ(; q = Q) # clock properties @test convert(Array, X^Q) ≈ I @@ -169,11 +163,11 @@ end @testset "Z_Q-symmetric Q-state Potts operators" for Q in 3:5 # array conversion _, _, W = weyl_heisenberg_matrices(Q, ComplexF64) - @test W * convert(Array, potts_X(; q=Q)) * W' ≈ convert(Array, potts_X(ZNIrrep{Q}; q=Q)) + @test W * convert(Array, potts_X(; q = Q)) * W' ≈ convert(Array, potts_X(ZNIrrep{Q}; q = Q)) # inferrability - X = @inferred potts_X(ZNIrrep{Q}; q=Q) - ZZ = @inferred potts_ZZ(ZNIrrep{Q}; q=Q) + X = @inferred potts_X(ZNIrrep{Q}; q = Q) + ZZ = @inferred potts_ZZ(ZNIrrep{Q}; q = Q) # unitarity @test X * X' ≈ X' * X diff --git a/test/tfim.jl b/test/tfim.jl index 88bc1bf..77b5f38 100644 --- a/test/tfim.jl +++ b/test/tfim.jl @@ -4,30 +4,30 @@ using TensorKit using Test E₀ = -1.273239 -alg = VUMPS(; maxiter=25, verbosity=0) +alg = VUMPS(; maxiter = 25, verbosity = 0) @testset "no symmetry" begin H = transverse_field_ising() ψ₀ = InfiniteMPS([ComplexSpace(2)], [ComplexSpace(16)]) - @test imag(expectation_value(ψ₀, H)) ≈ 0 atol = 1e-10 + @test imag(expectation_value(ψ₀, H)) ≈ 0 atol = 1.0e-10 ψ, envs, δ = find_groundstate(ψ₀, H, alg) - @test E₀ ≈ expectation_value(ψ, H, envs) atol = 1e-5 + @test E₀ ≈ expectation_value(ψ, H, envs) atol = 1.0e-5 end @testset "Z2 symmetry" begin H = transverse_field_ising(Z2Irrep) ψ₀ = InfiniteMPS([Rep[ℤ₂](0 => 1, 1 => 1)], [Rep[ℤ₂](0 => 8, 1 => 8)]) - @test imag(expectation_value(ψ₀, H)) ≈ 0 atol = 1e-10 + @test imag(expectation_value(ψ₀, H)) ≈ 0 atol = 1.0e-10 ψ, envs, δ = find_groundstate(ψ₀, H, alg) - @test E₀ ≈ expectation_value(ψ, H, envs) atol = 1e-5 + @test E₀ ≈ expectation_value(ψ, H, envs) atol = 1.0e-5 end @testset "fZ2 symmetry" begin H = transverse_field_ising(fℤ₂) ψ₀ = InfiniteMPS([Vect[fℤ₂](0 => 1, 1 => 1)], [Vect[fℤ₂](0 => 8, 1 => 8)]) - @test imag(expectation_value(ψ₀, H)) ≈ 0 atol = 1e-10 + @test imag(expectation_value(ψ₀, H)) ≈ 0 atol = 1.0e-10 ψ, envs, δ = find_groundstate(ψ₀, H, alg) - @test E₀ ≈ expectation_value(ψ, H, envs) atol = 1e-5 + @test E₀ ≈ expectation_value(ψ, H, envs) atol = 1.0e-5 end @testset "illegal symmetry" begin diff --git a/test/tjoperators.jl b/test/tjoperators.jl index 9265a69..a092e1e 100644 --- a/test/tjoperators.jl +++ b/test/tjoperators.jl @@ -3,41 +3,43 @@ using TensorKit using LinearAlgebra: eigvals import MPSKitModels: TJOperators as tJ -implemented_symmetries = [(Trivial, Trivial), - (Trivial, U1Irrep), - (U1Irrep, Trivial), - (U1Irrep, U1Irrep)] +implemented_symmetries = [ + (Trivial, Trivial), + (Trivial, U1Irrep), + (U1Irrep, Trivial), + (U1Irrep, U1Irrep), +] @testset "basic properties" begin for slave_fermion in (false, true), - particle_symmetry in (Trivial, U1Irrep), - spin_symmetry in (Trivial, U1Irrep, SU2Irrep) + particle_symmetry in (Trivial, U1Irrep), + spin_symmetry in (Trivial, U1Irrep, SU2Irrep) if (particle_symmetry, spin_symmetry) in implemented_symmetries # test hermiticity @test tJ.e_plusmin(particle_symmetry, spin_symmetry; slave_fermion)' ≈ - tJ.e_minplus(particle_symmetry, spin_symmetry; slave_fermion) + tJ.e_minplus(particle_symmetry, spin_symmetry; slave_fermion) if spin_symmetry !== SU2Irrep @test tJ.e_plusmin_down(particle_symmetry, spin_symmetry; slave_fermion)' ≈ - tJ.e_minplus_down(particle_symmetry, spin_symmetry; slave_fermion) + tJ.e_minplus_down(particle_symmetry, spin_symmetry; slave_fermion) @test tJ.e_plusmin_up(particle_symmetry, spin_symmetry; slave_fermion)' ≈ - tJ.e_minplus_up(particle_symmetry, spin_symmetry; slave_fermion) + tJ.e_minplus_up(particle_symmetry, spin_symmetry; slave_fermion) @test tJ.e_plusmin_down(particle_symmetry, spin_symmetry; slave_fermion)' ≈ - tJ.e_minplus_down(particle_symmetry, spin_symmetry; slave_fermion) + tJ.e_minplus_down(particle_symmetry, spin_symmetry; slave_fermion) @test tJ.e_plusmin_up(particle_symmetry, spin_symmetry; slave_fermion)' ≈ - tJ.e_minplus_up(particle_symmetry, spin_symmetry; slave_fermion) + tJ.e_minplus_up(particle_symmetry, spin_symmetry; slave_fermion) end # test number operator if spin_symmetry !== SU2Irrep pspace = tJ.tj_space(particle_symmetry, spin_symmetry; slave_fermion) @test tJ.e_number(particle_symmetry, spin_symmetry; slave_fermion) ≈ - tJ.e_number_up(particle_symmetry, spin_symmetry; slave_fermion) + - tJ.e_number_down(particle_symmetry, spin_symmetry; slave_fermion) + tJ.e_number_up(particle_symmetry, spin_symmetry; slave_fermion) + + tJ.e_number_down(particle_symmetry, spin_symmetry; slave_fermion) @test zeros(pspace, pspace) ≈ - tJ.e_number_up(particle_symmetry, spin_symmetry; slave_fermion) * - tJ.e_number_down(particle_symmetry, spin_symmetry; slave_fermion) ≈ - tJ.e_number_down(particle_symmetry, spin_symmetry; slave_fermion) * - tJ.e_number_up(particle_symmetry, spin_symmetry; slave_fermion) + tJ.e_number_up(particle_symmetry, spin_symmetry; slave_fermion) * + tJ.e_number_down(particle_symmetry, spin_symmetry; slave_fermion) ≈ + tJ.e_number_down(particle_symmetry, spin_symmetry; slave_fermion) * + tJ.e_number_up(particle_symmetry, spin_symmetry; slave_fermion) end # test spin operator @@ -47,9 +49,11 @@ implemented_symmetries = [(Trivial, Trivial), ε[mod1(i, 3), mod1(i + 1, 3), mod1(i + 2, 3)] = 1 ε[mod1(i, 3), mod1(i - 1, 3), mod1(i - 2, 3)] = -1 end - Svec = [tJ.S_x(particle_symmetry, spin_symmetry; slave_fermion), - tJ.S_y(particle_symmetry, spin_symmetry; slave_fermion), - tJ.S_z(particle_symmetry, spin_symmetry; slave_fermion)] + Svec = [ + tJ.S_x(particle_symmetry, spin_symmetry; slave_fermion), + tJ.S_y(particle_symmetry, spin_symmetry; slave_fermion), + tJ.S_z(particle_symmetry, spin_symmetry; slave_fermion), + ] # Hermiticity for s in Svec @test s' ≈ s @@ -59,12 +63,12 @@ implemented_symmetries = [(Trivial, Trivial), @test sum(tr(Svec[i]^2) for i in 1:3) / (2S + 1) ≈ S * (S + 1) # test S_plus and S_min @test tJ.S_plusmin(particle_symmetry, spin_symmetry; slave_fermion) ≈ - tJ.S_plus(particle_symmetry, spin_symmetry; slave_fermion) ⊗ - tJ.S_min(particle_symmetry, spin_symmetry; slave_fermion) + tJ.S_plus(particle_symmetry, spin_symmetry; slave_fermion) ⊗ + tJ.S_min(particle_symmetry, spin_symmetry; slave_fermion) # commutation relations for i in 1:3, j in 1:3 @test Svec[i] * Svec[j] - Svec[j] * Svec[i] ≈ - sum(im * ε[i, j, k] * Svec[k] for k in 1:3) + sum(im * ε[i, j, k] * Svec[k] for k in 1:3) end end else @@ -76,18 +80,20 @@ end function tjhamiltonian(particle_symmetry, spin_symmetry; t, J, mu, L, slave_fermion) num = tJ.e_number(particle_symmetry, spin_symmetry; slave_fermion) - hop_heis = (-t) * (tJ.e_plusmin(particle_symmetry, spin_symmetry; slave_fermion) + - tJ.e_minplus(particle_symmetry, spin_symmetry; slave_fermion)) + - J * - (tJ.S_exchange(particle_symmetry, spin_symmetry; slave_fermion) - - (1 / 4) * (num ⊗ num)) + hop_heis = (-t) * ( + tJ.e_plusmin(particle_symmetry, spin_symmetry; slave_fermion) + + tJ.e_minplus(particle_symmetry, spin_symmetry; slave_fermion) + ) + J * ( + tJ.S_exchange(particle_symmetry, spin_symmetry; slave_fermion) - + (1 / 4) * (num ⊗ num) + ) chemical_potential = (-mu) * num I = id(tJ.tj_space(particle_symmetry, spin_symmetry; slave_fermion)) H = sum(1:(L - 1)) do i return reduce(⊗, insert!(collect(Any, fill(I, L - 2)), i, hop_heis)) end + sum(1:L) do i - return reduce(⊗, insert!(collect(Any, fill(I, L - 1)), i, chemical_potential)) - end + return reduce(⊗, insert!(collect(Any, fill(I, L - 1)), i, chemical_potential)) + end return H end @@ -108,8 +114,9 @@ end if (particle_symmetry, spin_symmetry) == (Trivial, Trivial) continue end - H_symm = tjhamiltonian(particle_symmetry, spin_symmetry; t, J, mu, L, - slave_fermion) + H_symm = tjhamiltonian( + particle_symmetry, spin_symmetry; t, J, mu, L, slave_fermion + ) vals_symm = mapreduce(vcat, eigvals(H_symm)) do (c, v) return repeat(real.(v), dim(c)) end