diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index b715dc2ec..c42c701ef 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -3,7 +3,7 @@ module PEPSKit using LinearAlgebra, Statistics, Base.Threads, Base.Iterators, Printf using Random using Compat -using Accessors: @set, @reset +using Accessors: @set, @reset, @insert using VectorInterface import VectorInterface as VI @@ -125,6 +125,7 @@ include("algorithms/contractions/local_patch/densitymatrix/bp.jl") include("algorithms/contractions/bondenv/benv_tools.jl") include("algorithms/contractions/bondenv/gaugefix.jl") include("algorithms/contractions/bondenv/als_solve.jl") +include("algorithms/contractions/bondenv/benv_ntu.jl") include("algorithms/contractions/bondenv/benv_ctm.jl") include("algorithms/contractions/correlator/peps.jl") include("algorithms/contractions/correlator/pepo_purified.jl") @@ -157,6 +158,8 @@ include("algorithms/time_evolution/trotter_gate.jl") include("algorithms/time_evolution/time_evolve.jl") include("algorithms/time_evolution/simpleupdate.jl") include("algorithms/time_evolution/simpleupdate3site.jl") +include("algorithms/time_evolution/ntupdate.jl") +include("algorithms/time_evolution/ntupdate3site.jl") include("algorithms/time_evolution/gaugefix_su.jl") include("algorithms/bp/beliefpropagation.jl") @@ -201,7 +204,8 @@ export compress export absorb_weight export ALSTruncation, FullEnvTruncation -export SimpleUpdate +export NNEnv, NNpEnv +export SimpleUpdate, NeighbourUpdate export TimeEvolver, timestep, time_evolve export InfiniteSquareNetwork diff --git a/src/algorithms/contractions/bondenv/benv_ntu.jl b/src/algorithms/contractions/bondenv/benv_ntu.jl new file mode 100644 index 000000000..b7e32b6a1 --- /dev/null +++ b/src/algorithms/contractions/bondenv/benv_ntu.jl @@ -0,0 +1,232 @@ +#= +The construction of bond environment for Neighborhood Tensor Update (NTU) +is adapted from YASTN (https://github.com/yastn/yastn). +Copyright 2024 The YASTN Authors. All Rights Reserved. +Licensed under the Apache License, Version 2.0 +=# + +""" +Algorithms to construct bond environment for Neighborhood Tensor Update (NTU). +""" +abstract type NeighbourEnv end + +""" +For a rank-(2,2) tensor `T_{ik;jl}`, approximately decompose it to the +product of two rank-2 tensors `A_{ik;} * B_{;jl}` using truncated SVD +that keeps only the largest singular value in the charge-neutral sector. + +The kept singular value is the largest in the entire SVD spectrum if: +- `permute(T, ((1, 3), (2, 4)))` is positive semi-definite, and +- `space(T, 1) == space(T, 2)'` and `space(T, 3) == space(T, 4)`, +""" +function _svd_cut!(t::AbstractTensorMap{<:Any, <:Any, 2, 2}) + A, B = left_orth!(t; trunc = truncspace(oneunit(spacetype(t)))) + return removeunit(A, numind(A)), removeunit(B, 1) +end + +""" +Algorithm struct for "NTU-NN" bond environment. +""" +struct NNEnv <: NeighbourEnv end +""" +Calculate the bond environment within "NTU-NN" approximation. +`(row, col)` locates the site of `X` in the original `state`. +`dir` is the direction from X to `Y` (`NORTH`, `EAST`, `SOUTH`, or `WEST`). +`X` and `Y` must be orthogonal factors in the east-facing bond frame, +while `state` stays in its original orientation. +`allocator` supplies temporary storage for the common environment contraction. +``` + -1 ●=======● + ║ ║ + 0 ●===X== ==Y===● + ║ ║ + 1 ●=======● + -1 0 1 2 +``` +""" +function bondenv_ntu( + row::Int, col::Int, dir::Int, X, Y, state::InfiniteState, alg::NNEnv; + allocator = DefaultAllocator() + ) + neighbors = [(-1, 0), (0, -1), (1, 0), (1, 1), (0, 2), (-1, 1)] + m = collect_neighbors(state, row, col, neighbors; direction = dir) + X, Y = _prepare_site_tensor(X), _prepare_site_tensor(Y) + return _contract_ntu_NNEnv( + X, Y, hair_w(m[0, -1]), hair_e(m[0, 2]), + cor_nw(m[-1, 0]), cor_ne(m[-1, 1]), cor_sw(m[1, 0]), cor_se(m[1, 1]); allocator + ) +end + +""" +Algorithm struct for "NTU-NN+" bond environment. +""" +struct NNpEnv <: NeighbourEnv end +""" +Calculate the bond environment within "NTU-NN+" approximation. +`(row, col)` locates the site of `X` in the original `state`. +`dir` is the direction from X to `Y` (`NORTH`, `EAST`, `SOUTH`, or `WEST`). +`X` and `Y` must be orthogonal factors in the east-facing bond frame, +while `state` stays in its original orientation. +`allocator` supplies temporary storage for the common environment contraction. +``` + -2 ●.......● + ║ ║ + -1 ○===●=======●===○ + ║ ║ ║ ║ + 0 ●===●===X== ==Y===●===● + ║ ║ ║ ║ + 1 ○===●=======●===○ + ║ ║ + 2 ●.......● + -2 -1 0 1 2 3 +``` +Dotted lines and ○ are splitted using SVD with `truncrank(1)`. +""" +function bondenv_ntu( + row::Int, col::Int, dir::Int, X, Y, state::InfiniteState, alg::NNpEnv; + allocator = DefaultAllocator() + ) + neighbors = [ + (-1, -1), (0, -1), (1, -1), (1, 0), (1, 1), (1, 2), (0, 2), (-1, 2), + (-1, 1), (-1, 0), (0, -2), (2, 0), (2, 1), (0, 3), (-2, 1), (-2, 0), + ] + ms = collect_neighbors(state, row, col, neighbors; direction = dir) + X, Y = _prepare_site_tensor(X), _prepare_site_tensor(Y) + + # ---- hairs (size D^2) with a 1D auxiliary leg ---- + + @tensor top[-1 -2; -3 -4] := cor_nw(ms[-2, 0])[1 2 -1 -2] * cor_ne(ms[-2, 1])[-3 -4 1 2] + tl, tr = _svd_cut!(top) + + @tensor bot[-1 -2; -3 -4] := cor_sw(ms[2, 0])[-1 -2 1 2] * cor_se(ms[2, 1])[-3 -4 1 2] + bl, br = _svd_cut!(bot) + + nw = permute(cor_nw(ms[-1, -1]), ((3, 4), (1, 2))) + nw1, nw2 = _svd_cut!(nw) + + ne = permute(cor_ne(ms[-1, 2]), ((3, 4), (1, 2))) + ne1, ne2 = _svd_cut!(ne) + + sw = permute(cor_sw(ms[1, -1]), ((1, 2), (3, 4))) + sw1, sw2 = _svd_cut!(sw) + + se = permute(cor_se(ms[1, 2]), ((3, 4), (1, 2))) + se1, se2 = _svd_cut!(se) + + @tensoropt hW[DXw1 DXw0] := + hair_w(ms[0, -2])[Dw21 Dw20] * + nw1[Dnw11 Dnw10] * sw1[Dsw11 Dsw10] * + twistdual(ms[0, -1], 1)[phW Dnw10 DXw0 Dsw10 Dw20] * + conj(ms[0, -1][phW Dnw11 DXw1 Dsw11 Dw21]) + @tensoropt hE[DYe1 DYe0] := + hair_e(ms[0, 3])[De21 De20] * + ne2[Dne21 Dne20] * se2[Dse21 Dse20] * + twistdual(ms[0, 2], 1)[phE Dne20 De20 Dse20 DYe0] * + conj(ms[0, 2][phE Dne21 De21 Dse21 DYe1]) + @tensoropt NW[Dn1 Dn0 DXn1 DXn0] := + tl[Dtl1 Dtl0] * nw2[Dnw21 Dnw20] * + twistdual(ms[-1, 0], 1)[pNW Dtl0 Dn0 DXn0 Dnw20] * + conj(ms[-1, 0][pNW Dtl1 Dn1 DXn1 Dnw21]) + @tensoropt NE[DYn1 DYn0 Dn1 Dn0] := + tr[Dtr1 Dtr0] * ne1[Dne11 Dne10] * + twistdual(ms[-1, 1], 1)[pNE Dtr0 Dne10 DYn0 Dn0] * + conj(ms[-1, 1][pNE Dtr1 Dne11 DYn1 Dn1]) + @tensoropt SW[DXs1 DXs0 Ds1 Ds0] := + bl[Dbl1 Dbl0] * sw2[Dsw21 Dsw20] * + twistdual(ms[1, 0], 1)[pSW DXs0 Ds0 Dbl0 Dsw20] * + conj(ms[1, 0][pSW DXs1 Ds1 Dbl1 Dsw21]) + @tensoropt SE[DYs1 DYs0 Ds1 Ds0] := + br[Dbr1 Dbr0] * se1[Dse11 Dse10] * + twistdual(ms[1, 1], 1)[pSE DYs0 Dse10 Dbr0 Ds0] * + conj(ms[1, 1][pSE DYs1 Dse11 Dbr1 Ds1]) + return _contract_ntu_NNEnv(X, Y, hW, hE, NW, NE, SW, SE; allocator) +end + +# Common NN/NN+ network; each virtual bond has a bra and a ket index. +# NW -- NE +# | | +# hW --- X Y --- hE +# | | +# SW -- SE +# X/Y identify the central site; n/e/s/w identify its virtual leg. +# N/S join the two upper/lower corners; pX/pY are physical indices. +# Suffix 0 denotes ket and 1 denotes bra. +# The four open QR indices are Xq1, Yq1, Xq0, Yq0. +const _NTU_ENV_NETWORK = [ + [:Xw1, :Xw0], + [:Ye1, :Ye0], + [:N1, :N0, :Xn1, :Xn0], + [:Yn1, :Yn0, :N1, :N0], + [:Xs1, :Xs0, :S1, :S0], + [:Ys1, :Ys0, :S1, :S0], + [:pX, :Xn1, :Xq1, :Xs1, :Xw1], + [:pX, :Xn0, :Xq0, :Xs0, :Xw0], + [:pY, :Yn1, :Ye1, :Ys1, :Yq1], + [:pY, :Yn0, :Ye0, :Ys0, :Yq0], +] + +const _NTU_ORDER_CACHE = Dict{Tuple, Tuple}() +const _NTU_ORDER_CACHE_LOCK = ReentrantLock() + +""" +Choose the common environment's contraction order using actual leg dimensions. +The fixed ten-tensor network is independent of the MPO path and its length. +Floating-point costs avoid integer overflow; dense dimensions approximate symmetry-block costs. +Orders are cached by leg dimensions. +""" +function _ntu_contraction_order(tensors::NamedTuple) + dims = Tuple(dim(space(t, axis)) for t in tensors for axis in 1:numind(t)) + return lock(_NTU_ORDER_CACHE_LOCK) do + order = get(_NTU_ORDER_CACHE, dims, nothing) + isnothing(order) || return order + # Bound the number of retained dimension signatures. + length(_NTU_ORDER_CACHE) >= 64 && empty!(_NTU_ORDER_CACHE) + costs = Dict( + label => Float64(dim(space(t, axis))) + for (t, labels) in zip(tensors, _NTU_ENV_NETWORK) + for (axis, label) in enumerate(labels) + ) + tree, _ = TensorOperations.optimaltree(_NTU_ENV_NETWORK, costs) + return _NTU_ORDER_CACHE[dims] = Tuple(first(TensorOperations.tree2indexorder(tree, _NTU_ENV_NETWORK))) + end +end + +""" +Compile the selected order to ordinary `@tensor` contractions with statically known ranks. +Specializations depend on the contraction order and tensor types, not on leg dimensions. +Only temporaries use `allocator`; the returned environment owns its storage. +""" +@generated function _contract_ntu_kernel( + ::Val{Order}, tensors::NamedTuple, allocator + ) where {Order} + order = Expr(:tuple, Order...) + return macroexpand( + @__MODULE__, :( + @tensor allocator = allocator order = $order benv[Xq1 Yq1; Xq0 Yq0] := + tensors.hW[Xw1 Xw0] * tensors.hE[Ye1 Ye0] * + tensors.NW[N1 N0 Xn1 Xn0] * tensors.NE[Yn1 Yn0 N1 N0] * + tensors.SW[Xs1 Xs0 S1 S0] * tensors.SE[Ys1 Ys0 S1 S0] * + conj(tensors.Xbra[pX Xn1 Xq1 Xs1 Xw1]) * + tensors.Xket[pX Xn0 Xq0 Xs0 Xw0] * + conj(tensors.Ybra[pY Yn1 Ye1 Ys1 Yq1]) * + tensors.Yket[pY Yn0 Ye0 Ys0 Yq0] + ) + ) +end + +"Contract an NN or NN+ boundary while keeping the central bra and ket layers separate." +function _contract_ntu_NNEnv( + X::PEPSTensor, Y::PEPSTensor, hW, hE, NW, NE, SW, SE; + allocator = DefaultAllocator() + ) + # Bra factors are conjugated in the kernel; ket factors carry the physical-leg twist. + tensors = (; + hW, hE, NW, NE, SW, SE, Xbra = X, Xket = twistdual(X, 1), + Ybra = Y, Yket = twistdual(Y, 1), + ) + order = _ntu_contraction_order(tensors) + # The runtime-selected kernel always returns the same concrete rank-(2,2) type. + T = tensormaptype(spacetype(X), 2, 2, TensorKit.promote_storagetype(tensors...)) + benv = _contract_ntu_kernel(Val(Tuple(order)), tensors, allocator)::T + return normalize!(benv, Inf) +end diff --git a/src/algorithms/contractions/bondenv/benv_tools.jl b/src/algorithms/contractions/bondenv/benv_tools.jl index 7cd5b368a..967fea06f 100644 --- a/src/algorithms/contractions/bondenv/benv_tools.jl +++ b/src/algorithms/contractions/bondenv/benv_tools.jl @@ -1,3 +1,10 @@ +#= +The construction of bond environment for Neighborhood Tensor Update (NTU) +is adapted from YASTN (https://github.com/yastn/yastn). +Copyright 2024 The YASTN Authors. All Rights Reserved. +Licensed under the Apache License, Version 2.0 +=# + const BondEnv{T, S} = AbstractTensorMap{T, S, 2, 2} where {T <: Number, S <: ElementarySpace} const BondEnv3site{T, S} = AbstractTensorMap{T, S, 4, 4} where {T <: Number, S <: ElementarySpace} const Hair{T, S} = AbstractTensor{T, S, 2} where {T <: Number, S <: ElementarySpace} @@ -11,3 +18,170 @@ _prepare_site_tensor(t::PEPSTensor) = t _prepare_site_tensor(t::PEPOTensor) = first(fuse_physicalspaces(t)) _prepare_site_tensor(t::PEPSOrth) = permute(insertleftunit(t, 1), ((1,), (2, 3, 4, 5))) _prepare_site_tensor(t::PEPOOrth) = permute(t, ((1,), (2, 3, 4, 5))) + +""" +Extract tensors in an InfinitePEPS or 1-layer InfinitePEPO relative to `(row, col)` +in the original `state`. Neighbor offsets and returned tensors use the frame where +the bond direction `direction` points east. Only distinct requested neighbors are rotated. +""" +function collect_neighbors( + state::InfiniteState, row::Int, col::Int, neighbors::Vector{Tuple{Int, Int}}; + direction::Int = EAST + ) + ucell = size(state)[1:2] + origin = _bond_rotation(CartesianIndex(row, col), direction, EAST, ucell) + rotated_ucell = direction in (NORTH, SOUTH) ? reverse(ucell) : ucell + sites = map(neighbors) do nb + site = _bond_rotation(origin + CartesianIndex(nb), EAST, direction, rotated_ucell) + return CartesianIndex(mod1(site[1], ucell[1]), mod1(site[2], ucell[2])) + end + # Periodic images share a tensor, so rotate and fuse each distinct neighbor only once. + tensors = Dict(site => _prepare_site_tensor(_bond_rotation(state[site], direction, EAST)) for site in unique(sites)) + return Dict(nb => tensors[site] for (nb, site) in zip(neighbors, sites)) +end + +function _hair_axes(open_axs::NTuple{N, Int}) where {N} + @assert all(Base.Fix2(in, 2:5), open_axs) && allunique(open_axs) + return ntuple(Val(4 - N)) do i + n = 0 + for ax in 2:5 + ax in open_axs && continue + n += 1 + n == i && return ax + end + throw(ArgumentError("invalid open axes")) + end +end + +""" + benv_tensor(ket::PEPSTensor, bra::PEPSTensor, open_axs::NTuple{N, Int}) where {N} + benv_tensor(ket::PEPSTensor, bra::PEPSTensor, open_axs::NTuple{N, Int}, hairs::NTuple{Nh, H}) where {N, Nh, H <: Union{Nothing, Hair}} + +Contract the physical axes and the virtual axes of `ket` with `bra` to obtain the tensor on the boundary of the bond environment. +Virtual axes specified by `open_axs` (in ascending order) are not contracted. +Hair tensors can be inserted on contracted legs between `ket` and `bra`. + +# Examples + +- West "hair" tensor (`open_axs = (EAST,)`) + ``` + ╱| + ┌-----ket----- 2 + | ╱ | | + | | | | + | | | ╱ + └---|-bra----- 1 + |╱ + ``` +- Northwest corner tensor (`open_axs = (EAST, SOUTH)`, `hairs = (h, nothing)`) + ``` + ╱| + ┌-----ket----- 2 + | ╱ | h + | 4 | | + | | | + | | ╱ + └-----bra----- 1 + ╱ + 3 + ``` +- West edge tensor (`open_axs = (NORTH, EAST, SOUTH)`) + ``` + 2 + ╱ + ┌-----ket----- 4 + | ╱ | + | 6 | 1 + | | ╱ + └-----bra----- 3 + ╱ + 5 + ``` +""" +function benv_tensor( + ket::PEPSTensor, bra::PEPSTensor, open_axs::NTuple{N, Int} + ) where {N} + # no hair tensors to be attached to virtual legs + return benv_tensor(ket, bra, open_axs, ntuple(Returns(nothing), 4 - N)) +end +function benv_tensor( + ket::PEPSTensor, bra::PEPSTensor, + open_vaxs::NTuple{N, Int}, hairs::NTuple{Nh, Union{Nothing, H}} + ) where {N, Nh, H <: Hair} + @assert 1 <= N <= 3 && Nh == 4 - N + open_axs = open_vaxs .+ 1 + # axes to be contracted + hair_axs = _hair_axes(open_axs) + # attach hairs to ket + ket = twistdual!(copy(ket), 1) + axes = ntuple(identity, Val(5)) + for (h, ax) in zip(hairs, hair_axs) + twistdual!(ket, ax) + h === nothing && continue + axes, biperm = _permute_to_first(axes, ax) + # TODO: choose a better axis convention for hairs + ket = permute(h, ((1,), (2,))) * permute(ket, biperm) + end + perm_back = invperm(axes) + # combine bra and ket + cont_axs = (1, hair_axs...) + pbra = (open_axs, cont_axs) + pket = ( + map(p -> perm_back[p], cont_axs), + map(p -> perm_back[p], open_axs), + ) + pbraket = (ntuple(j -> isodd(j) ? (j + 1) ÷ 2 : (j ÷ 2) + N, 2N), ()) + return tensorcontract(bra, pbra, true, ket, pket, false, pbraket) +end + +#= Free axes of different boundary tensors +(C/E/H mean corner/edge/hair) + + H_n + | + + C_nw - - E_n - - C_ne + | | | + + | | | + H_w - E_w - - ket - - E_e - H_e + | | | + + | | | + C_sw - - E_s - - C_se + + | + H_s +=# +# construction of hairs +for (dir, open_axs) in [:n => (SOUTH,), :e => (WEST,), :s => (NORTH,), :w => (EAST,)] + fname = Symbol("hair_", dir) + @eval begin + $(fname)(ket) = benv_tensor(ket, ket, $open_axs) + $(fname)(ket, h1, h2, h3) = benv_tensor(ket, ket, $open_axs, (h1, h2, h3)) + end +end + +# construction of corners +for (dir, open_axs) in [ + :nw => (EAST, SOUTH), :ne => (SOUTH, WEST), + :se => (NORTH, WEST), :sw => (NORTH, EAST), + ] + fname = Symbol("cor_", dir) + @eval begin + $(fname)(ket) = benv_tensor(ket, ket, $open_axs) + $(fname)(ket, h1, h2) = benv_tensor(ket, ket, $open_axs, (h1, h2)) + end +end + +# construction of edges +for (dir, open_axs) in [ + :n => (EAST, SOUTH, WEST), :e => (NORTH, SOUTH, WEST), + :s => (NORTH, EAST, WEST), :w => (NORTH, EAST, SOUTH), + ] + fname = Symbol("edge_", dir) + @eval begin + $(fname)(ket) = benv_tensor(ket, ket, $open_axs) + $(fname)(ket, h) = benv_tensor(ket, ket, $open_axs, (h,)) + end +end diff --git a/src/algorithms/contractions/bondenv/gaugefix.jl b/src/algorithms/contractions/bondenv/gaugefix.jl index 19c982b35..ec384a002 100644 --- a/src/algorithms/contractions/bondenv/gaugefix.jl +++ b/src/algorithms/contractions/bondenv/gaugefix.jl @@ -108,6 +108,15 @@ to `X`. For example, when `X` is a `PEPSTensor`, function _fixgauge_benvX(X::PEPSOrth, Rinv::MPSBondTensor) return @plansor X[-1 -2 -3 -4] := X[-1 1 -3 -4] * Rinv[1; -2] end +function _fixgauge_benvX(X::PEPSTensor, Rinv::MPSBondTensor) + return @plansor X[-1; -2 -3 -4 -5] := X[-1; -2 1 -4 -5] * Rinv[1; -3] +end +function _fixgauge_benvX(X::PEPOOrth, Rinv::MPSBondTensor) + return @plansor X[-1 -2 -3 -4 -5] := X[-1 -2 1 -4 -5] * Rinv[1; -3] +end +function _fixgauge_benvX(X::PEPOTensor, Rinv::MPSBondTensor) + return @plansor X[-1 -2; -3 -4 -5 -6] := X[-1 -2; -3 1 -5 -6] * Rinv[1; -4] +end """ Apply the gauge transformation `Linv` for `Z` @@ -128,3 +137,12 @@ to `Y`. For example, when `Y` is a `PEPSTensor`, function _fixgauge_benvY(Y::PEPSOrth, Linv::MPSBondTensor) return @plansor Y[-1 -2 -3 -4] := Y[-1 -2 -3 1] * Linv[1; -4] end +function _fixgauge_benvY(Y::PEPSTensor, Linv::MPSBondTensor) + return @plansor Y[-1; -2 -3 -4 -5] := Y[-1; -2 -3 -4 1] * Linv[1; -5] +end +function _fixgauge_benvY(Y::PEPOOrth, Linv::MPSBondTensor) + return @plansor Y[-1 -2 -3 -4 -5] := Y[-1 -2 -3 -4 1] * Linv[1; -5] +end +function _fixgauge_benvY(Y::PEPOTensor, Linv::MPSBondTensor) + return @plansor Y[-1 -2; -3 -4 -5 -6] := Y[-1 -2; -3 -4 -5 1] * Linv[1; -6] +end diff --git a/src/algorithms/time_evolution/apply_gate.jl b/src/algorithms/time_evolution/apply_gate.jl index 5ed1c4cf8..b2a0ecce2 100644 --- a/src/algorithms/time_evolution/apply_gate.jl +++ b/src/algorithms/time_evolution/apply_gate.jl @@ -49,56 +49,7 @@ function _orient_bond(state::InfiniteState, sites::Vector{CartesianIndex{2}}, ga end """ - _bond_qr(T, k, gate_ax, [allocator]) -> Q, R, p - -Split `T = Q * R`, where the reduced tensor `R : x ← p ⊗ v` holds the physical axis `gate_ax` -and the virtual axis `k` (the bond) of `T`, while `Q` holds all other axes. `p` is the -permutation that brings the product back to the leg order of `T`. For example, for a -PEPSTensor `T[p; n e s w]` with the bond `v = e` on its east axis (`k = 2`), -``` - n n - | | - w - T - v = w - Q ← x ← R - v - | ↘ | ↘ - s p s p -``` -For a PEPOTensor, the physical axis that is not `gate_ax` stays in `Q`. -The permuted copy of `T` that is factorized is a temporary allocated with `allocator`. -""" -function _bond_qr( - T::Union{PEPSTensor, PEPOTensor}, k::Int, gate_ax::Int, allocator = DefaultAllocator() - ) - Np, M = numout(T), numind(T) - r = (gate_ax, Np + k) - q = TupleTools.deleteat(ntuple(identity, M), r) - cp = allocator_checkpoint!(allocator) - T′ = tensoralloc_add(scalartype(T), T, (q, r), false, Val(true), allocator) - tensoradd!(T′, T, (q, r), false, One(), Zero(), DefaultBackend(), allocator) - Q, R = left_orth!(T′; positive = true) - tensorfree!(T′, allocator) - allocator_reset!(allocator, cp) - return Q, R, _inv_mpo_perm((q..., r...), Val(Np)) -end - -""" - _undo_bond_qr(Q, X, pX, p, [allocator]) - -Inverse of [`_bond_qr`](@ref): contract `Q` with `X`, whose legs selected by `pX` are -ordered as `x ← p ⊗ v`, into the leg order specified by `p`. -""" -function _undo_bond_qr( - Q::AbstractTensorMap, X::AbstractTensorMap, pX, p, allocator = DefaultAllocator() - ) - N = numind(Q) - pQ = (ntuple(identity, N - 1), (N,)) - cp = allocator_checkpoint!(allocator) - T = tensorcontract(Q, pQ, false, X, pX, false, p, One(), DefaultBackend(), allocator) - allocator_reset!(allocator, cp) - return T -end - -""" - _apply_bond_gate(R₁, R₂, gate, trunc, [allocator]) -> U, S, Vᴴ, ϵ + _apply_bond_gate(R₁, R₂, gate, trunc, [allocator]) -> R₁, S, R₂, ϵ Apply the 2-site `gate` to the reduced tensors `R₁ : x₁ ← p₁ ⊗ v` and `R₂ : x₂ ← p₂ ⊗ v'` ``` @@ -115,7 +66,7 @@ a PEPO), and split the result as `U S Vᴴ` with a truncated SVD. The square roo absorbed into `U : x₁ ⊗ p₁ ← v′` and `Vᴴ : v′ ← x₂ ⊗ p₂`. """ function _apply_bond_gate( - R₁::AbstractTensorMap, R₂::AbstractTensorMap, gate::NNGate, + R₁::ReducedBondTensor, R₂::ReducedBondTensor, gate::NNGate, trunc::TruncationStrategy, allocator = DefaultAllocator() ) θ = if isdual(space(R₁, 2)) @@ -127,5 +78,12 @@ function _apply_bond_gate( end U, S, Vᴴ, ϵ = svd_trunc!(θ; trunc) absorb_s!(U, S, Vᴴ) - return U, S, Vᴴ, ϵ + if isdual(space(R₂, 3)) + U, S, Vᴴ = flip(U, numind(U)), _fliptwist_s!(S), flip(Vᴴ, 1) + end + R₁ = permute(U, ((1,), (2, 3))) + # TODO: Avoid copying Vᴴ into the reduced-tensor layout when reconstruction can + # consume its SVD layout directly. + R₂ = permute(Vᴴ, ((2,), (3, 1))) + return R₁, S, R₂, ϵ end diff --git a/src/algorithms/time_evolution/get_cluster.jl b/src/algorithms/time_evolution/get_cluster.jl index da94c9449..55daf6d47 100644 --- a/src/algorithms/time_evolution/get_cluster.jl +++ b/src/algorithms/time_evolution/get_cluster.jl @@ -72,6 +72,52 @@ function _nn_bonddir(site1::CartesianIndex{2}, site2::CartesianIndex{2}) end end +""" +Apply the tensor rotation that maps a bond from orientation `dir_from` +to orientation `dir_to`, where the directions are one of +`NORTH`, `EAST`, `SOUTH`, `WEST`. +""" +function _bond_rotation(x, dir_from::Int, dir_to::Int) + delta = mod(dir_to - dir_from, 4) + return if delta == 0 + x + elseif delta == 1 + rotr90(x) + elseif delta == 2 + rot180(x) + else # delta == 3 + rotl90(x) + end +end + +# Endpoint orthogonal factors keep their remaining physical axes before the four virtual axes. +function _bond_rotation(X::Union{PEPSOrth, PEPOOrth}, dir_from::Int, dir_to::Int) + axes = allind(X) + physical = axes[1:(end - 4)] + virtual = TupleTools.circshift(axes[(end - 3):end], dir_to - dir_from) + return permute(X, (physical..., virtual...)) +end + +""" +Map a site coordinate from the lattice frame where the bond has orientation +`dir_from` to the frame where the bond has orientation `dir_to`. +""" +function _bond_rotation( + site::CartesianIndex{2}, dir_from::Int, dir_to::Int, + unitcell::NTuple{2, Int} + ) + delta = mod(dir_to - dir_from, 4) + return if delta == 0 + site + elseif delta == 1 + siterotr90(site, unitcell) + elseif delta == 2 + siterot180(site, unitcell) + else # delta == 3 + siterotl90(site, unitcell) + end +end + """ Obtain the cluster `Ms` along the (open) path `sites` in `state`. @@ -120,3 +166,13 @@ function _get_cluster( Ms = map(site -> state[site], sites) return Ms, open_vaxs, perms end + +""" +Given a vector `Ms` of `AbstractTensorMap`s, +apply permutations `perms` to each tensor. +""" +function _permute_cluster(Ms::Vector{<:AbstractTensorMap}, perms::Vector{<:Tuple}) + return map(Ms, perms) do M, perm + return permute(M, perm) + end +end diff --git a/src/algorithms/time_evolution/ntupdate.jl b/src/algorithms/time_evolution/ntupdate.jl new file mode 100644 index 000000000..e6f38736e --- /dev/null +++ b/src/algorithms/time_evolution/ntupdate.jl @@ -0,0 +1,289 @@ +""" +$(TYPEDEF) + +Algorithm struct for neighbourhood tensor update (NTU) of InfinitePEPS or InfinitePEPO. + +## Fields + +$(TYPEDFIELDS) + +Reference: +- Physical Review B 104, 094411 (2021) +- Physical Review B 106, 195105 (2022) +""" +@kwdef struct NeighbourUpdate{ + TR <: Union{ALSTruncation, FullEnvTruncation}, + BE <: NeighbourEnv, + } <: TimeEvolution + "Bond truncation algorithm after applying time evolution gate" + opt_alg::TR = ALSTruncation(; trunc = FixedSpaceTruncation()) + "When true (or false), the Trotter gate is `exp(-H dt)` (or `exp(-iH dt)`)" + imaginary_time::Bool = true + "Algorithm to construct NTU bond environment." + bondenv_alg::BE = NNEnv() + "When true, fix gauge of bond environment" + fixgauge::Bool = true + "When true, assume bipartite unit cell structure" + bipartite::Bool = false +end + +""" +Internal state of neighbourhood tensor update algorithm +""" +struct NTUState{S <: InfiniteState, N <: Number} + "number of performed iterations" + iter::Int + "evolved time" + t::N + "PEPS/PEPO" + psi::S + # TODO: include BPEnv when incorporating BP +end + +function TimeEvolver(alg::NeighbourUpdate, dt::Number, nstep::Int, circuit, state::NTUState) + allocator = default_allocator(state.psi, SerialScheduler()) + return TimeEvolver(alg, dt, nstep, circuit, state, allocator) +end + +""" + TimeEvolver( + psi0::InfiniteState, H::LocalOperator, dt::Number, nstep::Int, + alg::NeighbourUpdate; t0::Number = 0.0, symmetrize_gates::Bool = false + ) + +Initialize a TimeEvolver with Hamiltonian `H` and neighbourhood tensor update `alg`, +starting from the initial state `psi0`. + +- The initial time is specified by `t0`. +- Use `symmetrize_gates = true` for second-order Trotter decomposition. +- Contraction scratch is owned by the evolver and reused across timesteps and resets. +""" +function TimeEvolver( + psi0::InfiniteState, H::LocalOperator, dt::Number, nstep::Int, + alg::NeighbourUpdate; t0::Number = 0.0, symmetrize_gates::Bool = false + ) + _timeevol_sanity_check(psi0, physicalspace(H), alg) + dt′ = _get_dt(psi0, dt, alg.imaginary_time) + gate = trotterize(H, dt′; symmetrize_gates) + # convert FixedSpaceTruncation to site-dependent `truncspace`s + if alg.opt_alg.trunc isa FixedSpaceTruncation + trunc = _get_fixedspacetrunc(psi0) + @reset alg.opt_alg.trunc = trunc + end + state = NTUState(0, t0, psi0) + return TimeEvolver(alg, dt, nstep, gate, state) +end + +""" +Neighbourhood tensor update optimized for nearest neighbor gates +utilizing reduced bond tensors with the physical leg. +""" +function _ntu_iter!( + state::InfiniteState, gate::NNGate, wts::SUWeight, + sites::Vector{CartesianIndex{2}}, alg::NeighbourUpdate; + allocator = DefaultAllocator() + ) + @assert length(sites) == 2 + trunc = only(_get_cluster_trunc(alg.opt_alg.trunc, sites)) + alg′ = (@set alg.opt_alg.trunc = trunc) + return _bond_truncate!(state, wts, Tuple(sites), (:first, :last), alg′; gate, allocator) +end + +""" +One round of neighbourhood tensor update, preserving the input state. +""" +function ntu_iter( + state::InfiniteState, circuit::LocalCircuit, alg::NeighbourUpdate; + reverse_trunc::Bool = false, + allocator = default_allocator(state, SerialScheduler()) + ) + state2, wts = copy(state), SUWeight(state) + info = (; fid = 1.0) + for (sites, gate) in circuit.gates + if length(sites) == 1 + # 1-site gate + # TODO: special treatment for bipartite state + site = only(sites) + state2[site] = _apply_sitegate(state2[site], gate) + info′ = (; fid = 1.0) + elseif length(sites) == 2 + (dir, r, c) = _nn_bonddir(sites...) + d = dir in (EAST, WEST) ? 1 : 2 + alg.bipartite && iseven(r) && continue + state2, wts, info′ = _ntu_iter!(state2, gate, wts, sites, alg; allocator) + (!alg.bipartite) && continue + if d == 1 + state2[r + 1, c + 1] = copy(state2[r, c]) + state2[r + 1, c] = copy(state2[r, c + 1]) + wts[1, r + 1, c + 1] = copy(wts[1, r, c]) + else + state2[r - 1, c - 1] = copy(state2[r, c]) + state2[r, c - 1] = copy(state2[r - 1, c]) + wts[2, r - 1, c - 1] = copy(wts[2, r, c]) + end + else + # N-site MPO gate (N ≥ 2) + alg.bipartite && error("Multi-site MPO gates are not compatible with bipartite states.") + state2, wts, info′ = _ntu_iter!(state2, gate, wts, sites, alg; reverse_trunc, allocator) + end + # record the worst fidelity + (info′.fid < info.fid) && (info = info′) + end + return state2, wts, info +end + +function Base.iterate(it::TimeEvolver{<:NeighbourUpdate}, state = it.state) + iter, t = state.iter, state.t + (iter == it.nstep) && return nothing + # reverse bond truncation order at even steps + # (isodd is used because iter is updated after applying the gate) + reverse_trunc = isodd(iter) + psi, wts, info = ntu_iter(state.psi, it.circuit, it.alg; reverse_trunc, allocator = it.allocator) + iter, t = iter + 1, t + it.dt + # update internal state + it.state = NTUState(iter, t, psi) + info = (; (; t, wts)..., info...) + return (psi, info), it.state +end + +""" + timestep( + it::TimeEvolver{<:NeighbourUpdate}, psi::InfiniteState; + iter::Int = it.state.iter, t::Number = it.state.t + ) -> (psi, info) + +Given the TimeEvolver iterator `it`, perform one step of NTU time evolution +on the input state `psi`. + +- Using `iter` and `t` to reset the current iteration number and evolved time + respectively of the TimeEvolver `it`. +""" +function MPSKit.timestep( + it::TimeEvolver{<:NeighbourUpdate}, psi::InfiniteState; + iter::Int = it.state.iter, t::Number = it.state.t + ) + _timeevol_sanity_check(psi, physicalspace(it.state.psi), it.alg) + state = NTUState(iter, t, psi) + result = iterate(it, state) + if result === nothing + @warn "TimeEvolver `it` has already reached the end." + return nothing + else + return first(result) + end +end + +""" + time_evolve(it; verbosity = 2, check_interval = 10) -> (psi, info) + time_evolve(it, H, env, ctm_alg; tol = 1.0e-7, verbosity = 2, check_interval = 10) -> (psi, info) + +Perform time evolution to the end of `NeighbourUpdate` TimeEvolver `it`, +or until convergence of energy set by a positive `tol`. + +To enable convergence check (for imaginary time evolution of InfinitePEPS only), +provide the Hamiltonian `H`, CTMRG environment `env`, CTMRG algorithm `ctm_alg` +and setting `tol > 0`. + +`check_interval` sets the number of iterations between energy checks +(for ground state search) and outputs of information. +""" +function MPSKit.time_evolve( + it::TimeEvolver{<:NeighbourUpdate}; + verbosity::Int = 2, check_interval::Int = 10 + ) + return with_verbosity(; verbosity) do + time_start = time0 = time() + @log_initialization "--- Time evolution (neighbourhood tensor update), dt = $(it.dt) ---" + info0 = nothing + for (psi, info) in it + iter = it.state.iter + stop = (iter == it.nstep) + showinfo = (iter == 1) || (iter % check_interval == 0) || stop + time1 = time() + if showinfo + Δλ = (info0 === nothing) ? NaN : compare_weights(info.wts, info0.wts) + @log_convergence @sprintf( + "NTU iter %d: t = %.2e, |Δλ| = %.3e. Time: %.2f s", + it.state.iter, it.state.t, Δλ, time1 - time0 + ) + end + if stop + time_end = time() + @log_convergence @sprintf("Time evolution finished in %.2f s", time_end - time_start) + return psi, info + end + info0, time0 = info, time() + end + end +end + +function MPSKit.time_evolve( + it::TimeEvolver{<:NeighbourUpdate, G, S}, + H::LocalOperator, env::CTMRGEnv, ctm_alg::CTMRGAlgorithm; + tol::Float64 = 1.0e-7, verbosity::Int = 2, check_interval::Int = 10 + ) where {G, S <: NTUState{<:InfinitePEPS}} + return with_verbosity(; verbosity) do + @log_initialization "--- Time evolution (neighbourhood tensor update), dt = $(it.dt) ---" + time_start = time0 = time() + psi0 = copy(it.state.psi) + @assert it.alg.imaginary_time "Only imaginary time evolution of InfinitePEPS allows convergence checking." + # initial energy + env, = leading_boundary(env, psi0, ctm_alg) + energy = real(expectation_value(psi0, H, env)) / prod(size(psi0)) + @log_convergence @sprintf("NTU iter 0: E = %.4e", energy) + info0 = (; energy, env) + # start evolving + energy0, ΔE = energy, 0.0 + iter0, t0 = it.state.iter, it.state.t + for (psi, info) in it + iter = it.state.iter + showinfo = (iter == 1) || (iter % check_interval == 0) || (iter == it.nstep) + !showinfo && continue + # bond weight change + # TODO: change to BPEnv change when incorporating BP + Δλ = hasproperty(info0, :wts) ? compare_weights(info.wts, info0.wts) : NaN + # reconverge environment + if !all(space(t) == space(t0) for (t, t0) in zip(psi.A, psi0.A)) + # recreate `env` from bond weights if psi virtual space changed + env = CTMRGEnv(info.wts) + end + env, = leading_boundary(env, psi, ctm_alg) + # measure energy + energy = real(expectation_value(psi, H, env)) / prod(size(psi)) + ΔE = energy - energy0 + info = @insert info.energy = energy + info = @insert info.env = env + # show information + time1 = time() + @log_convergence @sprintf( + "NTU iter %-6d: E = %.5f, ΔE = %.3e, |Δλ| = %.3e. Time: %.2f s", + it.state.iter, energy, ΔE, Δλ, time1 - time0 + ) + # determine whether to stop evolution + stop = false + if (ΔE <= 0 && abs(ΔE) < tol) + stop = true + @log_convergence "NTU: energy has converged." + end + if ΔE > 0 + stop = true + @warn "NTU: energy has increased. Abort evolution and return results from last check." + psi, info, energy = psi0, info0, energy0 + it.state = NTUState(iter0, t0, psi0) + end + if iter == it.nstep + stop = true + @warn "NTU: reached maximum iteration." + end + if stop + time_end = time() + @log_convergence @sprintf("Time evolution finished in %.2f s", time_end - time_start) + return psi, info + else + iter0, t0 = it.state.iter, it.state.t + psi0, energy0, info0 = psi, energy, info + end + time0 = time() + end + end +end diff --git a/src/algorithms/time_evolution/ntupdate3site.jl b/src/algorithms/time_evolution/ntupdate3site.jl new file mode 100644 index 000000000..5824088ff --- /dev/null +++ b/src/algorithms/time_evolution/ntupdate3site.jl @@ -0,0 +1,133 @@ +function _get_cluster_permute( + state::InfiniteState, sites::Vector{CartesianIndex{2}} + ) + Ms, _, perms = _get_cluster(state, sites) + Np = (state isa InfinitePEPS) ? Val(1) : Val(2) + invperms = map(p -> _inv_mpo_perm(p, Np), perms) + return _permute_cluster(Ms, perms), invperms +end + +""" +Apply an N-site MPO `gate` (N ≥ 2) to the owned NTU working state and weights. +""" +function _ntu_iter!( + state::InfiniteState, gate::Vector{T}, wts::SUWeight, + sites::Vector{CartesianIndex{2}}, alg::NeighbourUpdate; + reverse_trunc::Bool = false, allocator = DefaultAllocator() + ) where {T <: AbstractTensorMap} + truncs = _get_cluster_trunc(alg.opt_alg.trunc, sites) + + Ms, invperms = _get_cluster_permute(state, sites) + flips = [isdual(space(M, 1)) for M in Iterators.drop(Ms, 1)] + _flip_virtuals!(Ms, flips) # flip virtual arrows in `Ms` to ← + + # apply gate MPO without truncation + _apply_gatempo!(Ms, gate; allocator) + _flip_virtuals!(Ms, flips) # restore virtual arrows in `Ms` + for (M, s, invperm) in zip(Ms, sites, invperms) + state[s] = permute(M, invperm) + end + + # truncate each bond sequentially along the path + # if reverse_trunc, truncate in the reverse order + info = (; fid = 1.0) + nbond = length(sites) - 1 + for i in 1:nbond + stype1 = (i == 1) ? :first : :middle + stype2 = (i == nbond) ? :last : :middle + # the bond to be updated + bond_id = reverse_trunc ? (nbond + 1 - i) : i + # the sites involved + bondsites = if reverse_trunc + (sites[bond_id + 1], sites[bond_id]) + else + (sites[bond_id], sites[bond_id + 1]) + end + # truncate the bond + trunc = truncs[bond_id] + alg′ = (@set alg.opt_alg.trunc = trunc) + state, wts, info′ = _bond_truncate!(state, wts, bondsites, (stype1, stype2), alg′; allocator) + # record the worst fidelity + (info′.fid < info.fid) && (info = info′) + end + return state, wts, info +end + +""" +Truncate a nearest neighbor bond between `site1` and `site2` in place, +using the standard x direction `A ← B` for the local environment contractions. +""" +function _bond_truncate!( + state::InfiniteState, wts::SUWeight, + (site1, site2)::NTuple{2, CartesianIndex{2}}, + (stype1, stype2)::NTuple{2, Symbol}, + alg::NeighbourUpdate; gate::Union{NNGate, Nothing} = nothing, + allocator = DefaultAllocator() + ) + # Reduce in the original frame; orient only the orthogonal factors for the environment. + (dir, r, c) = _nn_bonddir(site1, site2) + A = state[site1] + B = state[site2] + dir_opposite = mod1(dir + 2, 4) + + # create bond environment + # MPO-expanded tensors need rank compression; NN gates are applied after reduction. + orth_kwargs = isnothing(gate) ? (; trunc = trunctol(; rtol = 1.0e-12)) : (; positive = true) + X, a = if stype1 == :first + bond_tensor_end(A, dir; orth_kwargs...) + else + @assert stype1 == :middle + bond_tensor_mid(A, dir; orth_kwargs...) + end + Y, b = if stype2 == :last + bond_tensor_end(B, dir_opposite; orth_kwargs...) + else + @assert stype2 == :middle + bond_tensor_mid(B, dir_opposite; orth_kwargs...) + end + X, Y = _bond_rotation(X, dir, EAST), _bond_rotation(Y, dir, EAST) + # Gauge fixing and bond truncation use the MPS-chain layouts on both sides. + a_chain = permute(a, ((1, 2), (3,))) + b_chain = permute(b, ((3, 2), (1,))) + benv = bondenv_ntu(site1[1], site1[2], dir, X, Y, state, alg.bondenv_alg; allocator) + @debug "cond(benv) before gauge fix: $(LinearAlgebra.cond(benv))" + if alg.fixgauge + Z = positive_approx(benv) + Z, a_chain, b_chain, (Linv, Rinv) = fixgauge_benv(Z, a_chain, b_chain) + X = _fixgauge_benvX(X, Rinv) + Y = _fixgauge_benvY(Y, Linv) + benv = Z' * Z + @debug "cond(L) = $(LinearAlgebra.cond(Linv)); cond(R): $(LinearAlgebra.cond(Rinv))" + @debug "cond(benv) after gauge fix: $(LinearAlgebra.cond(benv))" + end + + # (optional) apply the NN gate without truncation + if !(gate === nothing) + a = permute(a_chain, ((1,), (2, 3))) + b = permute(b_chain, ((3,), (2, 1))) + a, _, b, = _apply_bond_gate(a, b, gate, truncerror(; atol = 1.0e-15)) + a_chain = permute(a, ((1, 2), (3,))) + b_chain = permute(b, ((3, 2), (1,))) + end + a_chain, s, b_chain, info = bond_truncate(a_chain, b_chain, benv, alg.opt_alg) + + X, Y = _bond_rotation(X, EAST, dir), _bond_rotation(Y, EAST, dir) + a = permute(a_chain, ((1,), (2, 3))) + A = if stype1 == :first + undo_bond_tensor_end(a, X, dir) + else + undo_bond_tensor_mid(a, X, dir) + end + b = permute(b_chain, ((3,), (2, 1))) + B = if stype2 == :last + undo_bond_tensor_end(b, Y, dir_opposite) + else + undo_bond_tensor_mid(b, Y, dir_opposite) + end + + state[site1] = normalize!(A, Inf) + state[site2] = normalize!(B, Inf) + d = dir in (EAST, WEST) ? 1 : 2 + wts[d, r, c] = normalize!(dir in (WEST, SOUTH) ? transpose(s) : s, Inf) + return state, wts, info +end diff --git a/src/algorithms/time_evolution/simpleupdate.jl b/src/algorithms/time_evolution/simpleupdate.jl index 0e18e43b4..13108c8e8 100644 --- a/src/algorithms/time_evolution/simpleupdate.jl +++ b/src/algorithms/time_evolution/simpleupdate.jl @@ -100,12 +100,12 @@ function _su_iter!( ϵ = 0.0 local S for gate_ax in 1:numout(T₁) - Q₁, R₁, p₁ = _bond_qr(T₁, k₁, gate_ax, allocator) - Q₂, R₂, p₂ = _bond_qr(T₂, k₂, gate_ax, allocator) - U, S, Vᴴ, ϵ′ = _apply_bond_gate(R₁, R₂, gate, trunc, allocator) + X, a = bond_tensor_end(T₁, k₁; gate_ax, allocator, positive = true) + Y, b = bond_tensor_end(T₂, k₂; gate_ax, allocator, positive = true) + a, S, b, ϵ′ = _apply_bond_gate(a, b, gate, trunc, allocator) ϵ = max(ϵ, ϵ′) - T₁ = _undo_bond_qr(Q₁, U, ((1,), (2, 3)), p₁, allocator) - T₂ = _undo_bond_qr(Q₂, Vᴴ, ((2,), (3, 1)), p₂, allocator) + T₁ = undo_bond_tensor_end(a, X, k₁; gate_ax, allocator) + T₂ = undo_bond_tensor_end(b, Y, k₂; gate_ax, allocator) alg.purified && break # only apply gate to 1st physical leg end # remove environment weights diff --git a/src/algorithms/time_evolution/time_evolve.jl b/src/algorithms/time_evolution/time_evolve.jl index c55a04a88..d7bf17315 100644 --- a/src/algorithms/time_evolution/time_evolve.jl +++ b/src/algorithms/time_evolution/time_evolve.jl @@ -15,7 +15,7 @@ Iterator for Trotter-based time evolution of InfinitePEPS or InfinitePEPO. $(TYPEDFIELDS) """ mutable struct TimeEvolver{TE <: TimeEvolution, C, S, N <: Number, A} - "Time evolution algorithm (currently supported: `SimpleUpdate`)" + "Time evolution algorithm" alg::TE "Trotter time step" dt::N diff --git a/src/algorithms/truncation/bond_tensor.jl b/src/algorithms/truncation/bond_tensor.jl index 33df63d79..8be9f9f5a 100644 --- a/src/algorithms/truncation/bond_tensor.jl +++ b/src/algorithms/truncation/bond_tensor.jl @@ -1,115 +1,97 @@ +"Reduced bond factor `R[q; physical bond]`, with one outgoing and two incoming legs." +const ReducedBondTensor{T, S} = AbstractTensorMap{T, S, 1, 2} + +function _bond_tensor_axes(axes::Tuple, dir, gate_ax) + Np = length(axes) - 4 + @assert 1 <= dir <= 4 + @assert 1 <= gate_ax <= Np + reduced = (gate_ax, Np + dir) + return TupleTools.deleteat(axes, reduced), reduced +end + """ -Given the first tensor `A` in the cluster acted on by a gate, -obtain reduced tensor on its next bond. + bond_tensor_end(A, dir; gate_ax = 1, allocator = DefaultAllocator(), kwargs...) -For PEPSTensor, -``` - 1 - | - 4 - X ← 2 1 ← a - 3 - | ↘ - 3 2 -``` -For PEPOTensor, -``` - gate_ax = 1 gate_ax = 2 +Use `left_orth` to factor out the physical axis `gate_ax` and +virtual leg `dir` of `A`, returning `(X, R)` with `R[q; physical bond]`. - 1 2 2 2 - ↘ | | ↘ - 5 - X ← 3 1 ← a - 3 5 - X ← 3 1 ← a - 3 - | ↘ | ↘ - 4 2 4 1 +For a PEPS tensor `A` with `dir = EAST`, `gate_ax = 1`: +``` + n n + │ │ + w ─ A ─ e = w ─ X ← q ← R ─ e + │ ╲ │ ╲ + s p s p ``` """ -function bond_tensor_first(A::PEPSTensor; gate_ax::Integer = 1, kwargs...) - @assert gate_ax == 1 - X, a = left_orth!(permute(A, ((2, 4, 5), (1, 3)); copy = true); kwargs...) - X = permute(X, (1, 4, 2, 3)) - a = permute(a, ((1, 2), (3,))) - return a, X -end -function bond_tensor_first(A::PEPOTensor; gate_ax::Integer = 1, kwargs...) - @assert 1 <= gate_ax <= 2 - X, a = if gate_ax == 1 - left_orth!(permute(A, ((2, 3, 5, 6), (1, 4)); copy = true); kwargs...) - else - left_orth!(permute(A, ((1, 3, 5, 6), (2, 4)); copy = true); kwargs...) - end - X = permute(X, (1, 2, 5, 3, 4)) - a = permute(a, ((1, 2), (3,))) - return a, X +function bond_tensor_end(A::AbstractTensorMap{T, S, Np, 4}, dir::Int; gate_ax::Integer = 1, allocator = DefaultAllocator(), kwargs...) where {T, S, Np} + open, reduced = _bond_tensor_axes(allind(A), dir, gate_ax) + partition = (open, reduced) + cp = allocator_checkpoint!(allocator) + A′ = tensoralloc_add(scalartype(A), A, partition, false, Val(true), allocator) + tensoradd!(A′, A, partition, false, One(), Zero(), DefaultBackend(), allocator) + X, R = left_orth!(A′; kwargs...) + tensorfree!(A′, allocator) + allocator_reset!(allocator, cp) + order = TupleTools.insertafter(codomainind(X), Np + dir - 2, domainind(X)) + # TODO: can we cleanly avoid this? + X = permute(X, order) + return X, R end """ -Undo the decomposition in `bond_tensor_first`. + undo_bond_tensor_end(R, X, dir; gate_ax = 1, allocator = DefaultAllocator()) + +Undo `bond_tensor_end`, reconstructing directly in the original frame. +`R` has layout `R[q; physical bond]`; `X` and the result use the original frame. """ -function undo_bond_tensor_first(a::MPSTensor, X::PEPSOrth; gate_ax::Integer = 1) - @assert gate_ax == 1 - return @tensor A[-1; -2 -3 -4 -5] := X[-2 1 -4 -5] * a[1 -1 -3] -end -function undo_bond_tensor_first(a::MPSTensor, X::PEPOOrth; gate_ax::Integer = 1) - @assert 1 <= gate_ax <= 2 - if gate_ax == 1 - return @tensor A[-1 -2; -3 -4 -5 -6] := X[-2 -3 1 -5 -6] * a[1 -1 -4] - else - return @tensor A[-1 -2; -3 -4 -5 -6] := X[-1 -3 1 -5 -6] * a[1 -2 -4] - end +function undo_bond_tensor_end(R::ReducedBondTensor, X::AbstractTensorMap, dir::Int; gate_ax::Integer = 1, allocator = DefaultAllocator()) + Np = numind(X) - 3 + open, reduced = _bond_tensor_axes(ntuple(identity, Np + 4), dir, gate_ax) + contracted = Np - 1 + dir + free = TupleTools.deleteat(allind(X), contracted) + # Restore the original axes from the contraction's (open..., physical, bond) order. + order = invperm((open..., reduced...)) + output = (order[1:Np], order[(Np + 1):end]) + cp = allocator_checkpoint!(allocator) + A = tensorcontract( + X, (free, (contracted,)), false, R, ((1,), (2, 3)), false, + output, One(), DefaultBackend(), allocator + ) + allocator_reset!(allocator, cp) + return A end """ -Given the last tensor `A` in the cluster acted on by a gate, -obtain reduced tensor on its previous bond. + bond_tensor_mid(A, dir; allocator = DefaultAllocator(), kwargs...) -For PEPSTensor, -``` - 1 - | - 1 - b → 3 4 → Y - 2 - ↘ | - 2 3 -``` -For PEPOTensor, -``` - gate_ax = 1 gate_ax = 2 +Use QR decomposition to factor out virtual leg `dir` of `A`, +returning `(X, R)` with the physical axes kept in `X`. - 1 2 2 2 - ↘ | ↘ | - 1 - b → 3 5 → Y - 3 1 - b → 3 5 → Y - 3 - ↘ | | ↘ - 2 4 4 1 +For a PEPS tensor `A` with `dir = EAST`: ``` + n n + │ │ + w ─ A ─ e = w ─ X ← q ← R ─ e + │ ╲ │ ╲ ╲ + s p s p (1) +``` +The `(1)` leg of `R[q; auxiliary bond]` is a trivial physical leg. """ -function bond_tensor_last(A::PEPSTensor; gate_ax::Integer = 1, kwargs...) - @assert gate_ax == 1 - Y, b = left_orth!(permute(A, ((2, 3, 4), (1, 5)); copy = true); kwargs...) - Y = permute(Y, (1, 2, 3, 4)) - b = permute(b, ((3, 2), (1,))) - return b, Y -end -function bond_tensor_last(A::PEPOTensor; gate_ax::Integer = 1, kwargs...) - @assert 1 <= gate_ax <= 2 - Y, b = if gate_ax == 1 - left_orth!(permute(A, ((2, 3, 4, 5), (1, 6)); copy = true); kwargs...) - else - left_orth!(permute(A, ((1, 3, 4, 5), (2, 6)); copy = true); kwargs...) - end - Y = permute(Y, (1, 2, 3, 4, 5)) - b = permute(b, ((3, 2), (1,))) - return b, Y +function bond_tensor_mid(A::Union{PEPSTensor, PEPOTensor}, dir::Int; allocator = DefaultAllocator(), kwargs...) + # insert the auxiliary trivial physical leg in front of existing ones + # so it can be factored into R using `bond_tensor_end` with `gate_ax = 1` + X, R = bond_tensor_end(insertleftunit(A, Val(1)), dir; allocator, kwargs...) + return permute(X, (codomainind(A), domainind(A))), R end """ -Undo the decomposition in `bond_tensor_last`. + undo_bond_tensor_mid(R, X, dir; allocator = DefaultAllocator()) + +Reconstruct a middle tensor with `undo_bond_tensor_end`, then remove the trivial +physical leg introduced by `bond_tensor_mid`. """ -function undo_bond_tensor_last(b::MPSTensor, Y::PEPSOrth; gate_ax::Integer = 1) - @assert gate_ax == 1 - return @tensor A[-1; -2 -3 -4 -5] := b[-5 -1 1] * Y[-2 -3 -4 1] -end -function undo_bond_tensor_last(b::MPSTensor, Y::PEPOOrth; gate_ax::Integer = 1) - @assert 1 <= gate_ax <= 2 - if gate_ax == 1 - return @tensor A[-1 -2; -3 -4 -5 -6] := b[-6 -1 1] * Y[-2 -3 -4 -5 1] - else - return @tensor A[-1 -2; -3 -4 -5 -6] := b[-6 -2 1] * Y[-1 -3 -4 -5 1] - end +function undo_bond_tensor_mid(R::ReducedBondTensor, X::Union{PEPSTensor, PEPOTensor}, dir::Int; allocator = DefaultAllocator()) + A = undo_bond_tensor_end(R, X, dir; allocator) + return removeunit(A, Val(1)) end diff --git a/src/utility/util.jl b/src/utility/util.jl index 663c4ae50..e7f6b0502 100644 --- a/src/utility/util.jl +++ b/src/utility/util.jl @@ -173,3 +173,12 @@ function random_dual!(Vs::AbstractMatrix{E}; p = 0.7) where {E <: ElementarySpac end return Vs end + +"Move `ax` to the first codomain position, returning the new axis order and permutation." +function _permute_to_first(axes::NTuple{N, Int}, ax::Int) where {N} + domain_axes = TupleTools.deleteat(ntuple(identity, N), ax) + q = invperm(axes) + biperm = ((q[ax],), map(i -> q[i], domain_axes)) + new_axes = (ax, ntuple(i -> axes[biperm[2][i]], N - 1)...) + return new_axes, biperm +end diff --git a/test/Project.toml b/test/Project.toml index 7e019cf9c..ee51f984e 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -18,6 +18,7 @@ NVML = "611af6d1-644e-4c5d-bd58-854d7d1254b9" OptimKit = "77e91f04-9b3b-57a6-a776-40b61faaebe0" PEPSKit = "52969e89-939e-4361-9b68-9bc7cde4bdeb" ParallelTestRunner = "d3525ed8-44d0-4b2c-a655-542cee43accc" +Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" diff --git a/test/bondenv/benv_ntu.jl b/test/bondenv/benv_ntu.jl new file mode 100644 index 000000000..d576cc350 --- /dev/null +++ b/test/bondenv/benv_ntu.jl @@ -0,0 +1,21 @@ +using PEPSKit, CUDA, AMDGPU + +@isdefined(TestSuite) || include("../testsuite/TestSuite.jl") +using .TestSuite + +is_buildkite = get(ENV, "BUILDKITE", "false") == "true" + +if !is_buildkite + TestSuite.bondenv_ntu_neighbors(Vector) + TestSuite.bondenv_ntu(Vector) +end + +if CUDA.functional() + TestSuite.bondenv_ntu_neighbors(CuArray) + TestSuite.bondenv_ntu(CuArray) +end + +if AMDGPU.functional() + TestSuite.bondenv_ntu_neighbors(ROCArray) + TestSuite.bondenv_ntu(ROCArray) +end diff --git a/test/testsuite/TestSuite.jl b/test/testsuite/TestSuite.jl index 3e576386d..5e2e4d44d 100644 --- a/test/testsuite/TestSuite.jl +++ b/test/testsuite/TestSuite.jl @@ -25,6 +25,12 @@ module BondEnvCTM end using .BondEnvCTM +module BondEnvNTU + include("bondenv/benv_ntu.jl") + export bondenv_ntu, bondenv_ntu_neighbors +end +using .BondEnvNTU + module BondEnvGaugeFix include("bondenv/benv_gaugefix.jl") export bondenv_gaugefix @@ -229,10 +235,22 @@ using .TimeEvolClusterProjectors module TimeEvolSiteDepTruncation include("timeevol/sitedep_truncation.jl") - export timeevol_sitedep_rotation, timeevol_sitedep_su + export timeevol_sitedep_rotation, timeevol_sitedep_su, timeevol_sitedep_ntu end using .TimeEvolSiteDepTruncation +module TimeEvolHeisGroundstate + include("timeevol/heis_groundstate.jl") + export timeevol_heis_groundstate +end +using .TimeEvolHeisGroundstate + +module TimeEvolIsingRealtime + include("timeevol/tf_ising_realtime.jl") + export timeevol_ising_realtime +end +using .TimeEvolIsingRealtime + module TimeEvolJ1J2FiniteT include("timeevol/j1j2_finiteT.jl") export timeevol_j1j2_finiteT @@ -241,7 +259,7 @@ using .TimeEvolJ1J2FiniteT module TimeEvolTimestep include("timeevol/timestep.jl") - export timeevol_timestep + export timeevol_timestep, timeevol_timestep_ntu end using .TimeEvolTimestep diff --git a/test/testsuite/bondenv/benv_ctm.jl b/test/testsuite/bondenv/benv_ctm.jl index 8332d276b..2c1289dac 100644 --- a/test/testsuite/bondenv/benv_ctm.jl +++ b/test/testsuite/bondenv/benv_ctm.jl @@ -4,6 +4,7 @@ using PEPSKit using LinearAlgebra using Random using Adapt +using PEPSKit: EAST, WEST Nr, Nc = 2, 2 Envspace = Vect[FermionParity ⊠ U1Irrep]( @@ -41,8 +42,10 @@ function test_benv_ctm(AT, state::Union{InfinitePEPS, InfinitePEPO}) for row in 1:Nr, col in 1:Nc cp1 = col + 1 A, B = state[row, col], state[row, cp1] - a, X = PEPSKit.bond_tensor_first(A) - b, Y = PEPSKit.bond_tensor_last(B) + X, a = PEPSKit.bond_tensor_end(A, EAST) + Y, b = PEPSKit.bond_tensor_end(B, WEST) + a = permute(a, ((1, 2), (3,))) + b = permute(b, ((3, 2), (1,))) benv = PEPSKit.bondenv_ctm(row, col, X, Y, env) Z = PEPSKit.positive_approx(benv) # verify that gauge fixing can greatly reduce diff --git a/test/testsuite/bondenv/benv_ntu.jl b/test/testsuite/bondenv/benv_ntu.jl new file mode 100644 index 000000000..47c5c48f0 --- /dev/null +++ b/test/testsuite/bondenv/benv_ntu.jl @@ -0,0 +1,105 @@ +using Test +using TensorKit +using PEPSKit +using LinearAlgebra +using KrylovKit +using Random +using Adapt +using PEPSKit: NORTH, EAST, SOUTH, WEST + +Nr, Nc = 2, 2 + +"Check positivity and gauge fixing of a neighbourhood bond environment." +function test_ntu_env( + state::Union{InfinitePEPS, InfinitePEPO}, env_alg::Alg + ) where {Alg <: PEPSKit.NeighbourEnv} + @info "Testing $(typeof(env_alg))" + for row in 1:Nr, col in 1:Nc + A, B = state[row, col], state[row, col + 1] + X, a = PEPSKit.bond_tensor_end(A, EAST) + Y, b = PEPSKit.bond_tensor_end(B, WEST) + a = permute(a, ((1, 2), (3,))) + b = permute(b, ((3, 2), (1,))) + @tensor ab[DX DY; da db] := a[DX da D] * b[D db DY] + benv = PEPSKit.bondenv_ntu(row, col, EAST, X, Y, state, env_alg) + # this is a result of `bond_tensor_...` + @assert [isdual(space(benv, ax)) for ax in 1:numind(benv)] == [0, 0, 1, 1] + # NTU bond environments are exact and should be positive definite + @test benv' ≈ benv + benv = project_hermitian(benv) + nrm1 = PEPSKit.inner_prod(benv, ab, ab) + @test nrm1 ≈ real(nrm1) + D, U = eigh_full(benv) + @test minimum(D.data) >= 0 + # gauge fixing + Z = PEPSKit.sdiag_pow(D, 0.5) * U' + @assert benv ≈ Z' * Z + Z2, a2, b2, (Linv, Rinv) = PEPSKit.fixgauge_benv(Z, a, b) + benv2 = Z2' * Z2 + # gauge fixing should reduce condition number + cond = LinearAlgebra.cond(benv) + cond2 = LinearAlgebra.cond(benv2) + @test cond2 <= cond + @info "benv cond number: (gauge-fixed) $(cond2) ≤ $(cond) (initial)" + # verify gauge transformation of X, Y + @tensor a2b2[DX DY; da db] := a2[DX da D] * b2[D db DY] + nrm2 = PEPSKit.inner_prod(benv2, a2b2, a2b2) + X2 = PEPSKit._fixgauge_benvX(X, Rinv) + Y2 = PEPSKit._fixgauge_benvY(Y, Linv) + benv3 = PEPSKit.bondenv_ntu(row, col, PEPSKit.EAST, X2, Y2, state, env_alg) + benv3 *= norm(benv2, Inf) + nrm3 = PEPSKit.inner_prod(benv3, a2b2, a2b2) + @test benv2 ≈ benv3 + @test nrm1 ≈ nrm2 ≈ nrm3 + end + return +end + +"Check neighbor orientation and periodic-image reuse on a rectangular unit cell." +function bondenv_ntu_neighbors(AT) + Random.seed!(3719) + return @testset "Rectangular NTU neighbor mapping ($AT)" begin + for state in ( + InfinitePEPS(randn, ComplexF64, ℂ^2, ℂ^2; unitcell = (2, 3)), + InfinitePEPO(randn, ComplexF64, ℂ^2, ℂ^2; unitcell = (2, 3, 1)), + ) + state = adapt(AT, state) + for direction in (NORTH, EAST, SOUTH, WEST) + rotated = PEPSKit._bond_rotation(state, direction, PEPSKit.EAST) + origin = PEPSKit._bond_rotation(CartesianIndex(2, 3), direction, PEPSKit.EAST, size(state)[1:2]) + nr, nc = size(rotated)[1:2] + neighbors = [(0, 0), (nr, 0), (0, nc), (-2, 3), (1, -1)] + actual = PEPSKit.collect_neighbors(state, 2, 3, neighbors; direction) + expected = Dict(nb => PEPSKit._prepare_site_tensor(rotated[origin + CartesianIndex(nb)]) for nb in neighbors) + @test actual == expected + @test actual[0, 0] === actual[nr, 0] === actual[0, nc] + end + end + end +end + +"Check neighbourhood bond environments and their gauge transformations." +function bondenv_ntu(AT) + Random.seed!(20) + return @testset "NTU bond environment ($(sym), $AT)" for sym in [U1Irrep, FermionParity] + Pspace = Vect[sym](0 => 1, 1 => 1) + V2 = Vect[sym](0 => 4, 1 => 1) + V3 = Vect[sym](0 => 3, 1 => 2) + V4 = Vect[sym](0 => 2, 1 => 2) + V5 = Vect[sym](0 => 2, 1 => 3) + Pspaces = fill(Pspace, (Nr, Nc)) + Nspaces = [V2 V2; V4 V4] + Espaces = [V3 V5; V5 V3] + + for state in [ + InfinitePEPS(randn, ComplexF64, Pspaces, Nspaces, Espaces), + InfinitePEPO(randn, ComplexF64, Pspaces, Nspaces, Espaces), + ] + state = adapt(AT, state) + normalize!.(state.A, Inf) + for env_alg in (NNEnv(), NNpEnv()) + test_ntu_env(state, env_alg) + end + end + end +end diff --git a/test/testsuite/timeevol/heis_groundstate.jl b/test/testsuite/timeevol/heis_groundstate.jl new file mode 100644 index 000000000..3e8594d06 --- /dev/null +++ b/test/testsuite/timeevol/heis_groundstate.jl @@ -0,0 +1,48 @@ +using Test +using Random +using TensorKit +using PEPSKit +using Adapt + +function converge_env(AT, state, ctm_alg) + env0 = adapt(AT, CTMRGEnv(ones, scalartype(state), state, oneunit(spacetype(state)))) + env, = leading_boundary(env0, state, ctm_alg) + return env +end + +function timeevol_heis_groundstate(AT) + return @testset "Heisenberg ground state ($AT)" begin + Random.seed!(1457860) + Nr, Nc = 2, 2 + H = adapt( + AT, j1_j2_model( + Float64, U1Irrep, InfiniteSquare(Nr, Nc); + J1 = 1.0, J2 = 0.0, sublattice = false + ) + ) + Pspace = U1Space(1 // 2 => 1, -1 // 2 => 1) + Vspace = U1Space(0 => 2, 1 // 2 => 1, -1 // 2 => 1) + ψ0 = adapt(AT, InfinitePEPS(rand, Float64, Pspace, Vspace; unitcell = (Nr, Nc))) + trunc = truncerror(; atol = 1.0e-10) & truncrank(4) + ctm_alg = SequentialCTMRG(; tol = 1.0e-10, verbosity = 2, trunc = truncerror(; atol = 1.0e-10) & truncrank(16)) + + # prepare simple update state + su_alg = SimpleUpdate(; trunc) + evolver = TimeEvolver(ψ0, H, 0.01, 5000, su_alg, SUWeight(ψ0)) + ψ0, = time_evolve(evolver, H; tol = 1.0e-8, check_interval = 500) + env0 = converge_env(AT, ψ0, ctm_alg) + e0 = expectation_value(ψ0, H, env0) / (Nr * Nc) + @info "Simple update energy = $(e0)." + + # continue with NTU + @testset "NTU for ground state" begin + ntu_alg = NeighbourUpdate(; bondenv_alg = NNpEnv()) + evolver = TimeEvolver(ψ0, H, 0.01, 100, ntu_alg) + ψ, info = time_evolve(evolver, H, env0, ctm_alg) + env = info.env + e = expectation_value(ψ, H, env) / (Nr * Nc) + @info "NTU energy = $(e)" + @test e < e0 + end + end +end diff --git a/test/testsuite/timeevol/j1j2_finiteT.jl b/test/testsuite/timeevol/j1j2_finiteT.jl index fd32a7924..e2cb336b2 100644 --- a/test/testsuite/timeevol/j1j2_finiteT.jl +++ b/test/testsuite/timeevol/j1j2_finiteT.jl @@ -4,11 +4,7 @@ using TensorKit using PEPSKit using Adapt -# Benchmark energy from high-temperature expansion -# at β = 0.3, 0.6 -# Physical Review B 67, 014416 (2003) -bm = [-0.1235, -0.213] - +"Converge the CTMRG environment for a finite-temperature energy measurement." function converge_env(state, χ::Int) env0 = initialize_ctmrg_environment(state, ProductStateInitialization()) trunc1 = truncrank(χ) & truncerror(; atol = 1.0e-12) @@ -16,50 +12,55 @@ function converge_env(state, χ::Int) return env end +"Check purified simple and neighbourhood updates against high-temperature energies." function timeevol_j1j2_finiteT(AT) - return @testset "J1J2 finiteT ($AT)" begin - Nr, Nc = 2, 2 - ham = adapt( - AT, j1_j2_model( - Float64, SU2Irrep, InfiniteSquare(Nr, Nc); - J1 = 1.0, J2 = 0.5, sublattice = false - ) + # Benchmark energy from high-temperature expansion + # Physical Review B 67, 014416 (2003) + βs = [0.2, 0.4, 0.6] + bm = [-0.08624893, -0.15688984, -0.21300888] + + Nr, Nc = 2, 2 + ham = adapt( + AT, j1_j2_model( + Float64, SU2Irrep, InfiniteSquare(Nr, Nc); + J1 = 1.0, J2 = 0.5, sublattice = false ) - pepo0 = PEPSKit.infinite_temperature_density_matrix(ham) - wts0 = SUWeight(pepo0) + ) + pepo0 = PEPSKit.infinite_temperature_density_matrix(ham) + wts0 = SUWeight(pepo0) + dt, nstep, check_interval = 5.0e-3, 40, 40 + + @testset "Simple update" begin # 7 = 1 (spin-0) + 2 x 3 (spin-1) trunc_pepo = truncrank(7) & truncerror(; atol = 1.0e-12) - check_interval = 100 - dt, nstep = 1.0e-3, 600 - - # PEPO approach - alg = SimpleUpdate(; trunc = trunc_pepo, purified = false) - evolver = TimeEvolver(pepo0, ham, dt, nstep, alg, wts0) - pepo, wts, info = time_evolve(evolver; check_interval) - env = converge_env(InfinitePartitionFunction(pepo), 16) - energy = expectation_value(pepo, ham, env) / (Nr * Nc) - @info "β = $(dt * nstep): tr(ρH) = $(energy)" - @test dt * nstep ≈ info.t - @test energy ≈ bm[2] atol = 5.0e-3 - - # PEPS (purified PEPO) approach alg = SimpleUpdate(; trunc = trunc_pepo, purified = true) - evolver = TimeEvolver(pepo0, ham, dt, nstep, alg, wts0) - pepo, wts, info = time_evolve(evolver; check_interval) - env = converge_env(InfinitePartitionFunction(pepo), 16) - energy = expectation_value(pepo, ham, env) / (Nr * Nc) - @info "β = $(dt * nstep) / 2: tr(ρH) = $(energy)" - @test energy ≈ bm[1] atol = 5.0e-3 - - # test BP gauge fixing for purified iPEPO - bp_alg = BeliefPropagation(; maxiter = 100, tol = 1.0e-9) - bp_env, = leading_boundary(BPEnv(ones, Float64, pepo), pepo, bp_alg) - pepo, = gauge_fix(pepo, BPGauge(), bp_env) + pepo, wts = deepcopy(pepo0), deepcopy(wts0) + for (β, bme) in zip(βs, bm) + t0 = β - βs[1] + pepo, wts, info = time_evolve(pepo, ham, dt, nstep, alg, wts; t0, check_interval) + # measure energy + env = converge_env(InfinitePEPS(pepo), 16) + energy = expectation_value(pepo, ham, pepo, env) / (Nr * Nc) + @info "β = $(info.t): ⟨ρ|H|ρ⟩ = $(energy)" + @test energy ≈ bme atol = 5.0e-3 + end + end - env = converge_env(InfinitePEPS(pepo), 16) - energy = expectation_value(pepo, ham, pepo, env) / (Nr * Nc) - @info "β = $(dt * nstep): ⟨ρ|H|ρ⟩ = $(energy)" - @test dt * nstep ≈ info.t - @test energy ≈ bm[2] atol = 5.0e-3 + @testset "Neighbourhood tensor update" begin + trunc_pepo = truncrank(7) & truncerror(; atol = 1.0e-12) + opt_alg = ALSTruncation(; trunc = trunc_pepo, tol = 1.0e-10) + alg = NeighbourUpdate(; opt_alg, bondenv_alg = NNEnv()) + pepo = deepcopy(pepo0) + for (β, bme) in zip(βs, bm) + t0 = β - βs[1] + evolver = TimeEvolver(pepo, ham, dt, nstep, alg; t0) + pepo, info = time_evolve(evolver; check_interval) + # measure energy + env = converge_env(InfinitePEPS(pepo), 16) + energy = expectation_value(pepo, ham, pepo, env) / (Nr * Nc) + @info "β = $(info.t): ⟨ρ|H|ρ⟩ = $(energy)" + @test energy ≈ bme atol = 5.0e-3 + end end + return nothing end diff --git a/test/testsuite/timeevol/sitedep_truncation.jl b/test/testsuite/timeevol/sitedep_truncation.jl index 7addc0b86..3037b3a9b 100644 --- a/test/testsuite/timeevol/sitedep_truncation.jl +++ b/test/testsuite/timeevol/sitedep_truncation.jl @@ -2,7 +2,7 @@ using Test using Random using TensorKit using PEPSKit -using PEPSKit: _is_bipartite, _get_fixedspacetrunc +using PEPSKit: _is_bipartite, _get_fixedspacetrunc, NORTH, EAST using Adapt elt = Float64 @@ -67,3 +67,38 @@ function timeevol_sitedep_su(AT) end end end + +"Check that neighbourhood updates preserve site-dependent virtual spaces." +function timeevol_sitedep_ntu(AT) + Random.seed!(48736) + states = ( + adapt(AT, InfinitePEPS(randn, elt, Vps, Vns, Ves1)), + adapt(AT, InfinitePEPO(randn, elt, Vps, Vns, Ves2)), + ) + @testset "NTU on $(typeof(state0).name.wrapper), bipartite = $(bipartite)" for + (state0, bipartite) in Iterators.product(states, (false, true)) + J2 = 0.5 + if bipartite + state0[2, 1] = copy(state0[1, 2]) + state0[2, 2] = copy(state0[1, 1]) + J2 = 0.0 + end + ham = adapt(AT, j1_j2_model(elt, U1Irrep, InfiniteSquare(Nr, Nc); J1 = 1.0, J2, sublattice = false)) + # converted internally to SiteDependentTruncation + opt_alg = ALSTruncation(; trunc = FixedSpaceTruncation()) + alg = NeighbourUpdate(; opt_alg, bondenv_alg = NNEnv(), bipartite) + state, info = time_evolve(TimeEvolver(state0, ham, 0.1, 1, alg)) + for (t, t0) in zip(state.A, state0.A) + @test space(t) == space(t0) + end + for (r, c) in Iterators.product(1:Nr, 1:Nc) + @test space(info.wts[1, r, c], 1) == domain(state[r, c], EAST) + @test space(info.wts[2, r, c], 1) == domain(state[r, c], NORTH) + end + if bipartite + @test _is_bipartite(state) + @test _is_bipartite(info.wts) + end + end + return nothing +end diff --git a/test/testsuite/timeevol/tf_ising_realtime.jl b/test/testsuite/timeevol/tf_ising_realtime.jl new file mode 100644 index 000000000..7ddc53a0e --- /dev/null +++ b/test/testsuite/timeevol/tf_ising_realtime.jl @@ -0,0 +1,96 @@ +using Test +using TensorKit +import MPSKitModels: S_zz, σˣ +using PEPSKit +using Adapt +using Printf +using Accessors: @set + +const hc = 3.044382 +const formatter = Printf.Format("t = %.2f, ⟨σˣ⟩ = %.7e + %.7e im.") +# real time evolution of ⟨σx⟩ +# benchmark data from Physical Review B 104, 094411 (2021) +# Figure 6(a) calculated with D = 8 and χ = 32 +const data = [ + # 0.01 9.9920027e-1 + 0.06 9.7274912e-1 + 0.11 9.1973182e-1 + 0.16 8.6230618e-1 + 0.21 8.1894325e-1 + 0.26 8.0003708e-1 + 0.31 8.0081082e-1 + 0.36 8.0979257e-1 + # 0.41 8.1559623e-1 + # 0.46 8.1541661e-1 + # 0.51 8.1274128e-1 +] + +function timeevol_ising_realtime(AT) + return @testset "TF Ising real time ($AT)" begin + # the fully polarized state + peps0 = InfinitePEPS(zeros, ComplexF64, ℂ^2, ℂ^1; unitcell = (2, 2)) + for t in peps0.A + t[1, 1, 1, 1, 1] = 1.0 + t[2, 1, 1, 1, 1] = 1.0 + end + peps0 = adapt(AT, peps0) + lattice = collect(space(t, 1) for t in peps0.A) + + # Hamiltonian + op = adapt(AT, LocalOperator(lattice, [(1, 1)] => σˣ())) + ham = adapt(AT, transverse_field_ising(ComplexF64, Trivial, InfiniteSquare(2, 2); J = 1.0, g = hc)) + + # truncation strategy + Dcut, chi = 4, 16 + trunc_peps = truncerror(; atol = 1.0e-10) & truncrank(Dcut) + trunc_env = truncerror(; atol = 1.0e-10) & truncrank(chi) + + ctm_alg = SequentialCTMRG(; + tol = 1.0e-8, maxiter = 50, verbosity = 2, + trunc = trunc_env, projector_alg = :FullInfiniteProjector + ) + + interval = 5 + ntu_alg = NeighbourUpdate(; + opt_alg = FullEnvTruncation(; trunc = trunc_peps, tol = 1.0e-10), + bondenv_alg = NNEnv(), imaginary_time = false + ) + + # do one step of NTU to match benchmark data + peps0, = time_evolve(TimeEvolver(peps0, ham, 0.01, 6, ntu_alg)) + @info "Space of `peps0[1, 1]` = $(space(peps0[1, 1]))." + env0 = adapt(AT, CTMRGEnv(ones, ComplexF64, peps0, ℂ^1)) + env0, = leading_boundary(env0, peps0, ctm_alg) + # measure magnetization + magx = expectation_value(peps0, op, env0) + @info Printf.format(formatter, 0.06, real(magx), imag(magx)) + @test isapprox(magx, data[1, 2]; atol = 0.005) + + @testset "Neigborhood tensor update" begin + peps, env = deepcopy(peps0), deepcopy(env0) + count = 2 + evolver = TimeEvolver(peps, ham, 0.01, 30, ntu_alg; t0 = 0.06) + spaces0 = collect(space(t) for t in peps.A) + for (peps, info) in evolver + !(evolver.state.iter % interval == 0) && continue + spaces = collect(space(t) for t in peps.A) + if spaces0 == spaces + env, = leading_boundary(env, peps, ctm_alg) + else + env = complex(CTMRGEnv(info.wts)) + env, = leading_boundary(env, peps, ctm_alg) + end + # monitor the growth of env dimension + corner = env.corners[1, 1, 1] + corner_dim = dim.(space(corner, ax) for ax in 1:numind(corner)) + @info "Dimension of env.corner[1, 1, 1] = $(corner_dim)." + # measure magnetization + magx = expectation_value(peps, op, env) + @info Printf.format(formatter, info.t, real(magx), imag(magx)) + @test isapprox(magx, data[count, 2]; atol = 0.005) + count += 1 + spaces0 = spaces + end + end + end +end diff --git a/test/testsuite/timeevol/timestep.jl b/test/testsuite/timeevol/timestep.jl index f765e8562..fcd85aebc 100644 --- a/test/testsuite/timeevol/timestep.jl +++ b/test/testsuite/timeevol/timestep.jl @@ -39,3 +39,37 @@ function timeevol_timestep(AT) @test saved == snapshot # later scratch reuse must not overwrite earlier outputs end end + +"Check that all neighbourhood-update iteration interfaces agree exactly." +function timeevol_timestep_ntu(AT) + Nr, Nc = 2, 2 + H = adapt(AT, real(heisenberg_XYZ(ComplexF64, Trivial, InfiniteSquare(Nr, Nc); Jx = 1, Jy = 1, Jz = 1))) + Pspace, Vspace = ℂ^2, ℂ^4 + ψ0 = adapt(AT, InfinitePEPS(rand, Float64, Pspace, Vspace; unitcell = (Nr, Nc))) + dt, nstep = 0.1, 20 + trunc = truncerror(; atol = 1.0e-10) & truncrank(4) + + @testset "NeighbourUpdate timestep ($AT)" begin + alg = NeighbourUpdate(; opt_alg = ALSTruncation(; trunc)) + # manual timestep + evolver = TimeEvolver(ψ0, H, dt, nstep, alg) + ψ1, info1 = deepcopy(ψ0), nothing + for iter in 0:(nstep - 1) + ψ1, info1 = timestep(evolver, ψ1) + end + # time_evolve + evolver = TimeEvolver(ψ0, H, dt, nstep, alg) + ψ2, info2 = time_evolve(evolver) + # for-loop syntax + ## manually reset internal state of evolver + evolver.state = PEPSKit.NTUState(0, 0.0, ψ0) + ψ3, info3 = nothing, nothing, nothing + for state in evolver + ψ3, info3 = state + end + # results should be *exactly* the same + @test ψ1 == ψ2 == ψ3 + @test info1 == info2 == info3 + end + return nothing +end diff --git a/test/timeevol/heis_groundstate.jl b/test/timeevol/heis_groundstate.jl new file mode 100644 index 000000000..8a0919d26 --- /dev/null +++ b/test/timeevol/heis_groundstate.jl @@ -0,0 +1,18 @@ +using PEPSKit, CUDA, AMDGPU + +@isdefined(TestSuite) || include("../testsuite/TestSuite.jl") +using .TestSuite + +is_buildkite = get(ENV, "BUILDKITE", "false") == "true" + +if !is_buildkite + TestSuite.timeevol_heis_groundstate(Vector) +end + +if CUDA.functional() + TestSuite.timeevol_heis_groundstate(CuArray) +end + +if AMDGPU.functional() + TestSuite.timeevol_heis_groundstate(ROCArray) +end diff --git a/test/timeevol/sitedep_truncation.jl b/test/timeevol/sitedep_truncation.jl index f2b618155..b6ae056d9 100644 --- a/test/timeevol/sitedep_truncation.jl +++ b/test/timeevol/sitedep_truncation.jl @@ -8,6 +8,7 @@ is_buildkite = get(ENV, "BUILDKITE", "false") == "true" if !is_buildkite TestSuite.timeevol_sitedep_rotation(Vector) TestSuite.timeevol_sitedep_su(Vector) + TestSuite.timeevol_sitedep_ntu(Vector) end if CUDA.functional() diff --git a/test/timeevol/tf_ising_realtime.jl b/test/timeevol/tf_ising_realtime.jl new file mode 100644 index 000000000..42e3a7c4a --- /dev/null +++ b/test/timeevol/tf_ising_realtime.jl @@ -0,0 +1,18 @@ +using PEPSKit, CUDA, AMDGPU + +@isdefined(TestSuite) || include("../testsuite/TestSuite.jl") +using .TestSuite + +is_buildkite = get(ENV, "BUILDKITE", "false") == "true" + +if !is_buildkite + TestSuite.timeevol_ising_realtime(Vector) +end + +if CUDA.functional() + TestSuite.timeevol_ising_realtime(CuArray) +end + +if AMDGPU.functional() + TestSuite.timeevol_ising_realtime(ROCArray) +end diff --git a/test/timeevol/timestep.jl b/test/timeevol/timestep.jl index a3692114c..b6939a93d 100644 --- a/test/timeevol/timestep.jl +++ b/test/timeevol/timestep.jl @@ -7,12 +7,15 @@ is_buildkite = get(ENV, "BUILDKITE", "false") == "true" if !is_buildkite TestSuite.timeevol_timestep(Vector) + TestSuite.timeevol_timestep_ntu(Vector) end if CUDA.functional() TestSuite.timeevol_timestep(CuArray) + TestSuite.timeevol_timestep_ntu(CuArray) end if AMDGPU.functional() TestSuite.timeevol_timestep(ROCArray) + TestSuite.timeevol_timestep_ntu(ROCArray) end