From 53e53aac2d653725dab73de5ac4f0ee34779854c Mon Sep 17 00:00:00 2001 From: Lukas <37111893+lkdvos@users.noreply.github.com> Date: Tue, 9 Jul 2024 15:26:35 +0200 Subject: [PATCH] LocalOperator Refactor (#49) This is a major overhaul of the implementation of local operators. In particular, this PR enables the use of arbitrary interactions on a lattice through the `LocalOperator` construction. Additionally, this changes some conventions on how the sites of unit cells are labeled for CTMRG environments. Various other fixes and improvements are also incorporated --- .github/workflows/CI.yml | 28 +- Project.toml | 13 +- src/PEPSKit.jl | 7 +- src/algorithms/ctmrg.jl | 294 +++++++------- src/operators/infinitepepo.jl | 5 + src/operators/localoperator.jl | 694 ++++++++++++++++++++++++--------- src/operators/models.jl | 46 ++- src/operators/transferpepo.jl | 6 +- src/operators/transferpeps.jl | 6 +- src/states/abstractpeps.jl | 9 + src/states/infinitepeps.jl | 24 +- src/utility/autoopt.jl | 61 +++ src/utility/symmetrization.jl | 15 - test/boundarymps/vumps.jl | 7 +- test/ctmrg/gradients.jl | 5 +- test/ctmrg/unitcell.jl | 51 +++ test/heisenberg.jl | 14 +- test/pwave.jl | 6 +- test/runtests.jl | 18 +- test/tf_ising.jl | 46 +++ 20 files changed, 942 insertions(+), 413 deletions(-) create mode 100644 src/utility/autoopt.jl create mode 100644 test/ctmrg/unitcell.jl create mode 100644 test/tf_ising.jl diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 4d0a1f65..aad659d0 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -21,7 +21,7 @@ jobs: fail-fast: false matrix: version: - - '1.6' # LTS version + - '1.8' # LTS version - '1' # automatically expands to the latest stable 1.x release of Julia os: - ubuntu-latest @@ -47,29 +47,3 @@ jobs: CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }} 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@v2 - with: - version: ${{ matrix.version }} - arch: ${{ matrix.arch }} - - uses: julia-actions/cache@v2 - - 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/Project.toml b/Project.toml index f8b9f357..85623d71 100644 --- a/Project.toml +++ b/Project.toml @@ -15,6 +15,7 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec" +TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2" VectorInterface = "409d34a3-91d5-4945-b6ec-7529ddf182d8" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" @@ -24,20 +25,22 @@ ChainRulesCore = "1.0" Compat = "3.46, 4.2" KrylovKit = "0.6, 0.7, 0.8" LinearAlgebra = "1" -MPSKit = "0.10" +MPSKit = "0.11" OptimKit = "0.3" Printf = "1" +Random = "1" Statistics = "1" -TensorKit = "0.12" +TensorKit = "0.12.4" +TensorOperations = "4" VectorInterface = "0.4" Zygote = "0.6" -julia = "1.6" +julia = "1.8" [extras] -SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" ChainRulesTestUtils = "cdddcdb0-9152-4a09-a978-84456f9df70a" FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000" +SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] test = ["Test", "SafeTestsets", "ChainRulesTestUtils", "FiniteDifferences"] diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index 9fc0c239..46879a20 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -5,13 +5,14 @@ using Base: @kwdef using Compat using Accessors using VectorInterface -using TensorKit, KrylovKit, MPSKit, OptimKit +using TensorKit, KrylovKit, MPSKit, OptimKit, TensorOperations using ChainRulesCore, Zygote include("utility/util.jl") include("utility/eigsolve.jl") include("utility/rotations.jl") include("utility/hook_pullback.jl") +include("utility/autoopt.jl") include("states/abstractpeps.jl") include("states/infinitepeps.jl") @@ -59,7 +60,7 @@ module Defaults end export CTMRG, CTMRGEnv -export NLocalOperator, AnisotropicNNOperator, OnSite, NearestNeighbor +export LocalOperator export expectation_value, costfun export leading_boundary export PEPSOptimize, GeomSum, ManualIter, LinSolve @@ -69,6 +70,6 @@ export InfinitePEPO, InfiniteTransferPEPO export initializeMPS, initializePEPS export symmetrize, None, Depth, Full export showtypeofgrad -export square_lattice_heisenberg, square_lattice_pwave +export square_lattice_tf_ising, square_lattice_heisenberg, square_lattice_pwave end # module diff --git a/src/algorithms/ctmrg.jl b/src/algorithms/ctmrg.jl index 18c4d65f..b86ccfc8 100644 --- a/src/algorithms/ctmrg.jl +++ b/src/algorithms/ctmrg.jl @@ -140,18 +140,14 @@ function gauge_fix(envprev::CTMRGEnv{C,T}, envfinal::CTMRGEnv{C,T}) where {C,T} scalartype(T), MPSKit._lastspace(Tsfinal[end])' ← MPSKit._lastspace(M[end])', ) - - ρ_prev = transfermatrix_fixedpoint(Tsprev, M, ρinit) - ρ_final = transfermatrix_fixedpoint(Tsfinal, M, ρinit) + ρprev = transfermatrix_fixedpoint(Tsprev, M, ρinit) + ρfinal = transfermatrix_fixedpoint(Tsfinal, M, ρinit) # Decompose and multiply - Up, _, Vp = tsvd!(ρ_prev) - Uf, _, Vf = tsvd!(ρ_final) - Qprev = Up * Vp - Qfinal = Uf * Vf - σ = Qprev * Qfinal' + Qprev, = leftorth(ρprev) + Qfinal, = leftorth(ρfinal) - return σ + return Qprev * Qfinal' end cornersfix, edgesfix = fix_relative_phases(envfinal, signs) @@ -179,57 +175,66 @@ function transfermatrix_fixedpoint(tops, bottoms, ρinit) end # Explicit fixing of relative phases (doing this compactly in a loop is annoying) +function _contract_gauge_corner(corner, σ_in, σ_out) + @autoopt @tensor corner_fix[χ_in; χ_out] := + σ_in[χ_in; χ1] * corner[χ1; χ2] * conj(σ_out[χ_out; χ2]) +end +function _contract_gauge_edge(edge, σ_in, σ_out) + @autoopt @tensor edge_fix[χ_in D_above D_below; χ_out] := + σ_in[χ_in; χ1] * edge[χ1 D_above D_below; χ2] * conj(σ_out[χ_out; χ2]) +end function fix_relative_phases(envfinal::CTMRGEnv, signs) C1 = map(Iterators.product(axes(envfinal.corners)[2:3]...)) do (r, c) - @tensor Cfix[-1; -2] := - signs[WEST, _prev(r, end), c][-1 1] * - envfinal.corners[NORTHWEST, r, c][1; 2] * - conj(signs[NORTH, r, c][-2 2]) + _contract_gauge_corner( + envfinal.corners[NORTHWEST, r, c], + signs[WEST, r, c], + signs[NORTH, r, _next(c, end)], + ) end T1 = map(Iterators.product(axes(envfinal.edges)[2:3]...)) do (r, c) - @tensor Tfix[-1 -2 -3; -4] := - signs[NORTH, r, c][-1 1] * - envfinal.edges[NORTH, r, c][1 -2 -3; 2] * - conj(signs[NORTH, r, _next(c, end)][-4 2]) + _contract_gauge_edge( + envfinal.edges[NORTH, r, c], + signs[NORTH, r, c], + signs[NORTH, r, _next(c, end)], + ) end - C2 = map(Iterators.product(axes(envfinal.corners)[2:3]...)) do (r, c) - @tensor Cfix[-1; -2] := - signs[NORTH, r, _next(c, end)][-1 1] * - envfinal.corners[NORTHEAST, r, c][1; 2] * - conj(signs[EAST, r, c][-2 2]) + _contract_gauge_corner( + envfinal.corners[NORTHEAST, r, c], + signs[NORTH, r, c], + signs[EAST, _next(r, end), c], + ) end T2 = map(Iterators.product(axes(envfinal.edges)[2:3]...)) do (r, c) - @tensor Tfix[-1 -2 -3; -4] := - signs[EAST, r, c][-1 1] * - envfinal.edges[EAST, r, c][1 -2 -3; 2] * - conj(signs[EAST, _next(r, end), c][-4 2]) + _contract_gauge_edge( + envfinal.edges[EAST, r, c], signs[EAST, r, c], signs[EAST, _next(r, end), c] + ) end - C3 = map(Iterators.product(axes(envfinal.corners)[2:3]...)) do (r, c) - @tensor Cfix[-1; -2] := - signs[EAST, _next(r, end), c][-1 1] * - envfinal.corners[SOUTHEAST, r, c][1; 2] * - conj(signs[SOUTH, r, c][-2 2]) + _contract_gauge_corner( + envfinal.corners[SOUTHEAST, r, c], + signs[EAST, r, c], + signs[SOUTH, r, _prev(c, end)], + ) end T3 = map(Iterators.product(axes(envfinal.edges)[2:3]...)) do (r, c) - @tensor Tfix[-1 -2 -3; -4] := - signs[SOUTH, r, c][-1 1] * - envfinal.edges[SOUTH, r, c][1 -2 -3; 2] * - conj(signs[SOUTH, r, _prev(c, end)][-4 2]) + _contract_gauge_edge( + envfinal.edges[SOUTH, r, c], + signs[SOUTH, r, c], + signs[SOUTH, r, _prev(c, end)], + ) end - C4 = map(Iterators.product(axes(envfinal.corners)[2:3]...)) do (r, c) - @tensor Cfix[-1; -2] := - signs[SOUTH, r, _prev(c, end)][-1 1] * - envfinal.corners[SOUTHWEST, r, c][1; 2] * - conj(signs[WEST, r, c][-2 2]) + _contract_gauge_corner( + envfinal.corners[SOUTHWEST, r, c], + signs[SOUTH, r, c], + signs[WEST, _prev(r, end), c], + ) end T4 = map(Iterators.product(axes(envfinal.edges)[2:3]...)) do (r, c) - @tensor Tfix[-1 -2 -3; -4] := - signs[WEST, r, c][-1 1] * - envfinal.edges[WEST, r, c][1 -2 -3; 2] * - conj(signs[WEST, _prev(r, end), c][-4 2]) + _contract_gauge_edge( + envfinal.edges[WEST, r, c], signs[WEST, r, c], signs[WEST, _prev(r, end), c] + ) end return stack([C1, C2, C3, C4]; dims=1), stack([T1, T2, T3, T4]; dims=1) @@ -305,35 +310,36 @@ function left_move(state, env::CTMRGEnv{C,T}, alg::CTMRG) where {C,T} Pleft, Pright = Zygote.Buffer.(projector_type(T, size(state))) # Use Zygote.Buffer instead of @diffset to avoid ZeroTangent errors in _setindex for col in 1:size(state, 2) - cnext = _next(col, size(state, 2)) + cprev = _prev(col, size(state, 2)) # Compute projectors for row in 1:size(state, 1) + rprev = _prev(row, size(state, 1)) rnext = _next(row, size(state, 1)) - state_nw = state[row, col] - state_sw = rotate_north(state[rnext, col], WEST) # Enlarged corners - Q_sw = northwest_corner( - env.edges[SOUTH, _next(row, end), col], - env.corners[SOUTHWEST, _next(row, end), col], - env.edges[WEST, _next(row, end), col], - state_sw, + Q_sw = southwest_corner( + env.edges[SOUTH, _next(rnext, end), col], + env.corners[SOUTHWEST, _next(rnext, end), cprev], + env.edges[WEST, rnext, cprev], + state[rnext, col], ) Q_nw = northwest_corner( - env.edges[WEST, row, col], - env.corners[NORTHWEST, row, col], - env.edges[NORTH, row, col], - state_nw, + env.edges[WEST, row, cprev], + env.corners[NORTHWEST, rprev, cprev], + env.edges[NORTH, rprev, col], + state[row, col], ) # SVD half-infinite environment trscheme = if alg.fixedspace == true - truncspace(space(env.edges[WEST, row, cnext], 1)) + truncspace(space(env.edges[WEST, row, col], 1)) else alg.trscheme end - @tensor QQ[-1 -2 -3; -4 -5 -6] := Q_sw[-1 -2 -3; 1 2 3] * Q_nw[1 2 3; -4 -5 -6] + @autoopt @tensor QQ[χ_EB D_EBabove D_EBbelow; χ_ET D_ETabove D_ETbelow] := + Q_sw[χ_EB D_EBabove D_EBbelow; χ D1 D2] * + Q_nw[χ D1 D2; χ_ET D_ETabove D_ETbelow] U, S, V, ϵ_local = tsvd!(QQ; trunc=trscheme, alg=TensorKit.SVD()) ϵ = max(ϵ, ϵ_local / norm(S)) # TODO: check if we can just normalize enlarged corners s.t. trunc behaves a bit better @@ -355,81 +361,84 @@ function left_move(state, env::CTMRGEnv{C,T}, alg::CTMRG) where {C,T} # Use projectors to grow the corners & edges for row in 1:size(state, 1) rprev = _prev(row, size(state, 1)) - rnext = _next(row, size(state, 1)) C_sw, C_nw, T_w = grow_env_left( state[row, col], - Pleft[_prev(row, end), col], + Pleft[rprev, col], Pright[row, col], - env.corners[SOUTHWEST, rprev, col], - env.corners[NORTHWEST, rnext, col], - env.edges[SOUTH, rprev, col], - env.edges[WEST, row, col], - env.edges[NORTH, rnext, col], + env.corners[SOUTHWEST, row, cprev], + env.corners[NORTHWEST, row, cprev], + env.edges[SOUTH, row, col], + env.edges[WEST, row, cprev], + env.edges[NORTH, row, col], ) - @diffset corners[SOUTHWEST, rprev, cnext] = C_sw - @diffset corners[NORTHWEST, rnext, cnext] = C_nw - @diffset edges[WEST, row, cnext] = T_w + @diffset corners[SOUTHWEST, row, col] = C_sw / norm(C_sw) + @diffset corners[NORTHWEST, row, col] = C_nw / norm(C_nw) + @diffset edges[WEST, row, col] = T_w / norm(T_w) end - - @diffset corners[SOUTHWEST, :, cnext] ./= norm.(corners[SOUTHWEST, :, cnext]) - @diffset corners[NORTHWEST, :, cnext] ./= norm.(corners[NORTHWEST, :, cnext]) - @diffset edges[WEST, :, cnext] ./= norm.(edges[WEST, :, cnext]) end return CTMRGEnv(corners, edges), copy(Pleft), copy(Pright), ϵ end -# Compute enlarged NW corner -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]) +# Compute enlarged corners +function northwest_corner(edge_W, corner_NW, edge_N, peps_above, peps_below=peps_above) + @autoopt @tensor corner[χ_S D_Sabove D_Sbelow; χ_E D_Eabove D_Ebelow] := + edge_W[χ_S D1 D2; χ1] * + corner_NW[χ1; χ2] * + edge_N[χ2 D3 D4; χ_E] * + peps_above[d; D3 D_Eabove D_Sabove D1] * + conj(peps_below[d; D4 D_Ebelow D_Sbelow D2]) end - -# Compute enlarged NE corner -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]) +function northeast_corner(edge_N, corner_NE, edge_E, peps_above, peps_below=peps_above) + @autoopt @tensor corner[χ_W D_Wabove D_Wbelow; χ_S D_Sabove D_Sbelow] := + edge_N[χ_W D1 D2; χ1] * + corner_NE[χ1; χ2] * + edge_E[χ2 D3 D4; χ_S] * + peps_above[d; D1 D3 D_Sabove D_Wabove] * + conj(peps_below[d; D2 D4 D_Sbelow D_Wbelow]) end - -# Compute enlarged SE corner -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]) +function southeast_corner(edge_E, corner_SE, edge_S, peps_above, peps_below=peps_above) + @autoopt @tensor corner[χ_N D_Nabove D_Nbelow; χ_W D_Wabove D_Wbelow] := + edge_E[χ_N D1 D2; χ1] * + corner_SE[χ1; χ2] * + edge_S[χ2 D3 D4; χ_W] * + peps_above[d; D_Nabove D1 D3 D_Wabove] * + conj(peps_below[d; D_Nbelow D2 D4 D_Wbelow]) +end +function southwest_corner(edge_S, corner_SW, edge_W, peps_above, peps_below=peps_above) + @autoopt @tensor corner[χ_E D_Eabove D_Ebelow; χ_N D_Nabove D_Nbelow] := + edge_S[χ_E D1 D2; χ1] * + corner_SW[χ1; χ2] * + edge_W[χ2 D3 D4; χ_N] * + peps_above[d; D_Nabove D_Eabove D1 D3] * + conj(peps_below[d; D_Nbelow D_Ebelow D2 D4]) end # Build projectors from SVD and enlarged SW & NW corners function build_projectors( - U::AbstractTensorMap{E,3,1}, S, V::AbstractTensorMap{E,1,3}, Q_sw, Q_nw + U::AbstractTensorMap{E,3,1}, S, V::AbstractTensorMap{E,1,3}, Q_SW, Q_NW ) where {E<:ElementarySpace} isqS = sdiag_inv_sqrt(S) - Pl = Q_nw * V' * isqS - Pr = isqS * U' * Q_sw - return Pl, Pr + P_bottom = Q_NW * V' * isqS + P_top = isqS * U' * Q_SW + return P_bottom, P_top end # Apply projectors to entire left half-environment to grow SW & NW corners, and W edge -function grow_env_left(peps, Pl, Pr, C_sw, C_nw, T_s, T_w, T_n) - @tensor C_sw′[-1; -2] := C_sw[1; 4] * T_s[-1 2 3; 1] * Pl[4 2 3; -2] - @tensor C_nw′[-1; -2] := C_nw[1; 2] * T_n[2 3 4; -2] * Pr[-1; 1 3 4] - @tensor T_w′[-1 -2 -3; -4] := - T_w[1 2 3; 4] * - peps[9; 5 -2 7 2] * - conj(peps[9; 6 -3 8 3]) * - Pl[4 5 6; -4] * - Pr[-1; 1 7 8] - return C_sw′, C_nw′, T_w′ +function grow_env_left( + peps, P_bottom, P_top, corners_SW, corners_NW, edge_S, edge_W, edge_N +) + @autoopt @tensor corner_SW′[χ_E; χ_N] := + corners_SW[χ1; χ2] * edge_S[χ_E D1 D2; χ1] * P_bottom[χ2 D1 D2; χ_N] + @autoopt @tensor corner_NW′[χ_S; χ_E] := + corners_NW[χ1; χ2] * edge_N[χ2 D1 D2; χ_E] * P_top[χ_S; χ1 D1 D2] + @autoopt @tensor edge_W′[χ_S D_Eabove D_Ebelow; χ_N] := + edge_W[χ1 D1 D2; χ2] * + peps[d; D3 D_Eabove D5 D1] * + conj(peps[d; D4 D_Ebelow D6 D2]) * + P_bottom[χ2 D3 D4; χ_N] * + P_top[χ_S; χ1 D5 D6] + return corner_SW′, corner_NW′, edge_W′ end @doc """ @@ -442,35 +451,38 @@ function LinearAlgebra.norm(peps::InfinitePEPS, env::CTMRGEnv) total = one(scalartype(peps)) for r in 1:size(peps, 1), c in 1:size(peps, 2) - total *= @tensor env.edges[WEST, r, c][1 2 3; 4] * - env.corners[NORTHWEST, r, c][4; 5] * - env.edges[NORTH, r, c][5 6 7; 8] * - env.corners[NORTHEAST, r, c][8; 9] * - env.edges[EAST, r, c][9 10 11; 12] * - env.corners[SOUTHEAST, r, c][12; 13] * - env.edges[SOUTH, r, c][13 14 15; 16] * - env.corners[SOUTHWEST, r, c][16; 1] * - peps[r, c][17; 6 10 14 2] * - conj(peps[r, c][17; 7 11 15 3]) - - total *= @tensor env.corners[NORTHWEST, r, c][1; 2] * - env.corners[NORTHEAST, r, mod1(c - 1, end)][2; 3] * - env.corners[SOUTHEAST, mod1(r - 1, end), mod1(c - 1, end)][3; 4] * - env.corners[SOUTHWEST, mod1(r - 1, end), c][4; 1] - - total /= @tensor env.edges[WEST, r, c][1 10 11; 4] * - env.corners[NORTHWEST, r, c][4; 5] * - env.corners[NORTHEAST, r, mod1(c - 1, end)][5; 6] * - env.edges[EAST, r, mod1(c - 1, end)][6 10 11; 7] * - env.corners[SOUTHEAST, r, mod1(c - 1, end)][7; 8] * - env.corners[SOUTHWEST, r, c][8; 1] - - total /= @tensor env.corners[NORTHWEST, r, c][1; 2] * - env.edges[NORTH, r, c][2 10 11; 3] * - env.corners[NORTHEAST, r, c][3; 4] * - env.corners[SOUTHEAST, mod1(r - 1, end), c][4; 5] * - env.edges[SOUTH, mod1(r - 1, end), c][5 10 11; 6] * - env.corners[SOUTHWEST, mod1(r - 1, end), c][6; 1] + rprev = _prev(r, size(peps, 1)) + rnext = _next(r, size(peps, 1)) + cprev = _prev(c, size(peps, 2)) + cnext = _next(c, size(peps, 2)) + total *= @autoopt @tensor env.edges[WEST, r, cprev][χ1 D1 D2; χ2] * + env.corners[NORTHWEST, rprev, cprev][χ2; χ3] * + env.edges[NORTH, rprev, c][χ3 D3 D4; χ4] * + env.corners[NORTHEAST, rprev, cnext][χ4; χ5] * + env.edges[EAST, r, cnext][χ5 D5 D6; χ6] * + env.corners[SOUTHEAST, rnext, cnext][χ6; χ7] * + env.edges[SOUTH, rnext, c][χ7 D7 D8; χ8] * + env.corners[SOUTHWEST, rnext, cprev][χ8; χ1] * + peps[r, c][d; D3 D5 D7 D1] * + conj(peps[r, c][d; D4 D6 D8 D2]) + total *= tr( + env.corners[NORTHWEST, rprev, cprev] * + env.corners[NORTHEAST, rprev, c] * + env.corners[SOUTHEAST, r, c] * + env.corners[SOUTHWEST, r, cprev], + ) + total /= @autoopt @tensor env.edges[WEST, r, cprev][χ1 D1 D2; χ2] * + env.corners[NORTHWEST, rprev, cprev][χ2; χ3] * + env.corners[NORTHEAST, rprev, c][χ3; χ4] * + env.edges[EAST, r, c][χ4 D1 D2; χ5] * + env.corners[SOUTHEAST, rnext, c][χ5; χ6] * + env.corners[SOUTHWEST, rnext, cprev][χ6; χ1] + total /= @autoopt @tensor env.corners[NORTHWEST, rprev, cprev][χ1; χ2] * + env.edges[NORTH, rprev, c][χ2 D1 D2; χ3] * + env.corners[NORTHEAST, rprev, cnext][χ3; χ4] * + env.corners[SOUTHEAST, r, cnext][χ4; χ5] * + env.edges[SOUTH, r, c][χ5 D1 D2; χ6] * + env.corners[SOUTHWEST, r, cprev][χ6; χ1] end return total diff --git a/src/operators/infinitepepo.jl b/src/operators/infinitepepo.jl index d8e24e00..a33c7e2c 100644 --- a/src/operators/infinitepepo.jl +++ b/src/operators/infinitepepo.jl @@ -132,3 +132,8 @@ function initializePEPS( Espaces = repeat([vspace], size(T, 1), size(T, 2)) return InfinitePEPS(Pspaces, Nspaces, Espaces) end + +# Rotations +Base.rotl90(T::InfinitePEPO) = InfinitePEPO(stack(rotl90, eachslice(T.A; dims=3))) +Base.rotr90(T::InfinitePEPO) = InfinitePEPO(stack(rotr90, eachslice(T.A; dims=3))) +Base.rot180(T::InfinitePEPO) = InfinitePEPO(stack(rot180, eachslice(T.A; dims=3))) diff --git a/src/operators/localoperator.jl b/src/operators/localoperator.jl index 386869f8..929f82b0 100644 --- a/src/operators/localoperator.jl +++ b/src/operators/localoperator.jl @@ -1,219 +1,527 @@ -# TODO: change this implementation to a type-stable one +# Contraction of local operators on arbitrary lattice locations +# ------------------------------------------------------------- +import MPSKit: tensorexpr -abstract type AbstractInteraction end +# currently need this because MPSKit restricts tensor names to symbols +MPSKit.tensorexpr(ex::Expr, inds) = Expr(:ref, ex, inds...) +function MPSKit.tensorexpr(ex::Expr, indout, indin) + return Expr(:typed_vcat, ex, Expr(:row, indout...), Expr(:row, indin...)) +end """ - struct OnSite <: AbstractInteraction + contract_localoperator(inds, O, peps, env) -Trivial interaction representing terms that act on one isolated site. +Contract a local operator `O` on the PEPS `peps` at the indices `inds` using the environment `env`. """ -struct OnSite <: AbstractInteraction end +function contract_localoperator( + inds::NTuple{N,CartesianIndex{2}}, + O::AbstractTensorMap{S,N,N}, + ket::InfinitePEPS, + bra::InfinitePEPS, + env::CTMRGEnv, +) where {S,N} + static_inds = Val.(inds) + return _contract_localoperator(static_inds, O, ket, bra, env) +end +function contract_localoperator( + inds::NTuple{N,Tuple{Int,Int}}, + O::AbstractTensorMap{S,N,N}, + ket::InfinitePEPS, + bra::InfinitePEPS, + env::CTMRGEnv, +) where {S,N} + return contract_localoperator(CartesianIndex.(inds), O, ket, bra, env) +end -""" - struct NearestNeighbor <: AbstractInteraction +# This implements the contraction of an operator acting on sites `inds`. +# The generated function ensures that we can use @tensor to write dynamic contractions (and maximize performance). +@generated function _contract_localoperator( + inds::NTuple{N,Val}, + O::AbstractTensorMap{S,N,N}, + ket::InfinitePEPS, + bra::InfinitePEPS, + env::CTMRGEnv, +) where {S,N} + cartesian_inds = collect(CartesianIndex{2}, map(x -> x.parameters[1], inds.parameters)) # weird hack to extract information from Val + if !allunique(cartesian_inds) + throw(ArgumentError("Indices should not overlap: $cartesian_inds.")) + end -Interaction representing nearest neighbor terms that act on two adjacent sites. -""" -struct NearestNeighbor <: AbstractInteraction end + rmin, rmax = extrema(getindex.(cartesian_inds, 1)) + cmin, cmax = extrema(getindex.(cartesian_inds, 2)) -""" - struct NLocalOperator{I<:AbstractInteraction} - -Operator in form of a `AbstractTensorMap` which is parametrized by an interaction type. -Mostly, this is used to define Hamiltonian terms and observables. -""" -struct NLocalOperator{I<:AbstractInteraction} - op::AbstractTensorMap + gridsize = (rmax - rmin + 1, cmax - cmin + 1) + + corner_NW = tensorexpr( + :(env.corners[ + NORTHWEST, mod1($(rmin - 1), size(ket, 1)), mod1($(cmin - 1), size(ket, 2)) + ]), + (:C_NW_1,), + (:C_NW_2,), + ) + corner_NE = tensorexpr( + :(env.corners[ + NORTHEAST, mod1($(rmin - 1), size(ket, 1)), mod1($(cmax + 1), size(ket, 2)) + ]), + (:C_NE_1,), + (:C_NE_2,), + ) + corner_SE = tensorexpr( + :(env.corners[ + SOUTHEAST, mod1($(rmax + 1), size(ket, 1)), mod1($(cmax + 1), size(ket, 2)) + ]), + (:C_SE_1,), + (:C_SE_2,), + ) + corner_SW = tensorexpr( + :(env.corners[ + SOUTHWEST, mod1($(rmax + 1), size(ket, 1)), mod1($(cmin - 1), size(ket, 2)) + ]), + (:C_SW_1,), + (:C_SW_2,), + ) + + edges_N = map(1:gridsize[2]) do i + return tensorexpr( + :(env.edges[ + NORTH, + mod1($(rmin - 1), size(ket, 1)), + mod1($(cmin + i - 1), size(ket, 2)), + ]), + ( + (i == 1 ? :C_NW_2 : Symbol(:E_N_virtual, i - 1)), + Symbol(:E_N_top, i), + Symbol(:E_N_bot, i), + ), + ((i == gridsize[2] ? :C_NE_1 : Symbol(:E_N_virtual, i)),), + ) + end + + edges_E = map(1:gridsize[1]) do i + return tensorexpr( + :(env.edges[ + EAST, + mod1($(rmin + i - 1), size(ket, 1)), + mod1($(cmax + 1), size(ket, 2)), + ]), + ( + (i == 1 ? :C_NE_2 : Symbol(:E_E_virtual, i - 1)), + Symbol(:E_E_top, i), + Symbol(:E_E_bot, i), + ), + ((i == gridsize[1] ? :C_SE_1 : Symbol(:E_E_virtual, i)),), + ) + end + + edges_S = map(1:gridsize[2]) do i + return tensorexpr( + :(env.edges[ + SOUTH, + mod1($(rmax + 1), size(ket, 1)), + mod1($(cmin + i - 1), size(ket, 2)), + ]), + ( + (i == gridsize[2] ? :C_SE_2 : Symbol(:E_S_virtual, i)), + Symbol(:E_S_top, i), + Symbol(:E_S_bot, i), + ), + ((i == 1 ? :C_SW_1 : Symbol(:E_S_virtual, i - 1)),), + ) + end + + edges_W = map(1:gridsize[1]) do i + return tensorexpr( + :(env.edges[ + WEST, + mod1($(rmin + i - 1), size(ket, 1)), + mod1($(cmin - 1), size(ket, 2)), + ]), + ( + (i == gridsize[1] ? :C_SW_2 : Symbol(:E_W_virtual, i)), + Symbol(:E_W_top, i), + Symbol(:E_W_bot, i), + ), + ((i == 1 ? :C_NW_1 : Symbol(:E_W_virtual, i - 1)),), + ) + end + + operator = tensorexpr( + :O, ntuple(i -> Symbol(:O_out_, i), N), ntuple(i -> Symbol(:O_in_, i), N) + ) + + bra = map(Iterators.product(1:gridsize[1], 1:gridsize[2])) do (i, j) + inds_id = findfirst(==(CartesianIndex(rmin + i - 1, cmin + j - 1)), cartesian_inds) + physical_label = + isnothing(inds_id) ? Symbol(:physical, i, "_", j) : Symbol(:O_out_, inds_id) + return tensorexpr( + :(bra[ + mod1($(rmin + i - 1), size(bra, 1)), mod1($(cmin + j - 1), size(bra, 2)) + ]), + (physical_label,), + ( + (i == 1 ? Symbol(:E_N_bot, j) : Symbol(:bra_vertical, i - 1, "_", j)), + ( + if j == gridsize[2] + Symbol(:E_E_bot, i) + else + Symbol(:bra_horizontal, i, "_", j) + end + ), + ( + if i == gridsize[1] + Symbol(:E_S_bot, j) + else + Symbol(:bra_vertical, i, "_", j) + end + ), + (j == 1 ? Symbol(:E_W_bot, i) : Symbol(:bra_horizontal, i, "_", j - 1)), + ), + ) + end + + ket = map(Iterators.product(1:gridsize[1], 1:gridsize[2])) do (i, j) + inds_id = findfirst(==(CartesianIndex(rmin + i - 1, cmin + j - 1)), cartesian_inds) + physical_label = + isnothing(inds_id) ? Symbol(:physical, i, "_", j) : Symbol(:O_in_, inds_id) + return tensorexpr( + :(ket[ + mod1($(rmin + i - 1), size(ket, 1)), mod1($(cmin + j - 1), size(ket, 2)) + ]), + (physical_label,), + ( + (i == 1 ? Symbol(:E_N_top, j) : Symbol(:ket_vertical, i - 1, "_", j)), + ( + if j == gridsize[2] + Symbol(:E_E_top, i) + else + Symbol(:ket_horizontal, i, "_", j) + end + ), + ( + if i == gridsize[1] + Symbol(:E_S_top, j) + else + Symbol(:ket_vertical, i, "_", j) + end + ), + (j == 1 ? Symbol(:E_W_top, i) : Symbol(:ket_horizontal, i, "_", j - 1)), + ), + ) + end + + multiplication_ex = Expr( + :call, + :*, + corner_NW, + corner_NE, + corner_SE, + corner_SW, + edges_N..., + edges_E..., + edges_S..., + edges_W..., + ket..., + map(x -> Expr(:call, :conj, x), bra)..., + operator, + ) + + opt_ex = Expr(:tuple) + allinds = TensorOperations.getallindices(multiplication_ex) + for label in allinds + if startswith(String(label), "physical") || startswith(String(label), "O") + push!(opt_ex.args, :($label => $PEPS_PHYSICALDIM)) + elseif startswith(String(label), "ket") || startswith(String(label), "bra") + push!(opt_ex.args, :($label => $PEPS_BONDDIM)) + else + push!(opt_ex.args, :($label => $PEPS_ENVBONDDIM)) + end + end + + return quote + @tensor opt = $opt_ex $multiplication_ex + end end """ - struct AnisotropicNNOperator{I<:AbstractInteraction} - -Operator which includes an on-site term and two nearest-neighbor terms, vertical and horizontal. + contract_localnorm(inds, peps, env) + +Contract a local norm of the PEPS `peps` around indices `inds`. """ -struct AnisotropicNNOperator - h0::NLocalOperator{OnSite} - hx::NLocalOperator{NearestNeighbor} - hy::NLocalOperator{NearestNeighbor} -end -function AnisotropicNNOperator( - h0::AbstractTensorMap{S,1,1}, - hx::AbstractTensorMap{S,2,2}, - hy::AbstractTensorMap{S,2,2}=hx, -) where {S} - return AnisotropicNNOperator( - NLocalOperator{OnSite}(h0), - NLocalOperator{NearestNeighbor}(hx), - NLocalOperator{NearestNeighbor}(hy), - ) +function contract_localnorm( + inds::NTuple{N,CartesianIndex{2}}, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv +) where {N} + static_inds = Val.(inds) + return _contract_localnorm(static_inds, ket, bra, env) end -# TODO: include the on-site term in the two-site terms, to reduce number of contractions. +function contract_localnorm( + inds::NTuple{N,Tuple{Int,Int}}, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv +) where {N} + return contract_localnorm(CartesianIndex.(inds), ket, bra, env) +end +@generated function _contract_localnorm( + inds::NTuple{N,Val}, ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv +) where {N} + cartesian_inds = collect(CartesianIndex{2}, map(x -> x.parameters[1], inds.parameters)) # weird hack to extract information from Val + if !allunique(cartesian_inds) + throw(ArgumentError("Indices should not overlap: $cartesian_inds.")) + end -@doc """ - operator_env(peps::InfinitePEPS, env::CTMRGEnv, ::AbstractInteraction) + rmin, rmax = extrema(getindex.(cartesian_inds, 1)) + cmin, cmax = extrema(getindex.(cartesian_inds, 2)) -Contract a PEPS and a CTMRG environment to form an operator environment. -The open bonds correspond to the indices of an operator with the specified -`AbstractInteraction` type. -""" -operator_env - -function operator_env(peps::InfinitePEPS, env::CTMRGEnv, ::OnSite) - return map(Iterators.product(axes(env.corners, 2), axes(env.corners, 3))) do (r, c) - @tensor opt = true ρ[-1; -2] := - env.corners[NORTHWEST, r, c][1; 2] * - env.edges[NORTH, r, c][2 3 4; 5] * - env.corners[NORTHEAST, r, c][5; 6] * - env.edges[EAST, r, c][6 7 8; 9] * - env.corners[SOUTHEAST, r, c][9; 10] * - env.edges[SOUTH, r, c][10 11 12; 13] * - env.corners[SOUTHWEST, r, c][13; 14] * - env.edges[WEST, r, c][14 15 16; 1] * - peps[r, c][-1; 3 7 11 15] * - conj(peps[r, c][-2; 4 8 12 16]) - end -end - -function operator_env(peps::InfinitePEPS, env::CTMRGEnv, ::NearestNeighbor) - return map(Iterators.product(axes(env.corners, 2), axes(env.corners, 3))) do (r, c) - cnext = _next(c, size(peps, 2)) - @tensor opt = true ρ[-12 -18; -11 -20] := - env.corners[NORTHWEST, r, c][1; 3] * - env.edges[NORTH, r, c][3 5 8; 13] * - env.edges[NORTH, r, cnext][13 16 22; 23] * - env.corners[NORTHEAST, r, cnext][23; 24] * - env.edges[EAST, r, cnext][24 25 26; 27] * - env.corners[SOUTHEAST, r, cnext][27; 28] * - env.edges[SOUTH, r, cnext][28 17 21; 14] * - env.edges[SOUTH, r, c][14 6 10; 4] * - env.corners[SOUTHWEST, r, c][4; 2] * - env.edges[WEST, r, c][2 7 9; 1] * - peps[r, c][-12; 5 15 6 7] * - conj(peps[r, c][-11; 8 19 10 9]) * - peps[r, cnext][-18; 16 25 17 15] * - conj(peps[r, cnext][-20; 22 26 21 19]) - end -end - -@doc """ - MPSKit.expectation_value(peps::InfinitePEPS, env, O::NLocalOperator) - -Evaluate the expectation value of any `NLocalOperator` on each unit-cell entry -of `peps` and `env`. -""" MPSKit.expectation_value(::InfinitePEPS, ::Any, ::NLocalOperator) - -# Optimal contraction order is obtained by manually trying out some space sizes and using costcheck = warn -# in principle, we would like to use opt = true, but this does not give optimal results without also supplying costs -# However, due to a bug in tensoroperations this is currently not possible with integer labels. -# Thus, this is a workaround until the bug is fixed. (otherwise we'd need to rewrite all the labels to symbols...) - -# 1-site operator expectation values on unit cell -function MPSKit.expectation_value( - peps::InfinitePEPS, env::CTMRGEnv, O::NLocalOperator{OnSite} -) - return map(Iterators.product(axes(env.corners, 2), axes(env.corners, 3))) do (r, c) - o = @tensor order = (6, 2, 5, 10, 14, 13, 11, 15, 7, 9, 1, 3, 4, 8, 12, 16, 18, 17) begin - env.corners[NORTHWEST, r, c][1; 2] * - env.edges[NORTH, r, c][2 3 4; 5] * - env.corners[NORTHEAST, r, c][5; 6] * - env.edges[EAST, r, c][6 7 8; 9] * - env.corners[SOUTHEAST, r, c][9; 10] * - env.edges[SOUTH, r, c][10 11 12; 13] * - env.corners[SOUTHWEST, r, c][13; 14] * - env.edges[WEST, r, c][14 15 16; 1] * - peps[r, c][17; 3 7 11 15] * - conj(peps[r, c][18; 4 8 12 16]) * - O.op[18; 17] - end - n = @tensor order = (9, 13, 10, 5, 1, 2, 4, 16, 6, 8, 14, 12, 17, 3, 7, 11, 15) begin - env.corners[NORTHWEST, r, c][1; 2] * - env.edges[NORTH, r, c][2 3 4; 5] * - env.corners[NORTHEAST, r, c][5; 6] * - env.edges[EAST, r, c][6 7 8; 9] * - env.corners[SOUTHEAST, r, c][9; 10] * - env.edges[SOUTH, r, c][10 11 12; 13] * - env.corners[SOUTHWEST, r, c][13; 14] * - env.edges[WEST, r, c][14 15 16; 1] * - peps[r, c][17; 3 7 11 15] * - conj(peps[r, c][17; 4 8 12 16]) - end - o / n - end -end - -#! format: off -function MPSKit.expectation_value( - peps::InfinitePEPS, env, O::NLocalOperator{NearestNeighbor} -) - return map(Iterators.product(axes(env.corners, 2), axes(env.corners, 3))) do (r, c) - cnext = _next(c, size(peps, 2)) - o = @tensor order = ( - 28, 24, 23, 16, 25, 22, 26, 27, 17, 21, 4, 1, 3, 5, 7, 8, 9, 2, 6, 10, 14, 19, - 15, 13, 31, 32, 29, 30, - ) begin # physical spaces - env.corners[NORTHWEST, r, c][1; 3] * - env.edges[NORTH, r, c][3 5 8; 13] * - env.edges[NORTH, r, cnext][13 16 22; 23] * - env.corners[NORTHEAST, r, cnext][23; 24] * - env.edges[EAST, r, cnext][24 25 26; 27] * - env.corners[SOUTHEAST, r, cnext][27; 28] * - env.edges[SOUTH, r, cnext][28 17 21; 14] * - env.edges[SOUTH, r, c][14 6 10; 4] * - env.corners[SOUTHWEST, r, c][4; 2] * - env.edges[WEST, r, c][2 7 9; 1] * - peps[r, c][29; 5 15 6 7] * - conj(peps[r, c][31; 8 19 10 9]) * - peps[r, cnext][30; 16 25 17 15] * - conj(peps[r, cnext][32; 22 26 21 19]) * - O.op[31 32; 29 30] + gridsize = (rmax - rmin + 1, cmax - cmin + 1) + + corner_NW = tensorexpr( + :(env.corners[ + NORTHWEST, mod1($(rmin - 1), size(ket, 1)), mod1($(cmin - 1), size(ket, 2)) + ]), + (:C_NW_1,), + (:C_NW_2,), + ) + corner_NE = tensorexpr( + :(env.corners[ + NORTHEAST, mod1($(rmin - 1), size(ket, 1)), mod1($(cmax + 1), size(ket, 2)) + ]), + (:C_NE_1,), + (:C_NE_2,), + ) + corner_SE = tensorexpr( + :(env.corners[ + SOUTHEAST, mod1($(rmax + 1), size(ket, 1)), mod1($(cmax + 1), size(ket, 2)) + ]), + (:C_SE_1,), + (:C_SE_2,), + ) + corner_SW = tensorexpr( + :(env.corners[ + SOUTHWEST, mod1($(rmax + 1), size(ket, 1)), mod1($(cmin - 1), size(ket, 2)) + ]), + (:C_SW_1,), + (:C_SW_2,), + ) + + edges_N = map(1:gridsize[2]) do i + return tensorexpr( + :(env.edges[ + NORTH, + mod1($(rmin - 1), size(ket, 1)), + mod1($(cmin + i - 1), size(ket, 2)), + ]), + ( + (i == 1 ? :C_NW_2 : Symbol(:E_N_virtual, i - 1)), + Symbol(:E_N_top, i), + Symbol(:E_N_bot, i), + ), + ((i == gridsize[2] ? :C_NE_1 : Symbol(:E_N_virtual, i)),), + ) + end + + edges_E = map(1:gridsize[1]) do i + return tensorexpr( + :(env.edges[ + EAST, + mod1($(rmin + i - 1), size(ket, 1)), + mod1($(cmax + 1), size(ket, 2)), + ]), + ( + (i == 1 ? :C_NE_2 : Symbol(:E_E_virtual, i - 1)), + Symbol(:E_E_top, i), + Symbol(:E_E_bot, i), + ), + ((i == gridsize[1] ? :C_SE_1 : Symbol(:E_E_virtual, i)),), + ) + end + + edges_S = map(1:gridsize[2]) do i + return tensorexpr( + :(env.edges[ + SOUTH, + mod1($(rmax + 1), size(ket, 1)), + mod1($(cmin + i - 1), size(ket, 2)), + ]), + ( + (i == gridsize[2] ? :C_SE_2 : Symbol(:E_S_virtual, i)), + Symbol(:E_S_top, i), + Symbol(:E_S_bot, i), + ), + ((i == 1 ? :C_SW_1 : Symbol(:E_S_virtual, i - 1)),), + ) + end + + edges_W = map(1:gridsize[1]) do i + return tensorexpr( + :(env.edges[ + WEST, + mod1($(rmin + i - 1), size(ket, 1)), + mod1($(cmin - 1), size(ket, 2)), + ]), + ( + (i == gridsize[1] ? :C_SW_2 : Symbol(:E_W_virtual, i)), + Symbol(:E_W_top, i), + Symbol(:E_W_bot, i), + ), + ((i == 1 ? :C_NW_1 : Symbol(:E_W_virtual, i - 1)),), + ) + end + + bra = map(Iterators.product(1:gridsize[1], 1:gridsize[2])) do (i, j) + return tensorexpr( + :(bra[ + mod1($(rmin + i - 1), size(ket, 1)), mod1($(cmin + j - 1), size(ket, 2)) + ]), + (Symbol(:physical, i, "_", j),), + ( + (i == 1 ? Symbol(:E_N_bot, j) : Symbol(:bra_vertical, i - 1, "_", j)), + ( + if j == gridsize[2] + Symbol(:E_E_bot, i) + else + Symbol(:bra_horizontal, i, "_", j) + end + ), + ( + if i == gridsize[1] + Symbol(:E_S_bot, j) + else + Symbol(:bra_vertical, i, "_", j) + end + ), + (j == 1 ? Symbol(:E_W_bot, i) : Symbol(:bra_horizontal, i, "_", j - 1)), + ), + ) + end + + ket = map(Iterators.product(1:gridsize[1], 1:gridsize[2])) do (i, j) + return tensorexpr( + :(ket[ + mod1($(rmin + i - 1), size(ket, 1)), mod1($(cmin + j - 1), size(ket, 2)) + ]), + (Symbol(:physical, i, "_", j),), + ( + (i == 1 ? Symbol(:E_N_top, j) : Symbol(:ket_vertical, i - 1, "_", j)), + ( + if j == gridsize[2] + Symbol(:E_E_top, i) + else + Symbol(:ket_horizontal, i, "_", j) + end + ), + ( + if i == gridsize[1] + Symbol(:E_S_top, j) + else + Symbol(:ket_vertical, i, "_", j) + end + ), + (j == 1 ? Symbol(:E_W_top, i) : Symbol(:ket_horizontal, i, "_", j - 1)), + ), + ) + end + + multiplication_ex = Expr( + :call, + :*, + corner_NW, + corner_NE, + corner_SE, + corner_SW, + edges_N..., + edges_E..., + edges_S..., + edges_W..., + ket..., + map(x -> Expr(:call, :conj, x), bra)..., + ) + + opt_ex = Expr(:tuple) + allinds = TensorOperations.getallindices(multiplication_ex) + for label in allinds + if startswith(String(label), "physical") + push!(opt_ex.args, :($label => $PEPS_PHYSICALDIM)) + elseif startswith(String(label), "ket") || startswith(String(label), "bra") + push!(opt_ex.args, :($label => $PEPS_BONDDIM)) + else + push!(opt_ex.args, :($label => $PEPS_ENVBONDDIM)) end + end - n = @tensor order = ( - 2, 3, 1, 5, 7, 28, 24, 23, 16, 25, 30, 22, 26, 27, 17, 21, 14, 15, 6, 4, 13, 29, - 8, 19, 10, 9, - ) begin - env.corners[NORTHWEST, r, c][1; 3] * - env.edges[NORTH, r, c][3 5 8; 13] * - env.edges[NORTH, r, cnext][13 16 22; 23] * - env.corners[NORTHEAST, r, cnext][23; 24] * - env.edges[EAST, r, cnext][24 25 26; 27] * - env.corners[SOUTHEAST, r, cnext][27; 28] * - env.edges[SOUTH, r, cnext][28 17 21; 14] * - env.edges[SOUTH, r, c][14 6 10; 4] * - env.corners[SOUTHWEST, r, c][4; 2] * - env.edges[WEST, r, c][2 7 9; 1] * - peps[r, c][29; 5 15 6 7] * - conj(peps[r, c][29; 8 19 10 9]) * - peps[r, cnext][30; 16 25 17 15] * - conj(peps[r, cnext][30; 22 26 21 19]) + return quote + @tensor opt = $opt_ex $multiplication_ex + end +end + +# Hamiltonian consisting of local terms +# ------------------------------------- +struct LocalOperator{T<:Tuple,S} + lattice::Matrix{S} + terms::T +end +function LocalOperator(lattice::Matrix{S}, terms::Pair...) where {S} + lattice′ = PeriodicArray(lattice) + for (inds, operator) in terms + @assert operator isa AbstractTensorMap + @assert numout(operator) == numin(operator) == length(inds) + for i in 1:length(inds) + @assert space(operator, i) == lattice′[inds[i]] end - o / n end + return LocalOperator{typeof(terms),S}(lattice, terms) end -#! format: on """ - costfun(peps::InfinitePEPS, env, op::NLocalOperator{NearestNeighbor}) - -Compute the expectation value of a nearest-neighbor operator. -This is used to evaluate and differentiate the energy in ground-state PEPS optimizations. + checklattice(Bool, args...) + checklattice(args...) + +Helper function for checking lattice compatibility. The first version returns a boolean, +while the second version throws an error if the lattices do not match. """ -function costfun(peps::InfinitePEPS, env, op::NLocalOperator{NearestNeighbor}) - oh = sum(expectation_value(peps, env, op)) - ov = sum(expectation_value(rotl90(peps), rotl90(env), op)) - return real(oh + ov) +function checklattice(args...) + return checklattice(Bool, args...) || throw(ArgumentError("Lattice mismatch.")) end +function checklattice(::Type{Bool}, peps::InfinitePEPS, O::LocalOperator) + return size(peps) == size(O.lattice) +end +function checklattice(::Type{Bool}, H::LocalOperator, peps::InfinitePEPS) + return checklattice(Bool, peps, H) +end +@non_differentiable checklattice(args...) -""" - costfun(peps::InfinitePEPS, env, op::AnisotropicNNOperator) - -Compute the expectation value of an on-site and an anisotropic nearest-neighbor operator. -This is used to evaluate and differentiate the energy in ground-state PEPS optimizations. -""" -function costfun(peps::InfinitePEPS, env, op::AnisotropicNNOperator) - oos = sum(expectation_value(peps, env, op.h0)) - oh = sum(expectation_value(peps, env, op.hx)) - ov = sum(expectation_value(rotr90(peps), rotate_north(env, WEST), op.hy)) - #ov = sum(expectation_value(rotl90(peps), rotl90(env), op.hy)) - return real(oos + oh + ov) +function nearest_neighbour_hamiltonian( + lattice::Matrix{S}, h::AbstractTensorMap{S,2,2} +) where {S} + terms = [] + for I in eachindex(IndexCartesian(), lattice) + J1 = I + CartesianIndex(1, 0) + J2 = I + CartesianIndex(0, 1) + push!(terms, (I, J1) => h) + push!(terms, (I, J2) => h) + end + return LocalOperator(lattice, terms...) +end + +function Base.repeat(O::LocalOperator, m::Int, n::Int) + lattice = repeat(O.lattice, m, n) + terms = [] + for (inds, operator) in O.terms, i in 1:m, j in 1:n + offset = CartesianIndex((i - 1) * size(O.lattice, 1), (j - 1) * size(O.lattice, 2)) + push!(terms, (inds .+ Ref(offset)) => operator) + end + return LocalOperator(lattice, terms...) +end + +function MPSKit.expectation_value(peps::InfinitePEPS, O::LocalOperator, envs::CTMRGEnv) + checklattice(peps, O) + return sum(O.terms) do (inds, operator) + contract_localoperator(inds, operator, peps, peps, envs) / + contract_localnorm(inds, peps, peps, envs) + end +end + +function costfun(peps::InfinitePEPS, envs::CTMRGEnv, O::LocalOperator) + E = MPSKit.expectation_value(peps, O, envs) + ignore_derivatives() do + isapprox(imag(E), 0; atol=sqrt(eps(real(E)))) || + @warn "Expectation value is not real: $E." + end + return real(E) end diff --git a/src/operators/models.jl b/src/operators/models.jl index 0c36d7ad..0ea12807 100644 --- a/src/operators/models.jl +++ b/src/operators/models.jl @@ -2,31 +2,57 @@ # ------------------- """ - square_lattice_heisenberg(; Jx=-1, Jy=1, Jz=-1) + square_lattice_tf_ising(::Type{T}=ComplexF64; J=1, h=1, unitcell=(1, 1)) + +Square lattice transverse field Ising model. +""" +function square_lattice_tf_ising( + ::Type{T}=ComplexF64; J=1, h=1, unitcell::Tuple{Int,Int}=(1, 1) +) where {T<:Number} + physical_space = ComplexSpace(2) + lattice = fill(physical_space, 1, 1) + σx = TensorMap(T[0 1; 1 0], physical_space, physical_space) + σz = TensorMap(T[1 0; 0 -1], physical_space, physical_space) + hzz = nearest_neighbour_hamiltonian(lattice, -J / 4 * σz ⊗ σz) + return repeat( + LocalOperator(lattice, hzz.terms..., (CartesianIndex(1, 1),) => -J * h / 2 * σx), + unitcell..., + ) +end + +""" + square_lattice_heisenberg(::Type{T}=ComplexF64; Jx=-1, Jy=1, Jz=-1, unitcell=(1, 1)) Square lattice Heisenberg model. By default, this implements a single site unit cell via a sublattice rotation. """ function square_lattice_heisenberg( - ::Type{T}=ComplexF64; Jx=-1, Jy=1, Jz=-1 + ::Type{T}=ComplexF64; Jx=-1, Jy=1, Jz=-1, unitcell::Tuple{Int,Int}=(1, 1) ) where {T<:Number} physical_space = ComplexSpace(2) + lattice = fill(physical_space, 1, 1) σx = TensorMap(T[0 1; 1 0], physical_space, physical_space) σy = TensorMap(T[0 im; -im 0], physical_space, physical_space) σz = TensorMap(T[1 0; 0 -1], physical_space, physical_space) H = (Jx * σx ⊗ σx) + (Jy * σy ⊗ σy) + (Jz * σz ⊗ σz) - return NLocalOperator{NearestNeighbor}(H / 4) + return repeat(nearest_neighbour_hamiltonian(lattice, H / 4), unitcell...) end """ - square_lattice_pwave(; t=1, μ=2, Δ=1) + square_lattice_pwave(::Type{T}=ComplexF64; t=1, μ=2, Δ=1, unitcell=(1, 1)) Square lattice p-wave superconductor model. """ function square_lattice_pwave( - ::Type{T}=ComplexF64; t::Number=1, μ::Number=2, Δ::Number=1 + ::Type{T}=ComplexF64; + t::Number=1, + μ::Number=2, + Δ::Number=1, + unitcell::Tuple{Int,Int}=(1, 1), ) where {T<:Number} physical_space = Vect[FermionParity](0 => 1, 1 => 1) + lattice = fill(physical_space, 1, 1) + # on-site h0 = TensorMap(zeros, T, physical_space ← physical_space) block(h0, FermionParity(1)) .= -μ @@ -41,5 +67,13 @@ function square_lattice_pwave( block(hy, FermionParity(0)) .= [0 Δ*im; -Δ*im 0] block(hy, FermionParity(1)) .= [0 -t; -t 0] - return AnisotropicNNOperator(h0, hx, hy) + return repeat( + LocalOperator( + lattice, + (CartesianIndex(1, 1),) => h0, + (CartesianIndex(1, 1), CartesianIndex(1, 2)) => hx, + (CartesianIndex(1, 1), CartesianIndex(2, 1)) => hy, + ), + unitcell..., + ) end diff --git a/src/operators/transferpepo.jl b/src/operators/transferpepo.jl index eae23fd4..92a74fae 100644 --- a/src/operators/transferpepo.jl +++ b/src/operators/transferpepo.jl @@ -219,15 +219,13 @@ 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)) + return prod(product(1:size(st, 1), 1:size(st, 2))) do (i, j) 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( + return 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 f1597647..a191a0b7 100644 --- a/src/operators/transferpeps.jl +++ b/src/operators/transferpeps.jl @@ -154,17 +154,15 @@ 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)) + return prod(product(1:size(st, 1), 1:size(st, 2))) do (i, j) O_ij = opp[i, j] - retval[i, j] = @tensor leftenv(ca, i, j, st)[1 2 4; 7] * + return @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 @doc """ diff --git a/src/states/abstractpeps.jl b/src/states/abstractpeps.jl index c324e711..467b05c3 100644 --- a/src/states/abstractpeps.jl +++ b/src/states/abstractpeps.jl @@ -60,3 +60,12 @@ abstract type AbstractPEPS end Abstract supertype for a 2D projected entangled-pair operator. """ abstract type AbstractPEPO end + +# Rotations +Base.rotl90(t::PEPSTensor) = permute(t, ((1,), (3, 4, 5, 2))) +Base.rotr90(t::PEPSTensor) = permute(t, ((1,), (5, 2, 3, 4))) +Base.rot180(t::PEPSTensor) = permute(t, ((1,), (4, 5, 2, 3))) + +Base.rotl90(t::PEPOTensor) = permute(t, ((1, 2), (4, 5, 6, 3))) +Base.rotr90(t::PEPOTensor) = permute(t, ((1, 2), (6, 3, 4, 5))) +Base.rot180(t::PEPOTensor) = permute(t, ((1, 2), (5, 6, 3, 4))) diff --git a/src/states/infinitepeps.jl b/src/states/infinitepeps.jl index 48285bc7..4135e8f1 100644 --- a/src/states/infinitepeps.jl +++ b/src/states/infinitepeps.jl @@ -45,8 +45,8 @@ function InfinitePEPS( 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))) + Sspaces = adjoint.(circshift(Nspaces, (-1, 0))) + Wspaces = adjoint.(circshift(Espaces, (0, 1))) A = map(Pspaces, Nspaces, Espaces, Sspaces, Wspaces) do P, N, E, S, W return PEPSTensor(f, T, P, N, E, S, W) @@ -59,6 +59,21 @@ end InfinitePEPS(A; unitcell=(1, 1)) Create an `InfinitePEPS` by specifying a tensor and unit cell. + +The unit cell is labeled as a matrix which means that any tensor in the unit cell, +regardless if PEPS tensor or environment tensor, is obtained by shifting the row +and column index `[r, c]` by one, respectively: +``` + | | | +---C[r-1,c-1]---T[r-1,c]---T[r-1,c+1]--- + | || || +---T[r,c-1]=====AA[r,c]====AA[r,c+1]==== + | || || +---T[r+1,c-1]===AA[r+1,c]==AA[r+1,c+1]== + | || || +``` +The unit cell has periodic boundary conditions, so `[r, c]` is indexed modulo the +size of the unit cell. """ function InfinitePEPS(A::T; unitcell::Tuple{Int,Int}=(1, 1)) where {T<:PEPSTensor} return InfinitePEPS(fill(A, unitcell)) @@ -128,6 +143,11 @@ end # VectorInterface VectorInterface.zerovector(x::InfinitePEPS) = InfinitePEPS(zerovector(x.A)) +# Rotations +Base.rotl90(t::InfinitePEPS) = InfinitePEPS(rotl90(rotl90.(t.A))) +Base.rotr90(t::InfinitePEPS) = InfinitePEPS(rotr90(rotr90.(t.A))) +Base.rot180(t::InfinitePEPS) = InfinitePEPS(rot180(rot180.(t.A))) + # Chainrules function ChainRulesCore.rrule( ::typeof(Base.getindex), state::InfinitePEPS, row::Int, col::Int diff --git a/src/utility/autoopt.jl b/src/utility/autoopt.jl new file mode 100644 index 00000000..97985878 --- /dev/null +++ b/src/utility/autoopt.jl @@ -0,0 +1,61 @@ +# settings for determining contraction orders +const PEPS_PHYSICALDIM = TensorOperations.Power{:χ}(2, 0) # 2 +const PEPS_BONDDIM = TensorOperations.Power{:χ}(1, 1) # χ +const PEPS_ENVBONDDIM = TensorOperations.Power{:χ}(1, 2) # χ² + +""" + autoopt(ex) + +Preprocessor macro for `@tensor` which automatically inserts costs for all symbols that start with a pattern. +In particular, all labels that start with `d`, `D`, or `χ` are automatically inserted with the corresponding +costs. +""" +macro autoopt(ex) + @assert Meta.isexpr(ex, :macrocall) && ex.args[1] === Symbol("@tensor") "@autoopt expects a @tensor expression:\n$ex" + + # extract expression and kwargs + tensorexpr = ex.args[end] + kwargs = TensorOperations.parse_tensor_kwargs(ex.args[3:(end - 1)]) + + # extract pre-existing opt data + opt_id = findfirst(kwargs) do param + param[1] === :opt + end + if !isnothing(opt_id) + _, val = kwargs[opt_id] + if val isa Expr + optdict = TensorOperations.optdata(val, tensorexpr) + elseif val isa Bool && val + optdict = TensorOperations.optdata(tensorexpr) + else + throw(ArgumentError("Invalid use of `opt`")) + end + else + optdict = TensorOperations.optdata(tensorexpr) + end + + # insert costs for all labels starting with d, D, or χ + replace!(optdict) do (label, cost) + return if startswith(string(label), "d") + label => PEPS_PHYSICALDIM + elseif startswith(string(label), "D") + label => PEPS_BONDDIM + elseif startswith(string(label), "χ") + label => PEPS_ENVBONDDIM + else + label => cost + end + end + + # insert opt data into tensor expression + optexpr = Expr(:tuple, (Expr(:call, :(=>), k, v) for (k, v) in optdict)...) + if isnothing(opt_id) + insert!(ex.args, length(ex.args), :(opt = $optexpr)) + else + ex.args[opt_id + 2] = :(opt = $optexpr) + end + return esc(ex) +end + +# TODO: this definition should be in TensorOperations +TensorOperations.parsecost(x::TensorOperations.Power) = x diff --git a/src/utility/symmetrization.jl b/src/utility/symmetrization.jl index 55a97820..9b9a3b91 100644 --- a/src/utility/symmetrization.jl +++ b/src/utility/symmetrization.jl @@ -60,21 +60,6 @@ function herm_height_inv(x::Union{PEPSTensor,PEPOTensor}) end # rotation invariance -Base.rotl90(t::PEPSTensor) = permute(t, ((1,), (3, 4, 5, 2))) -Base.rotr90(t::PEPSTensor) = permute(t, ((1,), (5, 2, 3, 4))) -Base.rot180(t::PEPSTensor) = permute(t, ((1,), (4, 5, 2, 3))) - -Base.rotl90(t::PEPOTensor) = permute(t, ((1, 2), (4, 5, 6, 3))) -Base.rotr90(t::PEPOTensor) = permute(t, ((1, 2), (6, 3, 4, 5))) -Base.rot180(t::PEPOTensor) = permute(t, ((1, 2), (5, 6, 3, 4))) - -Base.rotl90(t::InfinitePEPS) = InfinitePEPS(rotl90(rotl90.(t.A))) -Base.rotr90(t::InfinitePEPS) = InfinitePEPS(rotr90(rotr90.(t.A))) -Base.rot180(t::InfinitePEPS) = InfinitePEPS(rot180(rot180.(t.A))) - -Base.rotl90(T::InfinitePEPO) = InfinitePEPO(stack(rotl90, eachslice(T.A; dims=3))) -Base.rotr90(T::InfinitePEPO) = InfinitePEPO(stack(rotr90, eachslice(T.A; dims=3))) -Base.rot180(T::InfinitePEPO) = InfinitePEPO(stack(rot180, eachslice(T.A; dims=3))) function rot_inv(x) return 0.25 * ( diff --git a/test/boundarymps/vumps.jl b/test/boundarymps/vumps.jl index 260f40b3..bbf69006 100644 --- a/test/boundarymps/vumps.jl +++ b/test/boundarymps/vumps.jl @@ -7,13 +7,14 @@ using LinearAlgebra Random.seed!(29384293742893) +const vumps_alg = VUMPS(; alg_eigsolve=MPSKit.Defaults.alg_eigsolve(; ishermitian=false)) @testset "(1, 1) PEPS" begin psi = InfinitePEPS(ComplexSpace(2), ComplexSpace(2)) T = PEPSKit.InfiniteTransferPEPS(psi, 1, 1) mps = PEPSKit.initializeMPS(T, [ComplexSpace(20)]) - mps, envs, ϵ = leading_boundary(mps, T, VUMPS()) + mps, envs, ϵ = leading_boundary(mps, T, vumps_alg) N = abs(sum(expectation_value(mps, T))) ctm = leading_boundary( @@ -29,7 +30,7 @@ end T = PEPSKit.TransferPEPSMultiline(psi, 1) mps = PEPSKit.initializeMPS(T, fill(ComplexSpace(20), 2, 2)) - mps, envs, ϵ = leading_boundary(mps, T, VUMPS()) + mps, envs, ϵ = leading_boundary(mps, T, vumps_alg) N = abs(prod(expectation_value(mps, T))) ctm = leading_boundary( @@ -67,6 +68,6 @@ end T = InfiniteTransferPEPO(psi, O, 1, 1) mps = PEPSKit.initializeMPS(T, [ComplexSpace(10)]) - mps, envs, ϵ = leading_boundary(mps, T, VUMPS()) + mps, envs, ϵ = leading_boundary(mps, T, vumps_alg) f = abs(prod(expectation_value(mps, T))) end diff --git a/test/ctmrg/gradients.jl b/test/ctmrg/gradients.jl index 2c9b86f9..1c81063a 100644 --- a/test/ctmrg/gradients.jl +++ b/test/ctmrg/gradients.jl @@ -27,7 +27,10 @@ steps = -0.01:0.005:0.01 ## Tests # ------ -@testset "AD CTMRG energy gradients for $(names[i]) model" for i in eachindex(models) +@testset "AD CTMRG energy gradients for $(names[i]) model" verbose = true for i in + eachindex( + models +) Pspace = Pspaces[i] Vspace = Pspaces[i] Espace = Espaces[i] diff --git a/test/ctmrg/unitcell.jl b/test/ctmrg/unitcell.jl new file mode 100644 index 00000000..9d8cad8a --- /dev/null +++ b/test/ctmrg/unitcell.jl @@ -0,0 +1,51 @@ +using Test +using Random +using PEPSKit +using PEPSKit: _prev, _next, ctmrg_iter +using TensorKit + +# initialize (3, 3) PEPS with unique N and E bond spaces +Random.seed!(91283219347) +unitcell = (3, 3) +Pspace = ℂ^2 +Nspace = ℂ^2 +Espace = ℂ^3 +stype = ComplexF64 +peps = InfinitePEPS(randn, stype, Pspace, Nspace, Espace; unitcell) + +# initialize (3, 3) environment with unique environment spaces +C_type = tensormaptype(spacetype(peps[1]), 1, 1, storagetype(peps[1])) +T_type = tensormaptype(spacetype(peps[1]), 3, 1, storagetype(peps[1])) +corners = Array{C_type}(undef, 4, unitcell...) +edges = Array{T_type}(undef, 4, unitcell...) +for r in 1:unitcell[1], c in 1:unitcell[2] + corners[1, _prev(r, end), _prev(c, end)] = TensorMap(randn, stype, ℂ^5, ℂ^2) + edges[1, _prev(r, end), c] = TensorMap( + randn, stype, ℂ^2 * space(peps[r, c], 1 + 1)' * space(peps[r, c], 1 + 1), ℂ^2 + ) + corners[2, _prev(r, end), _next(c, end)] = TensorMap(randn, stype, ℂ^2, ℂ^3) + edges[2, r, _next(c, end)] = TensorMap( + randn, stype, ℂ^3 * space(peps[r, c], 2 + 1)' * space(peps[r, c], 2 + 1), ℂ^3 + ) + corners[3, _next(r, end), _next(c, end)] = TensorMap(randn, stype, ℂ^3, ℂ^4) + edges[3, _next(r, end), c] = TensorMap( + randn, stype, ℂ^4 * space(peps[r, c], 3 + 1)' * space(peps[r, c], 3 + 1), ℂ^4 + ) + corners[4, _next(r, end), _prev(c, end)] = TensorMap(randn, stype, ℂ^4, ℂ^5) + edges[4, r, _prev(c, end)] = TensorMap( + randn, stype, ℂ^5 * space(peps[r, c], 4 + 1)' * space(peps[r, c], 4 + 1), ℂ^5 + ) +end +env = CTMRGEnv(corners, edges) + +# apply one CTMRG iteration with fixeds +ctm_alg = CTMRG(; fixedspace=true) +env′, = ctmrg_iter(peps, env, ctm_alg) + +# compute random expecation value to test matching bonds +random_tensor = TensorMap(randn, scalartype(peps), Pspace, Pspace) +random_op = repeat( + LocalOperator(fill(ℂ^2, 1, 1), (CartesianIndex(1, 1),) => random_tensor), unitcell... +) +@test expectation_value(peps, random_op, env) isa Number +@test expectation_value(peps, random_op, env′) isa Number diff --git a/test/heisenberg.jl b/test/heisenberg.jl index 35e12ef7..b72268ad 100644 --- a/test/heisenberg.jl +++ b/test/heisenberg.jl @@ -5,8 +5,7 @@ using TensorKit using KrylovKit using OptimKit -# Initialize parameters -H = square_lattice_heisenberg() +# initialize parameters χbond = 2 χenv = 16 ctm_alg = CTMRG(; trscheme=truncdim(χenv), tol=1e-10, miniter=4, maxiter=100, verbosity=1) @@ -20,6 +19,7 @@ opt_alg = PEPSOptimize(; # initialize states Random.seed!(91283219347) +H = square_lattice_heisenberg() psi_init = InfinitePEPS(2, χbond) env_init = leading_boundary(CTMRGEnv(psi_init; Venv=ComplexSpace(χenv)), psi_init, ctm_alg) @@ -27,3 +27,13 @@ env_init = leading_boundary(CTMRGEnv(psi_init; Venv=ComplexSpace(χenv)), psi_in result = fixedpoint(psi_init, H, opt_alg, env_init) @test result.E ≈ -0.6694421 atol = 1e-2 + +# same test but for 2x2 unit cell +H_2x2 = square_lattice_heisenberg(; unitcell=(2, 2)) +psi_init_2x2 = InfinitePEPS(2, χbond; unitcell=(2, 2)) +env_init_2x2 = leading_boundary( + CTMRGEnv(psi_init_2x2; Venv=ComplexSpace(χenv)), psi_init_2x2, ctm_alg +) +result_2x2 = fixedpoint(psi_init_2x2, H_2x2, opt_alg, env_init_2x2) + +@test result_2x2.E ≈ 4 * result.E atol = 1e-2 diff --git a/test/pwave.jl b/test/pwave.jl index 4ebb9c85..0c08c93b 100644 --- a/test/pwave.jl +++ b/test/pwave.jl @@ -6,7 +6,8 @@ using KrylovKit using OptimKit # Initialize parameters -H = square_lattice_pwave() +unitcell = (2, 2) +H = square_lattice_pwave(; unitcell) χbond = 2 χenv = 16 ctm_alg = CTMRG(; @@ -21,11 +22,10 @@ opt_alg = PEPSOptimize(; ) # initialize states -Random.seed!(96678827397) Pspace = Vect[FermionParity](0 => 1, 1 => 1) Vspace = Vect[FermionParity](0 => χbond ÷ 2, 1 => χbond ÷ 2) Envspace = Vect[FermionParity](0 => χenv ÷ 2, 1 => χenv ÷ 2) -psi_init = InfinitePEPS(Pspace, Vspace, Vspace; unitcell=(2, 2)) +psi_init = InfinitePEPS(Pspace, Vspace, Vspace; unitcell) env_init = leading_boundary(CTMRGEnv(psi_init; Venv=Envspace), psi_init, ctm_alg); # find fixedpoint diff --git a/test/runtests.jl b/test/runtests.jl index af2023e1..799726b8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,9 +1,16 @@ using SafeTestsets -const GROUP = get(ENV, "GROUP", "All") +# check if user supplied args +pat = r"(?:--group=)(\w+)" +arg_id = findfirst(contains(pat), ARGS) +const GROUP = if isnothing(arg_id) + uppercase(get(ENV, "GROUP", "ALL")) +else + uppercase(only(match(pat, ARGS[arg_id]).captures)) +end @time begin - if GROUP == "All" || GROUP == "CTMRG" + if GROUP == "ALL" || GROUP == "CTMRG" @time @safetestset "Gauge Fixing" begin include("ctmrg/gaugefix.jl") end @@ -11,12 +18,15 @@ const GROUP = get(ENV, "GROUP", "All") include("ctmrg/gradients.jl") end end - if GROUP == "All" || GROUP == "MPS" + if GROUP == "ALL" || GROUP == "MPS" @time @safetestset "VUMPS" begin include("boundarymps/vumps.jl") end end - if GROUP == "All" || GROUP == "EXAMPLES" + if GROUP == "ALL" || GROUP == "EXAMPLES" + @time @safetestset "Transverse Field Ising model" begin + include("tf_ising.jl") + end @time @safetestset "Heisenberg model" begin include("heisenberg.jl") end diff --git a/test/tf_ising.jl b/test/tf_ising.jl new file mode 100644 index 00000000..881bdd13 --- /dev/null +++ b/test/tf_ising.jl @@ -0,0 +1,46 @@ +using Test +using Random +using PEPSKit +using TensorKit +using KrylovKit +using OptimKit + +# References +# ---------- +# Classical Simulation of Infinite-Size Quantum Lattice Systems in Two Spatial Dimensions +# J. Jordan, R. Orús, G. Vidal, F. Verstraete, and J. I. Cirac +# Phys. Rev. Lett. 101, 250602 – Published 18 December 2008 +# (values estimated from plots) +# (factor of 2 in the energy and magnetisation due to convention differences) +h = 0.5 +e = -0.525 +mᶻ = 0.98 + +# initialize parameters +χbond = 2 +χenv = 16 +ctm_alg = CTMRG(; trscheme=truncdim(χenv), tol=1e-10, miniter=4, maxiter=100, verbosity=1) +opt_alg = PEPSOptimize(; + boundary_alg=ctm_alg, + optimizer=LBFGS(4; maxiter=100, gradtol=1e-3, verbosity=2), + gradient_alg=GMRES(; tol=1e-6, maxiter=100), + reuse_env=true, + verbosity=2, +) + +# initialize states +H = square_lattice_tf_ising(; h) +psi_init = InfinitePEPS(2, χbond) +env_init = leading_boundary(CTMRGEnv(psi_init; Venv=ComplexSpace(χenv)), psi_init, ctm_alg) + +# find fixedpoint +result = fixedpoint(psi_init, H, opt_alg, env_init) + +# compute magnetization +σz = TensorMap(scalartype(psi_init)[1 0; 0 -1], ℂ^2, ℂ^2) +M = LocalOperator(H.lattice, (CartesianIndex(1, 1),) => σz) +magn = expectation_value(result.peps, M, result.env) + +@test result.E ≈ e atol = 5e-2 +@test imag(magn) ≈ 0 atol = 1e-6 +@test abs(magn) ≈ mᶻ atol = 5e-2