From f87847fd7b1eda4ff7a951d96162698dd0ea6dd1 Mon Sep 17 00:00:00 2001 From: leburgel Date: Wed, 23 Aug 2023 11:52:54 +0200 Subject: [PATCH 1/9] Update for new TensorKit version --- examples/heisenberg.jl | 2 +- examples/test_vumps.jl | 2 +- src/operators/periodicpepo.jl | 4 ++-- src/states/abstractpeps.jl | 4 ++-- src/states/infinitepeps.jl | 4 ++-- 5 files changed, 8 insertions(+), 8 deletions(-) diff --git a/examples/heisenberg.jl b/examples/heisenberg.jl index fdb9cb32..2c9fa560 100644 --- a/examples/heisenberg.jl +++ b/examples/heisenberg.jl @@ -1,4 +1,4 @@ -using Revise, PEPSKit, TensorKit, TensorKitAD, Zygote, MPSKit +using Revise, PEPSKit, TensorKit, Zygote, MPSKit using MPSKitModels, LinearAlgebra, OptimKit using PEPSKit: NORTH,SOUTH,WEST,EAST,NORTHWEST,NORTHEAST,SOUTHEAST,SOUTHWEST,@diffset using JLD2,ChainRulesCore diff --git a/examples/test_vumps.jl b/examples/test_vumps.jl index 6bafb25c..fb92a213 100644 --- a/examples/test_vumps.jl +++ b/examples/test_vumps.jl @@ -1,4 +1,4 @@ -using Revise, PEPSKit, TensorKit, TensorKitAD, Zygote, MPSKit +using Revise, PEPSKit, TensorKit, Zygote, MPSKit p = InfinitePEPS(fill(ℂ^2,1,1),fill(ℂ^2,1,1)); diff --git a/src/operators/periodicpepo.jl b/src/operators/periodicpepo.jl index b3eb613b..153946c3 100644 --- a/src/operators/periodicpepo.jl +++ b/src/operators/periodicpepo.jl @@ -41,7 +41,7 @@ function PeriodicPEPO( Pspaces::AbstractArray{S,2}, Nspaces::AbstractArray{S,2}, Espaces::AbstractArray{S,2}=Nspaces -) where {S<:EuclideanSpace} +) where {S<:ElementarySpace} size(Pspaces) == size(Nspaces) == size(Espaces) || throw(ArgumentError("Input spaces should have equal sizes.")) @@ -61,7 +61,7 @@ end Allow users to pass in single space. """ -function PeriodicPEPO(Pspace::S, Nspace::S, Espace::S=Nspace) where {S<:EuclideanSpace} +function PeriodicPEPO(Pspace::S, Nspace::S, Espace::S=Nspace) where {S<:ElementarySpace} Pspaces = Array{S,2}(undef, (1, 1)) Pspaces[1, 1] = Pspace Nspaces = Array{S,2}(undef, (1, 1)) diff --git a/src/states/abstractpeps.jl b/src/states/abstractpeps.jl index 0d8f98d5..80f58beb 100644 --- a/src/states/abstractpeps.jl +++ b/src/states/abstractpeps.jl @@ -10,7 +10,7 @@ Default type for PEPS tensors with a single physical index, and 4 virtual indices, conventionally ordered as: T : P ← N ⊗ E ⊗ S ⊗ W. """ -const PEPSTensor{S} = AbstractTensorMap{S, 1, 4} where {S<:EuclideanSpace} +const PEPSTensor{S} = AbstractTensorMap{S, 1, 4} where {S<:ElementarySpace} """ @@ -19,7 +19,7 @@ const PEPSTensor{S} = AbstractTensorMap{S, 1, 4} where {S<:EuclideanSpace} Default type for PEPO tensors with a single incoming and outgoing physical index, and 4 virtual indices, conventionally ordered as: O : P ⊗ P' ← N ⊗ E ⊗ S ⊗ W. """ -const PEPOTensor{S} = AbstractTensorMap{S, 2, 4} where {S<:EuclideanSpace} +const PEPOTensor{S} = AbstractTensorMap{S, 2, 4} where {S<:ElementarySpace} ########################## ## Abstract state types ## diff --git a/src/states/infinitepeps.jl b/src/states/infinitepeps.jl index be31ac1f..49e50409 100644 --- a/src/states/infinitepeps.jl +++ b/src/states/infinitepeps.jl @@ -31,7 +31,7 @@ function InfinitePEPS( Pspaces::AbstractArray{S,2}, Nspaces::AbstractArray{S,2}, Espaces::AbstractArray{S,2}=Nspaces -) where {S<:EuclideanSpace} +) where {S<:ElementarySpace} size(Pspaces) == size(Nspaces) == size(Espaces) || throw(ArgumentError("Input spaces should have equal sizes.")) @@ -50,7 +50,7 @@ end Allow users to pass in single space. """ -function InfinitePEPS(Pspace::S, Nspace::S, Espace::S=Nspace) where {S<:EuclideanSpace} +function InfinitePEPS(Pspace::S, Nspace::S, Espace::S=Nspace) where {S<:ElementarySpace} Pspaces = Array{S,2}(undef, (1, 1)) Pspaces[1, 1] = Pspace Nspaces = Array{S,2}(undef, (1, 1)) From 7fb9ead2353ad6a7d5de87d634cd0d5131d2edda Mon Sep 17 00:00:00 2001 From: leburgel Date: Tue, 5 Sep 2023 11:04:52 +0200 Subject: [PATCH 2/9] Ready for progress --- examples/heisenberg.jl | 20 +++++++++----------- src/algorithms/ctmrg.jl | 12 ++++++------ src/environments/ctmrgenv.jl | 26 +++++++++++++------------- src/operators/transferpeps.jl | 19 ++++++++++++++++++- src/states/abstractpeps.jl | 4 ++-- src/utility/util.jl | 2 +- 6 files changed, 49 insertions(+), 34 deletions(-) diff --git a/examples/heisenberg.jl b/examples/heisenberg.jl index 2c9fa560..37ec0cc9 100644 --- a/examples/heisenberg.jl +++ b/examples/heisenberg.jl @@ -29,10 +29,10 @@ function iCtmGsEh(ψ::InfinitePEPS, env::PEPSKit.CTMRGEnv, H::AbstractTensorMap{ E = 0.0 for r in 1:size(ψ,1), c in 1:size(ψ,2) ρ = two_site_rho(r, c, ψ, env) - nn = @tensor ρ[1,2;1,2] - Eh = @tensor H[1,2;3,4]*ρ[1,2;3,4] + @tensor nn = ρ[1 2; 1 2] + @tensor Eh = H[1 2; 3 4] * ρ[1 2; 3 4] Eh = Eh / nn - E = E + Eh + E = E + Eh #@diffset Es[r,c] = Eh; end return real(E) @@ -40,13 +40,12 @@ end function H_expectation_value(ψ::InfinitePEPS, env::PEPSKit.CTMRGEnv, H::AbstractTensorMap{S,2,2}) where S Eh = iCtmGsEh(ψ, env, H) - ψ1 = rotl90(ψ) - env1 = PEPSKit.rotate_north(env,EAST); + env1 = PEPSKit.rotate_north(env, EAST) Ev = iCtmGsEh(ψ1, env1, H) E = real(Eh + Ev) return E -end +end function SqLatHeisenberg() Sx,Sy,Sz,_ = spinmatrices(1//2) @@ -72,13 +71,12 @@ function cfun(x) function fun(peps) env = leading_boundary(peps, alg_ctm, env) - x = H_expectation_value(peps, env, H) + x = H_expectation_value(peps, env, H) return x end - - ∂E = fun'(ψ) env = leading_boundary(ψ, alg_ctm, env) E = H_expectation_value(ψ, env, H) + ∂E = fun'(ψ) @assert !isnan(norm(∂E)) return E,∂E @@ -124,7 +122,7 @@ end alg_ctm = CTMRG( - verbose=10000, + verbose=1, tol=1e-10, trscheme=truncdim(10), miniter=4, @@ -137,7 +135,7 @@ function main(;d=2,D=2,Lx=1,Ly=1) optimize( cfun, (ψ,env), - ConjugateGradient(verbosity=3); + ConjugateGradient(verbosity=2); inner=my_inner, retract=my_retract, scale! = my_scale!, diff --git a/src/algorithms/ctmrg.jl b/src/algorithms/ctmrg.jl index 9561efda..e7edc4d6 100644 --- a/src/algorithms/ctmrg.jl +++ b/src/algorithms/ctmrg.jl @@ -29,7 +29,7 @@ function MPSKit.leading_boundary(peps_above::InfinitePEPS,peps_below::InfinitePE while (err>alg.tol&&iter<=alg.maxiter) || iter<=alg.miniter ϵ = 0.0 for i in 1:4 - envs,ϵ₀ = left_move(peps_above, peps_below,alg,envs); + envs,ϵ₀ = left_move(peps_above, peps_below, alg, envs); ϵ = max(ϵ,ϵ₀) envs = rotate_north(envs,EAST); peps_above = envs.peps_above; @@ -40,7 +40,7 @@ function MPSKit.leading_boundary(peps_above::InfinitePEPS,peps_below::InfinitePE err = abs(old_norm-new_norm) dϵ = abs((ϵ₁-ϵ)/ϵ₁) - @ignore_derivatives alg.verbose > 1 && @printf("%4d %.2e %.10e %.2e %.2e\n", + @ignore_derivatives alg.verbose > 1 && @printf("CTMRG: \titeration: %4d\t\terror: %.2e\t\tnorm: %.10e\t\tϵ: %.2e\t\tdϵ: %.2e\n", iter,err,abs(new_norm),ϵ,dϵ) old_norm = new_norm @@ -58,8 +58,8 @@ function left_move(peps_above::InfinitePEPS{PType},peps_below::InfinitePEPS{PTyp corners::typeof(envs.corners) = copy(envs.corners); edges::typeof(envs.edges) = copy(envs.edges); - above_projector_type = tensormaptype(spacetype(PType),1,3,storagetype(PType)); - below_projector_type = tensormaptype(spacetype(PType),3,1,storagetype(PType)); + above_projector_type = tensormaptype(spacetype(PType),1,3,storagetype(PType)) + below_projector_type = tensormaptype(spacetype(PType),3,1,storagetype(PType)) ϵ = 0.0 n0 = 1.0 n1 = 1.0 @@ -67,8 +67,8 @@ function left_move(peps_above::InfinitePEPS{PType},peps_below::InfinitePEPS{PTyp cop = mod1(col+1,size(peps_above,2)) com = mod1(col-1,size(peps_above,2)) - above_projs = Vector{above_projector_type}(undef,size(peps_above,1)); - below_projs = Vector{below_projector_type}(undef,size(peps_above,1)); + above_projs = Vector{above_projector_type}(undef,size(peps_above,1)) + below_projs = Vector{below_projector_type}(undef,size(peps_above,1)) # find all projectors for row in 1:size(peps_above,1) diff --git a/src/environments/ctmrgenv.jl b/src/environments/ctmrgenv.jl index ca204c36..fd7f7155 100644 --- a/src/environments/ctmrgenv.jl +++ b/src/environments/ctmrgenv.jl @@ -9,36 +9,36 @@ end CTMRGEnv(peps::InfinitePEPS) = CTMRGEnv(peps,peps); function CTMRGEnv(peps_above::InfinitePEPS{P},peps_below::InfinitePEPS{P}) where P - ou = oneunit(space(peps_above,1,1)); # the bogus space + ou = oneunit(space(peps_above,1,1)) # the bogus space - C_type = tensormaptype(spacetype(P),1,1,storagetype(P)); - T_type = tensormaptype(spacetype(P),3,1,storagetype(P)); # debatable how we should do the legs? + C_type = tensormaptype(spacetype(P),1,1,storagetype(P)) + T_type = tensormaptype(spacetype(P),3,1,storagetype(P)) # debatable how we should do the legs? #first index is de - corners = Array{C_type}(undef,4,size(peps_above)...); - edges = Array{T_type}(undef,4,size(peps_above)...); + corners = Array{C_type}(undef,4,size(peps_above)...) + edges = Array{T_type}(undef,4,size(peps_above)...) for dir in 1:4, i in 1:size(peps_above,1),j in 1:size(peps_above,2) @diffset corners[dir,i,j] = TensorMap(randn,eltype(P),ou,ou) @diffset edges[dir,i,j] = TensorMap(randn,eltype(P),ou*space(peps_above[i,j],dir+1)'*space(peps_below[i,j],dir+1),ou) end - @diffset corners[:,:,:]./=norm.(corners[:,:,:]); - @diffset edges[:,:,:]./=norm.(edges[:,:,:]); + @diffset corners[:,:,:]./=norm.(corners[:,:,:]) + @diffset edges[:,:,:]./=norm.(edges[:,:,:]) CTMRGEnv(peps_above,peps_below,corners,edges) end function Base.rotl90(envs::CTMRGEnv{P,C,T}) where {P,C,T} - n_peps_above = rotl90(envs.peps_above); - n_peps_below = rotl90(envs.peps_below); - n_corners = Array{C,3}(undef,size(envs.corners)...); - n_edges = Array{T,3}(undef,size(envs.edges)...); + n_peps_above = rotl90(envs.peps_above) + n_peps_below = rotl90(envs.peps_below) + n_corners = Array{C,3}(undef,size(envs.corners)...) + n_edges = Array{T,3}(undef,size(envs.edges)...) for dir in 1:4 dirm = mod1(dir-1,4) - @diffset n_corners[dirm,:,:] .= rotl90(envs.corners[dir,:,:]); - @diffset n_edges[dirm,:,:] .= rotl90(envs.edges[dir,:,:]); + @diffset n_corners[dirm,:,:] .= rotl90(envs.corners[dir,:,:]) + @diffset n_edges[dirm,:,:] .= rotl90(envs.edges[dir,:,:]) end CTMRGEnv(n_peps_above,n_peps_below,n_corners,n_edges) diff --git a/src/operators/transferpeps.jl b/src/operators/transferpeps.jl index ac972bb7..f41f3505 100644 --- a/src/operators/transferpeps.jl +++ b/src/operators/transferpeps.jl @@ -61,4 +61,21 @@ end function MPSKit.transfer_right(GR::GenericMPSTensor{S,3}, O::NTuple{2, PEPSTensor}, A::GenericMPSTensor{S,3}, Ā::GenericMPSTensor{S,3}) where S return @tensor GR′[-1 -2 -3; -4] := GR[7 4 2; 1] * conj(Ā[-4 6 3; 1]) * O[1][5; 9 4 6 -2] * conj(O[2][5; 8 2 3 -3]) * A[-1 9 8 7] -end \ No newline at end of file +end + +function MPSKit.expectation_value(st::MPSMultiline, O::TransferPEPSMultiline) + +end + +function MPSKit.expectation_value(st::MPSMultiline, ca::MPSKit.PerMPOInfEnv{H,V,S,A}) where {H<:TransferPEPSMultiline,V,S,A} + opp = ca.opp + retval = PeriodicArray{eltype(st.AC[1, 1]),2}(undef, size(st, 1), size(st, 2)) + for (i, j) in product(1:size(st, 1), 1:size(st, 2)) + retval[i, j] = @plansor leftenv(ca, i, j, st)[1 2; 3] * + opp[i, j][2 4; 6 5] * + st.AC[i, j][3 6; 7] * + rightenv(ca, i, j, st)[7 5; 8] * + conj(st.AC[i+1, j][1 4; 8]) + end + return retval +end diff --git a/src/states/abstractpeps.jl b/src/states/abstractpeps.jl index 80f58beb..12e5a762 100644 --- a/src/states/abstractpeps.jl +++ b/src/states/abstractpeps.jl @@ -40,5 +40,5 @@ Abstract supertype for a 2D projected entangled pairs operator. abstract type AbstractPEPO end -Base.rotl90(t::PEPSTensor) = permute(t,(1,),(3,4,5,2)); -Base.rotl90(t::PEPOTensor) = permute(t,(1,2),(4,5,6,3)); +Base.rotl90(t::PEPSTensor) = permute(t, ((1,), (3, 4, 5, 2))) +Base.rotl90(t::PEPOTensor) = permute(t, ((1, 2), (4, 5, 6, 3))); diff --git a/src/utility/util.jl b/src/utility/util.jl index 91ea4968..a78e8b16 100644 --- a/src/utility/util.jl +++ b/src/utility/util.jl @@ -47,7 +47,7 @@ function _setindex(a::AbstractArray,v,args...) b end -function ChainRulesCore.rrule(::typeof(_setindex),a::AbstractArray,tv,args...) +function ChainRulesCore.rrule(::typeof(_setindex), a::AbstractArray, tv, args...) t = _setindex(a,tv,args...); function toret(v) From c66b3a4de5a01b11ada587cd41c02776adedc80e Mon Sep 17 00:00:00 2001 From: leburgel Date: Fri, 29 Sep 2023 19:18:16 +0200 Subject: [PATCH 3/9] Update - `eltype` -> `scalartype` for `TensorMap`s - Clean up README - Set up tests stub - Update CI - Set up docs stub - Bump TensorKit compat - Format --- .JuliaFormatter.toml | 1 + .github/workflows/ci.yml | 58 +++- .github/workflows/docs.yml | 35 +++ .github/workflows/formatter.yml | 67 +++++ Project.toml | 5 + README.md | 30 +- docs/Project.toml | 5 + docs/make.jl | 28 ++ docs/src/assets/logo-dark.svg | 111 ++++++++ docs/src/assets/logo.svg | 105 +++++++ docs/src/examples/index.md | 1 + docs/src/index.md | 5 + docs/src/lib/lib.md | 1 + docs/src/man/intro.md | 1 + examples/heisenberg.jl | 103 ++++--- examples/test_vumps.jl | 9 +- src/PEPSKit.jl | 72 ++--- src/algorithms/ctmrg.jl | 280 ++++++++++++------- src/algorithms/expval.jl | 56 ++-- src/environments/ctmrgenv.jl | 45 +-- src/mpskit_glue/transferpeps_environments.jl | 72 +++-- src/operators/derivatives.jl | 60 ++-- src/operators/periodicpepo.jl | 20 +- src/operators/transferpeps.jl | 81 +++--- src/states/abstractpeps.jl | 6 +- src/states/infinitepeps.jl | 7 +- src/utility/rotations.jl | 2 +- src/utility/util.jl | 122 ++++---- test/Project.toml | 3 +- test/runtests.jl | 11 +- 30 files changed, 976 insertions(+), 426 deletions(-) create mode 100644 .JuliaFormatter.toml create mode 100644 .github/workflows/docs.yml create mode 100644 .github/workflows/formatter.yml create mode 100644 docs/Project.toml create mode 100644 docs/make.jl create mode 100644 docs/src/assets/logo-dark.svg create mode 100644 docs/src/assets/logo.svg create mode 100644 docs/src/examples/index.md create mode 100644 docs/src/index.md create mode 100644 docs/src/lib/lib.md create mode 100644 docs/src/man/intro.md diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml new file mode 100644 index 00000000..c7439503 --- /dev/null +++ b/.JuliaFormatter.toml @@ -0,0 +1 @@ +style = "blue" \ No newline at end of file diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 94aadd8a..17148087 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -1,17 +1,28 @@ name: CI on: - - push - - pull_request + push: + branches: + - 'master' + - 'main' + - 'release-' + tags: '*' + pull_request: + workflow_dispatch: + +concurrency: + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: true + jobs: test: name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} runs-on: ${{ matrix.os }} - continue-on-error: ${{ matrix.version == 'nightly' }} strategy: + fail-fast: false matrix: version: - - '1.5' - - 'nightly' + - '1.6' # LTS version + - '1' # automatically expands to the latest stable 1.x release of Julia os: - ubuntu-latest - macOS-latest @@ -19,14 +30,43 @@ jobs: arch: - x64 steps: - - uses: actions/checkout@v2 - - uses: julia-actions/setup-julia@latest + - uses: actions/checkout@v4 + - uses: julia-actions/setup-julia@v1 with: version: ${{ matrix.version }} arch: ${{ matrix.arch }} + - uses: julia-actions/cache@v1 - uses: julia-actions/julia-buildpkg@latest - uses: julia-actions/julia-runtest@latest - - uses: julia-actions/julia-uploadcodecov@latest env: JULIA_NUM_THREADS: 4 - CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }} + - uses: julia-actions/julia-processcoverage@v1 + - uses: codecov/codecov-action@v3 + with: + file: lcov.info + test-nightly: + needs: test + name: Julia nightly - ${{ matrix.os }} - ${{ matrix.arch }} + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + version: + - 'nightly' + os: + - ubuntu-latest + - macOS-latest + - windows-latest + arch: + - x64 + steps: + - uses: actions/checkout@v4 + - uses: julia-actions/setup-julia@v1 + with: + version: ${{ matrix.version }} + arch: ${{ matrix.arch }} + - uses: julia-actions/cache@v1 + - uses: julia-actions/julia-buildpkg@latest + - uses: julia-actions/julia-runtest@latest + env: + JULIA_NUM_THREADS: 4 \ No newline at end of file diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml new file mode 100644 index 00000000..f3269903 --- /dev/null +++ b/.github/workflows/docs.yml @@ -0,0 +1,35 @@ +name: Documentation + +on: + push: + branches: + - 'master' + - 'main' + - 'release-' + tags: '*' + pull_request: + +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: actions/checkout@v4 + - uses: julia-actions/setup-julia@latest + with: + version: ${{ matrix.version }} + arch: ${{ matrix.arch }} + - name: Install dependencies + run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()' + - name: Build and deploy + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # For authentication with GitHub Actions token + DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # For authentication with SSH deploy key + run: julia --project=docs/ docs/make.jl \ No newline at end of file diff --git a/.github/workflows/formatter.yml b/.github/workflows/formatter.yml new file mode 100644 index 00000000..a57dbf68 --- /dev/null +++ b/.github/workflows/formatter.yml @@ -0,0 +1,67 @@ +name: format-check + +on: + push: + branches: + - 'main' + - 'master' + - 'release-' + tags: '*' + pull_request: + +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@v4 + - name: Install JuliaFormatter and format + # This will use the latest version by default but you can set the version like so: + # + # julia -e 'using Pkg; Pkg.add(PackageSpec(name="JuliaFormatter", version="0.13.0"))' + run: | + julia -e 'using Pkg; Pkg.add(PackageSpec(name="JuliaFormatter"))' + 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@v3 + 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 + + - 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 }}" diff --git a/Project.toml b/Project.toml index bf745393..49957658 100644 --- a/Project.toml +++ b/Project.toml @@ -13,3 +13,8 @@ OptimKit = "77e91f04-9b3b-57a6-a776-40b61faaebe0" Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec" + +[compat] +julia = "1.6" +MPSKit = "0.9.1" +TensorKit = "0.12" diff --git a/README.md b/README.md index 35955e17..32f36fab 100644 --- a/README.md +++ b/README.md @@ -1,24 +1,20 @@ -# pepskit + + + PEPSKit.jl logo + -This is some semi-public code (I use it - it doesn't work - there is no documentation - it's still public) +# PepsKit.jl -In this package we +[![docs][docs-dev-img]][docs-dev-url] ![CI][ci-url] [![codecov][codecov-img]][codecov-url] - - define states +[docs-dev-img]: https://img.shields.io/badge/docs-dev-blue.svg +[docs-dev-url]: https://quantumghent.github.io/PEPSKit.jl/dev/ - - define operators +[codecov-img]: https://codecov.io/gh/quantumghent/PEPSKit.jl/branch/master/graph/badge.svg +[codecov-url]: https://codecov.io/gh/quantumghent/PEPSKit.jl - - define environments +[ci-url]: https://github.com/quantumghent/PEPSKit.jl/workflows/CI/badge.svg - - some are independent of the operator (boundary mpses,channels,...) +**Tools for working with projected entangled-pair states** - - some are dependent on the operator (hamiltonian channels) - - - define algorithms - - -To minimize the amount of complexity/code we typically - - - define operations acting on the north (finding the north boundary, finding the north hamiltonian channels) - - - define how to rotate the peps such that other directions now become north +It contracts, it optimizes, it may be broken at any point. diff --git a/docs/Project.toml b/docs/Project.toml new file mode 100644 index 00000000..3a52a5db --- /dev/null +++ b/docs/Project.toml @@ -0,0 +1,5 @@ +[deps] +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" + +[compat] +Documenter = "0.27" diff --git a/docs/make.jl b/docs/make.jl new file mode 100644 index 00000000..143b0c4b --- /dev/null +++ b/docs/make.jl @@ -0,0 +1,28 @@ +using Documenter +using PEPSKit + +makedocs(; + modules=[PEPSKit], + sitename="PEPSKit.jl", + format=Documenter.HTML(; + prettyurls=get(ENV, "CI", nothing) == "true", + mathengine=MathJax3( + Dict( + :loader => Dict("load" => ["[tex]/physics"]), + :tex => Dict( + "inlineMath" => [["\$", "\$"], ["\\(", "\\)"]], + "tags" => "ams", + "packages" => ["base", "ams", "autoload", "physics"], + ), + ), + ), + ), + pages=[ + "Home" => "index.md", + "Manual" => "man/intro.md", + "Examples" => "examples/index.md", + "Library" => "lib/lib.md", + ], +) + +deploydocs(; repo="https://github.com/quantumghent/PEPSKit.jl.git") diff --git a/docs/src/assets/logo-dark.svg b/docs/src/assets/logo-dark.svg new file mode 100644 index 00000000..9d8f3844 --- /dev/null +++ b/docs/src/assets/logo-dark.svg @@ -0,0 +1,111 @@ + + + + + + image/svg+xml + + + + + + + + + + + + + + + + diff --git a/docs/src/assets/logo.svg b/docs/src/assets/logo.svg new file mode 100644 index 00000000..ca13918e --- /dev/null +++ b/docs/src/assets/logo.svg @@ -0,0 +1,105 @@ + + + + + + image/svg+xml + + + + + + + + + + + + + + + + + diff --git a/docs/src/examples/index.md b/docs/src/examples/index.md new file mode 100644 index 00000000..fcc39fe6 --- /dev/null +++ b/docs/src/examples/index.md @@ -0,0 +1 @@ +Coming soon. \ No newline at end of file diff --git a/docs/src/index.md b/docs/src/index.md new file mode 100644 index 00000000..9da83b7a --- /dev/null +++ b/docs/src/index.md @@ -0,0 +1,5 @@ +# PEPSKit.jl + +**Tools for working with projected entangled-pair states** + +It contracts, it optimizes, it may be broken at any point. diff --git a/docs/src/lib/lib.md b/docs/src/lib/lib.md new file mode 100644 index 00000000..fcc39fe6 --- /dev/null +++ b/docs/src/lib/lib.md @@ -0,0 +1 @@ +Coming soon. \ No newline at end of file diff --git a/docs/src/man/intro.md b/docs/src/man/intro.md new file mode 100644 index 00000000..fcc39fe6 --- /dev/null +++ b/docs/src/man/intro.md @@ -0,0 +1 @@ +Coming soon. \ No newline at end of file diff --git a/examples/heisenberg.jl b/examples/heisenberg.jl index 37ec0cc9..7e8581d7 100644 --- a/examples/heisenberg.jl +++ b/examples/heisenberg.jl @@ -1,33 +1,35 @@ using Revise, PEPSKit, TensorKit, Zygote, MPSKit using MPSKitModels, LinearAlgebra, OptimKit -using PEPSKit: NORTH,SOUTH,WEST,EAST,NORTHWEST,NORTHEAST,SOUTHEAST,SOUTHWEST,@diffset -using JLD2,ChainRulesCore - +using PEPSKit: + NORTH, SOUTH, WEST, EAST, NORTHWEST, NORTHEAST, SOUTHEAST, SOUTHWEST, @diffset +using JLD2, ChainRulesCore function two_site_rho(r::Int, c::Int, ψ::InfinitePEPS, env::PEPSKit.CTMRGEnv) - cp = mod1(c+1,size(ψ,2)) - @tensor ρ[-11,-20;-12,-18] := env.corners[NORTHWEST,r,c][1,3] * - env.edges[WEST,r,c][2,7,9,1] * - env.corners[SOUTHWEST,r,c][4,2] * - env.edges[NORTH,r,c][3,5,8,13] * - env.edges[SOUTH,r,c][14,6,10,4] * - ψ[r,c][-12,5,15,6,7] * - conj(ψ[r,c][-11,8,19,10,9]) * - env.edges[NORTH,r,cp][13,16,22,23] * - env.edges[SOUTH,r,cp][28,17,21,14] * - ψ[r,cp][-18,16,25,17,15] * - conj(ψ[r,cp][-20,22,26,21,19]) * - env.corners[NORTHEAST,r,cp][23,24] * - env.edges[EAST,r,cp][24,25,26,27] * - env.corners[SOUTHEAST,r,cp][27,28] + cp = mod1(c + 1, size(ψ, 2)) + @tensor ρ[-11, -20; -12, -18] := + env.corners[NORTHWEST, r, c][1, 3] * + env.edges[WEST, r, c][2, 7, 9, 1] * + env.corners[SOUTHWEST, r, c][4, 2] * + env.edges[NORTH, r, c][3, 5, 8, 13] * + env.edges[SOUTH, r, c][14, 6, 10, 4] * + ψ[r, c][-12, 5, 15, 6, 7] * + conj(ψ[r, c][-11, 8, 19, 10, 9]) * + env.edges[NORTH, r, cp][13, 16, 22, 23] * + env.edges[SOUTH, r, cp][28, 17, 21, 14] * + ψ[r, cp][-18, 16, 25, 17, 15] * + conj(ψ[r, cp][-20, 22, 26, 21, 19]) * + env.corners[NORTHEAST, r, cp][23, 24] * + env.edges[EAST, r, cp][24, 25, 26, 27] * + env.corners[SOUTHEAST, r, cp][27, 28] return ρ end - -function iCtmGsEh(ψ::InfinitePEPS, env::PEPSKit.CTMRGEnv, H::AbstractTensorMap{S,2,2}) where S +function iCtmGsEh( + ψ::InfinitePEPS, env::PEPSKit.CTMRGEnv, H::AbstractTensorMap{S,2,2} +) where {S} #Es = Matrix{eltype(H)}(undef,size(ψ,1),size(ψ,2)) E = 0.0 - for r in 1:size(ψ,1), c in 1:size(ψ,2) + for r in 1:size(ψ, 1), c in 1:size(ψ, 2) ρ = two_site_rho(r, c, ψ, env) @tensor nn = ρ[1 2; 1 2] @tensor Eh = H[1 2; 3 4] * ρ[1 2; 3 4] @@ -38,36 +40,35 @@ function iCtmGsEh(ψ::InfinitePEPS, env::PEPSKit.CTMRGEnv, H::AbstractTensorMap{ return real(E) end -function H_expectation_value(ψ::InfinitePEPS, env::PEPSKit.CTMRGEnv, H::AbstractTensorMap{S,2,2}) where S +function H_expectation_value( + ψ::InfinitePEPS, env::PEPSKit.CTMRGEnv, H::AbstractTensorMap{S,2,2} +) where {S} Eh = iCtmGsEh(ψ, env, H) ψ1 = rotl90(ψ) env1 = PEPSKit.rotate_north(env, EAST) Ev = iCtmGsEh(ψ1, env1, H) E = real(Eh + Ev) return E -end +end function SqLatHeisenberg() - Sx,Sy,Sz,_ = spinmatrices(1//2) + Sx, Sy, Sz, _ = spinmatrices(1//2) Dphys = ComplexSpace(2) σx = TensorMap(Sx, Dphys, Dphys) σy = TensorMap(Sy, Dphys, Dphys) σz = TensorMap(Sz, Dphys, Dphys) - @tensor H[-1 -3; -2 -4] := - -σx[-1,-2] * σx[-3,-4] + - σy[-1,-2] * σy[-3,-4] + - -σz[-1,-2] * σz[-3,-4] + @tensor H[-1 -3; -2 -4] := + -σx[-1, -2] * σx[-3, -4] + σy[-1, -2] * σy[-3, -4] + -σz[-1, -2] * σz[-3, -4] return H end H = SqLatHeisenberg() - function cfun(x) - (ψ,env) = x + (ψ, env) = x function fun(peps) env = leading_boundary(peps, alg_ctm, env) @@ -79,20 +80,20 @@ function cfun(x) ∂E = fun'(ψ) @assert !isnan(norm(∂E)) - return E,∂E + return E, ∂E end # my_retract is not an in place function which should not change x -function my_retract(x,dx,α::Number) - (ϕ,env0) = x +function my_retract(x, dx, α::Number) + (ϕ, env0) = x ψ = deepcopy(ϕ) env = deepcopy(env0) ψ.A .+= dx.A .* α #env = leading_boundary(ψ, alg_ctm,env) - return (ψ,env),dx + return (ψ, env), dx end -my_inner(x,dx1,dx2) = real(dot(dx1,dx2)) +my_inner(x, dx1, dx2) = real(dot(dx1, dx2)) function my_add!(Y, X, a) Y.A .+= X.A .* a @@ -104,11 +105,10 @@ function my_scale!(η, β) return η end - function init_psi(d::Int, D::Int, Lx::Int, Ly::Int) - Pspaces = fill(ℂ^d,Lx,Ly) - Nspaces = fill(ℂ^D,Lx,Ly) - Espaces = fill(ℂ^D,Lx,Ly) + Pspaces = fill(ℂ^d, Lx, Ly) + Nspaces = fill(ℂ^D, Lx, Ly) + Espaces = fill(ℂ^D, Lx, Ly) Sspaces = adjoint.(circshift(Nspaces, (1, 0))) Wspaces = adjoint.(circshift(Espaces, (0, -1))) @@ -120,26 +120,19 @@ function init_psi(d::Int, D::Int, Lx::Int, Ly::Int) return InfinitePEPS(A) end +alg_ctm = CTMRG(; verbose=1, tol=1e-4, trscheme=truncdim(10), miniter=4, maxiter=200) -alg_ctm = CTMRG( - verbose=1, - tol=1e-10, - trscheme=truncdim(10), - miniter=4, - maxiter=200 - ) - -function main(;d=2,D=2,Lx=1,Ly=1) - ψ = init_psi(d,D,Lx,Ly) - env = leading_boundary(ψ, alg_ctm) +function main(; d=2, D=2, Lx=1, Ly=1) + ψ = init_psi(d, D, Lx, Ly) + env = leading_boundary(ψ, alg_ctm) optimize( - cfun, - (ψ,env), - ConjugateGradient(verbosity=2); + cfun, + (ψ, env), + ConjugateGradient(; verbosity=2); inner=my_inner, retract=my_retract, - scale! = my_scale!, - add! = my_add! + (scale!)=my_scale!, + (add!)=my_add!, ) return ψ end diff --git a/examples/test_vumps.jl b/examples/test_vumps.jl index fb92a213..b3492c2c 100644 --- a/examples/test_vumps.jl +++ b/examples/test_vumps.jl @@ -1,9 +1,8 @@ using Revise, PEPSKit, TensorKit, Zygote, MPSKit -p = InfinitePEPS(fill(ℂ^2,1,1),fill(ℂ^2,1,1)); +p = InfinitePEPS(fill(ℂ^2, 1, 1), fill(ℂ^2, 1, 1)); -trans = PEPSKit.InfiniteTransferPEPS(p,1,1); -mps = PEPSKit.initializeMPS(trans,[ℂ^5]); - -leading_boundary(mps,trans,VUMPS()) +trans = PEPSKit.InfiniteTransferPEPS(p, 1, 1); +mps = PEPSKit.initializeMPS(trans, [ℂ^5]); +leading_boundary(mps, trans, VUMPS()) diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index 8e5f43dc..b8f0eb22 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -1,38 +1,38 @@ module PEPSKit - using TensorKit, KrylovKit, MPSKit, OptimKit, Base.Threads, Base.Iterators, Parameters,Printf - using ChainRulesCore; - - import LinearAlgebra - - export CTMRG,CTMRG2 - export leading_boundary - - include("utility/util.jl") - - include("states/abstractpeps.jl") - include("states/infinitepeps.jl") - - include("operators/transferpeps.jl") - include("operators/derivatives.jl") - include("operators/periodicpepo.jl") - - include("environments/ctmrgenv.jl") - - include("mpskit_glue/transferpeps_environments.jl") - - include("algorithms/ctmrg.jl") - include("algorithms/expval.jl") - - include("utility/rotations.jl") - - - #default settings - module Defaults - const maxiter = 100 - const tol = 1e-12 - end - - export InfinitePEPS, InfiniteTransferPEPS - export PeriodicPEPO - export initializeMPS +using TensorKit, + KrylovKit, MPSKit, OptimKit, Base.Threads, Base.Iterators, Parameters, Printf +using ChainRulesCore + +using LinearAlgebra: LinearAlgebra + +export CTMRG, CTMRG2 +export leading_boundary + +include("utility/util.jl") + +include("states/abstractpeps.jl") +include("states/infinitepeps.jl") + +include("operators/transferpeps.jl") +include("operators/derivatives.jl") +include("operators/periodicpepo.jl") + +include("environments/ctmrgenv.jl") + +include("mpskit_glue/transferpeps_environments.jl") + +include("algorithms/ctmrg.jl") +include("algorithms/expval.jl") + +include("utility/rotations.jl") + +#default settings +module Defaults + const maxiter = 100 + const tol = 1e-12 +end + +export InfinitePEPS, InfiniteTransferPEPS +export PeriodicPEPO +export initializeMPS end # module diff --git a/src/algorithms/ctmrg.jl b/src/algorithms/ctmrg.jl index e7edc4d6..932ed866 100644 --- a/src/algorithms/ctmrg.jl +++ b/src/algorithms/ctmrg.jl @@ -15,9 +15,16 @@ end verbose::Integer = 0 end -MPSKit.leading_boundary(peps::InfinitePEPS,alg::CTMRG,envs = CTMRGEnv(peps)) = MPSKit.leading_boundary(peps,peps,alg,envs); - -function MPSKit.leading_boundary(peps_above::InfinitePEPS,peps_below::InfinitePEPS,alg::CTMRG,envs = CTMRGEnv(peps_above,peps_below)) +function MPSKit.leading_boundary(peps::InfinitePEPS, alg::CTMRG, envs=CTMRGEnv(peps)) + return MPSKit.leading_boundary(peps, peps, alg, envs) +end; + +function MPSKit.leading_boundary( + peps_above::InfinitePEPS, + peps_below::InfinitePEPS, + alg::CTMRG, + envs=CTMRGEnv(peps_above, peps_below), +) err = Inf iter = 1 @@ -26,22 +33,28 @@ function MPSKit.leading_boundary(peps_above::InfinitePEPS,peps_below::InfinitePE old_norm = 1.0 new_norm = old_norm ϵ₁ = 1.0 - while (err>alg.tol&&iter<=alg.maxiter) || iter<=alg.miniter + while (err > alg.tol && iter <= alg.maxiter) || iter <= alg.miniter ϵ = 0.0 for i in 1:4 - envs,ϵ₀ = left_move(peps_above, peps_below, alg, envs); - ϵ = max(ϵ,ϵ₀) - envs = rotate_north(envs,EAST); - peps_above = envs.peps_above; - peps_below = envs.peps_below; + envs, ϵ₀ = left_move(peps_above, peps_below, alg, envs) + ϵ = max(ϵ, ϵ₀) + envs = rotate_north(envs, EAST) + peps_above = envs.peps_above + peps_below = envs.peps_below end - + new_norm = contract_ctrmg(envs) - err = abs(old_norm-new_norm) - dϵ = abs((ϵ₁-ϵ)/ϵ₁) - @ignore_derivatives alg.verbose > 1 && @printf("CTMRG: \titeration: %4d\t\terror: %.2e\t\tnorm: %.10e\t\tϵ: %.2e\t\tdϵ: %.2e\n", - iter,err,abs(new_norm),ϵ,dϵ) + err = abs(old_norm - new_norm) + dϵ = abs((ϵ₁ - ϵ) / ϵ₁) + @ignore_derivatives alg.verbose > 1 && @printf( + "CTMRG: \titeration: %4d\t\terror: %.2e\t\tnorm: %.10e\t\tϵ: %.2e\t\tdϵ: %.2e\n", + iter, + err, + abs(new_norm), + ϵ, + dϵ + ) old_norm = new_norm ϵ₁ = ϵ @@ -49,88 +62,135 @@ function MPSKit.leading_boundary(peps_above::InfinitePEPS,peps_below::InfinitePE end #@ignore_derivatives @show iter, new_norm, err - @ignore_derivatives iter > alg.maxiter && alg.verbose > 0 && @warn "maxiter $(alg.maxiter) reached: error was $(err)" + @ignore_derivatives iter > alg.maxiter && + alg.verbose > 0 && + @warn "maxiter $(alg.maxiter) reached: error was $(err)" return envs end -function left_move(peps_above::InfinitePEPS{PType},peps_below::InfinitePEPS{PType},alg::CTMRG,envs::CTMRGEnv) where PType - corners::typeof(envs.corners) = copy(envs.corners); - edges::typeof(envs.edges) = copy(envs.edges); - - above_projector_type = tensormaptype(spacetype(PType),1,3,storagetype(PType)) - below_projector_type = tensormaptype(spacetype(PType),3,1,storagetype(PType)) +function left_move( + peps_above::InfinitePEPS{PType}, + peps_below::InfinitePEPS{PType}, + alg::CTMRG, + envs::CTMRGEnv, +) where {PType} + corners::typeof(envs.corners) = copy(envs.corners) + edges::typeof(envs.edges) = copy(envs.edges) + + above_projector_type = tensormaptype(spacetype(PType), 1, 3, storagetype(PType)) + below_projector_type = tensormaptype(spacetype(PType), 3, 1, storagetype(PType)) ϵ = 0.0 n0 = 1.0 n1 = 1.0 - for col in 1:size(peps_above,2) - cop = mod1(col+1,size(peps_above,2)) - com = mod1(col-1,size(peps_above,2)) + for col in 1:size(peps_above, 2) + cop = mod1(col + 1, size(peps_above, 2)) + com = mod1(col - 1, size(peps_above, 2)) - above_projs = Vector{above_projector_type}(undef,size(peps_above,1)) - below_projs = Vector{below_projector_type}(undef,size(peps_above,1)) + above_projs = Vector{above_projector_type}(undef, size(peps_above, 1)) + below_projs = Vector{below_projector_type}(undef, size(peps_above, 1)) # find all projectors - for row in 1:size(peps_above,1) - rop = mod1(row+1, size(peps_above,1)) - peps_above_nw = peps_above[row,col]; - peps_above_sw = rotate_north(peps_above[rop,col],WEST); - peps_below_nw = peps_below[row,col]; - peps_below_sw = rotate_north(peps_below[rop,col],WEST); - - Q1 = northwest_corner(envs.edges[SOUTH,mod1(row+1,end),col],envs.corners[SOUTHWEST,mod1(row+1,end),col],envs.edges[WEST,mod1(row+1,end),col],peps_above_sw,peps_below_sw); - Q2 = northwest_corner(envs.edges[WEST,row,col],envs.corners[NORTHWEST,row,col],envs.edges[NORTH,row,col],peps_above_nw,peps_below_nw); + for row in 1:size(peps_above, 1) + rop = mod1(row + 1, size(peps_above, 1)) + peps_above_nw = peps_above[row, col] + peps_above_sw = rotate_north(peps_above[rop, col], WEST) + peps_below_nw = peps_below[row, col] + peps_below_sw = rotate_north(peps_below[rop, col], WEST) + + Q1 = northwest_corner( + envs.edges[SOUTH, mod1(row + 1, end), col], + envs.corners[SOUTHWEST, mod1(row + 1, end), col], + envs.edges[WEST, mod1(row + 1, end), col], + peps_above_sw, + peps_below_sw, + ) + Q2 = northwest_corner( + envs.edges[WEST, row, col], + envs.corners[NORTHWEST, row, col], + envs.edges[NORTH, row, col], + peps_above_nw, + peps_below_nw, + ) + + trscheme = if alg.fixedspace == true + truncspace(space(envs.edges[WEST, row, cop], 1)) + else + alg.trscheme + end + #@ignore_derivatives @show norm(Q1*Q2) + (U, S, V) = tsvd(Q1 * Q2; trunc=trscheme, alg=SVD()) - trscheme = alg.fixedspace == true ? truncspace(space(envs.edges[WEST,row,cop],1)) : alg.trscheme - #@ignore_derivatives @show norm(Q1*Q2) - - (U,S,V) = tsvd(Q1*Q2,trunc = trscheme,alg = SVD()) - - @ignore_derivatives n0 = norm(Q1*Q2)^2 - @ignore_derivatives n1 = norm(U*S*V)^2 - @ignore_derivatives ϵ = max(ϵ, (n0-n1)/n0) + @ignore_derivatives n0 = norm(Q1 * Q2)^2 + @ignore_derivatives n1 = norm(U * S * V)^2 + @ignore_derivatives ϵ = max(ϵ, (n0 - n1) / n0) - isqS = sdiag_inv_sqrt(S); + isqS = sdiag_inv_sqrt(S) #Q = isqS*U'*Q1; #P = Q2*V'*isqS; - @tensor Q[-1;-2 -3 -4] := isqS[-1;1] * conj(U[2 3 4;1]) * Q1[2 3 4;-2 -3 -4] - @tensor P[-1 -2 -3;-4] := Q2[-1 -2 -3;1 2 3] * conj(V[4;1 2 3]) * isqS[4;-4] + @tensor Q[-1; -2 -3 -4] := isqS[-1; 1] * conj(U[2 3 4; 1]) * Q1[2 3 4; -2 -3 -4] + @tensor P[-1 -2 -3; -4] := Q2[-1 -2 -3; 1 2 3] * conj(V[4; 1 2 3]) * isqS[4; -4] - @diffset above_projs[row] = Q; - @diffset below_projs[row] = P; + @diffset above_projs[row] = Q + @diffset below_projs[row] = P end - - + #use the projectors to grow the corners/edges - for row in 1:size(peps_above,1) - Q = above_projs[row]; - P = below_projs[mod1(row-1,end)]; - rop = mod1(row+1,size(peps_above,1)) - rom = mod1(row-1,size(peps_above,1)) - - @diffset @tensor corners[NORTHWEST,rop,cop][-1;-2] := envs.corners[NORTHWEST,rop,col][1,2] * envs.edges[NORTH,rop,col][2,3,4,-2]*Q[-1;1 3 4] - @diffset @tensor corners[SOUTHWEST,rom,cop][-1;-2] := envs.corners[SOUTHWEST,rom,col][1,4] * envs.edges[SOUTH,rom,col][-1,2,3,1]*P[4 2 3;-2] - @diffset @tensor edges[WEST,row,cop][-1 -2 -3;-4] := envs.edges[WEST,row,col][1 2 3;4]* - peps_above[row,col][9;5 -2 7 2]* - conj(peps_below[row,col][9;6 -3 8 3])* - P[4 5 6;-4]* - Q[-1;1 7 8] + for row in 1:size(peps_above, 1) + Q = above_projs[row] + P = below_projs[mod1(row - 1, end)] + rop = mod1(row + 1, size(peps_above, 1)) + rom = mod1(row - 1, size(peps_above, 1)) + + @diffset @tensor corners[NORTHWEST, rop, cop][-1; -2] := + envs.corners[NORTHWEST, rop, col][1, 2] * + envs.edges[NORTH, rop, col][2, 3, 4, -2] * + Q[-1; 1 3 4] + @diffset @tensor corners[SOUTHWEST, rom, cop][-1; -2] := + envs.corners[SOUTHWEST, rom, col][1, 4] * + envs.edges[SOUTH, rom, col][-1, 2, 3, 1] * + P[4 2 3; -2] + @diffset @tensor edges[WEST, row, cop][-1 -2 -3; -4] := + envs.edges[WEST, row, col][1 2 3; 4] * + peps_above[row, col][9; 5 -2 7 2] * + conj(peps_below[row, col][9; 6 -3 8 3]) * + P[4 5 6; -4] * + Q[-1; 1 7 8] end - - @diffset corners[NORTHWEST,:,cop]./=norm.(corners[NORTHWEST,:,cop]); - @diffset edges[WEST,:,cop]./=norm.(edges[WEST,:,cop]); - @diffset corners[SOUTHWEST,:,cop]./=norm.(corners[SOUTHWEST,:,cop]); + @diffset corners[NORTHWEST, :, cop] ./= norm.(corners[NORTHWEST, :, cop]) + @diffset edges[WEST, :, cop] ./= norm.(edges[WEST, :, cop]) + @diffset corners[SOUTHWEST, :, cop] ./= norm.(corners[SOUTHWEST, :, cop]) end - - - return CTMRGEnv(peps_above,peps_below,corners,edges), ϵ + + return CTMRGEnv(peps_above, peps_below, corners, edges), ϵ end -northwest_corner(E4,C1,E1,peps_above,peps_below=peps_above) = @tensor corner[-1 -2 -3;-4 -5 -6] := E4[-1 1 2;3]*C1[3;4]*E1[4 5 6;-4]*peps_above[7;5 -5 -2 1]*conj(peps_below[7;6 -6 -3 2]) -northeast_corner(E1,C2,E2,peps_above,peps_below=peps_above) = @tensor corner[-1 -2 -3;-4 -5 -6] := E1[-1 1 2;3]*C2[3;4]*E2[4 5 6;-4]*peps_above[7;1 5 -5 -2]*conj(peps_below[7;2 6 -6 -3]) -southeast_corner(E2,C3,E3,peps_above,peps_below=peps_above) = @tensor corner[-1 -2 -3;-4 -5 -6] := E2[-1 1 2;3]*C3[3;4]*E3[4 5 6;-4]*peps_above[7;-2 1 5 -5]*conj(peps_below[7;-3 2 6 -6]) +function northwest_corner(E4, C1, E1, peps_above, peps_below=peps_above) + @tensor corner[-1 -2 -3; -4 -5 -6] := + E4[-1 1 2; 3] * + C1[3; 4] * + E1[4 5 6; -4] * + peps_above[7; 5 -5 -2 1] * + conj(peps_below[7; 6 -6 -3 2]) +end +function northeast_corner(E1, C2, E2, peps_above, peps_below=peps_above) + @tensor corner[-1 -2 -3; -4 -5 -6] := + E1[-1 1 2; 3] * + C2[3; 4] * + E2[4 5 6; -4] * + peps_above[7; 1 5 -5 -2] * + conj(peps_below[7; 2 6 -6 -3]) +end +function southeast_corner(E2, C3, E3, peps_above, peps_below=peps_above) + @tensor corner[-1 -2 -3; -4 -5 -6] := + E2[-1 1 2; 3] * + C3[3; 4] * + E3[4 5 6; -4] * + peps_above[7; -2 1 5 -5] * + conj(peps_below[7; -3 2 6 -6]) +end #= @@ -152,7 +212,6 @@ function MPSKit.leading_boundary(peps::InfinitePEPS,alg::CTMRG2,envs = CTMRGEnv( #@show new_norm err = abs(old_norm-new_norm) @ignore_derivatives mod(alg.verbose,alg.miniter)==0 && mod(iter,alg.verbose+1)==0 && @info "$(iter) $(err) $(new_norm)" - old_norm = new_norm iter += 1 @@ -216,43 +275,48 @@ function left_move(peps::InfinitePEPS{PType},alg::CTMRG2,envs::CTMRGEnv) where P @diffset corners[SOUTHWEST,:,col+1]./=norm.(corners[SOUTHWEST,:,col+1]); @diffset edges[WEST,:,col+1]./=norm.(edges[WEST,:,col+1]); end - + CTMRGEnv(peps,corners,edges); end =# -function contract_ctrmg(envs::CTMRGEnv,peps_above = envs.peps_above, peps_below = envs.peps_below) - total = 1.0+0im; - - for r in 1:size(peps_above,1), c in 1:size(peps_above,2) - total*=@tensor envs.edges[WEST,r,c][1 2 3;4]* - envs.corners[NORTHWEST,r,c][4;5]* - envs.edges[NORTH,r,c][5 6 7;8]* - envs.corners[NORTHEAST,r,c][8;9]* - envs.edges[EAST,r,c][9 10 11;12]* - envs.corners[SOUTHEAST,r,c][12;13]* - envs.edges[SOUTH,r,c][13 14 15;16]* - envs.corners[SOUTHWEST,r,c][16;1]* - peps_above[r,c][17;6 10 14 2]* - conj(peps_below[r,c][17;7 11 15 3]) - total *= tr(envs.corners[NORTHWEST,r,c]*envs.corners[NORTHEAST,r,mod1(c-1,end)]*envs.corners[SOUTHEAST,mod1(r-1,end),mod1(c-1,end)]*envs.corners[SOUTHWEST,mod1(r-1,end),c]) - - total /= @tensor envs.edges[WEST,r,c][1 10 11;4]* - envs.corners[NORTHWEST,r,c][4;5]* - envs.corners[NORTHEAST,r,mod1(c-1,end)][5;6]* - envs.edges[EAST,r,mod1(c-1,end)][6 10 11;7]* - envs.corners[SOUTHEAST,r,mod1(c-1,end)][7;8]* - envs.corners[SOUTHWEST,r,c][8;1] - - total /= @tensor envs.corners[NORTHWEST,r,c][1;2]* - envs.edges[NORTH,r,c][2 10 11;3]* - envs.corners[NORTHEAST,r,c][3;4]* - envs.corners[SOUTHEAST,mod1(r-1,end),c][4;5]* - envs.edges[SOUTH,mod1(r-1,end),c][5 10 11;6]* - envs.corners[SOUTHWEST,mod1(r-1,end),c][6;1] - +function contract_ctrmg( + envs::CTMRGEnv, peps_above=envs.peps_above, peps_below=envs.peps_below +) + total = 1.0 + 0im + + for r in 1:size(peps_above, 1), c in 1:size(peps_above, 2) + total *= @tensor envs.edges[WEST, r, c][1 2 3; 4] * + envs.corners[NORTHWEST, r, c][4; 5] * + envs.edges[NORTH, r, c][5 6 7; 8] * + envs.corners[NORTHEAST, r, c][8; 9] * + envs.edges[EAST, r, c][9 10 11; 12] * + envs.corners[SOUTHEAST, r, c][12; 13] * + envs.edges[SOUTH, r, c][13 14 15; 16] * + envs.corners[SOUTHWEST, r, c][16; 1] * + peps_above[r, c][17; 6 10 14 2] * + conj(peps_below[r, c][17; 7 11 15 3]) + total *= tr( + envs.corners[NORTHWEST, r, c] * + envs.corners[NORTHEAST, r, mod1(c - 1, end)] * + envs.corners[SOUTHEAST, mod1(r - 1, end), mod1(c - 1, end)] * + envs.corners[SOUTHWEST, mod1(r - 1, end), c], + ) + + total /= @tensor envs.edges[WEST, r, c][1 10 11; 4] * + envs.corners[NORTHWEST, r, c][4; 5] * + envs.corners[NORTHEAST, r, mod1(c - 1, end)][5; 6] * + envs.edges[EAST, r, mod1(c - 1, end)][6 10 11; 7] * + envs.corners[SOUTHEAST, r, mod1(c - 1, end)][7; 8] * + envs.corners[SOUTHWEST, r, c][8; 1] + + total /= @tensor envs.corners[NORTHWEST, r, c][1; 2] * + envs.edges[NORTH, r, c][2 10 11; 3] * + envs.corners[NORTHEAST, r, c][3; 4] * + envs.corners[SOUTHEAST, mod1(r - 1, end), c][4; 5] * + envs.edges[SOUTH, mod1(r - 1, end), c][5 10 11; 6] * + envs.corners[SOUTHWEST, mod1(r - 1, end), c][6; 1] end - - total -end + return total +end diff --git a/src/algorithms/expval.jl b/src/algorithms/expval.jl index 175f5ba1..99a7b1a5 100644 --- a/src/algorithms/expval.jl +++ b/src/algorithms/expval.jl @@ -1,34 +1,34 @@ -function MPSKit.expectation_value(envs::CTMRGEnv,ham::AbstractTensorMap{S,1,1}) where S +function MPSKit.expectation_value(envs::CTMRGEnv, ham::AbstractTensorMap{S,1,1}) where {S} @assert envs.peps_above == envs.peps_below - peps = envs.peps_above; - result = Matrix{eltype(ham)}(undef,size(peps,1),size(peps,2)); + peps = envs.peps_above + result = Matrix{eltype(ham)}(undef, size(peps, 1), size(peps, 2)) - for r in 1:size(peps,1), c in 1:size(peps,2) - e = @tensor envs.edges[WEST,r,c][1 2 3;4]* - envs.corners[NORTHWEST,r,c][4;5]* - envs.edges[NORTH,r,c][5 6 7;8]* - envs.corners[NORTHEAST,r,c][8;9]* - envs.edges[EAST,r,c][9 10 11;12]* - envs.corners[SOUTHEAST,r,c][12;13]* - envs.edges[SOUTH,r,c][13 14 15;16]* - envs.corners[SOUTHWEST,r,c][16;1]* - peps[r,c][17;6 10 14 2]* - conj(peps[r,c][18;7 11 15 3])* - ham[18;17] + for r in 1:size(peps, 1), c in 1:size(peps, 2) + e = @tensor envs.edges[WEST, r, c][1 2 3; 4] * + envs.corners[NORTHWEST, r, c][4; 5] * + envs.edges[NORTH, r, c][5 6 7; 8] * + envs.corners[NORTHEAST, r, c][8; 9] * + envs.edges[EAST, r, c][9 10 11; 12] * + envs.corners[SOUTHEAST, r, c][12; 13] * + envs.edges[SOUTH, r, c][13 14 15; 16] * + envs.corners[SOUTHWEST, r, c][16; 1] * + peps[r, c][17; 6 10 14 2] * + conj(peps[r, c][18; 7 11 15 3]) * + ham[18; 17] - n = @tensor envs.edges[WEST,r,c][1 2 3;4]* - envs.corners[NORTHWEST,r,c][4;5]* - envs.edges[NORTH,r,c][5 6 7;8]* - envs.corners[NORTHEAST,r,c][8;9]* - envs.edges[EAST,r,c][9 10 11;12]* - envs.corners[SOUTHEAST,r,c][12;13]* - envs.edges[SOUTH,r,c][13 14 15;16]* - envs.corners[SOUTHWEST,r,c][16;1]* - peps[r,c][17;6 10 14 2]* - conj(peps[r,c][17;7 11 15 3]) + n = @tensor envs.edges[WEST, r, c][1 2 3; 4] * + envs.corners[NORTHWEST, r, c][4; 5] * + envs.edges[NORTH, r, c][5 6 7; 8] * + envs.corners[NORTHEAST, r, c][8; 9] * + envs.edges[EAST, r, c][9 10 11; 12] * + envs.corners[SOUTHEAST, r, c][12; 13] * + envs.edges[SOUTH, r, c][13 14 15; 16] * + envs.corners[SOUTHWEST, r, c][16; 1] * + peps[r, c][17; 6 10 14 2] * + conj(peps[r, c][17; 7 11 15 3]) - @diffset result[r,c] = e/n; + @diffset result[r, c] = e / n end - - result + + return result end diff --git a/src/environments/ctmrgenv.jl b/src/environments/ctmrgenv.jl index fd7f7155..327273a3 100644 --- a/src/environments/ctmrgenv.jl +++ b/src/environments/ctmrgenv.jl @@ -6,42 +6,47 @@ struct CTMRGEnv{P,C,T} end # initialize ctmrg environments with some random tensors -CTMRGEnv(peps::InfinitePEPS) = CTMRGEnv(peps,peps); +CTMRGEnv(peps::InfinitePEPS) = CTMRGEnv(peps, peps); -function CTMRGEnv(peps_above::InfinitePEPS{P},peps_below::InfinitePEPS{P}) where P - ou = oneunit(space(peps_above,1,1)) # the bogus space +function CTMRGEnv(peps_above::InfinitePEPS{P}, peps_below::InfinitePEPS{P}) where {P} + ou = oneunit(space(peps_above, 1, 1)) # the bogus space - C_type = tensormaptype(spacetype(P),1,1,storagetype(P)) - T_type = tensormaptype(spacetype(P),3,1,storagetype(P)) # debatable how we should do the legs? + C_type = tensormaptype(spacetype(P), 1, 1, storagetype(P)) + T_type = tensormaptype(spacetype(P), 3, 1, storagetype(P)) # debatable how we should do the legs? #first index is de - corners = Array{C_type}(undef,4,size(peps_above)...) - edges = Array{T_type}(undef,4,size(peps_above)...) - - for dir in 1:4, i in 1:size(peps_above,1),j in 1:size(peps_above,2) - @diffset corners[dir,i,j] = TensorMap(randn,eltype(P),ou,ou) - @diffset edges[dir,i,j] = TensorMap(randn,eltype(P),ou*space(peps_above[i,j],dir+1)'*space(peps_below[i,j],dir+1),ou) + corners = Array{C_type}(undef, 4, size(peps_above)...) + edges = Array{T_type}(undef, 4, size(peps_above)...) + + for dir in 1:4, i in 1:size(peps_above, 1), j in 1:size(peps_above, 2) + @diffset corners[dir, i, j] = TensorMap(randn, eltype(P), ou, ou) + @diffset edges[dir, i, j] = TensorMap( + randn, + eltype(P), + ou * space(peps_above[i, j], dir + 1)' * space(peps_below[i, j], dir + 1), + ou, + ) end - @diffset corners[:,:,:]./=norm.(corners[:,:,:]) - @diffset edges[:,:,:]./=norm.(edges[:,:,:]) + @diffset corners[:, :, :] ./= norm.(corners[:, :, :]) + @diffset edges[:, :, :] ./= norm.(edges[:, :, :]) - CTMRGEnv(peps_above,peps_below,corners,edges) + return CTMRGEnv(peps_above, peps_below, corners, edges) end function Base.rotl90(envs::CTMRGEnv{P,C,T}) where {P,C,T} n_peps_above = rotl90(envs.peps_above) n_peps_below = rotl90(envs.peps_below) - n_corners = Array{C,3}(undef,size(envs.corners)...) - n_edges = Array{T,3}(undef,size(envs.edges)...) + n_corners = Array{C,3}(undef, size(envs.corners)...) + n_edges = Array{T,3}(undef, size(envs.edges)...) for dir in 1:4 - dirm = mod1(dir-1,4) - @diffset n_corners[dirm,:,:] .= rotl90(envs.corners[dir,:,:]) - @diffset n_edges[dirm,:,:] .= rotl90(envs.edges[dir,:,:]) + dirm = mod1(dir - 1, 4) + @diffset n_corners[dirm, :, :] .= rotl90(envs.corners[dir, :, :]) + @diffset n_edges[dirm, :, :] .= rotl90(envs.edges[dir, :, :]) end - CTMRGEnv(n_peps_above,n_peps_below,n_corners,n_edges) + return CTMRGEnv(n_peps_above, n_peps_below, n_corners, n_edges) end Base.eltype(envs::CTMRGEnv) = eltype(envs.corners[1]) diff --git a/src/mpskit_glue/transferpeps_environments.jl b/src/mpskit_glue/transferpeps_environments.jl index ff6c47a9..fb4c12ad 100644 --- a/src/mpskit_glue/transferpeps_environments.jl +++ b/src/mpskit_glue/transferpeps_environments.jl @@ -1,14 +1,26 @@ -MPSKit.environments(state::InfiniteMPS, O::InfiniteTransferPEPS; kwargs...) = environments(convert(MPSMultiline, state), convert(TransferPEPSMultiline, O); kwargs...) +function MPSKit.environments(state::InfiniteMPS, O::InfiniteTransferPEPS; kwargs...) + return environments( + convert(MPSMultiline, state), convert(TransferPEPSMultiline, O); kwargs... + ) +end import MPSKit.MPSMultiline -function MPSKit.environments(state::MPSMultiline, O::TransferPEPSMultiline; solver=MPSKit.Defaults.eigsolver) +function MPSKit.environments( + state::MPSMultiline, O::TransferPEPSMultiline; solver=MPSKit.Defaults.eigsolver +) (lw, rw) = MPSKit.mixed_fixpoints(state, O, state; solver) return MPSKit.PerMPOInfEnv(nothing, O, state, solver, lw, rw, ReentrantLock()) end -function MPSKit.mixed_fixpoints(above::MPSMultiline, O::TransferPEPSMultiline, below::MPSMultiline, init=gen_init_fps(above, O, below); solver=MPSKit.Defaults.eigsolver) +function MPSKit.mixed_fixpoints( + above::MPSMultiline, + O::TransferPEPSMultiline, + below::MPSMultiline, + init=gen_init_fps(above, O, below); + solver=MPSKit.Defaults.eigsolver, +) T = eltype(above) (numrows, numcols) = size(above) @@ -19,46 +31,54 @@ function MPSKit.mixed_fixpoints(above::MPSMultiline, O::TransferPEPSMultiline, b lefties = PeriodicArray{envtype,2}(undef, numrows, numcols) righties = PeriodicArray{envtype,2}(undef, numrows, numcols) - @threads for cr = 1:numrows + @threads for cr in 1:numrows c_above = above[cr] - c_below = below[cr+1] - + c_below = below[cr + 1] + (L0, R0) = init[cr] @sync begin Threads.@spawn begin E_LL = TransferMatrix($c_above.AL, $O[cr], $c_below.AL) (_, Ls, convhist) = eigsolve(flip(E_LL), $L0, 1, :LM, $solver) - convhist.converged < 1 && @info "left eigenvalue failed to converge $(convhist.normres)" + convhist.converged < 1 && + @info "left eigenvalue failed to converge $(convhist.normres)" L0 = first(Ls) end - + Threads.@spawn begin E_RR = TransferMatrix($c_above.AR, $O[cr], $c_below.AR) (_, Rs, convhist) = eigsolve(E_RR, $R0, 1, :LM, $solver) - convhist.converged < 1 && @info "right eigenvalue failed to converge $(convhist.normres)" + convhist.converged < 1 && + @info "right eigenvalue failed to converge $(convhist.normres)" R0 = first(Rs) end end lefties[cr, 1] = L0 for loc in 2:numcols - lefties[cr, loc] = lefties[cr, loc-1] * TransferMatrix(c_above.AL[loc-1], O[cr, loc-1], c_below.AL[loc-1]) + lefties[cr, loc] = + lefties[cr, loc - 1] * + TransferMatrix(c_above.AL[loc - 1], O[cr, loc - 1], c_below.AL[loc - 1]) end - renormfact::eltype(T) = dot(c_below.CR[0], PEPS_∂∂C(L0, R0) * c_above.CR[0]) + renormfact::scalartype(T) = dot(c_below.CR[0], PEPS_∂∂C(L0, R0) * c_above.CR[0]) righties[cr, end] = R0 / sqrt(renormfact) lefties[cr, 1] /= sqrt(renormfact) - for loc in numcols-1:-1:1 - righties[cr, loc] = TransferMatrix(c_above.AR[loc+1], O[cr, loc+1], c_below.AR[loc+1]) * righties[cr, loc+1] + for loc in (numcols - 1):-1:1 + righties[cr, loc] = + TransferMatrix(c_above.AR[loc + 1], O[cr, loc + 1], c_below.AR[loc + 1]) * + righties[cr, loc + 1] - renormfact = dot(c_below.CR[loc], PEPS_∂∂C(lefties[cr, loc+1], righties[cr, loc]) * c_above.CR[loc]) + renormfact = dot( + c_below.CR[loc], + PEPS_∂∂C(lefties[cr, loc + 1], righties[cr, loc]) * c_above.CR[loc], + ) righties[cr, loc] /= sqrt(renormfact) - lefties[cr, loc+1] /= sqrt(renormfact) + lefties[cr, loc + 1] /= sqrt(renormfact) end - end return (lefties, righties) @@ -68,8 +88,22 @@ function gen_init_fps(above::MPSMultiline, O::TransferPEPSMultiline, below::MPSM T = eltype(above) map(1:size(O, 1)) do cr - L0::T = TensorMap(rand, eltype(T), left_virtualspace(below, cr + 1, 0) * space(O[cr].top[1], 5)' * space(O[cr].bot[1], 5), left_virtualspace(above, cr, 0)) - R0::T = TensorMap(rand, eltype(T), right_virtualspace(above, cr, 0) * space(O[cr].top[1], 3)' * space(O[cr].bot[1], 3), right_virtualspace(below, cr + 1, 0)) + L0::T = TensorMap( + rand, + scalartype(T), + left_virtualspace(below, cr + 1, 0) * + space(O[cr].top[1], 5)' * + space(O[cr].bot[1], 5), + left_virtualspace(above, cr, 0), + ) + R0::T = TensorMap( + rand, + scalartype(T), + right_virtualspace(above, cr, 0) * + space(O[cr].top[1], 3)' * + space(O[cr].bot[1], 3), + right_virtualspace(below, cr + 1, 0), + ) (L0, R0) end -end \ No newline at end of file +end diff --git a/src/operators/derivatives.jl b/src/operators/derivatives.jl index 838c7c45..8442178a 100644 --- a/src/operators/derivatives.jl +++ b/src/operators/derivatives.jl @@ -22,37 +22,65 @@ struct PEPS_∂∂AC{T,O} GR::T end -function MPSKit.∂C(C::MPSBondTensor{S}, GL::GenericMPSTensor{S,3}, GR::GenericMPSTensor{S,3}) where {S} +function MPSKit.∂C( + C::MPSBondTensor{S}, GL::GenericMPSTensor{S,3}, GR::GenericMPSTensor{S,3} +) where {S} return @tensor C′[-1; -2] := GL[-1 3 4; 1] * C[1; 2] * GR[2 3 4; -2] end -function MPSKit.∂AC(AC::GenericMPSTensor{S,3}, O::NTuple{2,T}, GL::GenericMPSTensor{S,3}, GR::GenericMPSTensor{S,3}) where {S, T<:PEPSTensor} - return @tensor AC′[-1 -2 -3; -4] := GL[-1 8 9; 7] * AC[7 4 2; 1] * GR[1 6 3; -4] * O[1][5; 4 6 -2 8] * conj(O[2][5; 2 3 -3 9]) +function MPSKit.∂AC( + AC::GenericMPSTensor{S,3}, + O::NTuple{2,T}, + GL::GenericMPSTensor{S,3}, + GR::GenericMPSTensor{S,3}, +) where {S,T<:PEPSTensor} + return @tensor AC′[-1 -2 -3; -4] := + GL[-1 8 9; 7] * + AC[7 4 2; 1] * + GR[1 6 3; -4] * + O[1][5; 4 6 -2 8] * + conj(O[2][5; 2 3 -3 9]) end (H::PEPS_∂∂C)(x) = MPSKit.∂C(x, H.GL, H.GR) (H::PEPS_∂∂AC)(x) = MPSKit.∂AC(x, (H.top, H.bot), H.GL, H.GR) function MPSKit.∂AC(x::RecursiveVec, O::Tuple, GL, GR) - return RecursiveVec( + return RecursiveVec( circshift( - map((v, O1, O2, l, r) -> ∂AC(v, (O1, O2), l, r), x.vecs, O[1], O[2], GL, GR), - 1 - ) + map((v, O1, O2, l, r) -> ∂AC(v, (O1, O2), l, r), x.vecs, O[1], O[2], GL, GR), 1 + ), ) end - Base.:*(H::Union{<:PEPS_∂∂AC,<:PEPS_∂∂C}, v) = H(v) - # operator constructors -MPSKit.∂∂C(pos::Int, mps, mpo::InfiniteTransferPEPS, cache) = PEPS_∂∂C(leftenv(cache, pos + 1, mps), rightenv(cache, pos, mps)) -MPSKit.∂∂C(col::Int, mps, mpo::TransferPEPSMultiline, cache) = PEPS_∂∂C(leftenv(cache, col + 1, mps), rightenv(cache, col, mps)) -MPSKit.∂∂C(row::Int, col::Int, mps, mpo::TransferPEPSMultiline, cache) = PEPS_∂∂C(leftenv(cache, row, col + 1, mps), rightenv(cache, row, col, mps)) +function MPSKit.∂∂C(pos::Int, mps, mpo::InfiniteTransferPEPS, cache) + return PEPS_∂∂C(leftenv(cache, pos + 1, mps), rightenv(cache, pos, mps)) +end +function MPSKit.∂∂C(col::Int, mps, mpo::TransferPEPSMultiline, cache) + return PEPS_∂∂C(leftenv(cache, col + 1, mps), rightenv(cache, col, mps)) +end +function MPSKit.∂∂C(row::Int, col::Int, mps, mpo::TransferPEPSMultiline, cache) + return PEPS_∂∂C(leftenv(cache, row, col + 1, mps), rightenv(cache, row, col, mps)) +end -MPSKit.∂∂AC(pos::Int, mps, mpo::InfiniteTransferPEPS, cache) = PEPS_∂∂AC(mpo.top[pos], mpo.bot[pos], lefenv(cache, pos, mps), rightenv(cache, pos, mps)) -MPSKit.∂∂AC(row::Int, col::Int, mps, mpo::TransferPEPSMultiline, cache) = PEPS_∂∂AC(mpo[row,col]..., leftenv(cache, row, col, mps), rightenv(cache, row, col, mps)) +function MPSKit.∂∂AC(pos::Int, mps, mpo::InfiniteTransferPEPS, cache) + return PEPS_∂∂AC( + mpo.top[pos], mpo.bot[pos], lefenv(cache, pos, mps), rightenv(cache, pos, mps) + ) +end +function MPSKit.∂∂AC(row::Int, col::Int, mps, mpo::TransferPEPSMultiline, cache) + return PEPS_∂∂AC( + mpo[row, col]..., leftenv(cache, row, col, mps), rightenv(cache, row, col, mps) + ) +end function MPSKit.∂∂AC(col::Int, mps, mpo::TransferPEPSMultiline, cache) - return PEPS_∂∂AC(first.(mpo[:, col]), last.(mpo[:, col]), leftenv(cache, col, mps), rightenv(cache, col, mps)) -end \ No newline at end of file + return PEPS_∂∂AC( + first.(mpo[:, col]), + last.(mpo[:, col]), + leftenv(cache, col, mps), + rightenv(cache, col, mps), + ) +end diff --git a/src/operators/periodicpepo.jl b/src/operators/periodicpepo.jl index 153946c3..b0d47b3e 100644 --- a/src/operators/periodicpepo.jl +++ b/src/operators/periodicpepo.jl @@ -10,18 +10,17 @@ struct PeriodicPEPO{T<:PEPOTensor} <: AbstractPEPO Ivertical = CartesianIndex(-1, 0) Ihorizontal = CartesianIndex(0, 1) for I in CartesianIndices(A) - space(A[I], 3) == space(A[I+Ivertical], 5)' || throw(SpaceMismatch( - "North virtual space at site $(Tuple(I)) does not match." - )) - space(A[I], 4) == space(A[I+Ihorizontal], 6)' || throw(SpaceMismatch( - "East virtual space at site $(Tuple(I)) does not match." - )) + space(A[I], 3) == space(A[I + Ivertical], 5)' || throw( + SpaceMismatch("North virtual space at site $(Tuple(I)) does not match.") + ) + space(A[I], 4) == space(A[I + Ihorizontal], 6)' || throw( + SpaceMismatch("East virtual space at site $(Tuple(I)) does not match.") + ) end return new{T}(A) end end - ## Constructors """ InfinitePEPO(A::AbstractArray{T, 2}) @@ -40,7 +39,7 @@ Allow users to pass in arrays of spaces. function PeriodicPEPO( Pspaces::AbstractArray{S,2}, Nspaces::AbstractArray{S,2}, - Espaces::AbstractArray{S,2}=Nspaces + Espaces::AbstractArray{S,2}=Nspaces, ) where {S<:ElementarySpace} size(Pspaces) == size(Nspaces) == size(Espaces) || throw(ArgumentError("Input spaces should have equal sizes.")) @@ -81,7 +80,7 @@ function PeriodicPEPO(d::Integer, D::Integer) return PeriodicPEPO(PeriodicArray(reshape(T, (1, 1)))) end -PeriodicPEPO(d::Integer, D::Integer, Ls::Tuple{Integer}) = repeat(PeriodicPEPO(d,D),Ls...) +PeriodicPEPO(d::Integer, D::Integer, Ls::Tuple{Integer}) = repeat(PeriodicPEPO(d, D), Ls...) ## Shape and size Base.size(T::PeriodicPEPO) = size(T.A) @@ -96,5 +95,4 @@ Base.repeat(T::PeriodicPEPO, counts...) = PeriodicPEPO(repeat(T.A, counts...)) Base.getindex(T::PeriodicPEPO, args...) = getindex(T.A, args...); TensorKit.space(t::PeriodicPEPO, i, j) = space(t[i, j], 1) - -Base.rotl90(t::PeriodicPEPO) = PeriodicPEPO(rotl90(rotl90.(t.A))); \ No newline at end of file +Base.rotl90(t::PeriodicPEPO) = PeriodicPEPO(rotl90(rotl90.(t.A))); diff --git a/src/operators/transferpeps.jl b/src/operators/transferpeps.jl index f41f3505..e40b8365 100644 --- a/src/operators/transferpeps.jl +++ b/src/operators/transferpeps.jl @@ -1,5 +1,4 @@ - struct InfiniteTransferPEPS{T} top::PeriodicArray{T,1} bot::PeriodicArray{T,1} @@ -20,29 +19,25 @@ Base.getindex(O::InfiniteTransferPEPS, i) = (O.top[i], O.bot[i]) Base.iterate(O::InfiniteTransferPEPS, i=1) = i > length(O) ? nothing : (O[i], i + 1) function initializeMPS(O::InfiniteTransferPEPS, virtualspaces::AbstractArray{S,1}) where {S} - return InfiniteMPS( - [ - TensorMap( - rand, - MPSKit.Defaults.eltype, - virtualspaces[mod1(i-1, end)] * space(O.top[i], 2)' * space(O.bot[i], 2), - virtualspaces[i] - ) for i in 1:length(O) - ] - ) + return InfiniteMPS([ + TensorMap( + rand, + MPSKit.Defaults.eltype, + virtualspaces[mod1(i - 1, end)] * space(O.top[i], 2)' * space(O.bot[i], 2), + virtualspaces[i], + ) for i in 1:length(O) + ]) end function initializeMPS(O::InfiniteTransferPEPS, χ::Int) - return InfiniteMPS( - [ - TensorMap( - rand, - MPSKit.Defaults.eltype, - ℂ^χ * space(O.top[i], 2)' * space(O.bot[i], 2), - ℂ^χ - ) for i in 1:length(O) - ] - ) + return InfiniteMPS([ + TensorMap( + rand, + MPSKit.Defaults.eltype, + ℂ^χ * space(O.top[i], 2)' * space(O.bot[i], 2), + ℂ^χ, + ) for i in 1:length(O) + ]) end import MPSKit.GenericMPSTensor @@ -52,30 +47,48 @@ Base.convert(::Type{TransferPEPSMultiline}, O::InfiniteTransferPEPS) = MPSKit.Mu Base.getindex(t::TransferPEPSMultiline, i::Colon, j::Int) = Base.getindex.(t.data[i], j) Base.getindex(t::TransferPEPSMultiline, i::Int, j) = Base.getindex(t.data[i], j) - -function MPSKit.transfer_left(GL::GenericMPSTensor{S,3}, O::NTuple{2, PEPSTensor}, A::GenericMPSTensor{S,3}, Ā::GenericMPSTensor{S,3}) where S - return @tensor GL′[-1 -2 -3; -4] := GL[1 4 2; 7] * conj(Ā[1 3 6; -1]) * - O[1][5; 8 -2 3 4] * conj(O[2][5; 9 -3 6 2]) * A[7 8 9; -4] +function MPSKit.transfer_left( + GL::GenericMPSTensor{S,3}, + O::NTuple{2,PEPSTensor}, + A::GenericMPSTensor{S,3}, + Ā::GenericMPSTensor{S,3}, +) where {S} + return @tensor GL′[-1 -2 -3; -4] := + GL[1 4 2; 7] * + conj(Ā[1 3 6; -1]) * + O[1][5; 8 -2 3 4] * + conj(O[2][5; 9 -3 6 2]) * + A[7 8 9; -4] end -function MPSKit.transfer_right(GR::GenericMPSTensor{S,3}, O::NTuple{2, PEPSTensor}, A::GenericMPSTensor{S,3}, Ā::GenericMPSTensor{S,3}) where S - return @tensor GR′[-1 -2 -3; -4] := GR[7 4 2; 1] * conj(Ā[-4 6 3; 1]) * - O[1][5; 9 4 6 -2] * conj(O[2][5; 8 2 3 -3]) * A[-1 9 8 7] +function MPSKit.transfer_right( + GR::GenericMPSTensor{S,3}, + O::NTuple{2,PEPSTensor}, + A::GenericMPSTensor{S,3}, + Ā::GenericMPSTensor{S,3}, +) where {S} + return @tensor GR′[-1 -2 -3; -4] := + GR[7 4 2; 1] * + conj(Ā[-4 6 3; 1]) * + O[1][5; 9 4 6 -2] * + conj(O[2][5; 8 2 3 -3]) * + A[-1 9 8 7] end function MPSKit.expectation_value(st::MPSMultiline, O::TransferPEPSMultiline) - end -function MPSKit.expectation_value(st::MPSMultiline, ca::MPSKit.PerMPOInfEnv{H,V,S,A}) where {H<:TransferPEPSMultiline,V,S,A} +function MPSKit.expectation_value( + st::MPSMultiline, ca::MPSKit.PerMPOInfEnv{H,V,S,A} +) where {H<:TransferPEPSMultiline,V,S,A} opp = ca.opp retval = PeriodicArray{eltype(st.AC[1, 1]),2}(undef, size(st, 1), size(st, 2)) for (i, j) in product(1:size(st, 1), 1:size(st, 2)) retval[i, j] = @plansor leftenv(ca, i, j, st)[1 2; 3] * - opp[i, j][2 4; 6 5] * - st.AC[i, j][3 6; 7] * - rightenv(ca, i, j, st)[7 5; 8] * - conj(st.AC[i+1, j][1 4; 8]) + opp[i, j][2 4; 6 5] * + st.AC[i, j][3 6; 7] * + rightenv(ca, i, j, st)[7 5; 8] * + conj(st.AC[i + 1, j][1 4; 8]) end return retval end diff --git a/src/states/abstractpeps.jl b/src/states/abstractpeps.jl index 12e5a762..c74ab89f 100644 --- a/src/states/abstractpeps.jl +++ b/src/states/abstractpeps.jl @@ -10,8 +10,7 @@ Default type for PEPS tensors with a single physical index, and 4 virtual indices, conventionally ordered as: T : P ← N ⊗ E ⊗ S ⊗ W. """ -const PEPSTensor{S} = AbstractTensorMap{S, 1, 4} where {S<:ElementarySpace} - +const PEPSTensor{S} = AbstractTensorMap{S,1,4} where {S<:ElementarySpace} """ const PEPOTensor{S} @@ -19,7 +18,7 @@ const PEPSTensor{S} = AbstractTensorMap{S, 1, 4} where {S<:ElementarySpace} Default type for PEPO tensors with a single incoming and outgoing physical index, and 4 virtual indices, conventionally ordered as: O : P ⊗ P' ← N ⊗ E ⊗ S ⊗ W. """ -const PEPOTensor{S} = AbstractTensorMap{S, 2, 4} where {S<:ElementarySpace} +const PEPOTensor{S} = AbstractTensorMap{S,2,4} where {S<:ElementarySpace} ########################## ## Abstract state types ## @@ -39,6 +38,5 @@ Abstract supertype for a 2D projected entangled pairs operator. """ abstract type AbstractPEPO end - Base.rotl90(t::PEPSTensor) = permute(t, ((1,), (3, 4, 5, 2))) Base.rotl90(t::PEPOTensor) = permute(t, ((1, 2), (4, 5, 6, 3))); diff --git a/src/states/infinitepeps.jl b/src/states/infinitepeps.jl index 49e50409..a4345baa 100644 --- a/src/states/infinitepeps.jl +++ b/src/states/infinitepeps.jl @@ -11,7 +11,6 @@ struct InfinitePEPS{T<:PEPSTensor} <: AbstractPEPS end end - ## Constructors """ InfinitePEPS(A::AbstractArray{T, 2}) @@ -30,7 +29,7 @@ Allow users to pass in arrays of spaces. function InfinitePEPS( Pspaces::AbstractArray{S,2}, Nspaces::AbstractArray{S,2}, - Espaces::AbstractArray{S,2}=Nspaces + Espaces::AbstractArray{S,2}=Nspaces, ) where {S<:ElementarySpace} size(Pspaces) == size(Nspaces) == size(Espaces) || throw(ArgumentError("Input spaces should have equal sizes.")) @@ -72,7 +71,7 @@ end function InfinitePEPS(d::Integer, D::Integer, L::Integer) T = [TensorMap(rand, ComplexF64, ℂ^d ← ℂ^D ⊗ ℂ^D ⊗ (ℂ^D)' ⊗ (ℂ^D)')] - return InfinitePEPS(Array(repeat(T, L,L))) + return InfinitePEPS(Array(repeat(T, L, L))) end function InfinitePEPS(d::Integer, D::Integer, Ls::Tuple{Integer,Integer}) @@ -80,7 +79,6 @@ function InfinitePEPS(d::Integer, D::Integer, Ls::Tuple{Integer,Integer}) return InfinitePEPS(Array(repeat(T, Ls...))) end - ## Shape and size Base.size(T::InfinitePEPS) = size(T.A) Base.size(T::InfinitePEPS, i) = size(T.A, i) @@ -95,5 +93,4 @@ Base.repeat(T::InfinitePEPS, counts...) = InfinitePEPS(repeat(T.A, counts...)) Base.getindex(T::InfinitePEPS, args...) = getindex(T.A, args...); TensorKit.space(t::InfinitePEPS, i, j) = space(t[i, j], 1) - Base.rotl90(t::InfinitePEPS) = InfinitePEPS(rotl90(rotl90.(t.A))); diff --git a/src/utility/rotations.jl b/src/utility/rotations.jl index df63f838..58d96c3d 100644 --- a/src/utility/rotations.jl +++ b/src/utility/rotations.jl @@ -8,6 +8,6 @@ const NORTHEAST = 2; const SOUTHEAST = 3; const SOUTHWEST = 4; -rotate_north(t,dir) = mod1(dir,4) == NORTH ? t : rotate_north(rotl90(t),dir-1) +rotate_north(t, dir) = mod1(dir, 4) == NORTH ? t : rotate_north(rotl90(t), dir - 1) #Base.rotl90(t::Tuple) = rotl90.(t); diff --git a/src/utility/util.jl b/src/utility/util.jl index a78e8b16..5ebf8d4e 100644 --- a/src/utility/util.jl +++ b/src/utility/util.jl @@ -1,109 +1,121 @@ function sdiag_inv_sqrt(S::AbstractTensorMap) - toret = similar(S); - + toret = similar(S) + if sectortype(S) == Trivial - copyto!(toret.data,LinearAlgebra.diagm(LinearAlgebra.diag(S.data).^(-1/2))); + copyto!(toret.data, LinearAlgebra.diagm(LinearAlgebra.diag(S.data) .^ (-1 / 2))) else - for (k,b) in blocks(S) - copyto!(blocks(toret)[k],LinearAlgebra.diagm(LinearAlgebra.diag(b).^(-1/2))); + for (k, b) in blocks(S) + copyto!( + blocks(toret)[k], LinearAlgebra.diagm(LinearAlgebra.diag(b) .^ (-1 / 2)) + ) end end - - toret + return toret end -function ChainRulesCore.rrule(::typeof(sdiag_inv_sqrt),S::AbstractTensorMap) - toret = sdiag_inv_sqrt(S); - toret,c̄ -> (ChainRulesCore.NoTangent(),-1/2*_elementwise_mult(c̄,toret'^3)) +function ChainRulesCore.rrule(::typeof(sdiag_inv_sqrt), S::AbstractTensorMap) + toret = sdiag_inv_sqrt(S) + return toret, + c̄ -> (ChainRulesCore.NoTangent(), -1 / 2 * _elementwise_mult(c̄, toret'^3)) end -function _elementwise_mult(a::AbstractTensorMap,b::AbstractTensorMap) - dst = similar(a); - for (k,block) in blocks(dst) - copyto!(block,blocks(a)[k].*blocks(b)[k]); +function _elementwise_mult(a::AbstractTensorMap, b::AbstractTensorMap) + dst = similar(a) + for (k, block) in blocks(dst) + copyto!(block, blocks(a)[k] .* blocks(b)[k]) end - dst + return dst end #rotl90 appeared to lose PeriodicArray'ness, which in turn caused zygote problems #Base.rotl90(a::Array) = Array(rotl90(a)); #Base.rotr90(a::Array) = Array(rotr90(a)); -function ChainRulesCore.rrule(::typeof(rotl90),a::AbstractMatrix) +function ChainRulesCore.rrule(::typeof(rotl90), a::AbstractMatrix) function pb(x) if !iszero(x) - x = x isa Tangent ? ChainRulesCore.construct(typeof(a),ChainRulesCore.backing(x)) : x; - x = rotr90(x); + x = if x isa Tangent + ChainRulesCore.construct(typeof(a), ChainRulesCore.backing(x)) + else + x + end + x = rotr90(x) end - (ZeroTangent(),x) + return (ZeroTangent(), x) end - rotl90(a), pb + return rotl90(a), pb end -structure(t) = codomain(t)←domain(t); +structure(t) = codomain(t) ← domain(t); -function _setindex(a::AbstractArray,v,args...) - b::typeof(a) = copy(a); +function _setindex(a::AbstractArray, v, args...) + b::typeof(a) = copy(a) b[args...] = v - b + return b end -function ChainRulesCore.rrule(::typeof(_setindex), a::AbstractArray, tv, args...) - t = _setindex(a,tv,args...); - +function ChainRulesCore.rrule(::typeof(_setindex), a::AbstractArray, tv, args...) + t = _setindex(a, tv, args...) + function toret(v) if iszero(v) - backwards_tv = ZeroTangent(); - backwards_a = ZeroTangent(); + backwards_tv = ZeroTangent() + backwards_a = ZeroTangent() else - v = v isa Tangent ? ChainRulesCore.construct(typeof(a),ChainRulesCore.backing(v)) : v; - v = typeof(v) != typeof(a) ? convert(typeof(a),v) : v + v = if v isa Tangent + ChainRulesCore.construct(typeof(a), ChainRulesCore.backing(v)) + else + v + end + v = typeof(v) != typeof(a) ? convert(typeof(a), v) : v #v = convert(typeof(a),v); - backwards_tv = v[args...]; - backwards_a = copy(v); + backwards_tv = v[args...] + backwards_a = copy(v) if typeof(backwards_tv) == eltype(a) backwards_a[args...] = zero(v[args...]) else backwards_a[args...] = zero.(v[args...]) end end - (NoTangent(),backwards_a,backwards_tv,fill(ZeroTangent(),length(args))...) + return ( + NoTangent(), backwards_a, backwards_tv, fill(ZeroTangent(), length(args))... + ) end - t,toret + return t, toret end macro diffset(ex) - esc(parse_ex(ex)); + return esc(parse_ex(ex)) end parse_ex(ex) = ex function parse_ex(ex::Expr) - oppheads = (:(./=),:(.*=),:(.+=),:(.-=)); - opprep = (:(./),:(.*),:(.+),:(.-)); + oppheads = (:(./=), :(.*=), :(.+=), :(.-=)) + opprep = (:(./), :(.*), :(.+), :(.-)) if ex.head == :macrocall - parse_ex(macroexpand(PEPSKit,ex)) - elseif ex.head in ( :(.=), :(=)) && length(ex.args)==2 && is_indexing(ex.args[1]) - lhs = ex.args[1]; - rhs = ex.args[2]; + parse_ex(macroexpand(PEPSKit, ex)) + elseif ex.head in (:(.=), :(=)) && length(ex.args) == 2 && is_indexing(ex.args[1]) + lhs = ex.args[1] + rhs = ex.args[2] - vname = lhs.args[1]; - args = lhs.args[2:end]; + vname = lhs.args[1] + args = lhs.args[2:end] quote - $vname = _setindex($vname,$rhs,$(args...)) + $vname = _setindex($vname, $rhs, $(args...)) end - elseif ex.head in oppheads && length(ex.args)==2 && is_indexing(ex.args[1]) - hit = findfirst(x->x==ex.head,oppheads); - rep = opprep[hit]; + elseif ex.head in oppheads && length(ex.args) == 2 && is_indexing(ex.args[1]) + hit = findfirst(x -> x == ex.head, oppheads) + rep = opprep[hit] + + lhs = ex.args[1] + rhs = ex.args[2] - lhs = ex.args[1]; - rhs = ex.args[2]; + vname = lhs.args[1] + args = lhs.args[2:end] - vname = lhs.args[1]; - args = lhs.args[2:end]; - quote - $vname = _setindex($vname,$(rep)($lhs,$rhs),$(args...)) + $vname = _setindex($vname, $(rep)($lhs, $rhs), $(args...)) end else - return Expr(ex.head,parse_ex.(ex.args)...); + return Expr(ex.head, parse_ex.(ex.args)...) end end diff --git a/test/Project.toml b/test/Project.toml index 1d272b8e..b9240c28 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,6 +1,5 @@ [deps] +MPSKit = "bb1c41ca-d63c-52ed-829e-0820dda26502" OptimKit = "77e91f04-9b3b-57a6-a776-40b61faaebe0" TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" -TestExtras = "5ed8adda-3752-4e41-b88a-e8b09835ee3a" -MPSKit = "bb1c41ca-d63c-52ed-829e-0820dda26502" diff --git a/test/runtests.jl b/test/runtests.jl index b249b216..2c366366 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1 +1,10 @@ -using PEPSKit,MPSKit,TensorKit,Test,TestExtras +using Test, PEPSKit, MPSKit, TensorKit + +@testset "boundary MPS" begin + peps = InfinitePEPS(2, 3) + tpeps = InfiniteTransferPEPS(peps, 1, 1) + + mps = initializeMPS(tpeps, 4) + + mps, _, _ = leading_boundary(mps, tpeps, VUMPS()) +end From b957d92acb217f460a763541de23bae67e6d5a2e Mon Sep 17 00:00:00 2001 From: leburgel Date: Wed, 11 Oct 2023 12:17:18 +0200 Subject: [PATCH 4/9] Somewhat working PEPO code --- Project.toml | 7 +- examples/test_pepo.jl | 123 +++++++++ src/PEPSKit.jl | 22 +- src/algorithms/pepo_opt.jl | 273 +++++++++++++++++++ src/environments/boundarympsenv.jl | 243 +++++++++++++++++ src/mpskit_glue/transferpepo_environments.jl | 103 +++++++ src/operators/derivatives.jl | 233 ++++++++++++++++ src/operators/infinitepepo.jl | 157 +++++++++++ src/operators/periodicpepo.jl | 98 ------- src/operators/transferpepo.jl | 204 ++++++++++++++ src/operators/transferpeps.jl | 67 +++-- src/states/infinitepeps.jl | 49 +++- src/utility/symmetrization.jl | 181 ++++++++++++ test/runtests.jl | 11 +- 14 files changed, 1637 insertions(+), 134 deletions(-) create mode 100644 examples/test_pepo.jl create mode 100644 src/algorithms/pepo_opt.jl create mode 100644 src/environments/boundarympsenv.jl create mode 100644 src/mpskit_glue/transferpepo_environments.jl create mode 100644 src/operators/infinitepepo.jl delete mode 100644 src/operators/periodicpepo.jl create mode 100644 src/operators/transferpepo.jl create mode 100644 src/utility/symmetrization.jl diff --git a/Project.toml b/Project.toml index 49957658..fd27b9ff 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ authors = ["maarten "] version = "0.1.0" [deps] +Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" @@ -13,8 +14,10 @@ OptimKit = "77e91f04-9b3b-57a6-a776-40b61faaebe0" Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec" +VectorInterface = "409d34a3-91d5-4945-b6ec-7529ddf182d8" [compat] -julia = "1.6" -MPSKit = "0.9.1" +MPSKit = "0.10" TensorKit = "0.12" +VectorInterface = '0.4" +julia = "1.6" diff --git a/examples/test_pepo.jl b/examples/test_pepo.jl new file mode 100644 index 00000000..7c818555 --- /dev/null +++ b/examples/test_pepo.jl @@ -0,0 +1,123 @@ +# this should run once the compatocalypse has been resolved + +using Revise +using TensorOperations +using TensorKit +using MPSKit +using PEPSKit +using OptimKit + + +## Set some options + +optim_method = ConjugateGradient +optim_tol = 1e-5 +optim_maxiter = 100 +verbosity = 2 # >= 5 for debugging +boundary_maxiter = 50 +tol_min = 1e-12 +tol_max = 1e-4 +tol_factor = 1e-3 +symm = Full() +hermitian = true +boundary_method = VUMPS(; tol_galerkin=tol_max, maxiter=boundary_maxiter, dynamical_tols=true, eigs_tolfactor=1e-3, envs_tolfactor=1e-3, gauge_tolfactor=1e-6, tol_max=1e-4, verbose=verbosity >= 5) +# boundary_method = GradientGrassmann + + +## Set model parameters + +# 3D Ising temperature +beta = 1/3 # should give f = -1.0219139... @ D = 3 + +# bond dimensions +D = 3 +χ = 30 + +# pick size of unit cell +depth = 1 +width = 1 +height = 1 + + +## PEPO defintion: 3D classical Ising model + +t = ComplexF64[exp(beta) exp(-beta); exp(-beta) exp(beta)] +q = sqrt(t) + +O = zeros(2, 2, 2, 2, 2, 2) +O[1, 1, 1, 1, 1, 1] = 1 +O[2, 2, 2, 2, 2, 2] = 1 +@tensor o[-1 -2; -3 -4 -5 -6] := O[1 2; 3 4 5 6] * q[-1; 1] * q[-2; 2] * q[-3; 3] * q[-4; 4] * q[-5; 5] * q[-6; 6] + +M = zeros(2, 2, 2, 2, 2, 2) +M[1, 1, 1, 1, 1, 1] = 1 +M[2, 2, 2, 2, 2, 2] = -1 +@tensor m[-1 -2; -3 -4 -5 -6] := M[1 2; 3 4 5 6] * q[-1; 1] * q[-2; 2] * q[-3; 3] * q[-4; 4] * q[-5; 5] * q[-6; 6] + +o = TensorMap(o, ℂ^2 ⊗ (ℂ^2)' ← ℂ^2 ⊗ ℂ^2 ⊗ (ℂ^2)' ⊗ (ℂ^2)') +m = TensorMap(m, ℂ^2 ⊗ (ℂ^2)' ← ℂ^2 ⊗ ℂ^2 ⊗ (ℂ^2)' ⊗ (ℂ^2)') + +O = InfinitePEPO(repeat([o], depth, width, height)) + + +## Optimization algorithm + +pepo_alg = PEPOOptimize(; + optim_method, + optim_tol, + optim_maxiter, + verbosity, + boundary_method, + boundary_maxiter, + tol_min, + tol_max, + tol_factor, + symm, + hermitian, +) + + +## Initialize PEPS and MPS fixed points + +peps = symmetrize(initializePEPS(O, ℂ^D), symm) +envs = pepo_opt_environments(peps, O, pepo_alg.boundary_method; vspaces=[ℂ^χ], hermitian=pepo_alg.hermitian) +normalize!(peps, envs.peps_boundary) + + +## Perform optimization + +x, f, normgrad = leading_boundary(peps, O, pepo_alg, envs) + + +## Unpack result and check magnetization + +peps = x.state +above = x.envs.pepo_boundary.boundaries[1] +below = x.envs.pepo_boundary.boundaries[2] +lw, rw = x.envs.pepo_boundary.envs[3] + +pos = (1, 1) # pick a site + +if height == 1 # no crazy contractions for now... + row, col = pos + fliprow = size(peps, 1) - row + 1 # below starts counting from below + m = @tensor lw[row, col][13 8 10 18; 1] * + above.AC[row, col][1 9 11 15; 12] * + peps[row, col][5; 9 3 4 8] * + m[14 5; 11 6 7 10] * + rw[row, col][12 3 6 16; 2] * + conj(below.AC[fliprow, col][13 4 7 17; 2]) * + conj(peps[row, col][14; 15 16 17 18]) + + λ = @tensor lw[row, col][13 8 10 18; 1] * + above.AC[row, col][1 9 11 15; 12] * + peps[row, col][5; 9 3 4 8] * + t[14 5; 11 6 7 10] * + rw[row, col][12 3 6 16; 2] * + conj(below.AC[fliprow, col][13 4 7 17; 2]) * + conj(peps[row, col][14; 15 16 17 18]) + + println("m = $(m / λ)") # should be ~0.945 @ beta = 1/3 and D = 3 +end + +nothing diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index b8f0eb22..b2dcebbe 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -1,4 +1,7 @@ module PEPSKit + +using Accessors +using VectorInterface using TensorKit, KrylovKit, MPSKit, OptimKit, Base.Threads, Base.Iterators, Parameters, Printf using ChainRulesCore @@ -14,16 +17,22 @@ include("states/abstractpeps.jl") include("states/infinitepeps.jl") include("operators/transferpeps.jl") +include("operators/infinitepepo.jl") +include("operators/transferpepo.jl") include("operators/derivatives.jl") -include("operators/periodicpepo.jl") - -include("environments/ctmrgenv.jl") include("mpskit_glue/transferpeps_environments.jl") +include("mpskit_glue/transferpepo_environments.jl") + +include("environments/ctmrgenv.jl") +include("environments/boundarympsenv.jl") include("algorithms/ctmrg.jl") include("algorithms/expval.jl") +include("utility/symmetrization.jl") +include("algorithms/pepo_opt.jl") + include("utility/rotations.jl") #default settings @@ -33,6 +42,9 @@ module Defaults end export InfinitePEPS, InfiniteTransferPEPS -export PeriodicPEPO -export initializeMPS +export InfinitePEPO, InfiniteTransferPEPO +export initializeMPS, initializePEPS +export PEPOOptimize, pepo_opt_environments +export symmetrize, None, Depth, Full + end # module diff --git a/src/algorithms/pepo_opt.jl b/src/algorithms/pepo_opt.jl new file mode 100644 index 00000000..69053cc7 --- /dev/null +++ b/src/algorithms/pepo_opt.jl @@ -0,0 +1,273 @@ +# First go at PEPO fixed point optimization using PEPSKit + +## Characterize PEPO optimization manifold + +# point on optimization manifold contains PEPS, its 'boundary' +# the norm of the gradient is for convenience for now, but should be removed +mutable struct PEPOOptPoint{T,E,F} + state::T + envs::E + normgrad::F +end +PEPOOptPoint(state, envs) = PEPOOptPoint(state, envs, Inf) + +# the gradient is just a 2D array of PEPS tensors, is this sufficient? +function _retract(x::PEPOOptPoint, dx, α::Number) + # move the state, keep the rest? + new_state = copy(x.state) + new_state.A .+= dx .* α + return PEPOOptPoint(new_state, x.envs, x.normgrad), dx +end + +_inner(x, dx1, dx2) = real(dot(dx1, dx2)) + +function _add!(Y, X, a) + return add!(Y, X, a) +end + +function _scale!(η, β) + return scale!(η, β) +end + +## Utility; TODO: remove all of this + +function _nthroot(x::Real, n::Integer) + return if isodd(n) || x ≥ 0 + copysign(abs(x)^(1//n), x) + else + throw( + DomainError( + "Exponentiation yielding a complex result requires a complex argument. Replace nthroot(x, n) with complex(x)^(1//n).", + ), + ) + end +end +function _hacked_root(x, n) + m = abs(x) + p = angle(x) + mY = _nthroot(m, n) + pYs = p / n .+ 2 * pi / n .* (0:(n - 1)) + I = argmin(minimum([abs.(pYs .- p); abs.(2 * pi .- pYs .+ p)]; dims=1)) # find phase with minimal angular distance + return mY * exp(pYs[I] * 1im) +end + +isverbose(alg::VUMPS) = alg.verbose +isverbose(alg::GradientGrassmann) = alg.method.verbosity >= 0 + + +## Characterize environments for PEPO optimization + +mutable struct PEPOOptEnv{T,O,F} <: Cache + peps_boundary::T + pepo_boundary::O + alg::F +end + +# TODO: think about name; just picked something to avoid conflict with MPSKit.environments +function pepo_opt_environments( + peps::InfinitePEPS, + pepo::InfinitePEPO, + boundary_alg; + vspaces=[oneunit(space(peps, 1, 1))], + hermitian=false, + kwargs..., +) + # boundary_alg handles everything + peps_boundary = environments(peps, boundary_alg; vspaces, hermitian) + pepo_boundary = environments(peps, pepo, boundary_alg; vspaces, hermitian) + + return PEPOOptEnv(peps_boundary, pepo_boundary, boundary_alg) +end + +# I did overload recalculate!, this seemed to make sense +function MPSKit.recalculate!( + envs::PEPOOptEnv, + peps::InfinitePEPS, + pepo::InfinitePEPO; + tol=algtol(envs.alg), + hermitian=false, + kwargs..., +) + recalculate!(envs.peps_boundary, peps; tol, hermitian, kwargs...) + recalculate!(envs.pepo_boundary, peps, pepo; tol, hermitian, kwargs...) + + return envs +end + + +## PEPO fixed point optimization algorithm + +# first attempt at a PEPO-fixed-point-optimization algorithm, bit of a mess... +struct PEPOOptimize{A} + optim_method::OptimKit.OptimizationAlgorithm + optim_finalize!::Function + boundary_method::A + tol_min::Float64 # some basic tolerance scaling + tol_max::Float64 + tol_factor::Float64 + symm::SymmetrizationStyle # symmetrization + hermitian::Bool + + function PEPOOptimize(; + optim_method=ConjugateGradient, + (optim_finalize!)=OptimKit._finalize!, + optim_tol=Defaults.tol, + optim_maxiter=Defaults.maxiter, + verbosity=2, + boundary_method=VUMPS, + boundary_maxiter=Defaults.maxiter, + boundary_finalize=MPSKit.Defaults._finalize, + tol_min=1e-12, + tol_max=1e-5, + tol_factor=1e-3, + symm=None(), + hermitian=false, + ) + if isa(optim_method, OptimKit.OptimizationAlgorithm) + # We were given an optimisation method, just use it. + m = optim_method + elseif optim_method <: OptimKit.OptimizationAlgorithm + # We were given an optimisation method type, construct an instance of it. + m = optim_method(; + gradtol=optim_tol, maxiter=optim_maxiter, verbosity=verbosity + ) + else + msg = "optim_method should be either an instance or a subtype of `OptimKit.OptimizationAlgorithm`." + throw(ArgumentError(msg)) + end + if isa(boundary_method, Union{<:VUMPS,CTMRG,GradientGrassmann}) + bm = boundary_method + elseif boundary_method <: Union{Type{CTMRG},MPSKit.Algorithm} + # total syntax clusterfuck, need to clean this up + if boundary_method <: VUMPS + bm = boundary_method(; + tol_galerkin=tol_max, + verbose=verbosity >= 5, + maxiter=boundary_maxiter, + finalize=boundary_finalize, + ) + elseif boundary_method <: GradientGrassmann + bm = boundary_method(; + tol=tol_max, verbosity=verbosity - 4, maxiter=boundary_maxiter + ) + elseif boundary_method <: CTMRG + bm = method(; tol=tol_max, verbose=verbosity >= 5, maxiter=boundary_maxiter) + else + msg = "Unknown boundary contraction method." + throw(ArgumentError(msg)) + end + else + msg = "boundary_method should be a valid boundary contraction algorithm." + throw(ArgumentError(msg)) + end + return new{typeof(bm)}( + m, optim_finalize!, bm, tol_min, tol_max, tol_factor, symm, hermitian + ) + end +end + +# default PEPO optimization cost function for given PEPO and optimization algorithm +function pepo_opt_costfun( + op::InfinitePEPO, alg::PEPOOptimize) + D, W, H = size(op) + nrm = D * W * H + function pepo_opt_fg(x::PEPOOptPoint) + # unpack state + peps = symmetrize(x.state, alg.symm) + envs = x.envs + ng = x.normgrad + + # recompute environment with scaled tolerance + boundary_tol = min(max(alg.tol_min, ng * alg.tol_factor), alg.tol_max) + recalculate!( + envs, + peps, + op; + tol=boundary_tol, + hermitian=alg.hermitian, + ) + + # compute cost function + lambdas_peps = expectation_value(peps, envs.peps_boundary) + lambdas_pepo = expectation_value(peps, op, envs.pepo_boundary) + f = -log(real(prod(lambdas_pepo) / prod(lambdas_peps))) / nrm + + # compute gradient + ∂p_peps = ∂∂peps(peps, envs.peps_boundary) + ∂p_pepo = ∂∂peps(peps, op, envs.pepo_boundary) + grad = - (2 / nrm) .* ∂p_pepo ./ lambdas_pepo .+ (2 / nrm) .* ∂p_peps ./ lambdas_peps + grad = symmetrize(grad, alg.symm) + # TODO: test if gradient is actually correct + + # TODO: decide whether we want to: + # Actually update the state in place after symmetrization? + # Update the environments and norm of gradient after each cost function + # evaluation, or move this to the optim_finalize to only update after line search + # has terminated? + x.state = peps # is this sensible? + x.envs = envs + x.normgrad = sum(norm.(grad)) + + # some temporary debugging verbosity, to be removed + if isverbose(alg.boundary_method) + lambda_peps = prod(lambdas_peps) + lambda_pepo = prod(lambdas_pepo) + + lambda_pepo_root = _hacked_root(lambda_pepo, D * W * H) + rel_im_pepo = abs(imag(lambda_pepo_root) / abs(lambda_pepo_root)) + @printf( + "\n\tCurrent lambda_pepo: %f + %fim;\tRelative im. part: %e;\n", + real(lambda_pepo_root), + imag(lambda_pepo_root), + rel_im_pepo, + ) + + lambda_peps_root = _hacked_root(lambda_peps, D * W) + rel_im_peps = abs(imag(lambda_peps_root) / abs(lambda_peps_root)) + @printf( + "\tCurrent lambda_peps: %f + %fim;\tRelative im. part: %e;\n\n", + real(lambda_peps_root), + imag(lambda_peps_root), + rel_im_peps, + ) + + @printf("\tCurrent f: %f\n\n", f) + + lambdas_peps_bis = map(∂p_peps, peps.A) do ∂p, p + @tensor ∂p[1; 2 3 4 5] * conj(p[1; 2 3 4 5]) + end + lambdas_pepo_bis = map(∂p_pepo, peps.A) do ∂p, p + @tensor ∂p[1; 2 3 4 5] * conj(p[1; 2 3 4 5]) + end + @show diff_peps = maximum(abs.(lambdas_peps .- lambdas_peps_bis)) + @show diff_pepo = maximum(abs.(lambdas_pepo .- lambdas_pepo_bis)) + end + + return f, grad + end + return pepo_opt_fg +end + + +## The actual leading boundary routine + +function MPSKit.leading_boundary( + state::InfinitePEPS, + op::InfinitePEPO, + alg::PEPOOptimize, + envs=pepo_opt_environments(state, op, alg.boundary_method; hermitian=alg.hermitian), + fg=pepo_opt_costfun(op, alg), +) + res = optimize( + fg, + PEPOOptPoint(state, envs), + alg.optim_method; + retract=_retract, + inner=_inner, + (scale!)=_scale!, + (add!)=_add!, + (finalize!)=alg.optim_finalize!, + ) + (x, fx, gx, numfg, normgradhistory) = res + return x, fx, normgradhistory[end] +end diff --git a/src/environments/boundarympsenv.jl b/src/environments/boundarympsenv.jl new file mode 100644 index 00000000..33a861c7 --- /dev/null +++ b/src/environments/boundarympsenv.jl @@ -0,0 +1,243 @@ +# Some form of boundary MPS environments for infinite PEPS and PEPO routines + + +## Utility + +algtol(alg::VUMPS) = alg.tol_galerkin +algtol(alg::GradientGrassmann) = alg.method.gradtol +update_tol(alg::VUMPS, tol) = @set alg.tol_galerkin=tol +function update_tol(alg::GradientGrassmann, tol) # annoying disparity between typedef and actual constructor... + m = alg.method + m = @set m.gradtol=tol + return GradientGrassmann(; method=m, (finalize!)=alg.finalize!) +end + + +## Boundary MPS environment manager + +mutable struct BoundaryMPSEnv{A,E,F} <: Cache + boundaries::A + envs::E + alg::F +end + + +## PEPS boundary MPS + +function MPSKit.environments( + peps::InfinitePEPS, + alg::A=VUMPS(); + vspaces=[oneunit(space(peps, 1, 1))], + hermitian=false, + kwargs..., +) where {A<:Union{VUMPS,GradientGrassmann}} + tr = TransferPEPSMultiline(peps, 1) + return environments(tr, alg; vspaces, hermitian, kwargs...) +end + +function MPSKit.recalculate!( + envs::BoundaryMPSEnv, peps::InfinitePEPS; tol=algtol(envs.alg), hermitian=false, kwargs... +) + tr = TransferPEPSMultiline(peps, 1) + return recalculate!(envs, tr; tol, hermitian, kwargs...) +end + +# this probably needs a different name? +# computes norm-per-site of PEPS-PEPS overlap using above, below and mixed environments +function MPSKit.expectation_value(peps::InfinitePEPS, ca::BoundaryMPSEnv) + retval = PeriodicArray{scalartype(peps),2}(undef, size(peps)...) + above = ca.boundaries[1] + below = ca.boundaries[2] + (lw, rw) = ca.envs[3] + for (row, col) in product(1:size(peps, 1), 1:size(peps, 2)) + fliprow = size(peps, 1) - row + 1 # below starts counting from below + retval[row, col] = @tensor lw[row, col][1 2 4; 7] * + conj(below.AC[fliprow, col][1 3 6; 13]) * + peps[row, col][5; 8 11 3 2] * + conj(peps[row, col][5; 9 12 6 4]) * + above.AC[row, col][7 8 9; 10] * + rw[row, col][10 11 12; 13] + end + return retval +end + +function MPSKit.normalize!(peps::InfinitePEPS, envs::BoundaryMPSEnv) + norm_per_site = expectation_value(peps, envs) + for (row, col) in product(1:size(peps, 1), 1:size(peps, 2)) + scale!(peps[row, col], 1 / sqrt(norm_per_site[row, col])) + end +end + +# naming convention? this does not return an effective operator, so naming is probably misleading +function ∂∂peps(peps::InfinitePEPS{T}, ca::BoundaryMPSEnv) where {T<:PEPSTensor} + retval = PeriodicArray{T,2}(undef, size(peps)...) + above = ca.boundaries[1] + below = ca.boundaries[2] + (lw, rw) = ca.envs[3] + for (row, col) in product(1:size(peps, 1), 1:size(peps, 2)) + fliprow = size(peps, 1) - row + 1 # below starts counting from below + retval[row, col] = ∂peps( + above.AC[row, col], + below.AC[fliprow, col], + peps[row, col], + lw[row, col], + rw[row, col], + ) + end + return retval +end + + +## PEPO boundary MPS + +function MPSKit.environments( + peps::InfinitePEPS, + pepo::InfinitePEPO, + alg::A=VUMPS(); + vspaces=[oneunit(space(peps, 1, 1))], + hermitian=false, + kwargs..., +) where {A<:Union{VUMPS,GradientGrassmann}} + tr = TransferPEPOMultiline(peps, pepo, 1) + return environments(tr, alg; vspaces, hermitian, kwargs...) +end + +function MPSKit.recalculate!( + envs::BoundaryMPSEnv, + peps::InfinitePEPS, + pepo::InfinitePEPO; + tol=algtol(envs.alg), + hermitian=false, + kwargs..., +) + tr = TransferPEPOMultiline(peps, pepo, 1) + return recalculate!(envs, tr; tol, hermitian, kwargs...) +end + +# this probably needs a different name? +# computes norm-per-site of PEPS-PEPO-PEPS sandwich using above, below and mixed environments +function MPSKit.expectation_value( + peps::InfinitePEPS, pepo::InfinitePEPO, ca::BoundaryMPSEnv +) + retval = PeriodicArray{scalartype(peps),2}(undef, size(peps)...) + opp = TransferPEPOMultiline(peps, pepo, 1) # for convenience... + N = height(opp[1]) + 4 + above = ca.boundaries[1] + below = ca.boundaries[2] + (lw, rw) = ca.envs[3] + for (row, col) in product(1:size(peps, 1), 1:size(peps, 2)) + fliprow = size(peps, 1) - row + 1 # below starts counting from below + O_rc = opp[row, col] + GL´ = transfer_left(lw[row, col], O_rc, above.AC[row, col], below.AC[fliprow, col]) + retval[row, col] = TensorKit.TensorOperations.tensorscalar( + ncon([GL´, rw[row, col]], [[N, (2:(N - 1))..., 1], [(1:N)...]]) + ) + end + return retval +end + +# naming convention? this does not return an effective operator, so naming is probably misleading +function ∂∂peps( + peps::InfinitePEPS{T}, pepo::InfinitePEPO, ca::BoundaryMPSEnv +) where {T<:PEPSTensor} + retval = PeriodicArray{T,2}(undef, size(peps)...) + opp = TransferPEPOMultiline(peps, pepo, 1) # just for convenience + above = ca.boundaries[1] + below = ca.boundaries[2] + (lw, rw) = ca.envs[3] + for (row, col) in product(1:size(peps, 1), 1:size(peps, 2)) + fliprow = size(peps, 1) - row + 1 # below starts counting from below + O_rc = opp[row, col] + retval[row, col] = ∂peps( + above.AC[row, col], + below.AC[fliprow, col], + (O_rc[1], O_rc[3]), + lw[row, col], + rw[row, col], + ) + end + return retval +end + + +## The actual routines + +# because below starts counting from below +_flippedy(M::MPSMultiline) = MPSKit.Multiline(circshift(reverse(M.data), 1)) + +function MPSKit.environments( + tr::Union{TransferPEPSMultiline,TransferPEPOMultiline}, + alg::A=VUMPS(); + vspaces=[oneunit(space(tr, 1, 1))], # defaults to trivial + hermitian=false, +) where {A<:Union{VUMPS,GradientGrassmann}} + # above boundary + above = initializeMPS(tr, vspaces) + envs_above = environments(above, tr) + + # maybe below boundary + if hermitian # literally the same + below = above + envs_below = envs_above + envs_mixed = (envs_above.lw, envs_above.rw) + else + error("not implemented yet") + tr_dag = dagger(tr) + below = initializeMPS(tr_dag, vspaces) + envs_below = envirlnments(below, tr_dag) + + # mixed environments + envs_mixed = MPSKit.mixed_fixpoints(above, tr, _flippedy(below)) + end + + # collect + boundaries = [above, below] + envs = [envs_above, envs_below, envs_mixed] + + return recalculate!(BoundaryMPSEnv(boundaries, envs, alg), tr; hermitian) +end + +function MPSKit.recalculate!( + envs::BoundaryMPSEnv, + tr::Union{TransferPEPSMultiline,TransferPEPOMultiline}; + tol=algtol(envs.alg), + hermitian=false, +) + alg = algtol(envs.alg) == tol ? envs.alg : update_tol(envs.alg, tol) + + # above boundary + envs.envs[1].opp = tr # because pre-supplied environments only ever use their own operator... + envs.boundaries[1], envs.envs[1], err = leading_boundary( + envs.boundaries[1], tr, alg, envs.envs[1] + ) + + # maybe below boundary + if hermitian # literally the same + envs.boundaries[2], envs.envs[2] = envs.boundaries[1], envs.envs[1] + envs.envs[3] = (envs.envs[1].lw, envs.envs[1].rw) + else + error("not implemented yet") + tr_dag = dagger(tr) # TODO: actually define this... + envs.envs[2].opp = tr_dag # because pre-supplied environments only ever use their own operator... + envs.boundaries[2], envs.envs[2], err = leading_boundary( + envs.boundaries[2], tr_dag, alg, envs.envs[2] + ) + # mixed environments + envs.envs[3] = MPSKit.mixed_fixpoints( + envs.boundaries[1], tr, _flippedy(envs.boundaries[2]), envs.envs[3] + ) + end + + envs.alg = alg + + return envs +end + +## Channel MPS environments + +# TODO: actually try this again at some point? +mutable struct ChannelMPSEnv{A,C,F} <: Cache + boundaries::A + corners::C + alg::F +end diff --git a/src/mpskit_glue/transferpepo_environments.jl b/src/mpskit_glue/transferpepo_environments.jl new file mode 100644 index 00000000..e6f00161 --- /dev/null +++ b/src/mpskit_glue/transferpepo_environments.jl @@ -0,0 +1,103 @@ + +function MPSKit.environments(state::InfiniteMPS, O::InfiniteTransferPEPO; kwargs...) + return environments( + convert(MPSMultiline, state), convert(TransferPEPOMultiline, O); kwargs... + ) +end + +function MPSKit.environments( + state::MPSMultiline, O::TransferPEPOMultiline; solver=MPSKit.Defaults.eigsolver +) + (lw, rw) = MPSKit.mixed_fixpoints(state, O, state; solver) + return MPSKit.PerMPOInfEnv(nothing, O, state, solver, lw, rw, ReentrantLock()) +end + +function MPSKit.mixed_fixpoints( + above::MPSMultiline, + O::TransferPEPOMultiline, + below::MPSMultiline, + init=gen_init_fps(above, O, below); + solver=MPSKit.Defaults.eigsolver, +) + T = eltype(above) + + (numrows, numcols) = size(above) + @assert size(above) == size(O) + @assert size(below) == size(O) + + envtype = eltype(init[1]) + lefties = PeriodicArray{envtype,2}(undef, numrows, numcols) + righties = PeriodicArray{envtype,2}(undef, numrows, numcols) + + @threads for cr in 1:numrows + c_above = above[cr] + c_below = below[cr + 1] + + (L0, R0) = init[cr] + + @sync begin + Threads.@spawn begin + E_LL = TransferMatrix($c_above.AL, $O[cr], $c_below.AL) + (_, Ls, convhist) = eigsolve(flip(E_LL), $L0, 1, :LM, $solver) + convhist.converged < 1 && + @info "left eigenvalue failed to converge $(convhist.normres)" + L0 = first(Ls) + end + + Threads.@spawn begin + E_RR = TransferMatrix($c_above.AR, $O[cr], $c_below.AR) + (_, Rs, convhist) = eigsolve(E_RR, $R0, 1, :LM, $solver) + convhist.converged < 1 && + @info "right eigenvalue failed to converge $(convhist.normres)" + R0 = first(Rs) + end + end + + lefties[cr, 1] = L0 + for loc in 2:numcols + lefties[cr, loc] = + lefties[cr, loc - 1] * + TransferMatrix(c_above.AL[loc - 1], O[cr, loc - 1], c_below.AL[loc - 1]) + end + + renormfact::scalartype(T) = dot(c_below.CR[0], PEPO_∂∂C(L0, R0) * c_above.CR[0]) + + righties[cr, end] = R0 / sqrt(renormfact) + lefties[cr, 1] /= sqrt(renormfact) + + for loc in (numcols - 1):-1:1 + righties[cr, loc] = + TransferMatrix(c_above.AR[loc + 1], O[cr, loc + 1], c_below.AR[loc + 1]) * + righties[cr, loc + 1] + + renormfact = dot( + c_below.CR[loc], + PEPO_∂∂C(lefties[cr, loc + 1], righties[cr, loc]) * c_above.CR[loc], + ) + righties[cr, loc] /= sqrt(renormfact) + lefties[cr, loc + 1] /= sqrt(renormfact) + end + end + + return (lefties, righties) +end + +function gen_init_fps(above::MPSMultiline, O::TransferPEPOMultiline, below::MPSMultiline) + T = eltype(above) + + map(1:size(O, 1)) do cr + L0::T = TensorMap( + rand, + scalartype(T), + left_virtualspace(below, cr + 1, 0) * prod(adjoint.(west_spaces(O[cr], 1))), + left_virtualspace(above, cr, 0), + ) + R0::T = TensorMap( + rand, + scalartype(T), + right_virtualspace(above, cr, 0) * prod(adjoint.(east_spaces(O[cr], 1))), + right_virtualspace(below, cr + 1, 0), + ) + (L0, R0) + end +end diff --git a/src/operators/derivatives.jl b/src/operators/derivatives.jl index 8442178a..fb82f161 100644 --- a/src/operators/derivatives.jl +++ b/src/operators/derivatives.jl @@ -1,5 +1,8 @@ import MPSKit.GenericMPSTensor, MPSKit.MPSBondTensor +# PEPS +# ---- + """ struct PEPS_∂∂C{T<:GenericMPSTensor{S,N}} @@ -84,3 +87,233 @@ function MPSKit.∂∂AC(col::Int, mps, mpo::TransferPEPSMultiline, cache) rightenv(cache, col, mps), ) end + +# PEPS derivative +function ∂peps( + AC::GenericMPSTensor{S,3}, + ĀC::GenericMPSTensor{S,3}, + O::T, + GL::GenericMPSTensor{S,3}, + GR::GenericMPSTensor{S,3}, +) where {S,T<:PEPSTensor} + return @tensor ∂p[-1; -2 -3 -4 -5] := + GL[8 5 -5; 1] * + AC[1 6 -2; 7] * + O[-1; 6 3 4 5] * + GR[7 3 -3; 2] * + conj(ĀC[8 4 -4; 2]) +end + + +# PEPO +# ---- + +""" + struct PEPO_∂∂C{T<:GenericMPSTensor{S,N}} + +Represents the effective Hamiltonian for the zero-site derivative of an MPS. +""" +struct PEPO_∂∂C{T} + GL::T + GR::T +end + +""" + struct PEPO_∂∂AC{T,O,P} + +Represents the effective Hamiltonian for the one-site derivative of an MPS. +""" +struct PEPO_∂∂AC{T,O,P} + top::O + bot::O + mid::P + GL::T + GR::T +end + +# specialize simple case +function MPSKit.∂C( + C::MPSBondTensor{S}, GL::GenericMPSTensor{S,4}, GR::GenericMPSTensor{S,4} +) where {S} + return @tensor C′[-1; -2] := GL[-1 3 4 5; 1] * C[1; 2] * GR[2 3 4 5; -2] +end + +function MPSKit.∂C( + C::MPSBondTensor{S}, GL::GenericMPSTensor{S,N}, GR::GenericMPSTensor{S,N} +) where {S,N} + C′ = ncon([GL, C, GR], [[-1, ((2:N) .+ 1)..., 1], [1, 2], [2, ((2:N) .+ 1)..., -2]]) + return permute(C′, ((1,), (2,))) +end + +# specialize simple case +function MPSKit.∂AC( + AC::GenericMPSTensor{S,4}, + O::Tuple{T,T,Tuple{P}}, + GL::GenericMPSTensor{S,4}, + GR::GenericMPSTensor{S,4}, +) where {S,T<:PEPSTensor,P<:PEPOTensor} + return @tensor AC′[-1 -2 -3 -4; -5] := + GL[-1 2 4 7; 1] * + AC[1 3 5 8; 10] * + GR[10 11 12 13; -5] * + O[1][6; 3 11 -2 2] * + O[3][1][9 6; 5 12 -3 4] * + conj(O[2][9; 8 13 -4 7]) +end + +function MPSKit.∂AC( + AC::GenericMPSTensor{S,N}, + O::Tuple{T,T,Tuple{Vararg{P,H}}}, + GL::GenericMPSTensor{S,N}, + GR::GenericMPSTensor{S,N}, +) where {S,T<:PEPSTensor,P<:PEPOTensor,N,H} + # sanity check + @assert H == N - 3 + + # collect tensors in convenient order: AC, GL, GR, top, mid, bot + tensors = [AC, GL, GR, O[1], O[3]..., O[2]] + + # contraction order: AC, GL, top, mid..., bot, GR + + # number of contracted legs for full top-mid-bot stack + nlegs_tmb = 5 + 3 * H + + # assign and collect all contraction indices + indicesAC = [1, 3, ((1:3:((H + 1) * 3)) .+ 4)..., 2 + nlegs_tmb] + indicesGL = [-1, 2, ((1:3:((H + 1) * 3)) .+ 3)..., 1] + indicesGR = [((1:N) .+ (1 + nlegs_tmb))..., -(N + 1)] + indicesTop = [6, 3, 3 + nlegs_tmb, -2, 2] + indicesBot = [1 + nlegs_tmb, nlegs_tmb, 4 + H + nlegs_tmb, -N, nlegs_tmb - 1] + indicesMid = Vector{Vector{Int}}(undef, H) + for h in 1:H + indicesMid[h] = [ + 3 + 3 * (h + 1), 3 + 3 * h, 2 + 3 * h, 3 + h + nlegs_tmb, -(2 + h), 1 + 3 * h + ] + end + indices = [indicesAC, indicesGL, indicesGR, indicesTop, indicesMid..., indicesBot] + + # record conjflags + conjlist = [false, false, false, false, repeat([false], H)..., true] + + # perform contraction, permute to restore partition + AC′ = permute(ncon(tensors, indices, conjlist), (Tuple(1:N), (N + 1,))) + + return AC′ +end + +(H::PEPO_∂∂C)(x) = MPSKit.∂C(x, H.GL, H.GR) +(H::PEPO_∂∂AC)(x) = MPSKit.∂AC(x, (H.top, H.bot, H.mid), H.GL, H.GR) + +function MPSKit.∂AC(x::RecursiveVec, O::Tuple{T,T,P}, GL, GR) where {T,P} + return RecursiveVec( + circshift( + map( + (v, O1, O2, O3, l, r) -> ∂AC(v, (O1, O2, O3), l, r), + x.vecs, + O[1], + O[2], + O[3], + GL, + GR, + ), + 1, + ), + ) +end + +Base.:*(H::Union{<:PEPO_∂∂AC,<:PEPO_∂∂C}, v) = H(v) + +# operator constructors +function MPSKit.∂∂C(pos::Int, mps, ::InfiniteTransferPEPO, cache) + return PEPO_∂∂C(leftenv(cache, pos + 1, mps), rightenv(cache, pos, mps)) +end +function MPSKit.∂∂C(col::Int, mps, ::TransferPEPOMultiline, cache) + return PEPO_∂∂C(leftenv(cache, col + 1, mps), rightenv(cache, col, mps)) +end +function MPSKit.∂∂C(row::Int, col::Int, mps, ::TransferPEPOMultiline, cache) + return PEPO_∂∂C(leftenv(cache, row, col + 1, mps), rightenv(cache, row, col, mps)) +end + +function MPSKit.∂∂AC(pos::Int, mps, mpo::InfiniteTransferPEPO, cache) + return PEPO_∂∂AC(mpo[pos]..., lefenv(cache, pos, mps), rightenv(cache, pos, mps)) +end +function MPSKit.∂∂AC(row::Int, col::Int, mps, mpo::TransferPEPOMultiline, cache) + return PEPO_∂∂AC( + mpo[row, col]..., leftenv(cache, row, col, mps), rightenv(cache, row, col, mps) + ) +end +function MPSKit.∂∂AC(col::Int, mps, mpo::TransferPEPOMultiline, cache) + return PEPO_∂∂AC( + map(x -> x[1], mpo[:, col]), + map(x -> x[2], mpo[:, col]), + map(x -> x[3], mpo[:, col]), + leftenv(cache, col, mps), + rightenv(cache, col, mps), + ) +end + +# PEPS derivative + +# specialize simple case +function ∂peps( + AC::GenericMPSTensor{S,4}, + ĀC::GenericMPSTensor{S,4}, + O::Tuple{T,Tuple{P}}, + GL::GenericMPSTensor{S,4}, + GR::GenericMPSTensor{S,4}, +) where {S,T<:PEPSTensor,P<:PEPOTensor} + return @tensor ∂p[-1; -2 -3 -4 -5] := + GL[13 8 10 -5; 1] * + AC[1 9 11 -2; 12] * + O[1][5; 9 3 4 8] * + O[2][1][-1 5; 11 6 7 10] * + GR[12 3 6 -3; 2] * + conj(ĀC[13 4 7 -4; 2]) +end + +function ∂peps( + AC::GenericMPSTensor{S,N}, + ĀC::GenericMPSTensor{S,N}, + O::Tuple{T,Tuple{Vararg{P,H}}}, + GL::GenericMPSTensor{S,N}, + GR::GenericMPSTensor{S,N}, +) where {S,T,P,N,H} + # sanity check + @assert H == N - 3 + + # collect tensors in convenient order: AC, GL, top, mid, GR, ĀC + tensors = [AC, ĀC, GL, GR, O[1], O[2]...] + + # contraction order: AC, GL, top, mid..., bot, GR + + # number of contracted legs for full top-mid stack with AC and GL + nlegs_tm = 2 + 3 * H + + # assign and collect all contraction indices + indicesAC = [1, 3, ((1:3:((H) * 3)) .+ 4)..., -2, 2 + nlegs_tm] + indicesGL = [2 + nlegs_tm + (N - 1), 2, ((1:3:((H) * 3)) .+ 3)..., -5, 1] + indicesTop = [6, 3, 3 + nlegs_tm, 3 + nlegs_tm + (N - 1), 2] + indicesMid = Vector{Vector{Int}}(undef, H) + for h in 1:H + indicesMid[h] = [ + 3 + 3 * (h + 1), + 3 + 3 * h, + 2 + 3 * h, + 3 + h + nlegs_tm, + 3 + h + nlegs_tm + (N - 1), + 1 + 3 * h, + ] + end + indicesMid[end][1] = -1 # bottom physical leg is open + indicesGR = [((1:(N - 1)) .+ (1 + nlegs_tm))..., -3, nlegs_tm + 2 * N] + indicesĀC = [((1:(N - 1)) .+ (nlegs_tm + N))..., -4, nlegs_tm + 2 * N] + indices = [indicesAC, indicesĀC, indicesGL, indicesGR, indicesTop, indicesMid...] + + # record conjflags + conjlist = [false, true, false, false, false, repeat([false], H)...] + + # perform contraction, permute to restore partition + ∂p = permute(ncon(tensors, indices, conjlist), ((1,), Tuple(2:5))) + + return ∂p +end diff --git a/src/operators/infinitepepo.jl b/src/operators/infinitepepo.jl new file mode 100644 index 00000000..c8ffc9f0 --- /dev/null +++ b/src/operators/infinitepepo.jl @@ -0,0 +1,157 @@ +""" + struct InfinitePEPO{T<:PEPOTensor} + +Represents an infinte projected entangled-pair operator (PEPO) on a 3D cubic lattice. +""" +struct InfinitePEPO{T<:PEPOTensor} <: AbstractPEPO + A::Array{T,3} + + function InfinitePEPO(A::Array{T,3}) where {T<:PEPOTensor} + # space checks + for (d, w, h) in Tuple.(CartesianIndices(A)) + space(A[d, w, h], 1) == space(A[d, w, _next(h, end)], 2)' || + throw(SpaceMismatch("Physical space at site $((d, w, h)) does not match.")) + space(A[d, w, h], 3) == space(A[_prev(d, end), w, h], 5)' || throw( + SpaceMismatch("North virtual space at site $((d, w, h)) does not match."), + ) + space(A[d, w, h], 4) == space(A[d, _next(w, end), h], 6)' || throw( + SpaceMismatch("East virtual space at site $((d, w, h)) does not match.") + ) + end + return new{T}(A) + end +end + +## Constructors +""" + InfinitePEPO(A::AbstractArray{T, 2}) + +Allow users to pass in an array of tensors. +""" +function InfinitePEPO(A::AbstractArray{T,3}) where {T<:PEPOTensor} + return InfinitePEPO(Array(deepcopy(A))) +end + +""" + InfinitePEPO(Pspaces, Nspaces, Espaces) + +Allow users to pass in arrays of spaces. +""" +function InfinitePEPO( + Pspaces::AbstractArray{S,3}, + Nspaces::AbstractArray{S,3}, + Espaces::AbstractArray{S,3}=Nspaces, +) where {S<:ElementarySpace} + size(Pspaces) == size(Nspaces) == size(Espaces) || + throw(ArgumentError("Input spaces should have equal sizes.")) + + Sspaces = adjoint.(circshift(Nspaces, (1, 0, 0))) + Wspaces = adjoint.(circshift(Espaces, (0, -1, 0))) + Ppspaces = adjoint.(circshift(Pspaces, (0, 0, -1))) + + A = map(Pspaces, Ppspaces, Nspaces, Espaces, Sspaces, Wspaces) do P, Pp, N, E, S, W + return TensorMap(rand, ComplexF64, P * Pp ← N * E * S * W) + end + + return InfinitePEPO(A) +end + +""" + InfinitePEPO(Pspaces, Nspaces, Espaces) + +Allow users to pass in arrays of spaces, single layer special case. +""" +function InfinitePEPO( + Pspaces::AbstractArray{S,2}, + Nspaces::AbstractArray{S,2}, + Espaces::AbstractArray{S,2}=Nspaces, +) where {S<:ElementarySpace} + size(Pspaces) == size(Nspaces) == size(Espaces) || + throw(ArgumentError("Input spaces should have equal sizes.")) + + Pspaces = reshape(Pspaces, (size(Pspaces)..., 1)) + Nspaces = reshape(Pspaces, (size(Nspaces)..., 1)) + Espaces = reshape(Pspaces, (size(Espaces)..., 1)) + + return InfinitePEPO(Pspaces, Nspaces, Espaces) +end + +""" + InfinitePEPO(A) + +Allow users to pass in single tensor. +""" +function InfinitePEPO(A::T) where {T<:PEPOTensor} + As = Array{T,3}(undef, (1, 1, 1)) + As[1, 1, 1] = A + return InfinitePEPO(As) +end + +""" + InfinitePEPO(Pspace, Nspace, Espace) + +Allow users to pass in single space. +""" +function InfinitePEPO(Pspace::S, Nspace::S, Espace::S=Nspace) where {S<:ElementarySpace} + Pspaces = Array{S,3}(undef, (1, 1, 1)) + Pspaces[1, 1] = Pspace + Nspaces = Array{S,3}(undef, (1, 1, 1)) + Nspaces[1, 1] = Nspace + Espaces = Array{S,3}(undef, (1, 1, 1)) + Espaces[1, 1] = Espace + return InfinitePEPO(Pspaces, Nspaces, Espaces) +end + +""" + InfinitePEPO(d, D) + +Allow users to pass in integers. +""" +function InfinitePEPO(d::Integer, D::Integer) + T = TensorMap(rand, ComplexF64, ℂ^d ⊗ (ℂ^d)' ← ℂ^D ⊗ ℂ^D ⊗ (ℂ^D)' ⊗ (ℂ^D)') + return InfinitePEPO(T) +end + +""" +InfinitePEPO(d, D, L) +InfinitePEPO(d, D, (Lx, Ly, Lz))) + +Allow users to pass in integers and specify unit cell. +""" +function InfinitePEPO(d::Integer, D::Integer, L::Integer) + return InfinitePEPO(d, D, (L, L, L)) +end +function InfinitePEPO(d::Integer, D::Integer, Ls::NTuple{3,Integer}) + T = [TensorMap(rand, ComplexF64, ℂ^d ⊗ (ℂ^d)' ← ℂ^D ⊗ ℂ^D ⊗ (ℂ^D)' ⊗ (ℂ^D)')] + return InfinitePEPO(Array(repeat(T, Ls...))) +end + +## Shape and size +Base.size(T::InfinitePEPO) = size(T.A) +Base.size(T::InfinitePEPO, i) = size(T.A, i) +Base.length(T::InfinitePEPO) = length(T.A) +Base.eltype(T::InfinitePEPO) = eltype(T.A) +VectorInterface.scalartype(T::InfinitePEPO) = scalartype(T.A) + +## Copy +Base.copy(T::InfinitePEPO) = InfinitePEPO(copy(T.A)) +Base.similar(T::InfinitePEPO) = InfinitePEPO(similar(T.A)) +Base.repeat(T::InfinitePEPO, counts...) = InfinitePEPO(repeat(T.A, counts...)) + +Base.getindex(T::InfinitePEPO, args...) = Base.getindex(T.A, args...) +Base.axes(T::InfinitePEPO, args...) = axes(T.A, args...) +TensorKit.space(T::InfinitePEPO, i, j) = space(T[i, j, end], 1) + +Base.rotl90(T::InfinitePEPO) = InfinitePEPO(rotl90(rotl90.(T.A))); + +function initializePEPS( + T::InfinitePEPO{<:PEPOTensor{S}}, vspace::S +) where {S<:ElementarySpace} + Pspaces = Array{S,2}(undef, size(T, 1), size(T, 2)) + for (i, j) in product(1:size(T, 1), 1:size(T, 2)) + Pspaces[i, j] = space(T, i, j) + end + Nspaces = repeat([vspace], size(T, 1), size(T, 2)) + Espaces = repeat([vspace], size(T, 1), size(T, 2)) + return InfinitePEPS(Pspaces, Nspaces, Espaces) +end diff --git a/src/operators/periodicpepo.jl b/src/operators/periodicpepo.jl deleted file mode 100644 index b0d47b3e..00000000 --- a/src/operators/periodicpepo.jl +++ /dev/null @@ -1,98 +0,0 @@ -""" - struct PeriodicPEPO{T<:PEPOTensor} - -Represents an periodic PEPO on a 2D square lattice. -""" -struct PeriodicPEPO{T<:PEPOTensor} <: AbstractPEPO - A::PeriodicArray{T,2} - - function PeriodicPEPO(A::PeriodicArray{T,2}) where {T<:PEPOTensor} - Ivertical = CartesianIndex(-1, 0) - Ihorizontal = CartesianIndex(0, 1) - for I in CartesianIndices(A) - space(A[I], 3) == space(A[I + Ivertical], 5)' || throw( - SpaceMismatch("North virtual space at site $(Tuple(I)) does not match.") - ) - space(A[I], 4) == space(A[I + Ihorizontal], 6)' || throw( - SpaceMismatch("East virtual space at site $(Tuple(I)) does not match.") - ) - end - return new{T}(A) - end -end - -## Constructors -""" - InfinitePEPO(A::AbstractArray{T, 2}) - -Allow users to pass in an array of tensors. -""" -function PeriodicPEPO(A::AbstractArray{T,2}) where {T<:PEPOTensor} - return PeriodicPEPO(PeriodicArray(deepcopy(A))) -end - -""" - InfinitePEPO(Pspaces, Nspaces, Espaces) - -Allow users to pass in arrays of spaces. -""" -function PeriodicPEPO( - Pspaces::AbstractArray{S,2}, - Nspaces::AbstractArray{S,2}, - Espaces::AbstractArray{S,2}=Nspaces, -) where {S<:ElementarySpace} - size(Pspaces) == size(Nspaces) == size(Espaces) || - throw(ArgumentError("Input spaces should have equal sizes.")) - - Sspaces = adjoint.(circshift(Nspaces, (1, 0))) - Wspaces = adjoint.(circshift(Espaces, (0, -1))) - Ppspaces = adjoint.(Pspaces) - - A = map(Pspaces, Ppspaces, Nspaces, Espaces, Sspaces, Wspaces) do P, Pp, N, E, S, W - return TensorMap(rand, ComplexF64, P * Pp ← N * E * S * W) - end - - return PeriodicPEPO(A) -end - -""" - PeriodicPEPO(Pspace, Nspace, Espace) - -Allow users to pass in single space. -""" -function PeriodicPEPO(Pspace::S, Nspace::S, Espace::S=Nspace) where {S<:ElementarySpace} - Pspaces = Array{S,2}(undef, (1, 1)) - Pspaces[1, 1] = Pspace - Nspaces = Array{S,2}(undef, (1, 1)) - Nspaces[1, 1] = Nspace - Espaces = Array{S,2}(undef, (1, 1)) - Espaces[1, 1] = Espace - return PeriodicPEPO(Pspaces, Nspaces, Espaces) -end - -""" - InfinitePEPO(d, D) - -Allow users to pass in integers. -""" -function PeriodicPEPO(d::Integer, D::Integer) - T = [TensorMap(rand, ComplexF64, ℂ^d ⊗ ℂ^d ← ℂ^D ⊗ ℂ^D ⊗ (ℂ^D)' ⊗ (ℂ^D)')] - return PeriodicPEPO(PeriodicArray(reshape(T, (1, 1)))) -end - -PeriodicPEPO(d::Integer, D::Integer, Ls::Tuple{Integer}) = repeat(PeriodicPEPO(d, D), Ls...) - -## Shape and size -Base.size(T::PeriodicPEPO) = size(T.A) -Base.size(T::PeriodicPEPO, i) = size(T.A, i) -Base.length(T::PeriodicPEPO) = length(T.A) - -## Copy -Base.copy(T::PeriodicPEPO) = PeriodicPEPO(copy(T.A)) -Base.similar(T::PeriodicPEPO) = PeriodicPEPO(similar(T.A)) -Base.repeat(T::PeriodicPEPO, counts...) = PeriodicPEPO(repeat(T.A, counts...)) - -Base.getindex(T::PeriodicPEPO, args...) = getindex(T.A, args...); -TensorKit.space(t::PeriodicPEPO, i, j) = space(t[i, j], 1) - -Base.rotl90(t::PeriodicPEPO) = PeriodicPEPO(rotl90(rotl90.(t.A))); diff --git a/src/operators/transferpepo.jl b/src/operators/transferpepo.jl new file mode 100644 index 00000000..a05143c7 --- /dev/null +++ b/src/operators/transferpepo.jl @@ -0,0 +1,204 @@ +struct InfiniteTransferPEPO{T,O} + top::PeriodicArray{T,1} + mid::PeriodicArray{O,2} + bot::PeriodicArray{T,1} # this is assumed to be conjed by default, but do we want this? +end + +InfiniteTransferPEPO(top, mid) = InfiniteTransferPEPO(top, mid, top) + +function InfiniteTransferPEPO(T::InfinitePEPS, O::InfinitePEPO, dir, row) + T = rotate_north(T, dir) + O = rotate_north(O, dir) + return InfiniteTransferPEPO(PeriodicArray(T[row, :]), PeriodicArray(O[row, :, :])) +end + +Base.size(transfer::InfiniteTransferPEPO) = size(transfer.top) +Base.size(transfer::InfiniteTransferPEPO, args...) = size(transfer.top, args...) +Base.length(transfer::InfiniteTransferPEPO) = size(transfer, 1) +height(transfer::InfiniteTransferPEPO) = size(transfer.mid, 2) # will I ever need this? +Base.getindex(O::InfiniteTransferPEPO, i) = (O.top[i], O.bot[i], Tuple(O.mid[i, :])) # TODO: not too sure about this + +Base.iterate(O::InfiniteTransferPEPO, i=1) = i > length(O) ? nothing : (O[i], i + 1) + +function virtual_spaces(O::InfiniteTransferPEPO, i, dir) + return [ + space(O.top[i], dir + 1), + space.(O.mid[i, :], Ref(dir + 2))..., + space(O.bot[i], dir + 1)', + ] +end +north_spaces(O::InfiniteTransferPEPO, i) = virtual_spaces(O, i, NORTH) +east_spaces(O::InfiniteTransferPEPO, i) = virtual_spaces(O, i, EAST) +south_spaces(O::InfiniteTransferPEPO, i) = virtual_spaces(O, i, SOUTH) +west_spaces(O::InfiniteTransferPEPO, i) = virtual_spaces(O, i, WEST) + +function initializeMPS(O::InfiniteTransferPEPO, virtualspaces::AbstractArray{S,1}) where {S} + return InfiniteMPS([ + TensorMap( + rand, + MPSKit.Defaults.eltype, # should be scalartype of transfer PEPO? + virtualspaces[_prev(i, end)] * prod(adjoint.(north_spaces(O, i))), + virtualspaces[mod1(i, end)], + ) for i in 1:length(O) + ]) +end + +function initializeMPS(O::InfiniteTransferPEPO, χ::Int) + return InfiniteMPS([ + TensorMap( + rand, MPSKit.Defaults.eltype, ℂ^χ * prod(adjoint.(north_spaces(O, i))), ℂ^χ + ) for i in 1:length(O) + ]) +end + +const TransferPEPOMultiline = MPSKit.Multiline{<:InfiniteTransferPEPO} +Base.convert(::Type{TransferPEPOMultiline}, O::InfiniteTransferPEPO) = MPSKit.Multiline([O]) +Base.getindex(t::TransferPEPOMultiline, i::Colon, j::Int) = Base.getindex.(t.data[i], j) +Base.getindex(t::TransferPEPOMultiline, i::Int, j) = Base.getindex(t.data[i], j) + +# multiline patch +function TransferPEPOMultiline(T::InfinitePEPS, O::InfinitePEPO, dir) + return MPSKit.Multiline(map(cr -> InfiniteTransferPEPO(T, O, dir, cr), 1:size(T, 1))) +end + +# specialize simple case +function MPSKit.transfer_left( + GL::GenericMPSTensor{S,4}, + O::Tuple{T,T,Tuple{P}}, + A::GenericMPSTensor{S,4}, + Ā::GenericMPSTensor{S,4}, +) where {S,T<:PEPSTensor,P<:PEPOTensor} + @tensor GL′[-1 -2 -3 -4; -5] := + GL[10 7 4 2; 1] * + conj(Ā[10 11 12 13; -1]) * + O[1][8; 9 -2 11 7] * + O[3][1][5 8; 6 -3 12 4] * + conj(O[2][5; 3 -4 13 2]) * + A[1 9 6 3; -5] +end + +# it actually works, somehow +function MPSKit.transfer_left( + GL::GenericMPSTensor{S,N}, + O::Tuple{T,T,Tuple{Vararg{P,H}}}, + A::GenericMPSTensor{S,N}, + Ā::GenericMPSTensor{S,N}, +) where {S,T<:PEPSTensor,P<:PEPOTensor,N,H} + # sanity check + @assert H == N - 3 + + # collect tensors in convenient order: env, above, below, top, mid, bot + tensors = [GL, A, Ā, O[1], O[3]..., O[2]] + + # contraction order: GL, A, top, mid..., bot, Ā + + # number of contracted legs for full top-mid-bot stack + nlegs_tmb = 5 + 3 * H + + # assign and collect all contraction indices + indicesGL = [2 + nlegs_tmb, 2, ((1:3:((H + 1) * 3)) .+ 3)..., 1] + indicesA = [1, 3, ((1:3:((H + 1) * 3)) .+ 4)..., -(N + 1)] + indicesĀ = [((1:N) .+ (1 + nlegs_tmb))..., -1] + indicesTop = [6, 3, -2, 3 + nlegs_tmb, 2] + indicesBot = [1 + nlegs_tmb, nlegs_tmb, -N, 4 + H + nlegs_tmb, nlegs_tmb - 1] + indicesMid = Vector{Vector{Int}}(undef, H) + for h in 1:H + indicesMid[h] = [ + 3 + 3 * (h + 1), 3 + 3 * h, 2 + 3 * h, -(2 + h), 3 + h + nlegs_tmb, 1 + 3 * h + ] + end + indices = [indicesGL, indicesA, indicesĀ, indicesTop, indicesMid..., indicesBot] + + # record conjflags + conjlist = [false, false, true, false, repeat([false], H)..., true] + + # perform contraction, permute to restore partition + GL′ = permute(ncon(tensors, indices, conjlist), (Tuple(1:N), (N + 1,))) + + return GL′ +end + +# specialize simple case +function MPSKit.transfer_right( + GR::GenericMPSTensor{S,4}, + O::Tuple{T,T,Tuple{P}}, + A::GenericMPSTensor{S,4}, + Ā::GenericMPSTensor{S,4}, +) where {S,T<:PEPSTensor,P<:PEPOTensor} + return @tensor GR′[-1 -2 -3 -4; -5] := + GR[10 7 4 2; 1] * + conj(Ā[-5 9 6 3; 1]) * + O[1][8; 11 7 9 -2] * + O[3][1][5 8; 12 4 6 -3] * + conj(O[2][5; 13 2 3 -4]) * + A[-1 11 12 13; 10] +end + +function MPSKit.transfer_right( + GR::GenericMPSTensor{S,N}, + O::Tuple{T,T,Tuple{Vararg{P,H}}}, + A::GenericMPSTensor{S,N}, + Ā::GenericMPSTensor{S,N}, +) where {S,T<:PEPSTensor,P<:PEPOTensor,N,H} + # sanity check + @assert H == N - 3 + + # collect tensors in convenient order: env, above, below, top, mid, bot + tensors = [GR, A, Ā, O[1], O[3]..., O[2]] + + # contraction order: GR, A, top, mid..., bot, Ā + + # number of contracted legs for full top-mid-bot stack + nlegs_tmb = 5 + 3 * H + + # assign and collect all contraction indices + indicesGR = [1, 2, ((1:3:((H + 1) * 3)) .+ 3)..., 2 + nlegs_tmb] + indicesA = [-1, 3, ((1:3:((H + 1) * 3)) .+ 4)..., 1] + indicesĀ = [-(N + 1), ((2:N) .+ (1 + nlegs_tmb))..., 2 + nlegs_tmb] + indicesTop = [6, 3, 2, 3 + nlegs_tmb, -2] + indicesBot = [1 + nlegs_tmb, nlegs_tmb, nlegs_tmb - 1, 4 + H + nlegs_tmb, -N] + indicesMid = Vector{Vector{Int}}(undef, H) + for h in 1:H + indicesMid[h] = [ + 3 + 3 * (h + 1), 3 + 3 * h, 2 + 3 * h, 1 + 3 * h, 3 + h + nlegs_tmb, -(2 + h) + ] + end + indices = [indicesGR, indicesA, indicesĀ, indicesTop, indicesMid..., indicesBot] + + # record conjflags + conjlist = [false, false, true, false, repeat([false], H)..., true] + + # perform contraction, permute to restore partition + GR′ = permute(ncon(tensors, indices, conjlist), (Tuple(1:N), (N + 1,))) + + return GR′ +end + +function MPSKit.expectation_value(st::InfiniteMPS, transfer::InfiniteTransferPEPO) + return expectation_value( + convert(MPSMultiline, st), convert(TransferPEPOMultiline, transfer) + ) +end +function MPSKit.expectation_value(st::MPSMultiline, mpo::TransferPEPOMultiline) + return expectation_value(st, environments(st, mpo)) +end +function MPSKit.expectation_value( + st::MPSMultiline, ca::MPSKit.PerMPOInfEnv{H,V,S,A} +) where {H<:TransferPEPOMultiline,V,S,A} + return expectation_value(st, ca.opp, ca) +end +function MPSKit.expectation_value( + st::MPSMultiline, opp::TransferPEPOMultiline, ca::MPSKit.PerMPOInfEnv +) + retval = PeriodicArray{scalartype(st.AC[1, 1]),2}(undef, size(st, 1), size(st, 2)) + for (i, j) in product(1:size(st, 1), 1:size(st, 2)) + O_ij = opp[i, j] + N = height(opp[1]) + 4 + # just reuse left environment contraction? + GL´ = transfer_left(leftenv(ca, i, j, st), O_ij, st.AC[i, j], st.AC[i + 1, j]) + retval[i, j] = TensorKit.TensorOperations.tensorscalar( + ncon([GL´, rightenv(ca, i, j, st)], [[N, (2:(N - 1))..., 1], [(1:N)...]]) + ) + end + return retval +end diff --git a/src/operators/transferpeps.jl b/src/operators/transferpeps.jl index e40b8365..e5fc4738 100644 --- a/src/operators/transferpeps.jl +++ b/src/operators/transferpeps.jl @@ -22,9 +22,9 @@ function initializeMPS(O::InfiniteTransferPEPS, virtualspaces::AbstractArray{S,1 return InfiniteMPS([ TensorMap( rand, - MPSKit.Defaults.eltype, - virtualspaces[mod1(i - 1, end)] * space(O.top[i], 2)' * space(O.bot[i], 2), - virtualspaces[i], + MPSKit.Defaults.eltype, # should be scalartype of transfer PEPS? + virtualspaces[_prev(i, end)] * space(O.top[i], 2)' * space(O.bot[i], 2), + virtualspaces[mod1(i, end)], ) for i in 1:length(O) ]) end @@ -47,6 +47,21 @@ Base.convert(::Type{TransferPEPSMultiline}, O::InfiniteTransferPEPS) = MPSKit.Mu Base.getindex(t::TransferPEPSMultiline, i::Colon, j::Int) = Base.getindex.(t.data[i], j) Base.getindex(t::TransferPEPSMultiline, i::Int, j) = Base.getindex(t.data[i], j) +# multiline patch +function TransferPEPSMultiline(T::InfinitePEPS, dir) + return MPSKit.Multiline(map(cr -> InfiniteTransferPEPS(T, dir, cr), 1:size(T, 1))) +end +function initializeMPS(O::MPSKit.Multiline, virtualspaces::AbstractArray{S,2}) where {S} + mpss = map(cr -> initializeMPS(O[cr], virtualspaces[cr, :]), 1:size(O, 1)) + return MPSKit.Multiline(mpss) +end +function initializeMPS(O::MPSKit.Multiline, virtualspaces::AbstractArray{S,1}) where {S} + return initializeMPS(O, repeat(virtualspaces', length(O), 1)) +end +function initializeMPS(O::MPSKit.Multiline, χ::Int) + return initializeMPS(O, repeat([ℂ^χ], length(O), length(O[1]))) +end + function MPSKit.transfer_left( GL::GenericMPSTensor{S,3}, O::NTuple{2,PEPSTensor}, @@ -54,10 +69,10 @@ function MPSKit.transfer_left( Ā::GenericMPSTensor{S,3}, ) where {S} return @tensor GL′[-1 -2 -3; -4] := - GL[1 4 2; 7] * + GL[1 2 4; 7] * conj(Ā[1 3 6; -1]) * - O[1][5; 8 -2 3 4] * - conj(O[2][5; 9 -3 6 2]) * + O[1][5; 8 -2 3 2] * + conj(O[2][5; 9 -3 6 4]) * A[7 8 9; -4] end @@ -68,27 +83,43 @@ function MPSKit.transfer_right( Ā::GenericMPSTensor{S,3}, ) where {S} return @tensor GR′[-1 -2 -3; -4] := - GR[7 4 2; 1] * - conj(Ā[-4 6 3; 1]) * - O[1][5; 9 4 6 -2] * + GR[7 6 2; 1] * + conj(Ā[-4 4 3; 1]) * + O[1][5; 9 6 4 -2] * conj(O[2][5; 8 2 3 -3]) * A[-1 9 8 7] end -function MPSKit.expectation_value(st::MPSMultiline, O::TransferPEPSMultiline) +function MPSKit.expectation_value(st::InfiniteMPS, transfer::InfiniteTransferPEPS) + return expectation_value( + convert(MPSMultiline, st), convert(TransferPEPSMultiline, transfer) + ) +end +function MPSKit.expectation_value(st::MPSMultiline, mpo::TransferPEPSMultiline) + return expectation_value(st, environments(st, mpo)) end - function MPSKit.expectation_value( st::MPSMultiline, ca::MPSKit.PerMPOInfEnv{H,V,S,A} ) where {H<:TransferPEPSMultiline,V,S,A} - opp = ca.opp - retval = PeriodicArray{eltype(st.AC[1, 1]),2}(undef, size(st, 1), size(st, 2)) + return expectation_value(st, ca.opp, ca) +end +function MPSKit.expectation_value( + st::MPSMultiline, opp::TransferPEPSMultiline, ca::MPSKit.PerMPOInfEnv +) + retval = PeriodicArray{scalartype(st.AC[1, 1]),2}(undef, size(st, 1), size(st, 2)) for (i, j) in product(1:size(st, 1), 1:size(st, 2)) - retval[i, j] = @plansor leftenv(ca, i, j, st)[1 2; 3] * - opp[i, j][2 4; 6 5] * - st.AC[i, j][3 6; 7] * - rightenv(ca, i, j, st)[7 5; 8] * - conj(st.AC[i + 1, j][1 4; 8]) + O_ij = opp[i, j] + retval[i, j] = @tensor leftenv(ca, i, j, st)[1 2 4; 7] * + conj(st.AC[i + 1, j][1 3 6; 13]) * + O_ij[1][5; 8 11 3 2] * + conj(O_ij[2][5; 9 12 6 4]) * + st.AC[i, j][7 8 9; 10] * + rightenv(ca, i, j, st)[10 11 12; 13] end return retval end + + +# PEPS derivatives +# ---------------- + diff --git a/src/states/infinitepeps.jl b/src/states/infinitepeps.jl index a4345baa..83f37f49 100644 --- a/src/states/infinitepeps.jl +++ b/src/states/infinitepeps.jl @@ -1,12 +1,23 @@ +# not everything is a PeriodicArray anymore +_next(i, total) = mod1(i + 1, total) +_prev(i, total) = mod1(i - 1, total) + """ struct InfinitePEPS{T<:PEPSTensor} -Represents an infinite projected entangled pairs state on a 2D square lattice. +Represents an infinite projected entangled-pair state on a 2D square lattice. """ struct InfinitePEPS{T<:PEPSTensor} <: AbstractPEPS - A::Array{T,2} + A::Array{T,2} # TODO: switch back to PeriodicArray? function InfinitePEPS(A::Array{T,2}) where {T<:PEPSTensor} + for (d, w) in Tuple.(CartesianIndices(A)) + space(A[d, w], 2) == space(A[_prev(d, end), w], 4)' || throw( + SpaceMismatch("North virtual space at site $((d, w)) does not match.") + ) + space(A[d, w], 3) == space(A[d, _next(w, end)], 5)' || + throw(SpaceMismatch("East virtual space at site $((d, w)) does not match.")) + end return new{T}(A) end end @@ -44,6 +55,17 @@ function InfinitePEPS( return InfinitePEPS(A) end +""" + InfinitePEPS(A) + +Allow users to pass in single tensor. +""" +function InfinitePEPS(A::T) where {T<:PEPSTensor} + As = Array{T,2}(undef, (1, 1)) + As[1, 1] = A + return InfinitePEPS(As) +end + """ InfinitePEPS(Pspace, Nspace, Espace) @@ -65,16 +87,20 @@ end Allow users to pass in integers. """ function InfinitePEPS(d::Integer, D::Integer) - T = [TensorMap(rand, ComplexF64, ℂ^d ← ℂ^D ⊗ ℂ^D ⊗ (ℂ^D)' ⊗ (ℂ^D)')] - return InfinitePEPS(Array(reshape(T, (1, 1)))) + T = TensorMap(rand, ComplexF64, ℂ^d ← ℂ^D ⊗ ℂ^D ⊗ (ℂ^D)' ⊗ (ℂ^D)') + return InfinitePEPS(T) end +""" + InfinitePEPS(d, D, L) + InfinitePEPS(d, D, (Lx, Ly))) + +Allow users to pass in integers and specify unit cell. +""" function InfinitePEPS(d::Integer, D::Integer, L::Integer) - T = [TensorMap(rand, ComplexF64, ℂ^d ← ℂ^D ⊗ ℂ^D ⊗ (ℂ^D)' ⊗ (ℂ^D)')] - return InfinitePEPS(Array(repeat(T, L, L))) + return InfinitePEPS(d, D, (L, L)) end - -function InfinitePEPS(d::Integer, D::Integer, Ls::Tuple{Integer,Integer}) +function InfinitePEPS(d::Integer, D::Integer, Ls::NTuple{2,Integer}) T = [TensorMap(rand, ComplexF64, ℂ^d ← ℂ^D ⊗ ℂ^D ⊗ (ℂ^D)' ⊗ (ℂ^D)')] return InfinitePEPS(Array(repeat(T, Ls...))) end @@ -83,14 +109,17 @@ end Base.size(T::InfinitePEPS) = size(T.A) Base.size(T::InfinitePEPS, i) = size(T.A, i) Base.length(T::InfinitePEPS) = length(T.A) -Base.eltype(T::InfinitePEPS) = eltype(eltype(T.A)) +Base.eltype(T::InfinitePEPS) = eltype(T.A) +VectorInterface.scalartype(T::InfinitePEPS) = scalartype(T.A) ## Copy Base.copy(T::InfinitePEPS) = InfinitePEPS(copy(T.A)) Base.similar(T::InfinitePEPS) = InfinitePEPS(similar(T.A)) Base.repeat(T::InfinitePEPS, counts...) = InfinitePEPS(repeat(T.A, counts...)) -Base.getindex(T::InfinitePEPS, args...) = getindex(T.A, args...); +Base.getindex(T::InfinitePEPS, args...) = Base.getindex(T.A, args...) +Base.setindex!(T::InfinitePEPS, args...) = (Base.setindex!(T.A, args...); T) +Base.axes(T::InfinitePEPS, args...) = axes(T.A, args...) TensorKit.space(t::InfinitePEPS, i, j) = space(t[i, j], 1) Base.rotl90(t::InfinitePEPS) = InfinitePEPS(rotl90(rotl90.(t.A))); diff --git a/src/utility/symmetrization.jl b/src/utility/symmetrization.jl new file mode 100644 index 00000000..2f3484e1 --- /dev/null +++ b/src/utility/symmetrization.jl @@ -0,0 +1,181 @@ +# an initial attempt at some basic symmetrization routines for PEPS +# it looks horrible, but it works, maybe + +abstract type SymmetrizationStyle end + +struct None <: SymmetrizationStyle end +struct Depth <: SymmetrizationStyle end +struct Width <: SymmetrizationStyle end +struct Rot <: SymmetrizationStyle end +struct Full <: SymmetrizationStyle end + +# some rather shady definitions for 'hermitian conjugate' at the level of a single tensor +function herm_depth(x::PEPSTensor) + return permute(x', ((5,), (3, 2, 1, 4))) +end +function herm_depth(x::PEPOTensor) + return permute(x', ((5, 6), (3, 2, 1, 4))) +end + +function herm_width(x::PEPSTensor) + x = Permute(Conj(x), [1, 4, 3, 2, 5:(x.legs)]) + return permute(x', ((5,), (1, 4, 3, 2))) +end +function herm_width(x::PEPOTensor) + return permute(x', ((5, 6), (1, 4, 3, 2))) +end + +function herm_height(x::PEPOTensor) + x = Permute(Conj(x), [1:4, 6, 5]) + return permute(x', ((6, 5), (1, 2, 3, 4))) +end + +# hermitian invariance + +# make two TensorMap's have the same spaces, by force if necessary +# this is definitely not what you would want to do, but it circumvents having to think +# about what hermiticity means at the level of transfer operators, which is something +function _make_it_fit( + y::AbstractTensorMap{S,N₁,N₂}, x::AbstractTensorMap{S,N₁,N₂} +) where {S<:IndexSpace,N₁,N₂} + for i in 1:(N₁ + N₂) + if space(x, i) ≠ space(y, i) + f = unitary(space(x, i) ← space(y, i)) + y = permute( + ncon([f, y], [[-i, 1], [-(1:(i - 1))..., 1, -((i + 1):(N₁ + N₂))...]]), + (Tuple(1:N₁), Tuple((N₁ + 1):(N₁ + N₂))), + ) + end + end + return y +end + +function herm_depth_inv(x::Union{PEPSTensor,PEPOTensor}) + return 0.5 * (x + _make_it_fit(herm_depth(x), x)) +end + +function herm_width_inv(x::Union{PEPSTensor,PEPOTensor}) + return 0.5 * (x + _make_it_fit(herm_width(x), x)) +end + +function herm_height_inv(x::Union{PEPSTensor,PEPOTensor}) + return 0.5 * (x + _make_it_fit(herm_height(x), x)) +end + +# rotation invariance + +Base.rotr90(t::PEPSTensor) = permute(t, ((1,), (5, 2, 3, 4))) +Base.rot180(t::PEPSTensor) = permute(t, ((1,), (4, 5, 2, 3))) + +function rot_inv(x) + return 0.25 * ( + x + + _make_it_fit(rotl90(x), x) + + _make_it_fit(rot180(x), x) + + _make_it_fit(rotr90(x), x) + ) +end + + +## PEPS unit cell symmetrization + +PEPSLike = Union{InfinitePEPS,AbstractArray{<:PEPSTensor,2}} + +symmetrize(p::PEPSLike, ::None) = p + +function symmetrize(p::PEPSLike, ::Depth) + depth, width = size(p) + if mod(depth, 2) == 1 + for w in 1:width + p[ceil(Int, depth / 2), w] = herm_depth_inv(p[ceil(Int, depth / 2), w]) + end + end + for d in 1:floor(Int, depth / 2) + for w in 1:width + p[depth - d + 1, w] = _make_it_fit(herm_depth(p[d, w]), p[depth - d + 1, w]) + end + end + return p +end + +function symmetrize(p::PEPSLike, ::Width) + depth, width = size(p) + if mod(width, 2) == 1 + for d in 1:depth + p[d, ceil(Int, width / 2)] = herm_width_inv(p[d, ceil(Int, width / 2), h]) + end + end + for w in 1:floor(Int, width / 2) + for d in 1:depth + p[d, width - w + 1] = _make_it_fit(herm_width(p[d, w]), p[d, width - w + 1]) + end + end + return p +end + +function symmetrize(p::PEPSLike, ::Rot) + return error("TODO") +end + +function symmetrize(p::PEPSLike, ::Full) + # TODO: clean up this mess... + + # some auxiliary transformations + function symmetrize_corner(x::PEPSTensor) + return 0.5 * (x + _make_it_fit(permute(x', ((5,), (4, 3, 2, 1))), x)) + end + symmetrize_center(x::PEPSTensor) = herm_depth_inv(rot_inv(x)) + function symmetrize_mid_depth(x::PEPSTensor) + return x + _make_it_fit(permute(x', ((5,), (3, 2, 1, 4))), x) + end + + depth, width = size(p) + depth == width || + error("This only works for square unit cells.") + + odd = mod(depth, 2) + if odd == 1 + p[ceil(Int, depth / 2), ceil(Int, width / 2)] = symmetrize_center( + p[ceil(Int, depth / 2), ceil(Int, width / 2)] + ) + end + for d in 1:ceil(Int, depth / 2) + for w in 1:floor(Int, width / 2) + if d == w + p[d, w] = symmetrize_corner(p[d, w]) + p[d, width - w + 1] = _make_it_fit(rotr90(p[d, w]), p[d, width - w + 1]) + p[depth - d + 1, w] = _make_it_fit(herm_depth(p[d, w]), p[depth - d + 1, w]) + p[depth - d + 1, width - w + 1] = _make_it_fit( + herm_depth(rotr90(p[d, w])), p[depth - d + 1, width - w + 1] + ) + + elseif odd == 1 && d == ceil(Int, depth / 2) + p[d, w] = symmetrize_mid_depth(p[d, w]) + p[w, d] = _make_it_fit(rotr90(p[d, w]), p[w, d]) + p[d, width - w + 1] = _make_it_fit(rot180(p[d, w]), p[d, width - w + 1]) + p[width - w + 1, d] = _make_it_fit( + herm_depth(rotr90(p[d, w])), p[width - w + 1, d] + ) + + else + p[depth - d + 1, w] = _make_it_fit(herm_depth(p[d, w]), p[depth - d + 1, w]) + p[w, depth - d + 1] = _make_it_fit(rotr90(p[d, w]), p[w, depth - d + 1]) + p[width - w + 1, depth - d + 1] = _make_it_fit( + herm_depth(rotr90(p[d, w])), [width - w + 1, depth - d + 1] + ) + p[w, d] = _make_it_fit(rotr90(herm_depth(p[d, w])), p[w, d]) + p[width - w + 1, d] = _make_it_fit( + herm_depth(rotr90(herm_depth(p[d, w]))), p[width - w + 1, d] + ) + p[d, width - w + 1] = _make_it_fit( + rotr90(rotr90(herm_depth(p[d, w]))), p[d, width - w + 1] + ) + p[depth - d + 1, width - w + 1] = _make_it_fit( + herm_depth(rotr90(rotr90(herm_depth(p[d, w])))), + p[depth - d + 1, width - w + 1], + ) + end + end + end + return p +end diff --git a/test/runtests.jl b/test/runtests.jl index 2c366366..3255ddb4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,6 +1,6 @@ using Test, PEPSKit, MPSKit, TensorKit -@testset "boundary MPS" begin +@testset "Boundary MPS" begin peps = InfinitePEPS(2, 3) tpeps = InfiniteTransferPEPS(peps, 1, 1) @@ -8,3 +8,12 @@ using Test, PEPSKit, MPSKit, TensorKit mps, _, _ = leading_boundary(mps, tpeps, VUMPS()) end + +@testset "CTMRG" begin + peps = InfinitePEPS(2, 3) + tpeps = InfiniteTransferPEPS(peps, 1, 1) + + mps = initializeMPS(tpeps, 4) + + env = leading_boundary(peps, CTMRG(; tol=1e-10)) +end From 9eb4c0bb9ae65d829708804ce83af4ca41831b6a Mon Sep 17 00:00:00 2001 From: Lukas <37111893+lkdvos@users.noreply.github.com> Date: Mon, 16 Oct 2023 22:26:09 +0200 Subject: [PATCH 5/9] Fix typo --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index fd27b9ff..8028adf6 100644 --- a/Project.toml +++ b/Project.toml @@ -19,5 +19,5 @@ VectorInterface = "409d34a3-91d5-4945-b6ec-7529ddf182d8" [compat] MPSKit = "0.10" TensorKit = "0.12" -VectorInterface = '0.4" +VectorInterface = "0.4" julia = "1.6" From f9611c5fcc959bb1b09a480e68bb60516a278183 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Mon, 16 Oct 2023 23:21:49 +0200 Subject: [PATCH 6/9] eltype -> scalartype --- src/environments/ctmrgenv.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/environments/ctmrgenv.jl b/src/environments/ctmrgenv.jl index 327273a3..a00627e4 100644 --- a/src/environments/ctmrgenv.jl +++ b/src/environments/ctmrgenv.jl @@ -19,10 +19,10 @@ function CTMRGEnv(peps_above::InfinitePEPS{P}, peps_below::InfinitePEPS{P}) wher edges = Array{T_type}(undef, 4, size(peps_above)...) for dir in 1:4, i in 1:size(peps_above, 1), j in 1:size(peps_above, 2) - @diffset corners[dir, i, j] = TensorMap(randn, eltype(P), ou, ou) + @diffset corners[dir, i, j] = TensorMap(randn, scalartype(P), ou, ou) @diffset edges[dir, i, j] = TensorMap( randn, - eltype(P), + scalartype(P), ou * space(peps_above[i, j], dir + 1)' * space(peps_below[i, j], dir + 1), ou, ) From 0ee4f4419efafc276786d19f6018097322ad6b96 Mon Sep 17 00:00:00 2001 From: leburgel Date: Tue, 24 Oct 2023 13:47:27 +0200 Subject: [PATCH 7/9] Fix example --- examples/test_pepo.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/test_pepo.jl b/examples/test_pepo.jl index 7c818555..61ee1713 100644 --- a/examples/test_pepo.jl +++ b/examples/test_pepo.jl @@ -101,7 +101,7 @@ pos = (1, 1) # pick a site if height == 1 # no crazy contractions for now... row, col = pos fliprow = size(peps, 1) - row + 1 # below starts counting from below - m = @tensor lw[row, col][13 8 10 18; 1] * + mag = @tensor lw[row, col][13 8 10 18; 1] * above.AC[row, col][1 9 11 15; 12] * peps[row, col][5; 9 3 4 8] * m[14 5; 11 6 7 10] * @@ -112,12 +112,12 @@ if height == 1 # no crazy contractions for now... λ = @tensor lw[row, col][13 8 10 18; 1] * above.AC[row, col][1 9 11 15; 12] * peps[row, col][5; 9 3 4 8] * - t[14 5; 11 6 7 10] * + o[14 5; 11 6 7 10] * rw[row, col][12 3 6 16; 2] * conj(below.AC[fliprow, col][13 4 7 17; 2]) * conj(peps[row, col][14; 15 16 17 18]) - println("m = $(m / λ)") # should be ~0.945 @ beta = 1/3 and D = 3 + println("m = $(mag / λ)") # should be ~0.945 @ beta = 1/3 and D = 3 end nothing From 7ba859d60b4c505fa9846ed27fb91f35064e0171 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Tue, 24 Oct 2023 14:43:21 +0200 Subject: [PATCH 8/9] Refactor InfinitePEPS constructors --- src/states/abstractpeps.jl | 23 +++++++++++++++ src/states/infinitepeps.jl | 57 ++++++++++---------------------------- 2 files changed, 38 insertions(+), 42 deletions(-) diff --git a/src/states/abstractpeps.jl b/src/states/abstractpeps.jl index c74ab89f..490f4695 100644 --- a/src/states/abstractpeps.jl +++ b/src/states/abstractpeps.jl @@ -12,6 +12,29 @@ conventionally ordered as: T : P ← N ⊗ E ⊗ S ⊗ W. """ const PEPSTensor{S} = AbstractTensorMap{S,1,4} where {S<:ElementarySpace} +function PEPSTensor( + f, + ::Type{T}, + Pspace::S, + Nspace::S, + Espace::S=Nspace, + Sspace::S=Nspace', + Wspace::S=Espace', +) where {T,S<:ElementarySpace} + return TensorMap(f, T, Pspace ← Nspace ⊗ Espace ⊗ Sspace ⊗ Wspace) +end +function PEPSTensor( + f, + ::Type{T}, + Pspace::Int, + Nspace::Int, + Espace::Int=Nspace, + Sspace::Int=Nspace, + Wspace::Int=Espace, +) where {T} + return TensorMap(f, T, ℂ^Pspace ← ℂ^Nspace ⊗ ℂ^Espace ⊗ (ℂ^Sspace)' ⊗ (ℂ^Wspace)') +end + """ const PEPOTensor{S} diff --git a/src/states/infinitepeps.jl b/src/states/infinitepeps.jl index 83f37f49..ae4bb31e 100644 --- a/src/states/infinitepeps.jl +++ b/src/states/infinitepeps.jl @@ -29,7 +29,7 @@ end Allow users to pass in an array of tensors. """ function InfinitePEPS(A::AbstractArray{T,2}) where {T<:PEPSTensor} - return InfinitePEPS(Array(deepcopy(A))) + return InfinitePEPS(Array(deepcopy(A))) # TODO: find better way to copy end """ @@ -49,60 +49,33 @@ function InfinitePEPS( Wspaces = adjoint.(circshift(Espaces, (0, -1))) A = map(Pspaces, Nspaces, Espaces, Sspaces, Wspaces) do P, N, E, S, W - return TensorMap(rand, ComplexF64, P ← N * E * S * W) + return PEPSTensor(randn, ComplexF64, P, N, E, S, W) end return InfinitePEPS(A) end """ - InfinitePEPS(A) + InfinitePEPS(A; unitcell=(1, 1)) -Allow users to pass in single tensor. +Create an `InfinitePEPS` by specifying a tensor and unit cell. """ -function InfinitePEPS(A::T) where {T<:PEPSTensor} - As = Array{T,2}(undef, (1, 1)) - As[1, 1] = A - return InfinitePEPS(As) +function InfinitePEPS(A::T; unitcell::Tuple{Int,Int}=(1, 1)) where {T<:PEPSTensor} + return InfinitePEPS(fill(A, unitcell)) end """ - InfinitePEPS(Pspace, Nspace, Espace) + InfinitePEPS(Pspace, Nspace, [Espace]; unitcell=(1,1)) -Allow users to pass in single space. +Create an InfinitePEPS by specifying its spaces and unit cell. Spaces can be specified +either via `Int` or via `ElementarySpace`. """ -function InfinitePEPS(Pspace::S, Nspace::S, Espace::S=Nspace) where {S<:ElementarySpace} - Pspaces = Array{S,2}(undef, (1, 1)) - Pspaces[1, 1] = Pspace - Nspaces = Array{S,2}(undef, (1, 1)) - Nspaces[1, 1] = Nspace - Espaces = Array{S,2}(undef, (1, 1)) - Espaces[1, 1] = Espace - return InfinitePEPS(Pspaces, Nspaces, Espaces) -end - -""" - InfinitePEPS(d, D) - -Allow users to pass in integers. -""" -function InfinitePEPS(d::Integer, D::Integer) - T = TensorMap(rand, ComplexF64, ℂ^d ← ℂ^D ⊗ ℂ^D ⊗ (ℂ^D)' ⊗ (ℂ^D)') - return InfinitePEPS(T) -end - -""" - InfinitePEPS(d, D, L) - InfinitePEPS(d, D, (Lx, Ly))) - -Allow users to pass in integers and specify unit cell. -""" -function InfinitePEPS(d::Integer, D::Integer, L::Integer) - return InfinitePEPS(d, D, (L, L)) -end -function InfinitePEPS(d::Integer, D::Integer, Ls::NTuple{2,Integer}) - T = [TensorMap(rand, ComplexF64, ℂ^d ← ℂ^D ⊗ ℂ^D ⊗ (ℂ^D)' ⊗ (ℂ^D)')] - return InfinitePEPS(Array(repeat(T, Ls...))) +function InfinitePEPS(Pspace::S, Nspace::S, Espace::S=Nspace; unitcell::Tuple{Int,Int}=(1,1)) where {S<:Union{ElementarySpace,Int}} + return InfinitePEPS( + fill(Pspace, unitcell), + fill(Nspace, unitcell), + fill(Espace, unitcell), + ) end ## Shape and size From a99072c9c16de68adaf7e0896ee078a0f24413e0 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Tue, 24 Oct 2023 16:00:19 +0200 Subject: [PATCH 9/9] formatter --- examples/test_pepo.jl | 34 ++++++++++++++++++------------ src/algorithms/pepo_opt.jl | 16 +++----------- src/environments/boundarympsenv.jl | 17 +++++++-------- src/operators/derivatives.jl | 1 - src/operators/transferpeps.jl | 2 -- src/states/infinitepeps.jl | 8 +++---- src/utility/symmetrization.jl | 4 +--- 7 files changed, 36 insertions(+), 46 deletions(-) diff --git a/examples/test_pepo.jl b/examples/test_pepo.jl index 61ee1713..46db0e95 100644 --- a/examples/test_pepo.jl +++ b/examples/test_pepo.jl @@ -7,7 +7,6 @@ using MPSKit using PEPSKit using OptimKit - ## Set some options optim_method = ConjugateGradient @@ -20,14 +19,22 @@ tol_max = 1e-4 tol_factor = 1e-3 symm = Full() hermitian = true -boundary_method = VUMPS(; tol_galerkin=tol_max, maxiter=boundary_maxiter, dynamical_tols=true, eigs_tolfactor=1e-3, envs_tolfactor=1e-3, gauge_tolfactor=1e-6, tol_max=1e-4, verbose=verbosity >= 5) +boundary_method = VUMPS(; + tol_galerkin=tol_max, + maxiter=boundary_maxiter, + dynamical_tols=true, + eigs_tolfactor=1e-3, + envs_tolfactor=1e-3, + gauge_tolfactor=1e-6, + tol_max=1e-4, + verbose=verbosity >= 5, +) # boundary_method = GradientGrassmann - ## Set model parameters # 3D Ising temperature -beta = 1/3 # should give f = -1.0219139... @ D = 3 +beta = 1 / 3 # should give f = -1.0219139... @ D = 3 # bond dimensions D = 3 @@ -38,7 +45,6 @@ depth = 1 width = 1 height = 1 - ## PEPO defintion: 3D classical Ising model t = ComplexF64[exp(beta) exp(-beta); exp(-beta) exp(beta)] @@ -47,19 +53,20 @@ q = sqrt(t) O = zeros(2, 2, 2, 2, 2, 2) O[1, 1, 1, 1, 1, 1] = 1 O[2, 2, 2, 2, 2, 2] = 1 -@tensor o[-1 -2; -3 -4 -5 -6] := O[1 2; 3 4 5 6] * q[-1; 1] * q[-2; 2] * q[-3; 3] * q[-4; 4] * q[-5; 5] * q[-6; 6] +@tensor o[-1 -2; -3 -4 -5 -6] := + O[1 2; 3 4 5 6] * q[-1; 1] * q[-2; 2] * q[-3; 3] * q[-4; 4] * q[-5; 5] * q[-6; 6] M = zeros(2, 2, 2, 2, 2, 2) M[1, 1, 1, 1, 1, 1] = 1 M[2, 2, 2, 2, 2, 2] = -1 -@tensor m[-1 -2; -3 -4 -5 -6] := M[1 2; 3 4 5 6] * q[-1; 1] * q[-2; 2] * q[-3; 3] * q[-4; 4] * q[-5; 5] * q[-6; 6] +@tensor m[-1 -2; -3 -4 -5 -6] := + M[1 2; 3 4 5 6] * q[-1; 1] * q[-2; 2] * q[-3; 3] * q[-4; 4] * q[-5; 5] * q[-6; 6] o = TensorMap(o, ℂ^2 ⊗ (ℂ^2)' ← ℂ^2 ⊗ ℂ^2 ⊗ (ℂ^2)' ⊗ (ℂ^2)') m = TensorMap(m, ℂ^2 ⊗ (ℂ^2)' ← ℂ^2 ⊗ ℂ^2 ⊗ (ℂ^2)' ⊗ (ℂ^2)') O = InfinitePEPO(repeat([o], depth, width, height)) - ## Optimization algorithm pepo_alg = PEPOOptimize(; @@ -76,19 +83,18 @@ pepo_alg = PEPOOptimize(; hermitian, ) - ## Initialize PEPS and MPS fixed points peps = symmetrize(initializePEPS(O, ℂ^D), symm) -envs = pepo_opt_environments(peps, O, pepo_alg.boundary_method; vspaces=[ℂ^χ], hermitian=pepo_alg.hermitian) +envs = pepo_opt_environments( + peps, O, pepo_alg.boundary_method; vspaces=[ℂ^χ], hermitian=pepo_alg.hermitian +) normalize!(peps, envs.peps_boundary) - ## Perform optimization x, f, normgrad = leading_boundary(peps, O, pepo_alg, envs) - ## Unpack result and check magnetization peps = x.state @@ -108,7 +114,7 @@ if height == 1 # no crazy contractions for now... rw[row, col][12 3 6 16; 2] * conj(below.AC[fliprow, col][13 4 7 17; 2]) * conj(peps[row, col][14; 15 16 17 18]) - + λ = @tensor lw[row, col][13 8 10 18; 1] * above.AC[row, col][1 9 11 15; 12] * peps[row, col][5; 9 3 4 8] * @@ -116,7 +122,7 @@ if height == 1 # no crazy contractions for now... rw[row, col][12 3 6 16; 2] * conj(below.AC[fliprow, col][13 4 7 17; 2]) * conj(peps[row, col][14; 15 16 17 18]) - + println("m = $(mag / λ)") # should be ~0.945 @ beta = 1/3 and D = 3 end diff --git a/src/algorithms/pepo_opt.jl b/src/algorithms/pepo_opt.jl index 69053cc7..7182ccd5 100644 --- a/src/algorithms/pepo_opt.jl +++ b/src/algorithms/pepo_opt.jl @@ -54,7 +54,6 @@ end isverbose(alg::VUMPS) = alg.verbose isverbose(alg::GradientGrassmann) = alg.method.verbosity >= 0 - ## Characterize environments for PEPO optimization mutable struct PEPOOptEnv{T,O,F} <: Cache @@ -94,7 +93,6 @@ function MPSKit.recalculate!( return envs end - ## PEPO fixed point optimization algorithm # first attempt at a PEPO-fixed-point-optimization algorithm, bit of a mess... @@ -167,8 +165,7 @@ struct PEPOOptimize{A} end # default PEPO optimization cost function for given PEPO and optimization algorithm -function pepo_opt_costfun( - op::InfinitePEPO, alg::PEPOOptimize) +function pepo_opt_costfun(op::InfinitePEPO, alg::PEPOOptimize) D, W, H = size(op) nrm = D * W * H function pepo_opt_fg(x::PEPOOptPoint) @@ -179,13 +176,7 @@ function pepo_opt_costfun( # recompute environment with scaled tolerance boundary_tol = min(max(alg.tol_min, ng * alg.tol_factor), alg.tol_max) - recalculate!( - envs, - peps, - op; - tol=boundary_tol, - hermitian=alg.hermitian, - ) + recalculate!(envs, peps, op; tol=boundary_tol, hermitian=alg.hermitian) # compute cost function lambdas_peps = expectation_value(peps, envs.peps_boundary) @@ -195,7 +186,7 @@ function pepo_opt_costfun( # compute gradient ∂p_peps = ∂∂peps(peps, envs.peps_boundary) ∂p_pepo = ∂∂peps(peps, op, envs.pepo_boundary) - grad = - (2 / nrm) .* ∂p_pepo ./ lambdas_pepo .+ (2 / nrm) .* ∂p_peps ./ lambdas_peps + grad = -(2 / nrm) .* ∂p_pepo ./ lambdas_pepo .+ (2 / nrm) .* ∂p_peps ./ lambdas_peps grad = symmetrize(grad, alg.symm) # TODO: test if gradient is actually correct @@ -248,7 +239,6 @@ function pepo_opt_costfun( return pepo_opt_fg end - ## The actual leading boundary routine function MPSKit.leading_boundary( diff --git a/src/environments/boundarympsenv.jl b/src/environments/boundarympsenv.jl index 33a861c7..245ce7bc 100644 --- a/src/environments/boundarympsenv.jl +++ b/src/environments/boundarympsenv.jl @@ -1,18 +1,16 @@ # Some form of boundary MPS environments for infinite PEPS and PEPO routines - ## Utility algtol(alg::VUMPS) = alg.tol_galerkin algtol(alg::GradientGrassmann) = alg.method.gradtol -update_tol(alg::VUMPS, tol) = @set alg.tol_galerkin=tol +update_tol(alg::VUMPS, tol) = @set alg.tol_galerkin = tol function update_tol(alg::GradientGrassmann, tol) # annoying disparity between typedef and actual constructor... m = alg.method - m = @set m.gradtol=tol + m = @set m.gradtol = tol return GradientGrassmann(; method=m, (finalize!)=alg.finalize!) end - ## Boundary MPS environment manager mutable struct BoundaryMPSEnv{A,E,F} <: Cache @@ -21,7 +19,6 @@ mutable struct BoundaryMPSEnv{A,E,F} <: Cache alg::F end - ## PEPS boundary MPS function MPSKit.environments( @@ -36,7 +33,11 @@ function MPSKit.environments( end function MPSKit.recalculate!( - envs::BoundaryMPSEnv, peps::InfinitePEPS; tol=algtol(envs.alg), hermitian=false, kwargs... + envs::BoundaryMPSEnv, + peps::InfinitePEPS; + tol=algtol(envs.alg), + hermitian=false, + kwargs..., ) tr = TransferPEPSMultiline(peps, 1) return recalculate!(envs, tr; tol, hermitian, kwargs...) @@ -87,7 +88,6 @@ function ∂∂peps(peps::InfinitePEPS{T}, ca::BoundaryMPSEnv) where {T<:PEPSTen return retval end - ## PEPO boundary MPS function MPSKit.environments( @@ -159,7 +159,6 @@ function ∂∂peps( return retval end - ## The actual routines # because below starts counting from below @@ -229,7 +228,7 @@ function MPSKit.recalculate!( end envs.alg = alg - + return envs end diff --git a/src/operators/derivatives.jl b/src/operators/derivatives.jl index fb82f161..36c0d149 100644 --- a/src/operators/derivatives.jl +++ b/src/operators/derivatives.jl @@ -104,7 +104,6 @@ function ∂peps( conj(ĀC[8 4 -4; 2]) end - # PEPO # ---- diff --git a/src/operators/transferpeps.jl b/src/operators/transferpeps.jl index e5fc4738..05859bbd 100644 --- a/src/operators/transferpeps.jl +++ b/src/operators/transferpeps.jl @@ -119,7 +119,5 @@ function MPSKit.expectation_value( return retval end - # PEPS derivatives # ---------------- - diff --git a/src/states/infinitepeps.jl b/src/states/infinitepeps.jl index ae4bb31e..dc717656 100644 --- a/src/states/infinitepeps.jl +++ b/src/states/infinitepeps.jl @@ -70,11 +70,11 @@ end Create an InfinitePEPS by specifying its spaces and unit cell. Spaces can be specified either via `Int` or via `ElementarySpace`. """ -function InfinitePEPS(Pspace::S, Nspace::S, Espace::S=Nspace; unitcell::Tuple{Int,Int}=(1,1)) where {S<:Union{ElementarySpace,Int}} +function InfinitePEPS( + Pspace::S, Nspace::S, Espace::S=Nspace; unitcell::Tuple{Int,Int}=(1, 1) +) where {S<:Union{ElementarySpace,Int}} return InfinitePEPS( - fill(Pspace, unitcell), - fill(Nspace, unitcell), - fill(Espace, unitcell), + fill(Pspace, unitcell), fill(Nspace, unitcell), fill(Espace, unitcell) ) end diff --git a/src/utility/symmetrization.jl b/src/utility/symmetrization.jl index 2f3484e1..88cc6be1 100644 --- a/src/utility/symmetrization.jl +++ b/src/utility/symmetrization.jl @@ -76,7 +76,6 @@ function rot_inv(x) ) end - ## PEPS unit cell symmetrization PEPSLike = Union{InfinitePEPS,AbstractArray{<:PEPSTensor,2}} @@ -130,8 +129,7 @@ function symmetrize(p::PEPSLike, ::Full) end depth, width = size(p) - depth == width || - error("This only works for square unit cells.") + depth == width || error("This only works for square unit cells.") odd = mod(depth, 2) if odd == 1