diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index f1c537ee6..cf2e89377 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 @@ -89,6 +89,7 @@ include("algorithms/contractions/bp_contractions.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_1layer.jl") @@ -103,14 +104,17 @@ include("algorithms/ctmrg/c4v.jl") include("algorithms/truncation/truncationschemes.jl") include("algorithms/truncation/fullenv_truncation.jl") +include("algorithms/truncation/bond_tensor.jl") include("algorithms/truncation/bond_truncation.jl") -include("algorithms/time_evolution/trotter_gate.jl") include("algorithms/time_evolution/apply_gate.jl") include("algorithms/time_evolution/apply_mpo.jl") +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") @@ -143,7 +147,8 @@ export fixedpoint export absorb_weight export ALSTruncation, FullEnvTruncation -export SimpleUpdate +export NNEnv, NNpEnv, NNNEnv +export SimpleUpdate, NeighbourUpdate export TimeEvolver, timestep, time_evolve export InfiniteSquareNetwork @@ -151,6 +156,7 @@ export InfinitePartitionFunction export InfinitePEPS, InfiniteTransferPEPS export SUWeight export InfinitePEPO, InfiniteTransferPEPO +export infinite_temperature_density_matrix export BPEnv, BeliefPropagation export BPGauge, SUGauge diff --git a/src/algorithms/bp/beliefpropagation.jl b/src/algorithms/bp/beliefpropagation.jl index 0246d7e20..0aa900912 100644 --- a/src/algorithms/bp/beliefpropagation.jl +++ b/src/algorithms/bp/beliefpropagation.jl @@ -59,7 +59,7 @@ function leading_boundary(env₀::BPEnv, network::InfiniteSquareNetwork, alg::Be end function leading_boundary(env₀::BPEnv, state::InfiniteState, alg::BeliefPropagation) if alg.bipartite - @assert _state_bipartite_check(state) + _is_bipartite(state) || error("Input state is not bipartite") end return leading_boundary(env₀, InfiniteSquareNetwork(state), alg) end diff --git a/src/algorithms/bp/gaugefix.jl b/src/algorithms/bp/gaugefix.jl index 561dd2183..5c98dc265 100644 --- a/src/algorithms/bp/gaugefix.jl +++ b/src/algorithms/bp/gaugefix.jl @@ -7,16 +7,6 @@ Algorithm for gauging PEPS with belief propagation fixed point messages. # TODO: add options end -function _bpenv_bipartite_check(env::BPEnv) - for (r, c) in Iterators.product(1:2, 1:2) - r′, c′ = _next(r, 2), _next(c, 2) - if !all(env[:, r, c] .== env[:, r′, c′]) - return false - end - end - return true -end - """ gauge_fix(psi::Union{InfinitePEPS, InfinitePEPO}, alg::BPGauge, env::BPEnv) @@ -25,7 +15,7 @@ an [`InfinitePEPO`](@ref) interpreted as purified state with two physical legs) using fixed point environment `env` of belief propagation. """ function gauge_fix(psi::InfinitePEPS, alg::BPGauge, env::BPEnv) - bipartite = _state_bipartite_check(psi) && _bpenv_bipartite_check(env) + bipartite = _is_bipartite(psi) && _is_bipartite(env) psi′ = copy(psi) XXinv = map(eachcoordinate(psi, 1:2)) do I _, X, Xinv = _bp_gauge_fix!(CartesianIndex(I), psi′, env) diff --git a/src/algorithms/contractions/bondenv/als_solve.jl b/src/algorithms/contractions/bondenv/als_solve.jl index af083f618..5d049f84c 100644 --- a/src/algorithms/contractions/bondenv/als_solve.jl +++ b/src/algorithms/contractions/bondenv/als_solve.jl @@ -2,139 +2,148 @@ In the following, the names `Ra`, `Sa` etc comes from the fast full update article Physical Review B 92, 035142 (2015) =# - """ -$(SIGNATURES) +Contract the virtual legs between +``` + -- DX --a-- D --b-- DY -- + ↓ ↓ + da db +``` +""" +function _combine_ket(a::MPSTensor, b::MPSTensor) + return @tensor ket[DX DY; da db] := a[DX da; D] * b[D db; DY] +end -Construct the tensor -``` - ┌-----------------------------------┐ - | ┌----┐ | - └---| |- DX0 Db0 - b -- DY0 -┘ - | | ↓ - |benv| db - | | ↓ - ┌---| |- DX1 Db1 - b† - DY1 -┐ - | └----┘ | - └-----------------------------------┘ -``` -""" -function _tensor_Ra(benv::BondEnv, b::MPSTensor) - return @autoopt @tensor Ra[DX1 Db1; DX0 Db0] := ( - benv[DX1 DY1; DX0 DY0] * b[Db0 db; DY0] * conj(b[Db1 db; DY1]) - ) +function _combine_ket_for_svd(a::MPSTensor, b::MPSTensor) + return @tensor ket[DX da; db DY] := a[DX da; D] * b[D db; DY] end """ -$(SIGNATURES) - -Construct the tensor +Construct the norm with bra bond tensors removed ``` - ┌-----------------------------------┐ - | ┌----┐ | - └---| |- DX0 -- (a2 b2) -- DY0 --┘ - | | ↓ ↓ - |benv| da db - | | ↓ - ┌---| |- DX1 Db1 -- b† - DY1 --┐ - | └----┘ | - └-----------------------------------┘ + ┌benv-------┐ + ├---a---b---┤ + | ↓ ↓ | + ├-- --┤ + └-----------┘ ``` """ -function _tensor_Sa( - benv::BondEnv, b::MPSTensor, a2b2::AbstractTensorMap{T, S, 2, 2} - ) where {T <: Number, S <: ElementarySpace} - return @autoopt @tensor Sa[DX1 da; Db1] := ( - benv[DX1 DY1; DX0 DY0] * conj(b[Db1 db; DY1]) * a2b2[DX0 DY0; da db] - ) +function _benv_ket(benv::BondEnv, ket::AbstractTensorMap{T, S, 2, 2}) where {T, S} + return benv * twistdual(ket, 1:2) end """ -$(SIGNATURES) + _als_tensor_R(benv::BondEnv, xs::Vector{<:MPSTensor}, i::Int) -Construct the tensor -``` - ┌-----------------------------------┐ - | ┌----┐ | - └---| |- DX0 - a -- Da0 DY0 -┘ - | | ↓ - |benv| da - | | ↓ - ┌---| |- DX1 - a† - Da1 DY1 -┐ - | └----┘ | - └-----------------------------------┘ -``` -""" -function _tensor_Rb(benv::BondEnv, a::MPSTensor) - return @autoopt @tensor Rb[Da1 DY1; Da0 DY0] := ( - benv[DX1 DY1; DX0 DY0] * a[DX0 da; Da0] * conj(a[DX1 da; Da1]) - ) +Construct the bond environment around the `i`th bond tensor +in two-site ALS optimization. +``` + i = 1 i = 2 + ┌benv-------┐ ┌benv-------┐ + ├-- --b---┤ ├---a-- --┤ + | ↓ | | ↓ | + ├-- --b̄---┤ ├---ā-- --┤ + └-----------┘ └-----------┘ +``` +""" +_als_tensor_R(benv, xs, i::Int) = _als_tensor_R(benv, xs, Val(i)) +function _als_tensor_R(benv::BondEnv, xs::Vector{<:MPSTensor}, ::Val{1}) + return @tensor Ra[DX1 D1; DX0 D0] := + benv[DX1 DY1; DX0 DY0] * xs[2][D0 db; DY0] * conj(xs[2][D1 db; DY1]) +end +function _als_tensor_R(benv::BondEnv, xs::Vector{<:MPSTensor}, ::Val{2}) + return @tensor Rb[D1 DY1; D0 DY0] := + benv[DX1 DY1; DX0 DY0] * xs[1][DX0 da; D0] * conj(xs[1][DX1 da; D1]) end """ -$(SIGNATURES) - -Construct the tensor +Calculate the 2-site norm ``` - ┌-----------------------------------┐ - | ┌----┐ | - └---| |- DX0 -- (a2 b2) -- DY0 --┘ - | | ↓ ↓ - |benv| da db - | | ↓ - ┌---| |- DX1 -- a† - Da1 DY1 --┐ - | └----┘ | - └-----------------------------------┘ + ┌benv-------┐ + ├---a---b---┤ + | ↓ ↓ | + ├---ā---b̄---┤ + └-----------┘ ``` +using pre-calcuated partial contraction results. """ -function _tensor_Sb( - benv::BondEnv, a::MPSTensor, a2b2::AbstractTensorMap{T, S, 2, 2} - ) where {T <: Number, S <: ElementarySpace} - return @autoopt @tensor Sb[Da1 db; DY1] := ( - benv[DX1 DY1; DX0 DY0] * conj(a[DX1 da; Da1]) * a2b2[DX0 DY0; da db] - ) +function _als_norm( + ket::AbstractTensorMap{T, S, 2, 2}, benv_ket::AbstractTensorMap{T, S, 2, 2} + ) where {T, S} + return @tensor benv_ket[DX1 DY1; da db] * conj(ket[DX1 DY1; da db]) +end +function _als_norm(a::MPSTensor, Ra::BondEnv) + return @tensor Ra[DX1 D1; DX0 D0] * a[DX0 da; D0] * conj(a[DX1 da; D1]) end """ -$(SIGNATURES) + _als_tensor_S( + benv_ket::AbstractTensorMap{T, S, 2, 2}, + xs::Vector{<:MPSTensor}, i::Int + ) where {T <: Number, S <: ElementarySpace} -Calculate the inner product +Construct the overlap but with one of the bra bond tensor removed. ``` - ┌--------------------------------┐ - | ┌----┐ | - └---| |- DX0 - (a2 b2) - DY0 -┘ - | | ↓ ↓ - |benv| da db - | | ↓ ↓ - ┌---| |- DX1 - (a1 b1)†- DY1 -┐ - | └----┘ | - └--------------------------------┘ + i = 1 i = 2 + ┌benv-------┐ ┌benv-------┐ + ├---a₂==b₂--┤ ├---a₂==b₂--┤ + | ↓ ↓ | | ↓ ↓ | + ├-- --b̄---┤ ├---ā-- --┤ + └-----------┘ └-----------┘ ``` +The ket part is provided by the partial contraction `benv_ket`. """ -function inner_prod( - benv::BondEnv, a1b1::AbstractTensorMap{T, S, 2, 2}, a2b2::AbstractTensorMap{T, S, 2, 2} +_als_tensor_S(benv_ket, xs, i::Int) = _als_tensor_S(benv_ket, xs, Val(i)) +function _als_tensor_S( + benv_ket::AbstractTensorMap{T, S, 2, 2}, + xs::Vector{<:MPSTensor}, ::Val{1} ) where {T <: Number, S <: ElementarySpace} - return @autoopt @tensor benv[DX1 DY1; DX0 DY0] * - conj(a1b1[DX1 DY1; da db]) * a2b2[DX0 DY0; da db] + return @tensor Sa[DX1 da; D1] := + benv_ket[DX1 DY1; da db] * conj(xs[2][D1 db; DY1]) +end +function _als_tensor_S( + benv_ket::AbstractTensorMap{T, S, 2, 2}, + xs::Vector{<:MPSTensor}, ::Val{2} + ) where {T <: Number, S <: ElementarySpace} + return @tensor contractcheck = true Sb[D1 db; DY1] := + benv_ket[DX1 DY1; da db] * conj(xs[1][DX1 da; D1]) end """ -$(SIGNATURES) +Calculate the inner product (overlap) +``` + ┌benv-------┐ + ├---a₂--b₂--┤ + | ↓ ↓ | + ├---ā---b̄---┤ + └-----------┘ +``` +using pre-calculated partial contraction results. +""" +function _als_overlap(a::MPSTensor, Sa::MPSTensor) + # applies to b, Sb as well + # @tensor Sb[D1 db; DY1] * conj(b[D1 db; DY1]) + return @tensor Sa[DX1 da; D1] * conj(a[DX1 da; D1]) +end -Contract the axis between `a` and `b` tensors +""" +Calculate the 2-site ALS inner product ⟨a₁,b₁|a₂,b₂⟩ ``` - -- DX - a - D - b - DY -- - ↓ ↓ - da db + ┌benv-------┐ + ├---a₂--b₂--┤ + | ↓ ↓ | + ├---ā₁--b̄₁--┤ + └-----------┘ ``` +where `|bra⟩ = |a₁,b₁⟩` and `|ket⟩ = |a₂,b₂⟩`, +with virtual leg between a, b contracted. """ -function _combine_ab( - a::MPSTensor, b::AbstractTensorMap{T, S, 1, 2} +function inner_prod( + benv::BondEnv, bra::AbstractTensorMap{T, S, 2, 2}, + ket::AbstractTensorMap{T, S, 2, 2} ) where {T <: Number, S <: ElementarySpace} - return @tensor ab[DX DY; da db] := a[DX da; D] * b[D; db DY] -end -function _combine_ab(a::MPSTensor, b::MPSTensor) - return @tensor ab[DX DY; da db] := a[DX da; D] * b[D db; DY] + return @autoopt @tensor benv[DX1 DY1; DX0 DY0] * + conj(bra[DX1 DY1; da db]) * ket[DX0 DY0; da db] end """ @@ -161,10 +170,30 @@ function cost_function_als(benv, ψ1, ψ2) return cost, fid end +# applies to Rb, Sb, b as well +# b22 is the pre-calculated untruncated norm +function cost_function_als(Ra::BondEnv, Sa::MPSTensor, a::MPSTensor, b22::Real) + b11 = real(_als_norm(a, Ra)) + b12 = _als_overlap(a, Sa) + cost = b11 + b22 - 2 * real(b12) + fid = abs2(b12) / abs(b11 * b22) + return cost, fid +end + """ $(SIGNATURES) Solve the equations `Rx x = Sx` with initial guess `x0`. + +In ALS over `a`, `b`, if we fix `b`, the cost function can +be expressed in the `Ra`, `Sa` tensors as +``` + f(a†,a) = a† Ra a - a† Sa - Sa† a + const +``` +Therefore `f` is minimized when +``` + ∂f/∂ā = Ra a - Sa = 0 +``` """ function _solve_als( Rx::AbstractTensorMap{T, S, N, N}, diff --git a/src/algorithms/contractions/bondenv/benv_ctm.jl b/src/algorithms/contractions/bondenv/benv_ctm.jl index 45262cc67..8c6b1a199 100644 --- a/src/algorithms/contractions/bondenv/benv_ctm.jl +++ b/src/algorithms/contractions/bondenv/benv_ctm.jl @@ -9,8 +9,8 @@ Construct the environment (norm) tensor -1 0 1 2 ``` with `X` at position `[row, col]`. -`X, Y` are unitary tensors produced when finding the reduced site tensors -with `_qr_bond`; and `XX = X' X` and `YY = Y' Y` (stacked together). +`X, Y` are unitary tensors produced when finding the reduced site tensors, +and `XX = X' X` and `YY = Y' Y` (stacked together). Axis order: `[DX1 DY1; DX0 DY0]`, as in ``` @@ -24,8 +24,8 @@ Axis order: `[DX1 DY1; DX0 DY0]`, as in ``` """ function bondenv_ctm( - row::Int, col::Int, X::T, Y::T, env::CTMRGEnv - ) where {T} + row::Int, col::Int, X::TX, Y::TY, env::CTMRGEnv + ) where {TX, TY} Nr, Nc = size(env.corners)[[2, 3]] cm1 = _prev(col, Nc) cp1 = _next(col, Nc) diff --git a/src/algorithms/contractions/bondenv/benv_ntu.jl b/src/algorithms/contractions/bondenv/benv_ntu.jl new file mode 100644 index 000000000..ad4327f7b --- /dev/null +++ b/src/algorithms/contractions/bondenv/benv_ntu.jl @@ -0,0 +1,239 @@ +#= +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 + +""" +SVD with `truncrank(1)`. +""" +function _svd_cut!(t::AbstractTensorMap) + t1, s, t2 = svd_trunc!(t; trunc = truncrank(1)) + return absorb_s(t1, s, t2) +end + +""" +Algorithm struct for "NTU-NN" bond environment. +""" +struct NNEnv <: NeighbourEnv end +""" +Calculate the bond environment within "NTU-NN" approximation. +``` + -1 ●=======● + ║ ║ + 0 ●===X== ==Y===● + ║ ║ + 1 ●=======● + -1 0 1 2 +``` +""" +function bondenv_ntu( + row::Int, col::Int, X::TX, Y::TY, state::S, alg::NNEnv + ) where {TX, TY, S <: InfiniteState} + neighbors = [(-1, 0), (0, -1), (1, 0), (1, 1), (0, 2), (-1, 1)] + m = collect_neighbors(state, row, col, neighbors) + X, Y = _prepare_site_tensor(X), _prepare_site_tensor(Y) + # southwest half + @autoopt @tensor benv_sw[Dse1 Dse0 Dw1 Dw0 Dnw1 Dnw0] := + cor_se(m[1, 1])[Dse1 Dse0 Ds1 Ds0] * + cor_sw(m[1, 0])[Dsw1 Dsw0 Ds1 Ds0] * + edge_w(X, hair_w(m[0, -1]))[Dnw1 Dnw0 Dw1 Dw0 Dsw1 Dsw0] + normalize!(benv_sw, Inf) + # northeast half + @autoopt @tensor benv_ne[Dnw1 Dnw0 De1 De0 Dse1 Dse0] := + cor_nw(m[-1, 0])[Dn1 Dn0 Dnw1 Dnw0] * + cor_ne(m[-1, 1])[Dne1 Dne0 Dn1 Dn0] * + edge_e(Y, hair_e(m[0, 2]))[Dne1 Dne0 Dse1 Dse0 De1 De0] + normalize!(benv_ne, Inf) + @tensor benv[Dw1 De1; Dw0 De0] := + benv_sw[Dse1 Dse0 Dw1 Dw0 Dnw1 Dnw0] * benv_ne[Dnw1 Dnw0 De1 De0 Dse1 Dse0] + normalize!(benv, Inf) + return benv +end + +""" +Algorithm struct for "NTU-NN+" bond environment. +""" +struct NNpEnv <: NeighbourEnv end +""" +Calculate the bond environment within "NTU-NN+" approximation. +``` + -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, X::TX, Y::TY, state::S, alg::NNpEnv + ) where {TX, TY, S <: InfiniteState} + 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) + 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) + + m = ms[0, -1] + @tensoropt hW[anw DXw1 DXw0 asw] := + hair_w(ms[0, -2])[Dw21 Dw20] * + nw1[Dnw11 Dnw10 anw] * sw1[Dsw11 Dsw10 asw] * + twistdual(m, 1)[p Dnw10 DXw0 Dsw10 Dw20] * + conj(m[p Dnw11 DXw1 Dsw11 Dw21]) + + m = ms[0, 2] + @tensoropt hE[ane DYe1 DYe0 ase] := + hair_e(ms[0, 3])[De21 De20] * + ne2[ane Dne21 Dne20] * se2[ase Dse21 Dse20] * + twistdual(m, 1)[p Dne20 De20 Dse20 DYe0] * + conj(m[p Dne21 De21 Dse21 DYe1]) + + # ---- corners (size D^4) with a 1D auxiliary leg ---- + + m = ms[-1, 0] + @tensoropt NW[at Dn1 Dn0 DXn1 DXn0 anw] := + tl[Dtl1 Dtl0 at] * nw2[anw Dnw21 Dnw20] * + twistdual(m, 1)[p Dtl0 Dn0 DXn0 Dnw20] * + conj(m[p Dtl1 Dn1 DXn1 Dnw21]) + + m = ms[-1, 1] + @tensoropt NE[at Dn1 Dn0 DYn1 DYn0 ane] := + tr[at Dtr1 Dtr0] * ne1[Dne11 Dne10 ane] * + twistdual(m, 1)[p Dtr0 Dne10 DYn0 Dn0] * + conj(m[p Dtr1 Dne11 DYn1 Dn1]) + + m = ms[1, 0] + @tensoropt SW[asw DXs1 DXs0 Ds1 Ds0 ab] := + bl[Dbl1 Dbl0 ab] * sw2[asw Dsw21 Dsw20] * + twistdual(m, 1)[p DXs0 Ds0 Dbl0 Dsw20] * + conj(m[p DXs1 Ds1 Dbl1 Dsw21]) + + m = ms[1, 1] + @tensoropt SE[ase DYs1 DYs0 Ds1 Ds0 ab] := + br[ab Dbr1 Dbr0] * se1[Dse11 Dse10 ase] * + twistdual(m, 1)[p DYs0 Dse10 Dbr0 Ds0] * + conj(m[p DYs1 Dse11 Dbr1 Ds1]) + + # ---- left half ---- + @tensoropt benvL[at Dn1 Dn0 Dw1 Dw0 Ds1 Ds0 ab] := + hW[anw DXw1 DXw0 asw] * + NW[at Dn1 Dn0 DXn1 DXn0 anw] * + SW[asw DXs1 DXs0 Ds1 Ds0 ab] * + twistdual(X, 1)[p DXn0 Dw0 DXs0 DXw0] * + conj(X[p DXn1 Dw1 DXs1 DXw1]) + normalize!(benvL, Inf) + + # ---- right half ---- + + @tensoropt benvR[at Dn1 Dn0 De1 De0 Ds1 Ds0 ab] := + hE[ane DYe1 DYe0 ase] * + NE[at Dn1 Dn0 DYn1 DYn0 ane] * + SE[ase DYs1 DYs0 Ds1 Ds0 ab] * + twistdual(Y, 1)[p DYn0 DYe0 DYs0 De0] * + conj(Y[p DYn1 DYe1 DYs1 De1]) + normalize!(benvR, Inf) + + # ---- the full NN+ environment ---- + + @tensor benv[Dw1, De1; Dw0, De0] := + benvL[at Dn1 Dn0 Dw1 Dw0 Ds1 Ds0 ab] * + benvR[at Dn1 Dn0 De1 De0 Ds1 Ds0 ab] + normalize!(benv, Inf) + return benv +end + +""" +Algorithm struct for "NTU-NNN" bond environment. +""" +struct NNNEnv <: NeighbourEnv end +""" +Calculates the bond environment within "NTU-NNN" approximation. +``` + -1 ●===●=======●===● + ║ ║ ║ ║ + 0 ●===X== ==Y===● + ║ ║ ║ ║ + 1 ●===●=======●===● + -1 0 1 2 +``` +""" +function bondenv_ntu( + row::Int, col::Int, X::TX, Y::TY, state::S, alg::NNNEnv + ) where {TX, TY, S <: InfiniteState} + neighbors = [ + (-1, -1), (0, -1), (1, -1), + (1, 0), (1, 1), (1, 2), (0, 2), + (-1, 2), (-1, 1), (-1, 0), + ] + m = collect_neighbors(state, row, col, neighbors) + X, Y = _prepare_site_tensor(X), _prepare_site_tensor(Y) + #= left half + -1 ●======●=== -1/-2 + ║ ║ + 0 ●======X=== -3/-4 + ║ ║ + ....D1.....D2... + ║ ║ + 1 ●==D3==●=== -5/-6 + -1 0 + =# + vecl = enlarge_corner_nw(cor_nw(m[-1, -1]), edge_n(m[-1, 0]), edge_w(m[0, -1]), X) + @tensor vecl[:] := + cor_sw(m[1, -1])[D11 D10 D31 D30] * + edge_s(m[1, 0])[D21 D20 -5 -6 D31 D30] * + vecl[D11 D10 D21 D20 -1 -2 -3 -4] + normalize!(vecl, Inf) + #= right half + -1 -1/-2===●==D1==● + ║ ║ + ........D2.....D3... + ║ ║ + 0 -3/-4===Y======● + ║ ║ + 1 -5/-6===●======● + 0 1 + =# + vecr = enlarge_corner_se(cor_se(m[1, 2]), edge_s(m[1, 1]), edge_e(m[0, 2]), Y) + @tensor vecr[:] := + edge_n(m[-1, 1])[D11 D10 D21 D20 -1 -2] * + cor_ne(m[-1, 2])[D31 D30 D11 D10] * + vecr[D21 D20 D31 D30 -3 -4 -5 -6] + normalize!(vecr, Inf) + # combine left and right part + @tensor benv[-1 -2; -3 -4] := vecl[1 2 -1 -3 3 4] * vecr[1 2 -2 -4 3 4] + normalize!(benv, Inf) + return benv +end diff --git a/src/algorithms/contractions/bondenv/benv_tools.jl b/src/algorithms/contractions/bondenv/benv_tools.jl index 38400b919..8b2ff362e 100644 --- a/src/algorithms/contractions/bondenv/benv_tools.jl +++ b/src/algorithms/contractions/bondenv/benv_tools.jl @@ -1,8 +1,15 @@ +#= +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} # Orthogonal tensors obtained PEPSTensor/PEPOTensor -# with one physical leg factored out by `_qr_bond` +# with one physical leg factored out by `bond_tensor_...` const PEPSOrth{T, S} = AbstractTensor{T, S, 4} where {T <: Number, S <: ElementarySpace} const PEPOOrth{T, S} = AbstractTensor{T, S, 5} where {T <: Number, S <: ElementarySpace} @@ -11,3 +18,198 @@ _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 +at positions `neighbors` relative to `(row, col)` +""" +function collect_neighbors( + state::InfiniteState, row::Int, col::Int, neighbors::Vector{Tuple{Int, Int}} + ) + Nr, Nc = size(state) + return Dict( + nb => _prepare_site_tensor(state[mod1(row + nb[1], Nr), mod1(col + nb[2], Nc)]) + for nb in neighbors + ) +end + +""" + benv_tensor(ket::PEPSTensor, bra::PEPSTensor, open_axs::Vector{Int}) + benv_tensor(ket::PEPSTensor, bra::PEPSTensor, open_axs::Vector{Int}, hairs::Vector{H}) where {T, 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. + +# 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 = [1, 2, 3]`) + ``` + 2 + ╱ + ┌-----ket----- 4 + | ╱ | + | 6 | 1 + | | ╱ + └-----bra----- 3 + ╱ + 5 + ``` +""" +function benv_tensor( + ket::PEPSTensor, bra::PEPSTensor, open_axs::Vector{Int} + ) + # no hair tensors to be attached to virtual legs + return benv_tensor(ket, bra, open_axs, fill(nothing, 4 - length(open_axs))) +end +function benv_tensor( + ket::PEPSTensor, bra::PEPSTensor, open_axs::Vector{Int}, hairs::Vector{H} + ) where {H <: Union{Nothing, Hair}} + @assert length(hairs) == 4 - length(open_axs) + ket2, nax = copy(ket), numind(ket) + axs, open_axs2 = (2:5), open_axs .+ 1 + # contract with hair tensors + hair_axs = Tuple(ax for ax in axs if ax ∉ open_axs2) + for (h, ax) in zip(hairs, hair_axs) + if h === nothing + twistdual!(ket2, ax) + continue + end + # ensure the hair doesn't change the virtual spaces + @assert space(h, 1) == space(h, 2)' + ket_indices = collect(-1:-1:-nax) + ket_indices[ax] = 1 + ket2 = ncon([h, ket2], [[-ax, 1], ket_indices]) + end + # combine bra and ket + indexlist = [-collect(1:2:(2 * nax)), -collect(2:2:(2 * nax))] + for ax in 1:nax + if ax ∉ open_axs2 + indexlist[1][ax] = indexlist[2][ax] = ax + end + end + twistdual!(ket2, 1) + return ncon([bra, ket2], indexlist, [true, false]) +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 +=# +const open_axs_hair = Dict(:n => [SOUTH], :e => [WEST], :s => [NORTH], :w => [EAST]) +const open_axs_cor = Dict( + :nw => [EAST, SOUTH], :ne => [SOUTH, WEST], :se => [NORTH, WEST], :sw => [NORTH, EAST] +) +const open_axs_edge = Dict( + :n => [EAST, SOUTH, WEST], + :e => [NORTH, SOUTH, WEST], + :s => [NORTH, EAST, WEST], + :w => [NORTH, EAST, SOUTH], +) + +# construction of hairs +for (dir, open_axs) in open_axs_hair + 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 open_axs_cor + 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 open_axs_edge + 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 + +""" +Enlarge the northwest corner +``` + ctl══ D1 ══ et ══ -5/-6 + ║ ║ + D2 D3 + ║ ║ + el ══ D4 ══ X ═══ -7/-8 + ║ ║ + -1/-2 -3/-4 +``` +""" +function enlarge_corner_nw( + ctl::AbstractTensor{E, S, 4}, + et::AbstractTensor{E, S, 6}, el::AbstractTensor{E, S, 6}, + ket::PEPSTensor, bra::PEPSTensor = ket + ) where {E, S} + return @tensoropt ctl2[:] := ctl[D11 D10 D21 D20] * + et[-5 -6 D31 D30 D11 D10] * el[D21 D20 D41 D40 -1 -2] * + conj(bra[d D31 -7 -3 D41]) * twistdual(ket, 1)[d D30 -8 -4 D40] +end + +""" +Enlarge the southeast corner +``` + -1/-2 -3/-4 + ║ ║ + -5/-6 ═════ Y ══ D1 ═══ er + ║ ║ + D2 D3 + ║ ║ + -7/-8 ═════ eb ═ D4 ══ cbr +``` +""" +function enlarge_corner_se( + cbr::AbstractTensor{E, S, 4}, + eb::AbstractTensor{E, S, 6}, er::AbstractTensor{E, S, 6}, + ket::PEPSTensor, bra::PEPSTensor = ket + ) where {E, S} + return @tensoropt cbr2[:] := cbr[D31 D30 D41 D40] * + eb[D21 D20 D41 D40 -7 -8] * er[-3 -4 D31 D30 D11 D10] * + conj(bra[d -1 D11 D21 -5]) * twistdual(ket, 1)[d -2 D10 D20 -6] +end diff --git a/src/algorithms/contractions/bondenv/gaugefix.jl b/src/algorithms/contractions/bondenv/gaugefix.jl index 0572d4852..8d964d395 100644 --- a/src/algorithms/contractions/bondenv/gaugefix.jl +++ b/src/algorithms/contractions/bondenv/gaugefix.jl @@ -44,9 +44,7 @@ Reference: - Physical Review B 92, 035142 (2015) """ function fixgauge_benv( - Z::AbstractTensorMap{T, S, 1, 2}, - a::AbstractTensorMap{T, S, 1, 2}, - b::AbstractTensorMap{T, S, 2, 1}, + Z::AbstractTensorMap{T, S, 1, 2}, a::MPSTensor, b::MPSTensor ) where {T <: Number, S <: ElementarySpace} @assert !isdual(space(Z, 1)) @assert !isdual(space(a, 2)) @@ -83,7 +81,7 @@ function fixgauge_benv( ↓ -1 =# - @plansor a[-1; -2 -3] := R[-1; 1] * a[1; -2 -3] + @plansor a[-1 -2; -3] := R[-1; 1] * a[1 -2; -3] @plansor b[-1 -2; -3] := b[-1 -2; 1] * L[-3; 1] @plansor Z[-1; -2 -3] := Z[-1; 1 2] * Rinv[1; -2] * Linv[2; -3] (isdual(space(R, 1)) == isdual(space(R, 2))) && twist!(a, 1) @@ -92,29 +90,59 @@ function fixgauge_benv( end """ -When the (half) bond environment `Z` consists of two `PEPSOrth` tensors `X`, `Y` as +Apply the gauge transformation `Rinv` for `Z` ``` - ┌---------------┐ ┌-------------------┐ - | | = | | , - └---Z-- --┘ └--Z0---X-- --Y--┘ - ↓ ↓ + ┌-----------------------┐ + └---Z--(X)--Rinv-- ---┘ + ↓ ``` -apply the gauge transformation `Linv`, `Rinv` for `Z` to `X`, `Y`: +to `X`. For example, when `X` is a `PEPSTensor`, ``` - -1 -1 - | | - -4 - X - 1 - Rinv - -2 -4 - Linv - 1 - Y - -2 - | | - -3 -3 + -2 + | + -5 - X - 1 - Rinv - -3 + | ╲ + -4 -1 ``` """ -function _fixgauge_benvXY( - X::PEPSOrth{T, S}, - Y::PEPSOrth{T, S}, - Linv::AbstractTensorMap{T, S, 1, 1}, - Rinv::AbstractTensorMap{T, S, 1, 1}, - ) where {T <: Number, S <: ElementarySpace} - @plansor X[-1 -2 -3 -4] := X[-1 1 -3 -4] * Rinv[1; -2] - @plansor Y[-1 -2 -3 -4] := Y[-1 -2 -3 1] * Linv[1; -4] - return X, Y +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` +``` + ┌-----------------------┐ + └---Z--- ---Linv--(Y)--┘ + ↓ +``` +to `Y`. For example, when `Y` is a `PEPSTensor`, +``` + -2 + | + -5 - Linv - 1 - Y - -3 + | ╲ + -4 -1 +``` +""" +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 363b23f8d..67d4c63cf 100644 --- a/src/algorithms/time_evolution/apply_gate.jl +++ b/src/algorithms/time_evolution/apply_gate.jl @@ -1,3 +1,5 @@ +const NNGate{T, S} = AbstractTensorMap{T, S, 2, 2} + """ Apply 1-site `gate` on the PEPS or PEPO tensor `a`. """ @@ -21,98 +23,7 @@ end """ $(SIGNATURES) -Use QR decomposition on two tensors `A`, `B` connected by a bond to get the reduced tensors. -When `A`, `B` are PEPSTensors, -``` - 2 1 1 - | | | - 5 -A/B- 3 ====> 4 - X ← 2 1 ← a - 3 1 - b → 3 4 → Y - 2 - | ↘ | ↘ ↘ | - 4 1 3 2 2 3 -``` -When `A`, `B` are PEPOTensors, -- If `gate_ax = 1` -``` - 2 3 1 2 1 2 - ↘ | ↘ | ↘ | - 6 -A/B- 4 ====> 5 - X ← 3 1 ← a - 3 1 - b → 3 5 → Y - 3 - | ↘ | ↘ ↘ | - 5 1 4 2 2 4 -``` -- If `gate_ax = 2` -``` - 2 3 2 2 2 2 - ↘ | | ↘ ↘ | - 6 -A/B- 4 ====> 5 - X ← 3 1 ← a - 3 1 - b → 3 5 → Y - 3 - | ↘ | ↘ | ↘ - 5 1 4 1 4 1 -``` -""" -function _qr_bond(A::PT, B::PT; gate_ax::Int = 1, kwargs...) where {PT <: Union{PEPSTensor, PEPOTensor}} - @assert 1 <= gate_ax <= numout(A) - permA, permB, permX, permY = if A isa PEPSTensor - ((2, 4, 5), (1, 3)), ((2, 3, 4), (1, 5)), (1, 4, 2, 3), Tuple(1:4) - else - if gate_ax == 1 - ((2, 3, 5, 6), (1, 4)), ((2, 3, 4, 5), (1, 6)), (1, 2, 5, 3, 4), Tuple(1:5) - else - ((1, 3, 5, 6), (2, 4)), ((1, 3, 4, 5), (2, 6)), (1, 2, 5, 3, 4), Tuple(1:5) - end - end - X, a = left_orth!(permute(A, permA; copy = true); kwargs...) - Y, b = left_orth!(permute(B, permB; copy = true); kwargs...) - X, Y = permute(X, permX), permute(Y, permY) - b = permute(b, ((3, 2), (1,))) - return X, a, b, Y -end - -""" -$(SIGNATURES) - -Reconstruct the tensors connected by a bond from their `_qr_bond` results. -For PEPSTensors, -``` - -2 -2 - | | - -5- X - 1 - a - -3 -5 - b - 1 - Y - -3 - | ↘ ↘ | - -4 -1 -1 -4 -``` -For PEPOTensors -``` - -2 -3 -2 -3 - ↘ | ↘ | - -6- X - 1 - a - -4 -6 - b - 1 - Y - -4 - | ↘ ↘ | - -5 -1 -1 -5 - - -3 -2 -2 -3 - | ↘ ↘ | - -6- X - 1 - a - -4 -6 - b - 1 - Y - -4 - | ↘ | ↘ - -5 -1 -5 -1 -``` -""" -function _qr_bond_undo(X::PEPSOrth, a::AbstractTensorMap, b::AbstractTensorMap, Y::PEPSOrth) - @tensor A[-1; -2 -3 -4 -5] := X[-2 1 -4 -5] * a[1 -1 -3] - @tensor B[-1; -2 -3 -4 -5] := b[-5 -1 1] * Y[-2 -3 -4 1] - return A, B -end -function _qr_bond_undo(X::PEPOOrth, a::AbstractTensorMap, b::AbstractTensorMap, Y::PEPOOrth) - if !isdual(space(a, 2)) - @tensor A[-1 -2; -3 -4 -5 -6] := X[-2 -3 1 -5 -6] * a[1 -1 -4] - @tensor B[-1 -2; -3 -4 -5 -6] := b[-6 -1 1] * Y[-2 -3 -4 -5 1] - else - @tensor A[-1 -2; -3 -4 -5 -6] := X[-1 -3 1 -5 -6] * a[1 -2 -4] - @tensor B[-1 -2; -3 -4 -5 -6] := b[-6 -2 1] * Y[-1 -3 -4 -5 1] - end - return A, B -end - -""" -$(SIGNATURES) - -Apply 2-site `gate` on the reduced matrices `a`, `b` +Apply 2-site `gate` on the reduced bond tensors `a`, `b` ``` -1← a --- 3 --- b ← -4 -2 -3 ↓ ↓ ↓ ↓ @@ -123,10 +34,7 @@ Apply 2-site `gate` on the reduced matrices `a`, `b` -2 -3 -1← a --- 3 --- b ← -4 ``` """ -function _apply_gate( - a::AbstractTensorMap, b::AbstractTensorMap, - gate::AbstractTensorMap{T, S, 2, 2}, trunc::TruncationStrategy - ) where {T <: Number, S <: ElementarySpace} +function _apply_gate(a::MPSTensor, b::MPSTensor, gate::NNGate, trunc::TruncationStrategy) V = space(b, 1) need_flip = isdual(V) if isdual(space(a, 2)) @@ -134,15 +42,11 @@ function _apply_gate( else @tensor a2b2[-1 -2; -3 -4] := gate[-2 -3; 1 2] * a[-1 1 3] * b[3 2 -4] end - trunc = if trunc isa FixedSpaceTruncation - need_flip ? truncspace(flip(V)) : truncspace(V) - else - trunc - end a, s, b, ϵ = svd_trunc!(a2b2; trunc, alg = LAPACK_QRIteration()) a, b = absorb_s(a, s, b) if need_flip a, s, b = flip(a, numind(a)), _fliptwist_s(s), flip(b, 1) end + b = permute(b, ((1, 2), (3,))) return a, s, b, ϵ end diff --git a/src/algorithms/time_evolution/apply_mpo.jl b/src/algorithms/time_evolution/apply_mpo.jl index b7ae6f42f..60bfc6db7 100644 --- a/src/algorithms/time_evolution/apply_mpo.jl +++ b/src/algorithms/time_evolution/apply_mpo.jl @@ -229,14 +229,8 @@ function _get_allprojs( N = length(Ms) Rs, Ls = _get_allRLs(Ms) @assert length(truncs) == N - 1 - projs_errs = map(1:(N - 1)) do i - trunc = if isa(truncs[i], FixedSpaceTruncation) - tspace = space(Ms[i + 1], 1) - isdual(tspace) ? truncspace(flip(tspace)) : truncspace(tspace) - else - truncs[i] - end - return _proj_from_RL(Rs[i], Ls[i]; trunc) + projs_errs = map(zip(Rs, Ls, truncs)) do (R, L, trunc) + return _proj_from_RL(R, L; trunc) end Pas = map(Base.Fix2(getindex, 1), projs_errs) wts = map(Base.Fix2(getindex, 2), projs_errs) diff --git a/src/algorithms/time_evolution/gaugefix_su.jl b/src/algorithms/time_evolution/gaugefix_su.jl index 324a8d604..13a2782d6 100644 --- a/src/algorithms/time_evolution/gaugefix_su.jl +++ b/src/algorithms/time_evolution/gaugefix_su.jl @@ -37,22 +37,24 @@ end Fix the gauge of `psi` using trivial simple update. """ function gauge_fix(psi::InfiniteState, alg::SUGauge) + time0 = time() gates = _trivial_gates(scalartype(psi), physicalspace(psi)) - su_alg = SimpleUpdate(; trunc = FixedSpaceTruncation(), bipartite = _state_bipartite_check(psi)) + trunc = _get_fixedspacetrunc(psi) + su_alg = SimpleUpdate(; trunc, bipartite = _is_bipartite(psi)) wts0 = SUWeight(psi) # use default constructor to avoid calculation of exp(-H * 0) evolver = TimeEvolver(su_alg, 0.0, alg.maxiter, gates, SUState(0, 0.0, psi, wts0)) for (i, (psi′, wts, info)) in enumerate(evolver) ϵ = compare_weights(wts, wts0) if i >= alg.miniter && ϵ < alg.tol - @info "Trivial SU conv $i: |Δλ| = $ϵ." + @info "Trivial SU conv $i: |Δλ| = $ϵ, time = $(time() - time0) s" return psi′, wts, ϵ end if i == alg.maxiter - @warn "Trivial SU cancel $i: |Δλ| = $ϵ." + @warn "Trivial SU cancel $i: |Δλ| = $ϵ, time = $(time() - time0) s" return psi′, wts, ϵ end - wts0 = deepcopy(wts) + wts0 = wts end return end diff --git a/src/algorithms/time_evolution/ntupdate.jl b/src/algorithms/time_evolution/ntupdate.jl new file mode 100644 index 000000000..7c4281762 --- /dev/null +++ b/src/algorithms/time_evolution/ntupdate.jl @@ -0,0 +1,277 @@ +""" +$(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 = truncerror(; atol = 1.0e-10)) + "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 +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. +""" +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 + ) + @assert length(sites) == 2 + Nr, Nc = size(state) + trunc = only(_get_cluster_trunc(alg.opt_alg.trunc, sites, (Nr, Nc))) + alg′ = (@set alg.opt_alg.trunc = trunc) + return _bond_truncate(state, wts, Tuple(sites), (:first, :last), alg′; gate) +end + +""" +One round of neighbourhood tensor update +""" +function ntu_iter( + state::InfiniteState, circuit::LocalCircuit, alg::NeighbourUpdate + ) + Nr, Nc, = size(state) + 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 = sites[1] + r, c = mod1(site[1], Nr), mod1(site[2], Nc) + state2[r, c] = _apply_sitegate(state2[r, c], gate) + info′ = (; fid = 1.0) + elseif length(sites) == 2 + (d, r, c), = _nn_bondrev(sites..., (Nr, Nc)) + alg.bipartite && r > 1 && continue + state2, wts, info′ = _ntu_iter(state2, gate, wts, sites, alg) + (!alg.bipartite) && continue + if d == 1 + rp1, cp1 = _next(r, Nr), _next(c, Nc) + state2[rp1, cp1] = copy(state2[r, c]) + state2[rp1, c] = copy(state2[r, cp1]) + wts[1, rp1, cp1] = copy(wts[1, r, c]) + else + rm1, cm1 = _prev(r, Nr), _prev(c, Nc) + state2[rm1, cm1] = copy(state2[r, c]) + state2[r, cm1] = copy(state2[rm1, c]) + wts[2, rm1, cm1] = 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) + 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 + psi, wts, info = ntu_iter(state.psi, it.gate, it.alg) + 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::Float64 = 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::Float64 = 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::TimeEvolver{<:NeighbourUpdate}, + [H::LocalOperator, env::CTMRGEnv, ctm_alg::CTMRGAlgorithm]; + tol::Float64 = 1.0e-7, check_interval::Int = 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}; check_interval::Int = 50) + time_start = time0 = time() + @info "--- 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) + @info @sprintf( + "NTU iter %d: t = %.2e, |Δλ| = %.3e. Time: %.2f s", + it.state.iter, it.state.t, Δλ, time1 - time0 + ) + end + if stop + time_end = time() + @info @sprintf("Time evolution finished in %.2f s", time_end - time_start) + return psi, info + end + info0, time0 = info, time() + end + return +end + +function MPSKit.time_evolve( + it::TimeEvolver{<:NeighbourUpdate, G, S}, + H::LocalOperator, env::CTMRGEnv, ctm_alg::CTMRGAlgorithm; + tol::Float64 = 1.0e-7, check_interval::Int = 10 + ) where {G, S <: NTUState{<:InfinitePEPS}} + @info "--- 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)) + @info @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 + Δλ = 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() + @info @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 + @info "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 + @info "NTU: reached maximum iteration." + end + if stop + time_end = time() + @info @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 + return +end diff --git a/src/algorithms/time_evolution/ntupdate3site.jl b/src/algorithms/time_evolution/ntupdate3site.jl new file mode 100644 index 000000000..fa50cdb21 --- /dev/null +++ b/src/algorithms/time_evolution/ntupdate3site.jl @@ -0,0 +1,116 @@ +""" +Neighbourhood tensor update with N-site MPO `gate` (N ≥ 2). +""" +function _ntu_iter( + state::InfiniteState, gate::Vector{T}, wts::SUWeight, + sites::Vector{CartesianIndex{2}}, alg::NeighbourUpdate + ) where {T <: AbstractTensorMap} + Nr, Nc = size(state) + truncs = _get_cluster_trunc(alg.opt_alg.trunc, sites, (Nr, Nc)) + state, wts = copy(state), deepcopy(wts) + + Ms, _, invperms = _get_cluster(state, sites) + flips = [isdual(space(M, 1)) for M in Ms[2:end]] + _flip_virtuals!(Ms, flips) # flip virtual arrows in `Ms` to ← + + # apply gate MPO without truncation + _apply_gatempo!(Ms, gate) + _flip_virtuals!(Ms, flips) # restore virtual arrows in `Ms` + for (M, s, invperm) in zip(Ms, sites, invperms) + s′ = CartesianIndex(mod1(s[1], Nr), mod1(s[2], Nc)) + state[s′] = permute(M, invperm) + end + + # truncate each bond sequentially along the path + info = (; fid = 1.0) + nbond = length(sites) - 1 + for (i, bondsites) in enumerate(zip(sites, Iterators.drop(sites, 1))) + trunc = truncs[i] + alg′ = (@set alg.opt_alg.trunc = trunc) + stype1 = (i == 1) ? :first : :middle + stype2 = (i == nbond) ? :last : :middle + state, wts, info′ = _bond_truncate(state, wts, bondsites, (stype1, stype2), alg′) + # 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` +after rotating the bond to standard x direction `A ← B`. + +`bondtype` takes values in (1, 2, 3), meaning that the current bond is +(the first, a middle, the last) bond in the updated cluster. +""" +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 + ) + # rotate bond to standard x direction `A ← B` + ucell = size(state)[1:2] + bond, rev = _nn_bondrev(site1, site2, ucell) + state2 = _bond_rotation(state, bond[1], rev; inv = false) + wts2 = _bond_rotation(wts, bond[1], rev; inv = false) + + # rotated bond tensors + siteA = _bond_rotation(site1, bond[1], rev, ucell) + row, col = mod1.(Tuple(siteA), size(state2)[1:2]) + cp1 = _next(col, size(state2, 2)) + A, B = state2[row, col], state2[row, cp1] + + # create bond environment + qrtrunc = trunctol(; rtol = 1.0e-12) + a, X = if stype1 == :first + bond_tensor_first(A; trunc = qrtrunc) + else + @assert stype1 == :middle + bond_tensor_midnext(A; trunc = qrtrunc) + end + b, Y = if stype2 == :last + bond_tensor_last(B; trunc = qrtrunc) + else + @assert stype2 == :middle + bond_tensor_midprev(B; trunc = qrtrunc) + end + benv = bondenv_ntu(row, col, X, Y, state2, alg.bondenv_alg) + @debug "cond(benv) before gauge fix: $(LinearAlgebra.cond(benv))" + if alg.fixgauge + Z = positive_approx(benv) + Z, a, b, (Linv, Rinv) = fixgauge_benv(Z, a, b) + 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, s, b, = _apply_gate(a, b, gate, truncerror(; atol = 1.0e-15)) + end + a, s, b, info = bond_truncate(a, b, benv, alg.opt_alg) + + A = if stype1 == :first + undo_bond_tensor_first(a, X) + else + undo_bond_tensor_midnext(a, X) + end + B = if stype2 == :last + undo_bond_tensor_last(b, Y) + else + undo_bond_tensor_midprev(b, Y) + end + + normalize!(A, Inf) + normalize!(B, Inf) + normalize!(s, Inf) + state2[row, col], state2[row, cp1], wts2[1, row, col] = A, B, s + + # rotate back tensors and bond weight + state2 = _bond_rotation(state2, bond[1], rev; inv = true) + wts2 = _bond_rotation(wts2, bond[1], rev; inv = true) + return state2, wts2, info +end diff --git a/src/algorithms/time_evolution/simpleupdate.jl b/src/algorithms/time_evolution/simpleupdate.jl index a28f86d8b..5c3d30a19 100644 --- a/src/algorithms/time_evolution/simpleupdate.jl +++ b/src/algorithms/time_evolution/simpleupdate.jl @@ -7,19 +7,19 @@ Algorithm struct for simple update (SU) of InfinitePEPS or InfinitePEPO. $(TYPEDFIELDS) """ -@kwdef struct SimpleUpdate <: TimeEvolution +@kwdef struct SimpleUpdate{T <: TruncationStrategy} <: TimeEvolution "Truncation strategy for bonds updated by Trotter gates" - trunc::TruncationStrategy + trunc::T "When true (or false), the Trotter gate is `exp(-H dt)` (or `exp(-iH dt)`)" imaginary_time::Bool = true "When true, force decomposition of nearest neighbor gates to MPOs." force_mpo::Bool = false "When true, assume bipartite unit cell structure" bipartite::Bool = false - "(Only applicable to InfinitePEPO) + "(Only applicable to InfinitePEPO) When true, the PEPO is regarded as a purified PEPS, and updated as `|ρ(t + dt)⟩ = exp(-H dt/2) |ρ(t)⟩`. - When false, the PEPO is updated as + When false, the PEPO is updated as `ρ(t + dt) = exp(-H dt/2) ρ(t) exp(-H dt/2)`." purified::Bool = true end @@ -62,6 +62,12 @@ function TimeEvolver( # create Trotter gates gate = trotterize(H, dt′; symmetrize_gates, force_mpo = alg.force_mpo) state = SUState(0, t0, psi0, env0) + # convert FixedSpaceTruncation to site-dependent `truncspace`s + if alg.trunc isa FixedSpaceTruncation + trunc = _get_fixedspacetrunc(psi0) + @reset alg.trunc = trunc + end + # TODO: bipartite check for alg.trunc after equality is defined for all kinds of truncation strategies # TODO: check gates for bipartite case return TimeEvolver(alg, dt, nstep, gate, state) end @@ -79,6 +85,13 @@ function _bond_rotation(x, bonddir::Int, rev::Bool; inv::Bool = false) error("`bonddir` must be 1 (for x-bonds) or 2 (for y-bonds).") end end +function _bond_rotation(x::CartesianIndex{2}, bonddir::Int, rev::Bool, unitcell::NTuple{2, Int}) + return if bonddir == 1 + rev ? siterot180(x, unitcell) : x + else + rev ? siterotl90(x, unitcell) : siterotr90(x, unitcell) + end +end """ Simple update optimized for nearest neighbor gates @@ -89,8 +102,8 @@ function _su_iter!( sites::Vector{CartesianIndex{2}}, alg::SimpleUpdate ) Nr, Nc = size(state) - truncs = _get_cluster_trunc(alg.trunc, sites, (Nr, Nc)) - @assert length(sites) == 2 && length(truncs) == 1 + @assert length(sites) == 2 + trunc = only(_get_cluster_trunc(alg.trunc, sites, (Nr, Nc))) Ms, open_vaxs, = _get_cluster(state, sites, env; permute = false) normalize!.(Ms, Inf) # rotate @@ -100,10 +113,12 @@ function _su_iter!( ϵ, s = 0.0, nothing gate_axs = alg.purified ? (1:1) : (1:2) for gate_ax in gate_axs - X, a, b, Y = _qr_bond(A, B; gate_ax, positive = true) - a, s, b, ϵ′ = _apply_gate(a, b, gate, truncs[1]) + a, X = bond_tensor_first(A; gate_ax, positive = true) + b, Y = bond_tensor_last(B; gate_ax, positive = true) + a, s, b, ϵ′ = _apply_gate(a, b, gate, trunc) ϵ = max(ϵ, ϵ′) - A, B = _qr_bond_undo(X, a, b, Y) + A = undo_bond_tensor_first(a, X; gate_ax) + B = undo_bond_tensor_last(b, Y; gate_ax) end # rotate back A = _bond_rotation(A, bond[1], rev; inv = true) @@ -148,14 +163,14 @@ function su_iter( (!alg.bipartite) && continue if d == 1 rp1, cp1 = _next(r, Nr), _next(c, Nc) - state2[rp1, cp1] = deepcopy(state2[r, c]) - state2[rp1, c] = deepcopy(state2[r, cp1]) - env2[1, rp1, cp1] = deepcopy(env2[1, r, c]) + state2[rp1, cp1] = copy(state2[r, c]) + state2[rp1, c] = copy(state2[r, cp1]) + env2[1, rp1, cp1] = copy(env2[1, r, c]) else rm1, cm1 = _prev(r, Nr), _prev(c, Nc) - state2[rm1, cm1] = deepcopy(state2[r, c]) - state2[r, cm1] = deepcopy(state2[rm1, c]) - env2[2, rm1, cm1] = deepcopy(env2[2, r, c]) + state2[rm1, cm1] = copy(state2[r, c]) + state2[r, cm1] = copy(state2[rm1, c]) + env2[2, rm1, cm1] = copy(env2[2, r, c]) end else # N-site MPO gate (N ≥ 2) @@ -202,10 +217,8 @@ function MPSKit.timestep( end """ - time_evolve( - it::TimeEvolver{<:SimpleUpdate}; - tol::Float64 = 0.0, check_interval::Int = 500 - ) -> (psi, env, info) + time_evolve(it; check_interval = 500) -> (psi, env, info) + time_evolve(it, H; tol = 1.0e-8, check_interval = 500) -> (psi, env, info) Perform time evolution to the end of `TimeEvolver` iterator `it`, or until convergence of `SUWeight` set by a positive `tol`. @@ -215,15 +228,41 @@ or until convergence of `SUWeight` set by a positive `tol`. - `check_interval` sets the number of iterations between outputs of information. """ function MPSKit.time_evolve( - it::TimeEvolver{<:SimpleUpdate}; - tol::Float64 = 0.0, check_interval::Int = 500 + it::TimeEvolver{<:SimpleUpdate}; check_interval::Int = 500 ) time_start = time() - check_convergence = (tol > 0) @info "--- Time evolution (simple update), dt = $(it.dt) ---" - if check_convergence - @assert (it.state.psi isa InfinitePEPS) && it.alg.imaginary_time "Only imaginary time evolution of InfinitePEPS allows convergence checking." + env0, time0 = it.state.env, time() + for (psi, env, info) in it + iter = it.state.iter + diff = compare_weights(env0, env) + stop = (iter == it.nstep) + showinfo = (check_interval > 0) && + ((iter % check_interval == 0) || (iter == 1) || stop) + time1 = time() + if showinfo + @info "Space of x-weight at [1, 1] = $(space(env[1, 1, 1], 1))" + @info @sprintf("SU iter %-7d: |Δλ| = %.3e. Time = %.3f s/it", iter, diff, time1 - time0) + end + if stop + time_end = time() + @info @sprintf("Time evolution finished in %.2f s", time_end - time_start) + return psi, env, info + else + env0 = env + end + time0 = time() end + return +end + +function MPSKit.time_evolve( + it::TimeEvolver{<:SimpleUpdate, G, S}, H::LocalOperator; + tol::Float64 = 1.0e-8, check_interval::Int = 500 + ) where {G, S <: SUState{<:InfinitePEPS}} + time_start = time() + @info "--- Time evolution (simple update), dt = $(it.dt) ---" + @assert it.alg.imaginary_time "Only imaginary time evolution of InfinitePEPS allows convergence checking." env0, time0 = it.state.env, time() for (psi, env, info) in it iter = it.state.iter @@ -233,16 +272,20 @@ function MPSKit.time_evolve( ((iter % check_interval == 0) || (iter == 1) || stop) time1 = time() if showinfo + # TODO: convert to BPEnv instead + ctmenv = CTMRGEnv(env) + energy = real(expectation_value(psi, H, ctmenv)) / prod(size(psi)) @info "Space of x-weight at [1, 1] = $(space(env[1, 1, 1], 1))" - @info @sprintf("SU iter %-7d: |Δλ| = %.3e. Time = %.3f s/it", iter, diff, time1 - time0) + @info @sprintf( + "SU iter %-7d: E ≈ %.5f, |Δλ| = %.3e. Time = %.3f s/it", + iter, energy, diff, time1 - time0 + ) end - if check_convergence - if (iter == it.nstep) && (diff >= tol) - @warn "SU: bond weights have not converged." - end - if diff < tol - @info "SU: bond weights have converged." - end + if (iter == it.nstep) && (diff >= tol) + @warn "SU: bond weights have not converged." + end + if diff < tol + @info "SU: bond weights have converged." end if stop time_end = time() @@ -269,7 +312,6 @@ algorithm `alg`, time step `dt` for `nstep` number of steps. - Set `symmetrize_gates = true` for second-order Trotter decomposition. - Set `tol > 0` to enable convergence check (for imaginary time evolution of iPEPS only). - For other usages it should not be changed. - Use `t0` to specify the initial time of the evolution. - `check_interval` sets the interval to output information. Output during the evolution can be turned off by setting `check_interval <= 0`. - `info` is a NamedTuple containing information of the evolution, @@ -281,5 +323,9 @@ function MPSKit.time_evolve( tol::Float64 = 0.0, t0::Number = 0.0, check_interval::Int = 500 ) it = TimeEvolver(psi0, H, dt, nstep, alg, env0; t0, symmetrize_gates) - return time_evolve(it; tol, check_interval) + return if tol == 0 + time_evolve(it; check_interval) + else + time_evolve(it, H; tol, check_interval) + end end diff --git a/src/algorithms/time_evolution/simpleupdate3site.jl b/src/algorithms/time_evolution/simpleupdate3site.jl index 4b091c830..2fc9d096f 100644 --- a/src/algorithms/time_evolution/simpleupdate3site.jl +++ b/src/algorithms/time_evolution/simpleupdate3site.jl @@ -217,24 +217,14 @@ updated by the Trotter evolution MPO. """ function _get_cluster_trunc( trunc::TruncationStrategy, sites::Vector{CartesianIndex{2}}, - (Nrow, Ncol)::NTuple{2, Int} + unitcell::NTuple{2, Int} ) return map(zip(sites, Iterators.drop(sites, 1))) do (site1, site2) - diff = site2 - site1 - if diff == CartesianIndex(0, 1) - r, c = mod1(site1[1], Nrow), mod1(site1[2], Ncol) - return truncation_strategy(trunc, 1, r, c) - elseif diff == CartesianIndex(0, -1) - r, c = mod1(site2[1], Nrow), mod1(site2[2], Ncol) - return truncation_strategy(trunc, 1, r, c) - elseif diff == CartesianIndex(1, 0) - r, c = mod1(site2[1], Nrow), mod1(site2[2], Ncol) - return truncation_strategy(trunc, 2, r, c) - elseif diff == CartesianIndex(-1, 0) - r, c = mod1(site1[1], Nrow), mod1(site1[2], Ncol) - return truncation_strategy(trunc, 2, r, c) - else - error("The path `sites` contains a long-range bond.") + (d, r, c), rev = _nn_bondrev(site1, site2, unitcell) + t = truncation_strategy(trunc, d, r, c) + if rev && isa(t, TruncationSpace) + t = truncspace(flip(t.space)') end + return t end end diff --git a/src/algorithms/time_evolution/time_evolve.jl b/src/algorithms/time_evolution/time_evolve.jl index ee5f236c1..dd51c2f54 100644 --- a/src/algorithms/time_evolution/time_evolve.jl +++ b/src/algorithms/time_evolution/time_evolve.jl @@ -30,22 +30,6 @@ end Base.iterate(it::TimeEvolver) = iterate(it, it.state) -function _state_bipartite_check(psi::InfiniteState) - if isa(psi, InfinitePEPO) - @assert size(psi, 3) == 1 "Input InfinitePEPO is expected to have only one layer." - end - if !(size(psi, 1) == size(psi, 2) == 2) - return false - end - for (r, c) in Iterators.product(1:2, 1:2) - r′, c′ = _next(r, 2), _next(c, 2) - if psi[r, c] != psi[r′, c′] - return false - end - end - return true -end - """ Process the Trotter time step `dt` according to the intended usage. """ @@ -79,7 +63,7 @@ function _timeevol_sanity_check( @assert ψ₀ isa InfinitePEPO "alg.purified = false is only applicable to PEPO." end if hasfield(typeof(alg), :bipartite) && alg.bipartite - @assert _state_bipartite_check(ψ₀) "Input state is not bipartite with 2 x 2 unit cell." + @assert _is_bipartite(ψ₀) "Input state is not bipartite with 2 x 2 unit cell." end return nothing end @@ -94,3 +78,21 @@ function MPSKit.infinite_temperature_density_matrix(H::LocalOperator) end return InfinitePEPO(cat(A; dims = 3)) end + +""" +Get the `SiteDependentTruncation` used by time evolution +that preserves virtual spaces of `state`. +""" +function _get_fixedspacetrunc(state::InfiniteState) + if state isa InfinitePEPO + size(state, 3) != 1 && error("Input InfinitePEPO is expect to have only one layer.") + end + Nr, Nc = size(state) + return SiteDependentTruncation( + map(Iterators.product(1:2, 1:Nr, 1:Nc)) do (d, r, c) + V = domain(state[r, c], (d == 1) ? EAST : NORTH) + isdual(V) && (V = flip(V)) + return truncspace(V) + end + ) +end diff --git a/src/algorithms/time_evolution/trotter_gate.jl b/src/algorithms/time_evolution/trotter_gate.jl index 24248ad9c..67ae916d3 100644 --- a/src/algorithms/time_evolution/trotter_gate.jl +++ b/src/algorithms/time_evolution/trotter_gate.jl @@ -1,5 +1,3 @@ -const NNGate{T, S} = AbstractTensorMap{T, S, 2, 2} - """ Convert an N-site gate (N ≥ 2) to MPO by SVD, in which the axes are ordered as diff --git a/src/algorithms/truncation/bond_tensor.jl b/src/algorithms/truncation/bond_tensor.jl new file mode 100644 index 000000000..6b11fea2b --- /dev/null +++ b/src/algorithms/truncation/bond_tensor.jl @@ -0,0 +1,213 @@ +""" +Given the first tensor `A` in the cluster acted on by a gate, +obtain reduced tensor on its next bond. + +For PEPSTensor, +``` + 1 + | + 4 - X ← 2 1 ← a - 3 + | ↘ + 3 2 +``` +For PEPOTensor, +``` + gate_ax = 1 gate_ax = 2 + + 1 2 2 2 + ↘ | | ↘ + 5 - X ← 3 1 ← a - 3 5 - X ← 3 1 ← a - 3 + | ↘ | ↘ + 4 2 4 1 +``` +""" +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 +end + +""" +Undo the decomposition in `bond_tensor_first`. +""" +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 +end + +""" +Given the last tensor `A` in the cluster acted on by a gate, +obtain reduced tensor on its previous bond. + +For PEPSTensor, +``` + 1 + | + 1 - b → 3 4 → Y - 2 + ↘ | + 2 3 +``` +For PEPOTensor, +``` + gate_ax = 1 gate_ax = 2 + + 1 2 2 2 + ↘ | ↘ | + 1 - b → 3 5 → Y - 3 1 - b → 3 5 → Y - 3 + ↘ | | ↘ + 2 4 4 1 +``` +""" +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 +end + +""" +Undo the decomposition in `bond_tensor_last`. +""" +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 +end + +""" +Given a middle tensor `A` in the cluster acted on by a gate, +obtain reduced tensor on its next bond. + +For PEPSTensor, +``` + 2 + | + 5 - X ← 3 1 ← a - 3 + | ↘ ↘ + 4 1 (2) +``` +For PEPOTensor, +``` + 2 3 + ↘ | + 6 - X ← 4 1 ← a - 3 + | ↘ ↘ + 5 1 (2) +``` + +Here the physical leg on `a` is an auxiliary trivial leg. +""" +function bond_tensor_midnext(A::PEPSTensor; kwargs...) + X, a = left_orth!(permute(A, ((1, 2, 4, 5), (3,)); copy = true); kwargs...) + X = permute(X, ((1,), (2, 5, 3, 4))) + a = insertrightunit(a, 1) + return a, X +end +function bond_tensor_midnext(A::PEPOTensor; kwargs...) + X, a = left_orth!(permute(A, ((1, 2, 3, 5, 6), (4,)); copy = true); kwargs...) + X = permute(X, ((1, 2), (3, 6, 4, 5))) + a = insertrightunit(a, 1) + return a, X +end + +""" +Undo the decomposition in `bond_tensor_midprev`. +""" +function undo_bond_tensor_midnext(a::MPSTensor, X::PEPSTensor) + a = removeunit(a, 2) + return @tensor A[-1; -2 -3 -4 -5] := X[-1; -2 1 -4 -5] * a[1; -3] +end +function undo_bond_tensor_midnext(a::MPSTensor, X::PEPOTensor) + a = removeunit(a, 2) + return @tensor A[-1 -2; -3 -4 -5 -6] := X[-1 -2; -3 1 -5 -6] * a[1; -4] +end + +""" +Given a middle tensor `A` in the cluster acted on by a gate, +obtain reduced tensor on its previous bond. + +For PEPSTensor, +``` + 2 + | + 1 - b → 3 5 → Y - 3 + ↘ | ↘ + (2) 4 1 +``` +For PEPOTensor, +``` + 2 3 + ↘ | + 1 - b → 3 6 → Y - 4 + ↘ | ↘ + (2) 5 1 +``` + +Here the physical leg on `a` is an auxiliary trivial leg. +""" +function bond_tensor_midprev(A::PEPSTensor; kwargs...) + Y, b = left_orth!(permute(A, ((1, 2, 3, 4), (5,)); copy = true); kwargs...) + Y = permute(Y, ((1,), (2, 3, 4, 5))) + b = permute(b, ((2,), (1,))) + b = insertrightunit(b, 1) + return b, Y +end +function bond_tensor_midprev(A::PEPOTensor; kwargs...) + Y, b = left_orth!(permute(A, ((1, 2, 3, 4, 5), (6,)); copy = true); kwargs...) + Y = permute(Y, ((1, 2), (3, 4, 5, 6))) + b = permute(b, ((2,), (1,))) + b = insertrightunit(b, 1) + return b, Y +end + +""" +Undo the decomposition in `bond_tensor_midprev`. +""" +function undo_bond_tensor_midprev(b::MPSTensor, Y::PEPSTensor) + b = removeunit(b, 2) + return @tensor A[-1; -2 -3 -4 -5] := b[-5; 1] * Y[-1; -2 -3 -4 1] +end +function undo_bond_tensor_midprev(b::MPSTensor, Y::PEPOTensor) + b = removeunit(b, 2) + return @tensor A[-1 -2; -3 -4 -5 -6] := b[-6; 1] * Y[-1 -2; -3 -4 -5 1] +end diff --git a/src/algorithms/truncation/bond_truncation.jl b/src/algorithms/truncation/bond_truncation.jl index 57fd68cc7..1a95b8d34 100644 --- a/src/algorithms/truncation/bond_truncation.jl +++ b/src/algorithms/truncation/bond_truncation.jl @@ -18,8 +18,8 @@ The truncation algorithm can be constructed from the following keyword arguments * `tol::Float64=1e-9` : ALS converges when the relative change in bond SVD spectrum between two iterations is smaller than `tol`. * `check_interval::Int=0` : Set number of iterations to print information. Output is suppressed when `check_interval <= 0`. """ -@kwdef struct ALSTruncation - trunc::TruncationStrategy +@kwdef struct ALSTruncation{T <: TruncationStrategy} + trunc::T maxiter::Int = 50 tol::Float64 = 1.0e-9 check_interval::Int = 0 @@ -35,7 +35,21 @@ function _als_message( end """ - bond_truncate(a::AbstractTensorMap{T,S,2,1}, b::AbstractTensorMap{T,S,1,2}, benv::BondEnv{T,S}, alg) -> U, S, V, info +Initialize truncated bond tensors for 2-site ALS +""" +function _als_init_truncate( + ket2::AbstractTensorMap{T, S, 2, 2}, trunc::TruncationStrategy + ) where {T, S} + a, s0, b = svd_trunc!(permute(ket2, ((1, 3), (4, 2)); copy = true); trunc) + a, b = absorb_s(a, s0, b) + # put b in MPS axis order + b = permute(b, ((1, 2), (3,))) + xs = [a, b] + return xs, s0 +end + +""" + bond_truncate(a::MPSTensor, b::MPSTensor, benv::BondEnv, alg) -> U, S, V, info After time-evolving the reduced tensors `a` and `b` connected by a bond, truncate the bond dimension using the bond environment tensor `benv`. @@ -49,19 +63,14 @@ truncate the bond dimension using the bond environment tensor `benv`. └-----------------------┘ ``` The truncation algorithm `alg` can be either `FullEnvTruncation` or `ALSTruncation`. -The index order of `a` or `b` is +The index order of `a` or `b` follows `MPSTensor` convention. ``` 1 -a/b- 3 - ↓ a[1 2; 3] - 2 b[1; 2 3] + ↓ [1 2; 3] + 2 ``` """ -function bond_truncate( - a::AbstractTensorMap{T, S, 2, 1}, - b::AbstractTensorMap{T, S, 1, 2}, - benv::BondEnv{T, S}, - alg::ALSTruncation, - ) where {T <: Number, S <: ElementarySpace} +function bond_truncate(a::MPSTensor, b::MPSTensor, benv::BondEnv, alg::ALSTruncation) # dual check of physical index @assert !isdual(space(a, 2)) @assert !isdual(space(b, 2)) @@ -69,40 +78,39 @@ function bond_truncate( need_flip = isdual(space(b, 1)) time00 = time() verbose = (alg.check_interval > 0) - a2b2 = _combine_ab(a, b) - # initialize truncated a, b - perm_ab = ((1, 3), (4, 2)) - a, s0, b = svd_trunc(permute(a2b2, perm_ab); trunc = alg.trunc) - a, b = absorb_s(a, s0, b) - # put b in MPS axis order - b = permute(b, ((1, 2), (3,))) - ab = _combine_ab(a, b) + + # untruncated things + ket2 = _combine_ket(a, b) + benv_ket2 = _benv_ket(benv, ket2) + b22 = _als_norm(ket2, benv_ket2) + + # initialize truncated bond tensors and bond weight + xs, s0 = _als_init_truncate(ket2, alg.trunc) + + # initialize ALS cache + Rs = [_als_tensor_R(benv, xs, i) for i in 1:2] + Ss = [_als_tensor_S(benv_ket2, xs, i) for i in 1:2] + # cost function will be normalized by initial value - cost00, fid = cost_function_als(benv, ab, a2b2) + cost00, fid = cost_function_als(Rs[1], Ss[1], xs[1], b22) cost0, fid0, Δcost, Δfid, Δs = cost00, fid, NaN, NaN, NaN verbose && @info "ALS init" * _als_message(0, cost0, fid, Δcost, Δfid, Δs, 0.0) + for iter in 1:(alg.maxiter) time0 = time() - #= - Fixing `b`, the cost function can be expressed in the R, S tensors as - ``` - f(a†,a) = a† Ra a - a† Sa - Sa† a + const - ``` - `f` is minimized when - ∂f/∂ā = Ra a - Sa = 0 - =# - Ra = _tensor_Ra(benv, b) - Sa = _tensor_Sa(benv, b, a2b2) - a, info_a = _solve_als(Ra, Sa, a) - # Fixing `a`, solve for `b` from `Rb b = Sb` - Rb = _tensor_Rb(benv, a) - Sb = _tensor_Sb(benv, a, a2b2) - b, info_b = _solve_als(Rb, Sb, b) - @debug "Bond truncation info" info_a info_b - ab = _combine_ab(a, b) - cost, fid = cost_function_als(benv, ab, a2b2) + for (i, (Rx, Sx, x)) in enumerate(zip(Rs, Ss, xs)) + # TODO: option to use pinv + xs[i], info_x = _solve_als(Rx, Sx, x) + @debug "Bond truncation info $(i):" info_x + # update R, S for the next site + i_next = _next(i, 2) + Rs[i_next] = _als_tensor_R(benv, xs, i_next) + Ss[i_next] = _als_tensor_S(benv_ket2, xs, i_next) + end + # cost function and local fidelity + cost, fid = cost_function_als(Rs[1], Ss[1], xs[1], b22) # TODO: replace with truncated svdvals (without calculating u, vh) - _, s, _ = svd_trunc!(permute(ab, perm_ab); trunc = alg.trunc) + _, s, _ = svd_trunc!(_combine_ket_for_svd(xs...); trunc = alg.trunc) # fidelity, cost and normalized bond-s change s_nrm = norm(s0, Inf) Δs = _singular_value_distance(s, s0) / s_nrm @@ -129,20 +137,16 @@ function bond_truncate( end converge && break end - a, s, b = svd_trunc!(permute(_combine_ab(a, b), perm_ab); trunc = alg.trunc) + a, s, b = svd_trunc!(_combine_ket_for_svd(xs...); trunc = alg.trunc) a, b = absorb_s(a, s, b) + b = permute(b, ((1, 2), (3,))) if need_flip a, s, b = flip(a, numind(a)), _fliptwist_s(s), flip(b, 1) end return a, s, b, (; fid, Δfid, Δs) end -function bond_truncate( - a::AbstractTensorMap{T, S, 2, 1}, - b::AbstractTensorMap{T, S, 1, 2}, - benv::BondEnv{T, S}, - alg::FullEnvTruncation, - ) where {T <: Number, S <: ElementarySpace} +function bond_truncate(a::MPSTensor, b::MPSTensor, benv::BondEnv, alg::FullEnvTruncation) # dual check of physical index @assert !isdual(space(a, 2)) @assert !isdual(space(b, 2)) @@ -150,14 +154,14 @@ function bond_truncate( need_flip = isdual(space(b, 1)) #= initialize bond matrix using QR as `Ra Lb` - --- a == b --- ==> - Qa ← Ra == Rb ← Qb - + --- a == b --- ==> - Qa ← Ra == Rb → Qb - ↓ ↓ ↓ ↓ =# Qa, Ra = left_orth(a; positive = true) - Rb, Qb = right_orth(b; positive = true) - @assert !isdual(space(Ra, 1)) && !isdual(space(Qb, 1)) - @tensor b0[-1; -2] := Ra[-1 1] * Rb[1 -2] - #= initialize bond environment around `Ra Lb` + b = permute(b, ((3, 2), (1,)); copy = true) + Qb, Rb = left_orth!(b; positive = true) + @tensor b0[-1; -2] := Ra[-1 1] * Rb[-2 1] + #= initialize bond environment around `Ra Rb` ┌--------------------------------------┐ | ┌----┐ | @@ -169,15 +173,14 @@ function bond_truncate( | └----┘ | └--------------------------------------┘ =# - @tensor benv2[-1 -2; -3 -4] := ( - benv[1 2; 3 4] * conj(Qa[1 5 -1]) * conj(Qb[-2 6 2]) * Qa[3 5 -3] * Qb[-4 6 4] - ) + @tensor benv2[-1 -2; -3 -4] := benv[1 2; 3 4] * + conj(Qa[1 5 -1]) * conj(Qb[2 6 -2]) * Qa[3 5 -3] * Qb[4 6 -4] # optimize bond matrix u, s, vh, info = fullenv_truncate(b0, benv2, alg) u, vh = absorb_s(u, s, vh) # truncate a, b tensors with u, s, vh @tensor a[-1 -2; -3] := Qa[-1 -2 3] * u[3 -3] - @tensor b[-1; -2 -3] := vh[-1 1] * Qb[1 -2 -3] + @tensor b[-1 -2; -3] := vh[-1 1] * Qb[-3 -2 1] if need_flip a, s, b = flip(a, numind(a)), _fliptwist_s(s), flip(b, 1) end diff --git a/src/algorithms/truncation/fullenv_truncation.jl b/src/algorithms/truncation/fullenv_truncation.jl index fd9c685bd..4ff6bac2b 100644 --- a/src/algorithms/truncation/fullenv_truncation.jl +++ b/src/algorithms/truncation/fullenv_truncation.jl @@ -23,8 +23,8 @@ The truncation algorithm can be constructed from the following keyword arguments * [Glen Evenbly, Phys. Rev. B 98, 085155 (2018)](@cite evenbly_gauge_2018). """ -@kwdef struct FullEnvTruncation - trunc::TruncationStrategy +@kwdef struct FullEnvTruncation{T <: TruncationStrategy} + trunc::T maxiter::Int = 50 tol::Float64 = 1.0e-9 trunc_init::Bool = true @@ -75,13 +75,13 @@ function _fet_message( end """ - fullenv_truncate(benv::BondEnv{T,S}, b0::AbstractTensorMap{T,S,1,1}, alg::FullEnvTruncation) -> U, S, V, info + fullenv_truncate(b0, benv::BondEnv, alg::FullEnvTruncation) -> U, S, V, info Perform full environment truncation algorithm from [Phys. Rev. B 98, 085155 (2018)](@cite evenbly_gauge_2018) on `benv`. -Given a fixed state `|b0⟩` with bond matrix `b0` -and the corresponding positive-definite bond environment `benv`, +Given a fixed state `|b0⟩` with bond matrix `b0` and the +corresponding positive-definite bond environment `benv`, find the state `|b⟩` with truncated bond matrix `b = u s v†` that maximizes the fidelity (not normalized by `⟨b0|b0⟩`) ``` @@ -215,11 +215,12 @@ function fullenv_truncate( b1 = similar(b0) s0 = deepcopy(s) Δfid, Δs, fid, fid0 = NaN, NaN, 0.0, 0.0 + @tensor benv_b0[-1 -2] := benv[-1 -2; 3 4] * b0[3; 4] for iter in 1:(alg.maxiter) time0 = time() # update `← r - = ← s ← v† -` @tensor r[-1 -2] := s[-1; 1] * vh[1; -2] - @tensor p[-1 -2] := conj(u[1; -1]) * benv[1 -2; 3 4] * b0[3; 4] + @tensor p[-1 -2] := conj(u[1; -1]) * benv_b0[1 -2] @tensor B[-1 -2; -3 -4] := conj(u[1; -1]) * benv[1 -2; 3 -4] * u[3; -3] _linearmap_twist!(p) _linearmap_twist!(B) @@ -228,7 +229,7 @@ function fullenv_truncate( u, s, vh = svd_trunc(b1; trunc = alg.trunc) # update `- l ← = - u ← s ←` @tensor l[-1 -2] := u[-1; 1] * s[1; -2] - @tensor p[-1 -2] := conj(vh[-2; 2]) * benv[-1 2; 3 4] * b0[3; 4] + @tensor p[-1 -2] := conj(vh[-2; 2]) * benv_b0[-1 2] @tensor B[-1 -2; -3 -4] := conj(vh[-2; 2]) * benv[-1 2; -3 4] * vh[-4; 4] _linearmap_twist!(p) _linearmap_twist!(B) diff --git a/src/algorithms/truncation/truncationschemes.jl b/src/algorithms/truncation/truncationschemes.jl index 71b5e9dbe..228f9065d 100644 --- a/src/algorithms/truncation/truncationschemes.jl +++ b/src/algorithms/truncation/truncationschemes.jl @@ -1,16 +1,41 @@ """ $(TYPEDEF) -CTMRG specific truncation strategy for `svd_trunc` which keeps the bond space on which the SVD -is performed fixed. Since different environment directions and unit cell entries might -have different spaces, this truncation style is different from `TruncationSpace`. +SVD truncation strategy which preserves the `CTMRGEnv` environment virtual spaces, +or `InfinitePEPS`, `InfinitePEPO` virtual spaces. """ struct FixedSpaceTruncation <: TruncationStrategy end +""" +$(TYPEDEF) + +SVD truncation strategy specified for each nearest neighbor bond +in `InfinitePEPS`, `InfinitePEPO`: +- `trunc[1, r, c]` applies to the x-bond between `[r, c]` and `[r, c+1]`. + If it is a `TruncationSpace`, the space refers to the east domain of + `[r, c]` or its `flip`, whichever is non-dual. +- `trunc[2, r, c]` applies to the y-bond between `[r, c]` and `[r-1, c]`. + If it is a `TruncationSpace`, the space refers to the north domain of + `[r, c]` or its `flip`, whichever is non-dual. +""" struct SiteDependentTruncation{T <: TruncationStrategy} <: TruncationStrategy truncs::Array{T, 3} + + function SiteDependentTruncation(truncs::Array{T, 3}) where {T} + # TODO: generalize it to CTMRGEnv + size(truncs, 1) != 2 && throw( + DimensionMismatch( + "The first dimension of `truncs` must have a size of 2. Got $(size(truncs, 1))." + ) + ) + return new{T}(truncs) + end end +Base.getindex(trunc::SiteDependentTruncation, args...) = Base.getindex(trunc.truncs, args...) + +# TODO: _is_bipartite(trunc::SiteDependentTruncation) + const TRUNCATION_STRATEGY_SYMBOLS = IdDict{Symbol, Type{<:TruncationStrategy}}( :notrunc => MatrixAlgebraKit.NoTruncation, :truncerror => MatrixAlgebraKit.TruncationByError, @@ -40,28 +65,5 @@ end function truncation_strategy( trunc::SiteDependentTruncation, direction::Int, row::Int, col::Int ) - return trunc.truncs[direction, row, col] -end - -# TODO: type piracy -Base.rotl90(trunc::TruncationStrategy) = trunc - -function Base.rotl90(trunc::SiteDependentTruncation) - directions, rows, cols = size(trunc.truncs) - truncs_rotated = similar(trunc.truncs, directions, cols, rows) - - if directions == 2 - truncs_rotated[NORTH, :, :] = circshift( - rotl90(trunc.truncs[EAST, :, :]), (0, -1) - ) - truncs_rotated[EAST, :, :] = rotl90(trunc.truncs[NORTH, :, :]) - elseif directions == 4 - for dir in 1:4 - dir′ = _prev(dir, 4) - truncs_rotated[dir′, :, :] = rotl90(trunc.truncs[dir, :, :]) - end - else - throw(ArgumentError("Unsupported number of directions for rotl90: $directions")) - end - return SiteDependentTruncation(truncs_rotated) + return trunc[direction, row, col] end diff --git a/src/environments/bp_environments.jl b/src/environments/bp_environments.jl index 765d31c0c..ea7ae9301 100644 --- a/src/environments/bp_environments.jl +++ b/src/environments/bp_environments.jl @@ -158,6 +158,15 @@ function eachcoordinate(x::BPEnv, dirs) return collect(Iterators.product(dirs, axes(x, 2), axes(x, 3))) end +## Bipartite check +function _is_bipartite(env::BPEnv) + (size(env, 2) == size(env, 3) == 2) || (return false) + for (d, c) in Iterators.product(axes(env, 1), axes(env, 3)) + (env[d, 1, c] == env[d, 2, _next(c, 2)]) || (return false) + end + return true +end + # conversion to CTMRGEnv """ CTMRGEnv(bp_env::BPEnv) diff --git a/src/environments/suweight.jl b/src/environments/suweight.jl index 2966d4155..cd24cb111 100644 --- a/src/environments/suweight.jl +++ b/src/environments/suweight.jl @@ -138,6 +138,15 @@ TensorKit.spacetype(::Type{T}) where {E, T <: SUWeight{E}} = spacetype(E) TensorKit.sectortype(w::SUWeight) = sectortype(typeof(w)) TensorKit.sectortype(::Type{<:SUWeight{T}}) where {T} = sectortype(spacetype(T)) +## Bipartite check +function _is_bipartite(wts::SUWeight) + (size(wts, 2) == size(wts, 3) == 2) || (return false) + for (d, c) in Iterators.product(1:2, 1:2) + (wts[d, 1, c] == wts[d, 2, _next(c, 2)]) || (return false) + end + return true +end + ## (Approximate) equality function Base.:(==)(wts1::SUWeight, wts2::SUWeight) return wts1.data == wts2.data diff --git a/src/operators/infinitepepo.jl b/src/operators/infinitepepo.jl index 172c56b05..0eb44444a 100644 --- a/src/operators/infinitepepo.jl +++ b/src/operators/infinitepepo.jl @@ -202,6 +202,14 @@ function InfinitePEPS(ρ::InfinitePEPO) ) end +## Bipartite check for PEPS/PEPO +function _is_bipartite(psi::InfiniteState) + (size(psi, 1) == size(psi, 2) == 2) || (return false) + for (c, h) in Iterators.product(1:2, 1:size(psi, 3)) + (psi[1, c, h] == psi[2, _next(c, 2), h]) || (return false) + end + return true +end ## Vector interface diff --git a/test/bondenv/benv_ctm.jl b/test/bondenv/benv_ctm.jl index 7676c5139..ee081520f 100644 --- a/test/bondenv/benv_ctm.jl +++ b/test/bondenv/benv_ctm.jl @@ -19,7 +19,7 @@ function get_hubbard_peps(t::Float64 = 1.0, U::Float64 = 8.0) wts = SUWeight(peps) alg = SimpleUpdate(; trunc = truncerror(; atol = 1.0e-10) & truncrank(4)) evolver = TimeEvolver(peps, H, 1.0e-2, 10000, alg, wts) - peps, = time_evolve(evolver; tol = 1.0e-8, check_interval = 2000) + peps, = time_evolve(evolver, H; tol = 1.0e-8, check_interval = 2000) normalize!.(peps.A, Inf) return peps end @@ -42,7 +42,8 @@ function test_benv_ctm(state::Union{InfinitePEPS, InfinitePEPO}) for row in 1:Nr, col in 1:Nc cp1 = PEPSKit._next(col, Nc) A, B = state.A[row, col], state.A[row, cp1] - X, a, b, Y = PEPSKit._qr_bond(A, B) + a, X = PEPSKit.bond_tensor_first(A) + b, Y = PEPSKit.bond_tensor_last(B) 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/bondenv/benv_gaugefix.jl b/test/bondenv/benv_gaugefix.jl index bdd5cecf7..4be536b0d 100644 --- a/test/bondenv/benv_gaugefix.jl +++ b/test/bondenv/benv_gaugefix.jl @@ -30,14 +30,15 @@ for V1 in Vs, V2 in Vs, V3 in Vs ↓ ↓ ↓ -1 -2 -3 =# - a = randn(ComplexF64, V1 ← Vphy' ⊗ V2) + a = randn(ComplexF64, V1 ⊗ Vphy ← V2) b = randn(ComplexF64, V2 ⊗ Vphy ← V3) - @tensor half[:] := Z[-1; 1 3] * a[1; -2 2] * b[2 -3; 3] + @tensor half[:] := Z[-1; 1 3] * a[1 -2; 2] * b[2 -3; 3] Z2, a2, b2, (Linv, Rinv) = PEPSKit.fixgauge_benv(Z, a, b) - @tensor half2[:] := Z2[-1; 1 3] * a2[1; -2 2] * b2[2 -3; 3] + @tensor half2[:] := Z2[-1; 1 3] * a2[1 -2; 2] * b2[2 -3; 3] @test half ≈ half2 # test gauge transformation of X, Y - X2, Y2 = PEPSKit._fixgauge_benvXY(X, Y, Linv, Rinv) + X2 = PEPSKit._fixgauge_benvX(X, Rinv) + Y2 = PEPSKit._fixgauge_benvY(Y, Linv) @tensor Z2_[p; Xe Yw] := Z0[p; Xn Xs Xw Yn Ye Ys] * X2[Xn Xe Xs Xw] * Y2[Yn Ye Ys Yw] @test Z2 ≈ Z2_ end diff --git a/test/bondenv/benv_ntu.jl b/test/bondenv/benv_ntu.jl new file mode 100644 index 000000000..24cb5ab0c --- /dev/null +++ b/test/bondenv/benv_ntu.jl @@ -0,0 +1,74 @@ +using Test +using TensorKit +using PEPSKit +using LinearAlgebra +using KrylovKit +using Random + +Nr, Nc = 2, 2 +Random.seed!(20) + +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 + cp1 = PEPSKit._next(col, Nc) + A, B = state.A[row, col], state.A[row, cp1] + a, X = PEPSKit.bond_tensor_first(A) + b, Y = PEPSKit.bond_tensor_last(B) + @tensor ab[DX DY; da db] := a[DX da D] * b[D db DY] + benv = PEPSKit.bondenv_ntu(row, col, 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, 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 + +@testset "NTU bond environment ($(sym))" 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), + ] + normalize!.(state.A, Inf) + for env_alg in (NNEnv(), NNpEnv(), NNNEnv()) + test_ntu_env(state, env_alg) + end + end +end diff --git a/test/bondenv/bond_truncate.jl b/test/bondenv/bond_truncate.jl index 25f56c200..415d7c0d4 100644 --- a/test/bondenv/bond_truncate.jl +++ b/test/bondenv/bond_truncate.jl @@ -5,40 +5,50 @@ using TensorKit using PEPSKit using LinearAlgebra using KrylovKit -using PEPSKit: cost_function_als +using PEPSKit: bond_truncate, cost_function_als +using PEPSKit: _combine_ket, _combine_ket_for_svd Random.seed!(0) maxiter = 600 -check_interval = 20 -trunc = truncerror(; atol = 1.0e-10) & truncrank(8) -Vext = Vect[FermionParity](0 => 100, 1 => 100) -Vint = Vect[FermionParity](0 => 6, 1 => 6) -Vphy = Vect[FermionParity](0 => 1, 1 => 2) -perm_ab = ((1, 3), (4, 2)) -for Vbondl in (Vint, Vint'), Vbondr in (Vint, Vint') - Vbond = Vbondl ⊗ Vbondr +check_interval = 30 +elt = Float64 +# simulating the situation of applying a 2-site gate +# to a bond with virtual dimension D, physical dimension d. +d, D = 2, 4 +trunc = truncerror(; atol = 1.0e-10) & truncrank(D) +Vphy = Vect[FermionParity](0 => div(d, 2), 1 => div(d, 2)) +Vqro = Vect[FermionParity](0 => div(d * D, 2), 1 => div(d * D, 2)) +# virtual dimension of gate MPO is d^2 +Vint = Vect[FermionParity](0 => div(d^2 * D, 2), 1 => div(d^2 * D, 2)) +for Vl in (Vqro, Vqro'), Vr in (Vqro, Vqro') # random positive-definite environment - Z = randn(Float64, Vext ← Vbond) + Vbond = Vl ⊗ Vr + Dext = dim(Vbond) + Vext = Vect[FermionParity](0 => div(Dext, 2) + 1, 1 => div(Dext, 2) + 1) + Z = randn(elt, Vext ← Vbond) + normalize!(Z, Inf) benv = Z' * Z - normalize!(benv, Inf) - # untruncated bond tensor - a2b2 = randn(Float64, Vbondl ⊗ Vbondr ← Vphy' ⊗ Vphy') - a2, s, b2 = svd_compact(permute(a2b2, perm_ab)) - a2, b2 = PEPSKit.absorb_s(a2, s, b2) + @info "Dimension of benv = $(Dext)" + # untruncated bond tensors + a2 = randn(elt, Vl ⊗ Vphy ← Vint) + b2 = randn(elt, Vint ⊗ Vphy ← Vr') # bond tensor (truncated SVD initialization) - a0, s, b0 = svd_trunc(permute(a2b2, perm_ab); trunc = trunc) + a2b2 = _combine_ket(a2, b2) + a0, s, b0 = svd_trunc(permute(a2b2, ((1, 3), (4, 2))); trunc = trunc) a0, b0 = PEPSKit.absorb_s(a0, s, b0) - fid0 = cost_function_als(benv, PEPSKit._combine_ab(a0, b0), a2b2)[2] + b0 = permute(b0, ((1, 2), (3,))) + fid0 = cost_function_als(benv, _combine_ket(a0, b0), a2b2)[2] @info "Fidelity of simple SVD truncation = $fid0.\n" ss = Dict{String, DiagonalTensorMap}() + # FET is slower when d is large for (label, alg) in ( ("ALS", ALSTruncation(; trunc, maxiter, check_interval)), ("FET", FullEnvTruncation(; trunc, maxiter, check_interval, trunc_init = false)), ) - a1, ss[label], b1, info = PEPSKit.bond_truncate(a2, b2, benv, alg) + a1, ss[label], b1, info = bond_truncate(a2, b2, benv, alg) @info "$label improved fidelity = $(info.fid)." # display(ss[label]) - @test info.fid ≈ cost_function_als(benv, PEPSKit._combine_ab(a1, b1), a2b2)[2] + @test info.fid ≈ cost_function_als(benv, _combine_ket(a1, b1), a2b2)[2] @test info.fid > fid0 end @test isapprox(ss["ALS"], ss["FET"], atol = 1.0e-3) diff --git a/test/bp/gaugefix.jl b/test/bp/gaugefix.jl index 941c8cbea..350a4f205 100644 --- a/test/bp/gaugefix.jl +++ b/test/bp/gaugefix.jl @@ -3,10 +3,13 @@ using Random using TensorKit using PEPSKit using PEPSKit: compare_weights, random_dual!, twistdual +using PEPSKit: _next, _is_bipartite -@testset "Compare BP and SU ($S, posdef msgs = $h)" for (S, h) in - Iterators.product([U1Irrep, FermionParity], [true, false]) - unitcell = (2, 3) +@testset "BP vs SU ($S, bipartite = $(bipartite), posdef msgs = $h)" for + (S, bipartite, h) in Iterators.product( + [U1Irrep, FermionParity], [true, false], [true, false] + ) + unitcell = bipartite ? (2, 2) : (2, 3) elt = ComplexF64 maxiter, tol = 100, 1.0e-9 Random.seed!(52840679) @@ -32,25 +35,48 @@ using PEPSKit: compare_weights, random_dual!, twistdual end end Nspaces, Espaces = random_dual!(Nspaces), random_dual!(Espaces) + if bipartite + for c in 1:2 + cp1 = _next(c, 2) + Pspaces[2, c] = Pspaces[1, cp1] + Nspaces[2, c] = Nspaces[1, cp1] + Espaces[2, c] = Espaces[1, cp1] + end + end peps0 = InfinitePEPS(randn, elt, Pspaces, Nspaces, Espaces) + if bipartite + for c in 1:2 + peps0[2, c] = copy(peps0[1, _next(c, 2)]) + end + end # start by gauging with SU peps1, wts1 = gauge_fix(peps0, SUGauge(; maxiter, tol)) for (a0, a1) in zip(peps0.A, peps1.A) @test space(a0) == space(a1) end + if bipartite + @test _is_bipartite(peps1) + @test _is_bipartite(wts1) + end normalize!.(wts1.data) # find BP fixed point and SUWeight - bp_alg = BeliefPropagation(; maxiter, tol, project_hermitian = h) + bp_alg = BeliefPropagation(; maxiter, tol, bipartite, project_hermitian = h) env = BPEnv(randn, elt, peps1; posdef = h) env, err = leading_boundary(env, peps1, bp_alg) + if bipartite + @test _is_bipartite(env) + end wts2 = SUWeight(env) normalize!.(wts2.data) @test compare_weights(wts1, wts2) < 1.0e-9 bpg_alg = BPGauge() peps2, XXinv = @constinferred gauge_fix(peps1, bpg_alg, env) + if bipartite + @test _is_bipartite(peps2) + end for (a1, a2) in zip(peps1.A, peps2.A) @test space(a1) == space(a2) end diff --git a/test/ctmrg/suweight.jl b/test/ctmrg/suweight.jl index 3e90e28cc..12ff64e09 100644 --- a/test/ctmrg/suweight.jl +++ b/test/ctmrg/suweight.jl @@ -2,7 +2,7 @@ using Test using Random using TensorKit using PEPSKit -using PEPSKit: str, twistdual, _prev, _next +using PEPSKit: str, twistdual, _prev, _next, unitcell Vps = Dict( Z2Irrep => Vect[Z2Irrep](0 => 1, 1 => 2), @@ -43,12 +43,13 @@ end normalize!.(wts.data, Inf) end env = CTMRGEnv(wts) - for (r, c) in Tuple.(CartesianIndices(peps.A)) + for idx in CartesianIndices(unitcell(peps)) + r, c = Tuple(idx) ρ1 = su_rdm_1x1(r, c, peps, wts) if init == :trivial @test ρ1 ≈ su_rdm_1x1(r, c, peps, nothing) end - ρ2 = reduced_densitymatrix(((r, c),), peps, env) + ρ2 = reduced_densitymatrix([idx], peps, env) @test ρ1 ≈ ρ2 end end diff --git a/test/runtests.jl b/test/runtests.jl index 69f670e35..fb3848f85 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -93,6 +93,9 @@ end @time @safetestset "Gauge fixing" begin include("bondenv/benv_gaugefix.jl") end + @time @safetestset "Exact NTU bond environments" begin + include("bondenv/benv_ntu.jl") + end @time @safetestset "Full update bond environment" begin include("bondenv/benv_ctm.jl") end @@ -104,7 +107,7 @@ end @time @safetestset "Cluster truncation with projectors" begin include("timeevol/cluster_projectors.jl") end - @time @safetestset "Time evolution with site-dependent truncation" begin + @time @safetestset "Fixed-space and site-dependent truncation" begin include("timeevol/sitedep_truncation.jl") end @time @safetestset "Transverse field Ising model at finite temperature" begin @@ -113,6 +116,9 @@ end @time @safetestset "J1-J2 model at finite temperature" begin include("timeevol/j1j2_finiteT.jl") end + @time @safetestset "Transverse Field Ising model real-time evolution" begin + include("timeevol/tf_ising_realtime.jl") + end end if GROUP == "ALL" || GROUP == "TOOLBOX" @time @safetestset "Density matrix from double-layer PEPO" begin diff --git a/test/timeevol/cluster_projectors.jl b/test/timeevol/cluster_projectors.jl index 53fad0a68..1276ad2f8 100644 --- a/test/timeevol/cluster_projectors.jl +++ b/test/timeevol/cluster_projectors.jl @@ -35,7 +35,8 @@ Vspaces = [ flips = [isdual(space(M, 1)) for M in Ms1[2:end]] # no truncation Ms2 = _flip_virtuals!(deepcopy(Ms1), flips) - wts2, ϵs, = _cluster_truncate!(Ms2, fill(FixedSpaceTruncation(), N - 1)) + truncs = [truncrank(dim(space(M, 1))) for M in Ms2[2:end]] + wts2, ϵs, = _cluster_truncate!(Ms2, truncs) @test all((ϵ == 0) for ϵ in ϵs) normalize!.(Ms2, Inf) @test fidelity_cluster(Ms1, Ms2) ≈ 1.0 diff --git a/test/timeevol/j1j2_finiteT.jl b/test/timeevol/j1j2_finiteT.jl index 63576e6a1..b834b3522 100644 --- a/test/timeevol/j1j2_finiteT.jl +++ b/test/timeevol/j1j2_finiteT.jl @@ -5,13 +5,12 @@ import MPSKitModels: σˣ, σᶻ using PEPSKit # Benchmark energy from high-temperature expansion -# at β = 0.3, 0.6 -# Physical Review B 86, 045139 (2012) Fig. 15-16 -bm = [-0.1235, -0.213] +const βs = [0.2, 0.4, 0.6] +const bm = [-0.08624893, -0.15688984, -0.21300888] function converge_env(state, χ::Int) trunc1 = truncrank(χ) & truncerror(; atol = 1.0e-12) - env0 = CTMRGEnv(ones, Float64, state, Vect[SU2Irrep](0 => 1)) + env0 = CTMRGEnv(ones, Float64, state, oneunit(spacetype(state))) env, = leading_boundary(env0, state; alg = :sequential, trunc = trunc1, tol = 1.0e-10) return env end @@ -23,37 +22,37 @@ ham = j1_j2_model( ) pepo0 = PEPSKit.infinite_temperature_density_matrix(ham) wts0 = SUWeight(pepo0) -# 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 +dt, nstep, check_interval = 5.0e-3, 40, 40 -# 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) +@testset "Simple update" begin + # 7 = 1 (spin-0) + 2 x 3 (spin-1) + trunc_pepo = truncrank(7) & truncerror(; atol = 1.0e-12) + alg = SimpleUpdate(; trunc = trunc_pepo, purified = true) + 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(4) & 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 = 2.0e-2 + end +end diff --git a/test/timeevol/sitedep_truncation.jl b/test/timeevol/sitedep_truncation.jl index a36f3de75..0efde3c5c 100644 --- a/test/timeevol/sitedep_truncation.jl +++ b/test/timeevol/sitedep_truncation.jl @@ -1,70 +1,73 @@ using Test -using LinearAlgebra using Random using TensorKit using PEPSKit -using PEPSKit: NORTH, EAST, _next +using PEPSKit: _is_bipartite -function get_bonddims(peps::InfinitePEPS) - xdims = collect(dim(domain(t, EAST)) for t in peps.A) - ydims = collect(dim(domain(t, NORTH)) for t in peps.A) - return stack([xdims, ydims]; dims = 1) -end - -function get_bonddims(wts::SUWeight) - xdims = collect(dim(space(wt, 1)) for wt in wts[1, :, :]) - ydims = collect(dim(space(wt, 1)) for wt in wts[2, :, :]) - return stack([xdims, ydims]; dims = 1) -end +elt = Float64 +Nr, Nc = 2, 2 +Vps = fill(U1Space(1 / 2 => 1, -1 / 2 => 1), (Nr, Nc)) +Vns = [ + U1Space(0 => 1, 1 => 2, -1 => 1) U1Space(0 => 1, 1 => 2, -1 => 1)'; + U1Space(0 => 1, 1 => 2, -1 => 1)' U1Space(0 => 1, 1 => 2, -1 => 1) +] +Ves1 = [ + U1Space(1 / 2 => 1, -1 / 2 => 2, -3 / 2 => 1)' U1Space(0 => 1, 1 => 1, -1 => 2); + U1Space(0 => 1, 1 => 1, -1 => 2) U1Space(1 / 2 => 1, -1 / 2 => 2, -3 / 2 => 1)' +] +Ves2 = [ + U1Space(0 => 1, 1 => 2, -1 => 1)' U1Space(0 => 1, 1 => 1, -1 => 2); + U1Space(0 => 1, 1 => 1, -1 => 2) U1Space(0 => 1, 1 => 2, -1 => 1)' +] +Venv = U1Space(0 => 2, 1 => 1, -1 => 1) +states = ( + InfinitePEPS(randn, elt, Vps, Vns, Ves1), + InfinitePEPO(randn, elt, Vps, Vns, Ves2), +) -@testset "Simple update: bipartite 2-site" begin - Nr, Nc = 2, 2 - ham = real(heisenberg_XYZ(InfiniteSquare(Nr, Nc); Jx = 1.0, Jy = 1.0, Jz = 1.0)) - Random.seed!(100) - peps0 = InfinitePEPS(rand, Float64, ℂ^2, ℂ^10; unitcell = (Nr, Nc)) - # make state bipartite - for r in 1:2 - peps0.A[_next(r, 2), 2] = copy(peps0.A[r, 1]) +@testset "Simple update on $(typeof(state0).name.wrapper), bipartite = $(bipartite)" for + (state0, bipartite) in Iterators.product(states, (true, false)) + J2 = 0.5 + if bipartite + state0[2, 1] = copy(state0[1, 2]) + state0[2, 2] = copy(state0[1, 1]) + J2 = 0.0 + end + ham = j1_j2_model(elt, U1Irrep, InfiniteSquare(Nr, Nc); J1 = 1.0, J2, sublattice = false) + # converted internally to SiteDependentTruncation + alg = SimpleUpdate(; trunc = FixedSpaceTruncation(), bipartite) + wts0 = SUWeight(state0) + state, wts, = time_evolve(state0, ham, 0.1, 1, alg, wts0) + for (t, t0) in zip(state.A, state0.A) + @test space(t) == space(t0) end - env0 = SUWeight(peps0) - normalize!.(peps0.A, Inf) - # set trunc to be compatible with bipartite structure - bonddims = stack([[6 4; 4 6], [5 7; 7 5]]; dims = 1) - trunc = SiteDependentTruncation(collect(truncrank(d) for d in bonddims)) - alg = SimpleUpdate(; trunc, bipartite = true) - peps, env, = time_evolve(peps0, ham, 1.0e-2, 4, alg, env0) - @test get_bonddims(peps) == bonddims - @test get_bonddims(env) == bonddims - # check bipartite structure is preserved - for col in 1:2 - cp1 = PEPSKit._next(col, 2) - @test ( - peps.A[1, col] == peps.A[2, cp1] && - env[1, 1, col] == env[1, 2, cp1] && - env[2, 1, col] == env[2, 2, cp1] - ) + for (wt, wt0) in zip(wts.data, wts0.data) + @test space(wt) == space(wt0) + end + if bipartite + @test _is_bipartite(state) + @test _is_bipartite(wts) end end -@testset "Simple update: generic 2-site and 3-site" begin - Nr, Nc = 3, 4 - Random.seed!(100) - peps0 = InfinitePEPS(rand, Float64, ℂ^2, ℂ^10; unitcell = (Nr, Nc)) - normalize!.(peps0.A, Inf) - env0 = SUWeight(peps0) - # Site dependent truncation - bonddims = rand(2:8, 2, Nr, Nc) - @show bonddims - trunc = SiteDependentTruncation(collect(truncrank(d) for d in bonddims)) - alg = SimpleUpdate(; trunc) - # 2-site SU - ham = j1_j2_model(Float64, Trivial, InfiniteSquare(Nr, Nc); J2 = 0.0, sublattice = false) - peps, env, = time_evolve(peps0, ham, 1.0e-2, 4, alg, env0) - @test get_bonddims(peps) == bonddims - @test get_bonddims(env) == bonddims - # 3-site SU - ham = j1_j2_model(Float64, Trivial, InfiniteSquare(Nr, Nc); J2 = 0.2, sublattice = false) - peps, env, = time_evolve(peps0, ham, 1.0e-2, 4, alg, env0) - @test get_bonddims(peps) == bonddims - @test get_bonddims(env) == bonddims +@testset "NTU on $(typeof(state0).name.wrapper), bipartite = $(bipartite)" for + (state0, bipartite) in Iterators.product(states, (true, false)) + J2 = 0.5 + if bipartite + state0[2, 1] = copy(state0[1, 2]) + state0[2, 2] = copy(state0[1, 1]) + J2 = 0.0 + end + ham = 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 + if bipartite + @test _is_bipartite(state) + @test _is_bipartite(info.wts) + end end diff --git a/test/timeevol/tf_ising_realtime.jl b/test/timeevol/tf_ising_realtime.jl new file mode 100644 index 000000000..b34318dcf --- /dev/null +++ b/test/timeevol/tf_ising_realtime.jl @@ -0,0 +1,90 @@ +using Test +using TensorKit +import MPSKitModels: S_zz, σˣ +using PEPSKit +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 +] + +# 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 +lattice = collect(space(t, 1) for t in peps0.A) + +# Hamiltonian +op = LocalOperator(lattice, ((1, 1),) => σˣ()) +ham = 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 = :fullinfinite +) + +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 = 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 diff --git a/test/timeevol/timestep.jl b/test/timeevol/timestep.jl index 39153a898..ce86a4569 100644 --- a/test/timeevol/timestep.jl +++ b/test/timeevol/timestep.jl @@ -3,14 +3,16 @@ using Random using TensorKit using PEPSKit +Nr, Nc = 2, 2 +H = real(heisenberg_XYZ(ComplexF64, Trivial, InfiniteSquare(Nr, Nc); Jx = 1, Jy = 1, Jz = 1)) +Pspace, Vspace = ℂ^2, ℂ^4 +ψ0 = InfinitePEPS(rand, Float64, Pspace, Vspace; unitcell = (Nr, Nc)) +dt, nstep = 0.1, 20 +trunc = truncerror(; atol = 1.0e-10) & truncrank(4) + @testset "SimpleUpdate timestep" begin - Nr, Nc = 2, 2 - H = real(heisenberg_XYZ(ComplexF64, Trivial, InfiniteSquare(Nr, Nc); Jx = 1, Jy = 1, Jz = 1)) - Pspace, Vspace = ℂ^2, ℂ^4 - ψ0 = InfinitePEPS(rand, Float64, Pspace, Vspace; unitcell = (Nr, Nc)) + alg = SimpleUpdate(; trunc) env0 = SUWeight(ψ0) - alg = SimpleUpdate(; trunc = truncerror(; atol = 1.0e-10) & truncrank(4)) - dt, nstep = 1.0e-2, 50 # manual timestep evolver = TimeEvolver(ψ0, H, dt, nstep, alg, env0) ψ1, env1, info1 = deepcopy(ψ0), deepcopy(env0), nothing @@ -31,3 +33,26 @@ using PEPSKit @test env1 == env2 == env3 @test info1 == info2 == info3 end + +@testset "NeighbourUpdate timestep" 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