diff --git a/docs/src/examples/hubbard_su/index.md b/docs/src/examples/hubbard_su/index.md index 6b347c15d..5af3f5045 100644 --- a/docs/src/examples/hubbard_su/index.md +++ b/docs/src/examples/hubbard_su/index.md @@ -130,7 +130,7 @@ is achieved by using `trunc=truncrank(χ)` with different `χ`s in the CTMRG run χenv₀, χenv = 6, 16 env_space = Vect[fℤ₂](0 => χenv₀ / 2, 1 => χenv₀ / 2) normalize!.(peps.A, Inf) -env = CTMRGEnv(wts, peps) +env = CTMRGEnv(wts) for χ in [χenv₀, χenv] global env, = leading_boundary( env, peps; alg = :sequential, tol = 1.0e-8, maxiter = 50, trunc = truncrank(χ) diff --git a/docs/src/examples/hubbard_su/main.ipynb b/docs/src/examples/hubbard_su/main.ipynb index a15f8199b..4dbda8971 100644 --- a/docs/src/examples/hubbard_su/main.ipynb +++ b/docs/src/examples/hubbard_su/main.ipynb @@ -140,7 +140,7 @@ "χenv₀, χenv = 6, 16\n", "env_space = Vect[fℤ₂](0 => χenv₀ / 2, 1 => χenv₀ / 2)\n", "normalize!.(peps.A, Inf)\n", - "env = CTMRGEnv(wts, peps)\n", + "env = CTMRGEnv(wts)\n", "for χ in [χenv₀, χenv]\n", " global env, = leading_boundary(\n", " env, peps; alg = :sequential, tol = 1.0e-8, maxiter = 50, trunc = truncrank(χ)\n", diff --git a/examples/hubbard_su/main.jl b/examples/hubbard_su/main.jl index 786dab680..2bb0f411c 100644 --- a/examples/hubbard_su/main.jl +++ b/examples/hubbard_su/main.jl @@ -83,7 +83,7 @@ is achieved by using `trunc=truncrank(χ)` with different `χ`s in the CTMRG run χenv₀, χenv = 6, 16 env_space = Vect[fℤ₂](0 => χenv₀ / 2, 1 => χenv₀ / 2) normalize!.(peps.A, Inf) -env = CTMRGEnv(wts, peps) +env = CTMRGEnv(wts) for χ in [χenv₀, χenv] global env, = leading_boundary( env, peps; alg = :sequential, tol = 1.0e-8, maxiter = 50, trunc = truncrank(χ) diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index 4a39d6bd3..57d359cd3 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -1,6 +1,7 @@ module PEPSKit using LinearAlgebra, Statistics, Base.Threads, Base.Iterators, Printf +using Random using Compat using Accessors: @set, @reset using VectorInterface diff --git a/src/algorithms/time_evolution/evoltools.jl b/src/algorithms/time_evolution/evoltools.jl index 40ceee032..146ede367 100644 --- a/src/algorithms/time_evolution/evoltools.jl +++ b/src/algorithms/time_evolution/evoltools.jl @@ -233,7 +233,7 @@ function _apply_gate( a, b = absorb_s(a, s, b) if need_flip - a, s, b = flip_svd(a, s, b) + a, s, b = flip(a, numind(a)), _fliptwist_s(s), flip(b, 1) end return a, s, b, ϵ end diff --git a/src/algorithms/time_evolution/simpleupdate3site.jl b/src/algorithms/time_evolution/simpleupdate3site.jl index d72d61b8a..895df8dca 100644 --- a/src/algorithms/time_evolution/simpleupdate3site.jl +++ b/src/algorithms/time_evolution/simpleupdate3site.jl @@ -479,7 +479,7 @@ function _su3site_se!( Ms = [first(_unfuse_physicalspace(M, Vphy)) for (M, Vphy) in zip(Ms, Vphys)] end for (wt, wt_idx, flip) in zip(wts, wt_idxs, flips) - env[CartesianIndex(wt_idx)] = normalize(flip ? _flip_s(wt) : wt, Inf) + env[CartesianIndex(wt_idx)] = normalize(flip ? _fliptwist_s(wt) : wt, Inf) end end # restore virtual arrows in `Ms` diff --git a/src/algorithms/truncation/bond_truncation.jl b/src/algorithms/truncation/bond_truncation.jl index d4d1fde07..9d70947b1 100644 --- a/src/algorithms/truncation/bond_truncation.jl +++ b/src/algorithms/truncation/bond_truncation.jl @@ -139,7 +139,7 @@ function bond_truncate( a, s, b = svd_trunc!(permute(_combine_ab(a, b), perm_ab); trunc = alg.trunc) a, b = absorb_s(a, s, b) if need_flip - a, s, b = flip_svd(a, s, b) + a, s, b = flip(a, numind(a)), _fliptwist_s(s), flip(b, 1) end return a, s, b, (; fid, Δfid, Δs) end @@ -186,7 +186,7 @@ function bond_truncate( @tensor a[-1 -2; -3] := Qa[-1 -2 3] * u[3 -3] @tensor b[-1; -2 -3] := vh[-1 1] * Qb[1 -2 -3] if need_flip - a, s, b = flip_svd(a, s, b) + a, s, b = flip(a, numind(a)), _fliptwist_s(s), flip(b, 1) end return a, s, b, info end diff --git a/src/environments/suweight.jl b/src/environments/suweight.jl index 75e3180bc..ea4491aec 100644 --- a/src/environments/suweight.jl +++ b/src/environments/suweight.jl @@ -8,32 +8,23 @@ const PEPSWeight{T, S} = AbstractTensorMap{T, S, 1, 1} """ struct SUWeight{E<:PEPSWeight} -Schmidt bond weights used in simple/cluster update. -Weight elements are always real and non-negative. -The domain and codomain of each weight matrix must be an un-dualed `ElementarySpace`. +Schmidt bond weights used in simple/cluster update. +Each weight is a real and semi-positive definite +`DiagonalTensorMap`, with the same codomain and domain. -For a square lattice InfinitePEPS, the weights are placed as -``` - | - -T[r-1,c]- - | - wt[2,r,c] - | | - --T[r,c]--wt[1,r,c]--T[r,c+1]-- - | | -``` +On the square lattice, +- `wt[1,r,c]` is on the x-bond between `[r,c]` and `[r,c+1]`; +- `wt[2,r,c]` is on the y-bond between `[r,c]` and `[r-1,c]`. Axis order of each weight matrix is ``` - x-weights: - 1 ← x ← 2 or 2 → x → 1 - - y-weights: - 2 1 - ↓ ↑ - y or y - ↓ ↑ - 1 2 + x-weights: y-weights: + + 1 - x - 2 2 + | + y + | + 1 ``` ## Fields @@ -54,37 +45,26 @@ function SUWeight(data::Array{E, 3}) where {E <: PEPSWeight} for wt in data isa(wt, DiagonalTensorMap) || error("Each weight matrix should be a DiagonalTensorMap") + # in case TensorKit drops this requirement domain(wt, 1) == codomain(wt, 1) || error("Domain and codomain of each weight matrix must be the same.") - !isdual(codomain(wt, 1)) || - error("Domain and codomain of each weight matrix cannot be a dual space.") all(wt.data .>= 0) || error("Weight elements must be non-negative.") end return SUWeight{E}(data) end function SUWeight(wts_mats::AbstractMatrix{E}...) where {E <: PEPSWeight} - n_mat = length(wts_mats) - Nr, Nc = size(wts_mats[1]) - @assert all((Nr, Nc) == size(wts_mat) for wts_mat in wts_mats) - weights = collect( - wts_mats[d][r, c] for (d, r, c) in Iterators.product(1:n_mat, 1:Nr, 1:Nc) - ) - return SUWeight(weights) + return SUWeight(stack(wts_mats; dims = 1)) end """ - SUWeight(Nspaces::M, [Espaces::M]) where {M<:AbstractMatrix{<:Union{Int,ElementarySpace}}} + SUWeight(Nspaces::M, [Espaces::M]) where {M<:AbstractMatrix{<:ElementarySpace}} -Create an SUWeight by specifying the vertical (north) or horizontal (east) virtual bond spaces. -Each individual space can be specified as either an `Int` or an `ElementarySpace`. -The weights are initialized as identity matrices of element type `Float64`. +Create a trivial `SUWeight` by specifying the vertical (north) or horizontal (east) virtual bond spaces. """ function SUWeight( Nspaces::M, Espaces::M = Nspaces - ) where {M <: AbstractMatrix{<:Union{Int, ElementarySpace}}} - @assert all(!isdual, Nspaces) - @assert all(!isdual, Espaces) + ) where {M <: AbstractMatrix{<:ElementarySpace}} @assert size(Nspaces) == size(Espaces) Nr, Nc = size(Nspaces) weights = map(Iterators.product(1:2, 1:Nr, 1:Nc)) do (d, r, c) @@ -97,9 +77,8 @@ end """ SUWeight(Nspace::S, Espace::S=Nspace; unitcell::Tuple{Int,Int}=(1, 1)) where {S<:ElementarySpace} -Create an SUWeight by specifying its vertical (north) and horizontal (east) +Create a trivial `SUWeight` by specifying its vertical (north) and horizontal (east) as `ElementarySpace`s) and unit cell size. -The weights are initialized as identity matrices of element type `Float64`. """ function SUWeight( Nspace::S, Espace::S = Nspace; unitcell::Tuple{Int, Int} = (1, 1) @@ -110,40 +89,36 @@ end """ SUWeight(peps::InfinitePEPS) -Create an SUWeight for a given InfinitePEPS. -The weights are initialized as identity matrices of element type `Float64`. +Create a trivial `SUWeight` for a given InfinitePEPS. """ function SUWeight(peps::InfinitePEPS) - Nspaces = map(peps.A) do t - V = domain(t, NORTH) - isdual(V) ? V' : V - end - Espaces = map(peps.A) do t - V = domain(t, EAST) - isdual(V) ? V' : V - end + Nspaces = map(Base.Fix2(domain, NORTH), unitcell(peps)) + Espaces = map(Base.Fix2(domain, EAST), unitcell(peps)) return SUWeight(Nspaces, Espaces) end """ SUWeight(pepo::InfinitePEPO) -Create an SUWeight for a given one-layer InfinitePEPO. -The weights are initialized as identity matrices of element type `Float64`. +Create a trivial `SUWeight` for a given one-layer InfinitePEPO. """ function SUWeight(pepo::InfinitePEPO) @assert size(pepo, 3) == 1 - Nspaces = map(@view(pepo.A[:, :, 1])) do t - V = north_virtualspace(t) - isdual(V) ? V' : V - end - Espaces = map(@view(pepo.A[:, :, 1])) do t - V = east_virtualspace(t) - isdual(V) ? V' : V - end + Nspaces = map(Base.Fix2(domain, NORTH), @view(unitcell(pepo)[:, :, 1])) + Espaces = map(Base.Fix2(domain, EAST), @view(unitcell(pepo)[:, :, 1])) return SUWeight(Nspaces, Espaces) end +Random.rand!(wts::SUWeight) = rand!(Random.default_rng(), wts) +function Random.rand!(rng::Random.AbstractRNG, wts::SUWeight) + foreach(wts.data) do wt + for (_, b) in blocks(wt) + sort!(rand!(rng, b.diag); rev = true) + end + end + return wts +end + ## Shape and size Base.size(W::SUWeight) = size(W.data) Base.size(W::SUWeight, i) = size(W.data, i) @@ -182,7 +157,8 @@ end function Base.show(io::IO, ::MIME"text/plain", wts::SUWeight) println(io, typeof(wts)) for idx in CartesianIndices(wts.data) - println(io, Tuple(idx), ":") + print(io, Tuple(idx), ": ") + println(space(wts.data[idx])) for (k, b) in blocks(wts.data[idx]) println(io, k, " = ", diag(b)) end @@ -194,9 +170,9 @@ end absorb_weight(t::Union{PEPSTensor, PEPOTensor}, weights::SUWeight, row::Int, col::Int, ax::Int; inv::Bool = false) absorb_weight(t::Union{PEPSTensor, PEPOTensor}, weights::SUWeight, row::Int, col::Int, ax::NTuple{N, Int}; inv::Bool = false) -Absorb or remove environment weight on an axis of tensor `t` known to be located at -position (`row`, `col`) in the unit cell of an InfinitePEPS or InfinitePEPO. -Weights around the tensor at `(row, col)` are +Absorb or remove (in a twist-free way) the square root of environment weight +on an axis of the PEPS/PEPO tensor `t` known to be at position (`row`, `col`) +in the unit cell of an InfinitePEPS/InfinitePEPO. The involved weights are ``` | [2,r,c] @@ -209,7 +185,7 @@ Weights around the tensor at `(row, col)` are ## Arguments -- `t::PT` : PEPSTensor or PEPOTensor to which the weight will be absorbed. +- `t::Union{PEPSTensor, PEPOTensor}` : PEPSTensor or PEPOTensor to which the weight will be absorbed. - `weights::SUWeight` : All simple update weights. - `row::Int` : The row index specifying the position in the tensor network. - `col::Int` : The column index specifying the position in the tensor network. @@ -252,11 +228,9 @@ function absorb_weight( ) t_idx = [(n - nout == ax) ? 1 : -n for n in 1:ntol] ax′ = ax + nout - wt_idx = if isdual(space(t, ax′)) - [1, -ax′] # t ← wt - else - [-ax′, 1] # t → wt - end + wt_idx = (ax == NORTH || ax == EAST) ? [1, -ax′] : [-ax′, 1] + # make absorption/removal twist-free + twistdual!(wt, 1) return permute(ncon((t, wt), (t_idx, wt_idx)), (Tuple(1:nout), Tuple((nout + 1):ntol))) end function absorb_weight( @@ -278,13 +252,13 @@ end y₁₁ y₁₂ y₁₃ | | | ..x₁₃...┼---x₁₁---┼---x₁₂---┼---x₁₃--- - | | | + | | | 2 y₂₁ y₂₂ y₂₃ | | | | y ..x₂₃...┼---x₂₁---┼---x₂₂---┼---x₂₃--- | - | | | - y₃₁ y₃₂ y₃₃ -- x -- - | | | + | | | 1 + y₃₁ y₃₂ y₃₃ + | | | 1 -- x -- 2 ..x₃₃...┼---x₃₁---┼---x₃₂---┼---x₃₃--- : : : y₁₁ y₁₂ y₁₃ @@ -297,19 +271,20 @@ end x₁₃ x₂₃ x₃₃ | | | --y₁₃---┼---y₂₃---┼---y₃₃---┼...y₁₃... - | | | + | | | 2 x₁₂ x₂₂ x₃₂ | | | | x --y₁₂---┼---y₂₂---┼---y₃₂---┼...y₁₂... | - | | | - x₁₁ x₂₁ x₃₁ -- y -- - | | | + | | | 1 + x₁₁ x₂₁ x₃₁ + | | | 2 -- y -- 1 --y₁₁---┼---y₂₁---┼---y₃₁---┼...y₁₁... : : : x₁₃ x₂₃ x₃₃ : : : ``` - x/y-weights are exchanged. + - need to further transpose x-weights. - need to further move 1st column of x-weights to the last column. - `rotr90`: @@ -319,18 +294,19 @@ end : : : ..y₁₁...┼---y₃₁---┼---y₂₁---┼---y₁₁--- | | | - x₃₁ x₂₁ x₁₁ -- y -- - | | | - ..y₁₂...┼---y₃₂---┼---y₂₂---┼---y₁₂--- | - | | | x - x₃₂ x₂₂ x₁₂ | + x₃₁ x₂₁ x₁₁ 1 -- y -- 2 | | | - ..y₁₃...┼---y₃₃---┼---y₂₃---┼---y₁₃--- + ..y₁₂...┼---y₃₂---┼---y₂₂---┼---y₁₂--- 1 + | | | | + x₃₂ x₂₂ x₁₂ x + | | | | + ..y₁₃...┼---y₃₃---┼---y₂₃---┼---y₁₃--- 2 | | | x₃₃ x₂₃ x₁₃ | | | ``` - x/y-weights are exchanged. + - need to further transpose y-weights. - need to further move last row of y-weights to the 1st row. - `rot180`: @@ -340,40 +316,87 @@ end : : : --x₃₃---┼---x₃₂---┼---x₃₁---┼...x₃₃... | | | - y₃₃ y₃₂ y₃₁ -- x -- + y₃₃ y₃₂ y₃₁ 2 -- x -- 1 | | | - --x₂₃---┼---x₂₂---┼---x₂₁---┼...x₂₃... | - | | | y - y₂₃ y₂₂ y₂₁ | - | | | - --x₁₃---┼---x₁₂---┼---x₁₁---┼...x₁₃... + --x₂₃---┼---x₂₂---┼---x₂₁---┼...x₂₃... 1 + | | | | + y₂₃ y₂₂ y₂₁ y + | | | | + --x₁₃---┼---x₁₂---┼---x₁₁---┼...x₁₃... 2 | | | y₁₃ y₁₂ y₁₁ | | | ``` + - need to transpose all weights. - need to move 1st column of x-weights to the last column. - need to move last row of y-weights to the 1st row. =# +function _rotl90_wts_x(wts_x::AbstractMatrix{<:PEPSWeight}) + wts_y = rotl90(wts_x) + return wts_y +end +function _rotr90_wts_x(wts_x::AbstractMatrix{<:PEPSWeight}) + wts_y = circshift(rotr90(wts_x), (1, 0)) + for (i, wt) in enumerate(wts_y) + wts_y[i] = DiagonalTensorMap(transpose(wt; copy = true)) + end + return wts_y +end +function _rot180_wts_x(wts_x::AbstractMatrix{<:PEPSWeight}) + wts_x_ = circshift(rot180(wts_x), (0, -1)) + for (i, wt) in enumerate(wts_x_) + wts_x_[i] = DiagonalTensorMap(transpose(wt; copy = true)) + end + return wts_x_ +end + +function _rotl90_wts_y(wts_y::AbstractMatrix{<:PEPSWeight}) + wts_x = circshift(rotl90(wts_y), (0, -1)) + for (i, wt) in enumerate(wts_x) + wts_x[i] = DiagonalTensorMap(transpose(wt; copy = true)) + end + return wts_x +end +function _rotr90_wts_y(wts_y::AbstractMatrix{<:PEPSWeight}) + wts_x = rotr90(wts_y) + return wts_x +end +function _rot180_wts_y(wts_y::AbstractMatrix{<:PEPSWeight}) + wts_y_ = circshift(rot180(wts_y), (1, 0)) + for (i, wt) in enumerate(wts_y_) + wts_y_[i] = DiagonalTensorMap(transpose(wt; copy = true)) + end + return wts_y_ +end + function Base.rotl90(wts::SUWeight) - wts_x = circshift(rotl90(wts[2, :, :]), (0, -1)) - wts_y = rotl90(wts[1, :, :]) + wts_y = _rotl90_wts_x(wts[1, :, :]) + wts_x = _rotl90_wts_y(wts[2, :, :]) return SUWeight(wts_x, wts_y) end function Base.rotr90(wts::SUWeight) - wts_x = rotr90(wts[2, :, :]) - wts_y = circshift(rotr90(wts[1, :, :]), (1, 0)) + wts_y = _rotr90_wts_x(wts[1, :, :]) + wts_x = _rotr90_wts_y(wts[2, :, :]) return SUWeight(wts_x, wts_y) end function Base.rot180(wts::SUWeight) - wts_x = circshift(rot180(wts[1, :, :]), (0, -1)) - wts_y = circshift(rot180(wts[2, :, :]), (1, 0)) + wts_x = _rot180_wts_x(wts[1, :, :]) + wts_y = _rot180_wts_y(wts[2, :, :]) return SUWeight(wts_x, wts_y) end -function _CTMRGEnv(wts::SUWeight, flips::Array{Bool, 3}) - @assert size(wts) == size(flips) +""" + CTMRGEnv(wts::SUWeight) + +Construct a CTMRG environment with a trivial environment space +(bond dimension χ = 1) from SUWeight `wts`, +which has the same real scalartype as ``wts`. +""" +function CTMRGEnv(wts::SUWeight) _, Nr, Nc = size(wts) + elt = scalartype(wts) + V_env = oneunit(spacetype(wts)) edges = map(Iterators.product(1:4, 1:Nr, 1:Nc)) do (d, r, c) wt_idx = if d == NORTH CartesianIndex(2, _next(r, Nr), c) @@ -384,35 +407,16 @@ function _CTMRGEnv(wts::SUWeight, flips::Array{Bool, 3}) else # WEST CartesianIndex(1, r, c) end - wt = deepcopy(wts[wt_idx]) - if (flips[wt_idx] && d in (NORTH, EAST)) || (!flips[wt_idx] && d in (SOUTH, WEST)) - wt = permute(wt, ((2,), (1,))) + wt = if d in (NORTH, EAST) + twist!(repartition(wts[wt_idx], 2, 0; copy = true), 1) + else + permute(wts[wt_idx], ((2, 1), ()); copy = true) end - twistdual!(wt, 1) - wt = insertrightunit(insertleftunit(wt, 1)) - return permute(wt, ((1, 2, 3), (4,))) + # attach identity on environment space + return insertleftunit(insertleftunit(wt), 1) end corners = map(CartesianIndices(edges)) do idx - return TensorKit.id(Float64, codomain(edges[idx], 1)) + return TensorKit.id(elt, V_env) end return CTMRGEnv(corners, edges) end - -""" - CTMRGEnv(wts::SUWeight, peps::InfinitePEPS) - -Construct a CTMRG environment with bond dimension χ = 1 -for an InfinitePEPS `peps` from the accompanying SUWeight `wts`. -The scalartype of the returned environment is always `Float64`. -""" -function CTMRGEnv(wts::SUWeight, peps::InfinitePEPS) - _, Nr, Nc = size(wts) - flips = map(Iterators.product(1:2, 1:Nr, 1:Nc)) do (d, r, c) - return if d == 1 - isdual(domain(peps[r, c], EAST)) - else - isdual(domain(peps[r, c], NORTH)) - end - end - return _CTMRGEnv(wts, flips) -end diff --git a/src/utility/util.jl b/src/utility/util.jl index f5c5c31f3..a63a2a97f 100644 --- a/src/utility/util.jl +++ b/src/utility/util.jl @@ -100,23 +100,7 @@ function absorb_s(U::AbstractTensorMap, S::DiagonalTensorMap, V::AbstractTensorM return U * sqrt_S, sqrt_S * V end -""" - flip_svd(u::AbstractTensorMap, s::DiagonalTensorMap, vh::AbstractTensorMap) - -Given SVD result `u ← s ← vh`, flip the arrow between the three tensors -to `u2 → s2 → vh2` such that -``` - u * s * vh = (@tensor t2[-1, ...; -2, ...] := u2[-1, ...; 2] * s2[1; 2] * vh2[1; -2, ...]) -``` -The axis orders for `s`, `s2` are -``` - 1 ← s ← 2 2 → s2 → 1 -``` -""" -function flip_svd(u::AbstractTensorMap, s::DiagonalTensorMap, vh::AbstractTensorMap) - return flip(u, numind(u)), _flip_s(s), flip(vh, 1) -end -_flip_s(s::DiagonalTensorMap) = permute(DiagonalTensorMap(flip(s, (1, 2))), ((2,), (1,))) +_fliptwist_s(s::DiagonalTensorMap) = twist!(DiagonalTensorMap(flip(s, 1:2)), 1) """ twistdual(t::AbstractTensorMap, i) diff --git a/test/ctmrg/suweight.jl b/test/ctmrg/suweight.jl index 1945bd3b2..3e90e28cc 100644 --- a/test/ctmrg/suweight.jl +++ b/test/ctmrg/suweight.jl @@ -2,38 +2,52 @@ using Test using Random using TensorKit using PEPSKit -using PEPSKit: str +using PEPSKit: str, twistdual, _prev, _next -function rand_wts(peps::InfinitePEPS) - wts = SUWeight(peps) - for idx in CartesianIndices(wts.data) - n = length(wts.data[idx].data) - wts.data[idx].data[:] = sort(rand(Float64, n); lt = !isless) - end - return wts -end +Vps = Dict( + Z2Irrep => Vect[Z2Irrep](0 => 1, 1 => 2), + U1Irrep => Vect[U1Irrep](0 => 2, 1 => 2, -1 => 1), + FermionParity => Vect[FermionParity](0 => 1, 1 => 2), +) +Vvs = Dict( + Z2Irrep => Vect[Z2Irrep](0 => 2, 1 => 2), + U1Irrep => Vect[U1Irrep](0 => 3, 1 => 1, -1 => 2), + FermionParity => Vect[FermionParity](0 => 2, 1 => 2), +) -function su_rdm_1x1(row::Int, col::Int, peps::InfinitePEPS, wts::SUWeight) +function su_rdm_1x1( + row::Int, col::Int, peps::InfinitePEPS, wts::Union{Nothing, SUWeight} = nothing + ) Nr, Nc = size(peps) @assert 1 <= row <= Nr && 1 <= col <= Nc t = peps.A[row, col] - t = absorb_weight(t, wts, row, col, Tuple(1:4)) - ρ = t * t' + if !(wts === nothing) + t = absorb_weight(t, wts, row, col, Tuple(1:4)) + end + # contract local ⟨t|t⟩ without virtual twists + @tensor ρ[k; b] := conj(t[b; n e s w]) * twistdual(t, 2:5)[k; n e s w] return ρ / str(ρ) end -Vps = [Vect[U1Irrep](0 => 2, 1 => 2, -1 => 1), Vect[FermionParity](0 => 1, 1 => 2)] -Vvs = [Vect[U1Irrep](0 => 3, 1 => 1, -1 => 2), Vect[FermionParity](0 => 2, 1 => 2)] +@testset "SUWeight ($(init) init, $(sect))" for (init, sect) in + Iterators.product([:trivial, :random], keys(Vps)) -for (Vp, Vv) in zip(Vps, Vvs) + Vp, Vv = Vps[sect], Vvs[sect] Nspaces = [Vv Vv' Vv; Vv' Vv Vv'] Espaces = [Vv Vv Vv'; Vv Vv' Vv'] Pspaces = fill(Vp, size(Nspaces)) peps = InfinitePEPS(randn, ComplexF64, Pspaces, Nspaces, Espaces) - wts = rand_wts(peps) - env = CTMRGEnv(wts, peps) + wts = SUWeight(peps) + if init != :trivial + rand!(wts) + normalize!.(wts.data, Inf) + end + env = CTMRGEnv(wts) for (r, c) in Tuple.(CartesianIndices(peps.A)) ρ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) @test ρ1 ≈ ρ2 end diff --git a/test/timeevol/cluster_projectors.jl b/test/timeevol/cluster_projectors.jl index f232c5df4..035ef3826 100644 --- a/test/timeevol/cluster_projectors.jl +++ b/test/timeevol/cluster_projectors.jl @@ -111,7 +111,7 @@ end Espace = Vect[FermionParity ⊠ U1Irrep]((0, 0) => 8, (1, 1 // 2) => 4, (1, -1 // 2) => 4) trunc_env0 = truncerror(; atol = 1.0e-12) & truncrank(4) trunc_env = truncerror(; atol = 1.0e-12) & truncrank(16) - peps = InfinitePEPS(rand, Float64, Pspace, Vspace; unitcell = (Nr, Nc)) + peps = InfinitePEPS(rand, Float64, Pspace, Vspace, Vspace'; unitcell = (Nr, Nc)) wts = SUWeight(peps) ham = real( hubbard_model( @@ -127,7 +127,7 @@ end peps, wts, = time_evolve(peps, ham, dt, 10000, alg, wts; tol, check_interval = 1000) end normalize!.(peps.A, Inf) - env = CTMRGEnv(wts, peps) + env = CTMRGEnv(wts) env, = leading_boundary(env, peps; tol = ctmrg_tol, trunc = trunc_env0) env, = leading_boundary(env, peps; tol = ctmrg_tol, trunc = trunc_env) e_site = cost_function(peps, env, ham) / (Nr * Nc) diff --git a/test/types/suweight.jl b/test/types/suweight.jl index e5b142485..cc7440920 100644 --- a/test/types/suweight.jl +++ b/test/types/suweight.jl @@ -18,43 +18,20 @@ function test_rotation(wts::SUWeight) return nothing end -function test_rotation(peps::InfinitePEPS, wts::SUWeight) +function test_rotation(state::Union{InfinitePEPS, InfinitePEPO}, wts::SUWeight) + InfiniteState = (state isa InfinitePEPS) ? InfinitePEPS : InfinitePEPO + A1 = InfiniteState( + map(CartesianIndices(state.A)) do idx + return absorb_weight(state.A[idx], wts, idx[1], idx[2], Tuple(1:4)) + end + ) for n in 1:4 rot = compose_n(rotl90, n) - A1 = InfinitePEPS( - collect( - absorb_weight(peps.A[idx], wts, idx[1], idx[2], Tuple(1:4)) for - idx in CartesianIndices(peps.A) - ), - ) - peps2, wts2 = rot(peps), rot(wts) - A2 = InfinitePEPS( - collect( - absorb_weight(peps2.A[idx], wts2, idx[1], idx[2], Tuple(1:4)) for - idx in CartesianIndices(peps2.A) - ), - ) - @test A2 ≈ rot(A1) - end - return -end - -function test_rotation(pepo::InfinitePEPO, wts::SUWeight) - @assert size(pepo, 3) == 1 - for n in 1:4 - rot = compose_n(rotl90, n) - A1 = InfinitePEPO( - collect( - absorb_weight(pepo.A[idx], wts, idx[1], idx[2], Tuple(1:4)) for - idx in CartesianIndices(pepo.A) - ), - ) - pepo2, wts2 = rot(pepo), rot(wts) - A2 = InfinitePEPO( - collect( - absorb_weight(pepo2.A[idx], wts2, idx[1], idx[2], Tuple(1:4)) for - idx in CartesianIndices(pepo2.A) - ), + state2, wts2 = rot(state), rot(wts) + A2 = InfiniteState( + map(CartesianIndices(state2.A)) do idx + return absorb_weight(state2.A[idx], wts2, idx[1], idx[2], Tuple(1:4)) + end ) @test A2 ≈ rot(A1) end @@ -71,14 +48,19 @@ Vs = ( Nr, Nc = 2, 3 peps = InfinitePEPS(rand, Float64, Vphy, Vs[2], Vs[1]'; unitcell = (Nr, Nc)) pepo = InfinitePEPO(rand, Float64, Vphy, Vs[2], Vs[1]'; unitcell = (Nr, Nc, 1)) -wts = collect( - svd_compact(rand(Float64, Vs[dir] ← Vs[dir]))[2] for dir in 1:2, r in 1:Nr, c in 1:Nc -) -wts = SUWeight(wts) +wts = SUWeight(peps) +rand!(wts) +normalize!.(wts.data, Inf) +# check that elements of wts are successfully randomized +@test !(wts ≈ SUWeight(peps)) @test sectortype(wts) === sectortype(Vs[1]) @test spacetype(wts) === spacetype(Vs[1]) test_rotation(wts) test_rotation(peps, wts) + +wts = SUWeight(pepo) +rand!(wts) +normalize!.(wts.data, Inf) test_rotation(pepo, wts)