diff --git a/.github/workflows/Tests.yml b/.github/workflows/Tests.yml index 621e19c4b..c88bccd70 100644 --- a/.github/workflows/Tests.yml +++ b/.github/workflows/Tests.yml @@ -26,6 +26,7 @@ jobs: - 'lts' # minimal supported version - '1' group: + - bp - types - ctmrg - boundarymps diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index 57d359cd3..8fd007a71 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -51,11 +51,13 @@ include("operators/models.jl") include("environments/ctmrg_environments.jl") include("environments/vumps_environments.jl") include("environments/suweight.jl") +include("environments/bp_environments.jl") include("algorithms/contractions/ctmrg_contractions.jl") include("algorithms/contractions/transfer.jl") include("algorithms/contractions/localoperator.jl") include("algorithms/contractions/vumps_contractions.jl") +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") @@ -78,6 +80,10 @@ include("algorithms/time_evolution/evoltools.jl") include("algorithms/time_evolution/time_evolve.jl") include("algorithms/time_evolution/simpleupdate.jl") include("algorithms/time_evolution/simpleupdate3site.jl") +include("algorithms/time_evolution/gaugefix_su.jl") + +include("algorithms/bp/beliefpropagation.jl") +include("algorithms/bp/gaugefix.jl") include("algorithms/transfermatrix.jl") include("algorithms/toolbox.jl") @@ -114,6 +120,9 @@ export InfinitePartitionFunction export InfinitePEPS, InfiniteTransferPEPS export SUWeight export InfinitePEPO, InfiniteTransferPEPO + +export BPEnv, BeliefPropagation, BPGauge + export initialize_mps, initializePEPS export ReflectDepth, ReflectWidth, Rotate, RotateReflect export symmetrize!, symmetrize_retract_and_finalize! diff --git a/src/algorithms/bp/beliefpropagation.jl b/src/algorithms/bp/beliefpropagation.jl new file mode 100644 index 000000000..aedc5deb8 --- /dev/null +++ b/src/algorithms/bp/beliefpropagation.jl @@ -0,0 +1,111 @@ +""" + struct BeliefPropagation + +Algorithm for computing the belief propagation fixed point messages. + +## Fields + +$(TYPEDFIELDS) +""" +@kwdef struct BeliefPropagation + "Stopping criterion for the BP iterations in relative trace norm difference" + tol::Float64 = 1.0e-6 + + "Minimal number of BP iterations" + miniter::Int = 2 + + "Maximal number of BP iterations" + maxiter::Int = 50 + + "Toggle for projecting messages onto the hermitian subspace immediately after update through BP equation" + project_hermitian::Bool = true + + "Output verbosity level" + verbosity::Int = 2 +end + +""" + leading_boundary(env₀::BPEnv, network, alg::BeliefPropagation) + +Contract `network` in the BP approximation and return the corresponding messages. +""" +function leading_boundary(env₀::BPEnv, network::InfiniteSquareNetwork, alg::BeliefPropagation) + return LoggingExtras.withlevel(; alg.verbosity) do + env = deepcopy(env₀) + log = MPSKit.IterLog("BP") + ϵ = Inf + @infov 1 loginit!(log, ϵ) + for iter in 1:(alg.maxiter) + env′ = bp_iteration(network, env, alg) + ϵ = oftype(ϵ, tr_distance(env, env′)) + env = env′ + + if ϵ <= alg.tol && iter >= alg.miniter + @infov 2 logfinish!(log, iter, ϵ) + break + end + if iter ≥ alg.maxiter + @warnv 1 logcancel!(log, iter, ϵ) + else + @infov 3 logiter!(log, iter, ϵ) + end + end + + return env, ϵ + end +end +function leading_boundary(env₀::BPEnv, state, args...) + return leading_boundary(env₀, InfiniteSquareNetwork(state), args...) +end + +""" +One iteration to update the BP environment. +""" +function bp_iteration(network::InfiniteSquareNetwork, env::BPEnv, alg::BeliefPropagation) + messages = map(eachindex(env)) do I + M = update_message(I, network, env) + normalize!(M) + alg.project_hermitian && (M = project_hermitian!!(M)) + return M + end + return BPEnv(messages) +end + +""" +Update the BP message in `env.messages[I]`. +""" +function update_message(I::CartesianIndex{3}, network::InfiniteSquareNetwork, env::BPEnv) + dir, row, col = Tuple(I) + (1 <= dir <= 4) || throw(ArgumentError("Invalid direction $dir")) + + A = network[row, col] + dir == SOUTH || (M_north = env[NORTH, _prev(row, end), col]) + dir == WEST || (M_east = env[EAST, row, _next(col, end)]) + dir == NORTH || (M_south = env[SOUTH, _next(row, end), col]) + dir == EAST || (M_west = env[WEST, row, _prev(col, end)]) + + return if dir == NORTH + contract_north_message(A, M_west, M_north, M_east) + elseif dir == EAST + contract_east_message(A, M_north, M_east, M_south) + elseif dir == SOUTH + contract_south_message(A, M_east, M_south, M_west) + else # dir == WEST + contract_west_message(A, M_south, M_west, M_north) + end +end + +function tr_distance(A::BPEnv, B::BPEnv) + return sum(zip(A.messages, B.messages)) do (a, b) + return trnorm(add(a, b, -inv(tr(b)), inv(tr(a)))) + end / length(A.messages) +end + +function trnorm(M::AbstractTensorMap, p::Real = 1) + return TensorKit._norm(svdvals(M), p, zero(real(scalartype(M)))) +end +function trnorm!(M::AbstractTensorMap, p::Real = 1) + return TensorKit._norm(svdvals!(M), p, zero(real(scalartype(M)))) +end + +project_hermitian!!(t) = add(t, t', 1 / 2, 1 / 2) diff --git a/src/algorithms/bp/gaugefix.jl b/src/algorithms/bp/gaugefix.jl new file mode 100644 index 000000000..126d66047 --- /dev/null +++ b/src/algorithms/bp/gaugefix.jl @@ -0,0 +1,127 @@ +""" + struct BPGauge + +Algorithm for gauging PEPS with belief propagation fixed point messages. +""" +@kwdef struct BPGauge + # TODO: add options +end + +""" +$(SIGNATURES) + +Fix the gauge of `psi` using fixed point environment `env` of belief propagation. +""" +function gauge_fix(psi::InfinitePEPS, alg::BPGauge, env::BPEnv) + psi′ = copy(psi) + XXinv = map(eachcoordinate(psi, 1:2)) do I + _, X, Xinv = _bp_gauge_fix!(CartesianIndex(I), psi′, env) + return X, Xinv + end + return psi′, XXinv +end + +function _sqrt_bp_messages(I::CartesianIndex{3}, env::BPEnv) + dir, row, col = Tuple(I) + @assert dir == NORTH || dir == EAST + M12 = env[dir, dir == NORTH ? _prev(row, end) : row, dir == EAST ? _next(col, end) : col] + sqrtM12, isqrtM12 = sqrt_invsqrt(twist(M12, 1)) + M21 = env[dir + 2, row, col] + sqrtM21, isqrtM21 = sqrt_invsqrt(M21) + return sqrtM12, isqrtM12, sqrtM21, isqrtM21 +end + +""" + _bp_gauge_fix!(I, psi::InfinitePEPS, env::BPEnv) -> psi, X, X⁻¹ + +For the bond at direction `I[1]` (which can be `NORTH` or `EAST`) +from site `I[2], I[3]`, we identify the following gauge matrices, +along the canonical direction of the PEPS arrows (`SOUTH ← NORTH` or `WEST ← EAST`): + +```math + I = √M₁₂⁻¹ √M₁₂ √M₂₁ √M₂₁⁻¹ + = √M₁₂⁻¹ (U Λ Vᴴ) √M₂₁⁻¹ + = (√M₁₂⁻¹ U √Λ) (√Λ Vᴴ √M₂₁⁻¹) + = X X⁻¹ +``` + +Which are then used to update the gauge of `psi`. Thus, by convention `X` is attached to the `SOUTH`/`WEST` directions +and `X⁻¹` is attached to the `NORTH`/`EAST` directions. +""" +function _bp_gauge_fix!(I::CartesianIndex{3}, psi::InfinitePEPS, env::BPEnv) + dir, row, col = Tuple(I) + @assert dir == NORTH || dir == EAST + + sqrtM12, isqrtM12, sqrtM21, isqrtM21 = _sqrt_bp_messages(I, env) + U, Λ, Vᴴ = svd_compact!(sqrtM12 * sqrtM21) + sqrtΛ = sdiag_pow(Λ, 1 / 2) + X = isqrtM12 * U * sqrtΛ + invX = sqrtΛ * Vᴴ * isqrtM21 + if isdual(space(sqrtM12, 1)) + X, invX = twist(flip(X, 2), 1), flip(invX, 1) + end + if dir == NORTH + psi[row, col] = absorb_north_message(psi[row, col], X) + psi[_prev(row, end), col] = absorb_south_message(psi[_prev(row, end), col], invX) + elseif dir == EAST + psi[row, col] = absorb_east_message(psi[row, col], X) + psi[row, _next(col, end)] = absorb_west_message(psi[row, _next(col, end)], invX) + end + return psi, X, invX +end + +""" + SUWeight(env::BPEnv) + +Construct `SUWeight` from belief propagation fixed point environment `env`. +""" +function SUWeight(env::BPEnv) + wts = map(Iterators.product(1:2, axes(env, 2), axes(env, 3))) do (dir′, row, col) + I = CartesianIndex(mod1(dir′ + 1, 2), row, col) + sqrtM12, _, sqrtM21, _ = _sqrt_bp_messages(I, env) + Λ = svd_vals!(sqrtM12 * sqrtM21) + return isdual(space(sqrtM12, 1)) ? _fliptwist_s(Λ) : Λ + end + return SUWeight(wts) +end + +""" + BPEnv(wts::SUWeight) + +Convert fixed point weights `wts` of trivial simple update +to a belief propagation environment. +""" +function BPEnv(wts::SUWeight) + messages = map(Iterators.product(1:4, axes(wts, 2), axes(wts, 3))) do (d, r, c) + wt = if d == NORTH + twist(wts[2, _next(r, end), c], 1) + elseif d == EAST + twist(wts[1, r, _prev(c, end)], 1) + elseif d == SOUTH + copy(wts[2, r, c]) + else # WEST + copy(wts[1, r, c]) + end + return TensorMap(wt) + end + return BPEnv(messages) +end + +function sqrt_invsqrt(A::PEPSMessage) + if isposdef(A) + D, V = eigh_full(A) + sqrtA = V * sdiag_pow(D, 1 / 2) * V' + isqrtA = V * sdiag_pow(D, -1 / 2) * V' + else + D, V = eig_full(A) + V⁻¹ = inv(V) + sqrtA = V * sdiag_pow(D, 1 / 2) * V⁻¹ + isqrtA = V * sdiag_pow(D, -1 / 2) * V⁻¹ + if scalartype(A) <: Real + # TODO: is this valid? + sqrtA = real(sqrtA) + isqrtA = real(isqrtA) + end + end + return sqrtA, isqrtA +end diff --git a/src/algorithms/contractions/bp_contractions.jl b/src/algorithms/contractions/bp_contractions.jl new file mode 100644 index 000000000..62e32177e --- /dev/null +++ b/src/algorithms/contractions/bp_contractions.jl @@ -0,0 +1,249 @@ +const PEPSMessage = AbstractTensorMap{<:Any, <:Any, 1, 1} + +# Belief Propagation Updates +# -------------------------- +function contract_north_message( + A::PEPSSandwich, M_west::PEPSMessage, M_north::PEPSMessage, M_east::PEPSMessage + ) + return @autoopt @tensor begin + M_north′[DSt; DSb] := + ket(A)[d; DNt DEt DSt DWt] * conj(bra(A)[d; DNb DEb DSb DWb]) * + M_west[DWb; DWt] * M_north[DNt; DNb] * M_east[DEt; DEb] + end +end +function contract_east_message( + A::PEPSSandwich, M_north::PEPSMessage, M_east::PEPSMessage, M_south::PEPSMessage + ) + return @autoopt @tensor begin + M_east′[DWt; DWb] := + ket(A)[d; DNt DEt DSt DWt] * conj(bra(A)[d; DNb DEb DSb DWb]) * + M_north[DNt; DNb] * M_east[DEt; DEb] * M_south[DSb; DSt] + end +end +function contract_south_message( + A::PEPSSandwich, M_east::PEPSMessage, M_south::PEPSMessage, M_west::PEPSMessage + ) + return @autoopt @tensor begin + M_south′[DNb; DNt] := + ket(A)[d; DNt DEt DSt DWt] * conj(bra(A)[d; DNb DEb DSb DWb]) * + M_east[DEt; DEb] * M_south[DSb; DSt] * M_west[DWb; DWt] + end +end +function contract_west_message( + A::PEPSSandwich, M_south::PEPSMessage, M_west::PEPSMessage, M_north::PEPSMessage + ) + return @autoopt @tensor begin + M_west′[DEb; DEt] := + ket(A)[d; DNt DEt DSt DWt] * conj(bra(A)[d; DNb DEb DSb DWb]) * + M_south[DSb; DSt] * M_west[DWb; DWt] * M_north[DNt; DNb] + end +end + +absorb_north_message(A::PEPSTensor, M::PEPSMessage) = + @tensor A′[d; N' E S W] := A[d; N E S W] * M[N; N'] +absorb_east_message(A::PEPSTensor, M::PEPSMessage) = + @tensor A′[d; N E' S W] := A[d; N E S W] * M[E; E'] +absorb_south_message(A::PEPSTensor, M::PEPSMessage) = + @tensor A′[d; N E S' W] := A[d; N E S W] * M[S'; S] +absorb_west_message(A::PEPSTensor, M::PEPSMessage) = + @tensor A′[d; N E S W'] := A[d; N E S W] * M[W'; W] + +# Belief Propagation Expectation values +# ------------------------------------- +function MPSKit.expectation_value(peps::InfinitePEPS, O::LocalOperator, env::BPEnv) + checklattice(peps, O) + term_vals = dtmap([O.terms...]) do (inds, operator) # OhMyThreads can't iterate over O.terms directly + contract_local_operator(inds, operator, peps, peps, env) / + contract_local_norm(inds, peps, peps, env) + end + return sum(term_vals) +end + +function contract_local_operator( + inds::NTuple{1, CartesianIndex{2}}, + O::AbstractTensorMap{<:Any, <:Any, 1, 1}, + ket::InfinitePEPS, + bra::InfinitePEPS, + env::BPEnv, + ) + row, col = Tuple(only(inds)) + M_north = env.messages[NORTH, _prev(row, end), mod1(col, end)] + M_east = env.messages[EAST, mod1(row, end), _next(col, end)] + M_south = env.messages[SOUTH, _next(row, end), mod1(col, end)] + M_west = env.messages[WEST, mod1(row, end), _prev(col, end)] + + return @autoopt @tensor begin + ket[mod1(row, end), mod1(col, end)][dt; DNt DEt DSt DWt] * + conj(bra[mod1(row, end), mod1(col, end)][db; DNb DEb DSb DWb]) * + O[db; dt] * + M_north[DNt; DNb] * + M_east[DEt; DEb] * + M_south[DSb; DSt] * + M_west[DWb; DWt] + end +end + +function contract_local_norm( + inds::NTuple{1, CartesianIndex{2}}, ket::InfinitePEPS, bra::InfinitePEPS, env::BPEnv + ) + row, col = Tuple(only(inds)) + M_north = env.messages[NORTH, _prev(row, end), mod1(col, end)] + M_east = env.messages[EAST, mod1(row, end), _next(col, end)] + M_south = env.messages[SOUTH, _next(row, end), mod1(col, end)] + M_west = env.messages[WEST, mod1(row, end), _prev(col, end)] + + return @autoopt @tensor begin + ket[mod1(row, end), mod1(col, end)][d; DNt DEt DSt DWt] * + conj(bra[mod1(row, end), mod1(col, end)][d; DNb DEb DSb DWb]) * + M_north[DNt; DNb] * + M_east[DEt; DEb] * + M_south[DSb; DSt] * + M_west[DWb; DWt] + end +end + +function contract_local_operator( + inds::NTuple{2, CartesianIndex{2}}, + O::AbstractTensorMap{<:Any, <:Any, 2, 2}, + ket::InfinitePEPS, + bra::InfinitePEPS, + env::BPEnv, + ) + ind_relative = inds[2] - inds[1] + return if ind_relative == CartesianIndex(1, 0) + contract_vertical_operator(inds[1], O, ket, bra, env) + elseif ind_relative == CartesianIndex(0, 1) + contract_horizontal_operator(inds[1], O, ket, bra, env) + else + error("Only contractions for nearest neighbor bonds are implemented.") + end +end + +function contract_vertical_operator( + coord::CartesianIndex{2}, + O::AbstractTensorMap{<:Any, <:Any, 2, 2}, + ket::InfinitePEPS, + bra::InfinitePEPS, + env::BPEnv, + ) + row, col = Tuple(coord) + M_north = env.messages[NORTH, _prev(row, end), mod1(col, end)] + M_northeast = env.messages[EAST, mod1(row, end), _next(col, end)] + M_southeast = env.messages[EAST, _next(row, end), _next(col, end)] + M_south = env.messages[SOUTH, mod1(row + 2, end), mod1(col, end)] + M_southwest = env.messages[WEST, _next(row, end), _prev(col, end)] + M_northwest = env.messages[WEST, mod1(row, end), _prev(col, end)] + + return @autoopt @tensor ket[mod1(row, end), mod1(col, end)][dNt; DNt DNEt DMt DNWt] * + ket[_next(row, end), mod1(col, end)][dSt; DMt DSEt DSt DSWt] * + conj(bra[mod1(row, end), mod1(col, end)][dNb; DNb DNEb DMb DNWb]) * + conj(bra[_next(row, end), mod1(col, end)][dSb; DMb DSEb DSb DSWb]) * + M_north[DNt; DNb] * + M_northeast[DNEt; DNEb] * + M_southeast[DSEt; DSEb] * + M_south[DSb; DSt] * + M_southwest[DSWb; DSWt] * + M_northwest[DNWb; DNWt] * + O[dNb dSb; dNt dSt] +end + +function contract_horizontal_operator( + coord::CartesianIndex{2}, + O::AbstractTensorMap{<:Any, <:Any, 2, 2}, + ket::InfinitePEPS, + bra::InfinitePEPS, + env::BPEnv, + ) + row, col = Tuple(coord) + M_west = env.messages[WEST, mod1(row, end), _prev(col, end)] + M_northwest = env.messages[NORTH, _prev(row, end), mod1(col, end)] + M_northeast = env.messages[NORTH, _prev(row, end), _next(col, end)] + M_east = env.messages[EAST, mod1(row, end), mod1(col + 2, end)] + M_southeast = env.messages[SOUTH, _next(row, end), _next(col, end)] + M_southwest = env.messages[SOUTH, _next(row, end), mod1(col, end)] + A_west = ket[mod1(row, end), mod1(col, end)] + Ā_west = bra[mod1(row, end), mod1(col, end)] + A_east = ket[mod1(row, end), _next(col, end)] + Ā_east = bra[mod1(row, end), _next(col, end)] + + return @autoopt @tensor begin + A_west[dWt; DNWt DMt DSWt DWt] * + A_east[dEt; DNEt DEt DSEt DMt] * + conj(Ā_west[dWb; DNWb DMb DSWb DWb]) * + conj(Ā_east[dEb; DNEb DEb DSEb DMb]) * + M_west[DWb; DWt] * + M_northwest[DNWt; DNWb] * + M_northeast[DNEt; DNEb] * + M_east[DEt; DEb] * + M_southeast[DSEb; DSEt] * + M_southwest[DSWb; DSWt] * + O[dWb dEb; dWt dEt] + end +end + +function contract_local_norm( + inds::NTuple{2, CartesianIndex{2}}, ket::InfinitePEPS, bra::InfinitePEPS, env::BPEnv + ) + ind_relative = inds[2] - inds[1] + return if ind_relative == CartesianIndex(1, 0) + contract_vertical_norm(inds[1], ket, bra, env) + elseif ind_relative == CartesianIndex(0, 1) + contract_horizontal_norm(inds[1], ket, bra, env) + else + error("Only contractions for nearest neighbor bonds are implemented.") + end +end + +function contract_vertical_norm( + coord::CartesianIndex{2}, ket::InfinitePEPS, bra::InfinitePEPS, env::BPEnv + ) + row, col = Tuple(coord) + M_north = env.messages[NORTH, _prev(row, end), mod1(col, end)] + M_northeast = env.messages[EAST, mod1(row, end), _next(col, end)] + M_southeast = env.messages[EAST, _next(row, end), _next(col, end)] + M_south = env.messages[SOUTH, mod1(row + 2, end), mod1(col, end)] + M_southwest = env.messages[WEST, _next(row, end), _prev(col, end)] + M_northwest = env.messages[WEST, mod1(row, end), _prev(col, end)] + + return @autoopt @tensor ket[mod1(row, end), mod1(col, end)][dN; DNt DNEt DMt DNWt] * + ket[_next(row, end), mod1(col, end)][dS; DMt DSEt DSt DSWt] * + conj(bra[mod1(row, end), mod1(col, end)][dN; DNb DNEb DMb DNWb]) * + conj(bra[_next(row, end), mod1(col, end)][dS; DMb DSEb DSb DSWb]) * + M_north[DNt; DNb] * + M_northeast[DNEt; DNEb] * + M_southeast[DSEt; DSEb] * + M_south[DSb; DSt] * + M_southwest[DSWb; DSWt] * + M_northwest[DNWb; DNWt] +end + +function contract_horizontal_norm( + coord::CartesianIndex{2}, ket::InfinitePEPS, bra::InfinitePEPS, env::BPEnv + ) + row, col = Tuple(coord) + + M_west = env.messages[WEST, mod1(row, end), _prev(col, end)] + M_northwest = env.messages[NORTH, _prev(row, end), mod1(col, end)] + M_northeast = env.messages[NORTH, _prev(row, end), _next(col, end)] + M_east = env.messages[EAST, mod1(row, end), mod1(col + 2, end)] + M_southeast = env.messages[SOUTH, _next(row, end), _next(col, end)] + M_southwest = env.messages[SOUTH, _next(row, end), mod1(col, end)] + + A_west = ket[mod1(row, end), mod1(col, end)] + Ā_west = bra[mod1(row, end), mod1(col, end)] + A_east = ket[mod1(row, end), _next(col, end)] + Ā_east = bra[mod1(row, end), _next(col, end)] + + return @autoopt @tensor begin + A_west[dW; DNWt DMt DSWt DWt] * + A_east[dE; DNEt DEt DSEt DMt] * + conj(Ā_west[dW; DNWb DMb DSWb DWb]) * + conj(Ā_east[dE; DNEb DEb DSEb DMb]) * + M_west[DWb; DWt] * + M_northwest[DNWt; DNWb] * + M_northeast[DNEt; DNEb] * + M_east[DEt; DEb] * + M_southeast[DSEb; DSEt] * + M_southwest[DSWb; DSWt] + end +end diff --git a/src/algorithms/time_evolution/gaugefix_su.jl b/src/algorithms/time_evolution/gaugefix_su.jl new file mode 100644 index 000000000..87c215dc6 --- /dev/null +++ b/src/algorithms/time_evolution/gaugefix_su.jl @@ -0,0 +1,62 @@ +""" + struct SUGauge + +Algorithm for fixing gauge of an iPEPS using trivial simple update +(with identity gates). + +## Fields + +$(TYPEDFIELDS) +""" +@kwdef struct SUGauge + "Stopping criterion for the trivial SU iterations in weight difference" + tol::Float64 = 1.0e-10 + "Minimal number of SU iterations" + miniter::Int = 2 + "Maximal number of SU iterations" + maxiter::Int = 100 +end + +""" +A LocalOperator consisting of identity gates on all nearest neighbor bonds. +""" +function _trivial_gates(elt::Type{<:Number}, lattice::Matrix{S}) where {S <: ElementarySpace} + terms = [] + for site1 in CartesianIndices(lattice) + r1, c1 = mod1.(Tuple(site1), size(lattice)) + for d in (CartesianIndex(1, 0), CartesianIndex(0, 1)) + site2 = site1 + d + r2, c2 = mod1.(Tuple(site2), size(lattice)) + V1, V2 = lattice[r1, c1], lattice[r2, c2] + h = TensorKit.id(elt, V1 ⊗ V2) + push!(terms, (site1, site2) => h) + end + end + return LocalOperator(lattice, terms...) +end + +""" +$(SIGNATURES) + +Fix the gauge of `psi` using trivial simple update. +""" +function gauge_fix(psi::InfinitePEPS, alg::SUGauge) + gates = _trivial_gates(scalartype(psi), physicalspace(psi)) + su_alg = SimpleUpdate(; trunc = FixedSpaceTruncation()) + 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: |Δλ| = $ϵ." + return psi′, wts, ϵ + end + if i == alg.maxiter + @warn "Trivial SU cancel $i: |Δλ| = $ϵ." + return psi′, wts, ϵ + end + wts0 = deepcopy(wts) + end + return +end diff --git a/src/environments/bp_environments.jl b/src/environments/bp_environments.jl new file mode 100644 index 000000000..9e868736d --- /dev/null +++ b/src/environments/bp_environments.jl @@ -0,0 +1,225 @@ +""" +$(TYPEDEF) + +Belief propagation (BP) environment for a square lattice norm network, +containing a 4 x rows x cols array of message tensors, defined for +each *oriented* nearest neighbor bond in the network. + +The message tensors connect to the network tensors +`P` at site `[r,c]` in the unit cell as: +``` + m[1,r-1,c] + | + m[4,r,c-1]------P[r,c]------m[2,r,c+1] + | + m[3,r+1,c] +``` +- `[1,r-1,c]`: message from `P[r-1,c]` to `P[r,c]` (axis order: ket ← bra) +- `[2,r,c+1]`: message from `P[r,c+1]` to `P[r,c]` (axis order: ket ← bra) +- `[3,r+1,c]`: message from `P[r+1,c]` to `P[r,c]` (axis order: bra ← ket) +- `[4,r,c-1]`: message from `P[r,c-1]` to `P[r,c]` (axis order: bra ← ket) + +## Fields + +$(TYPEDFIELDS) +""" +struct BPEnv{T} + "4 x rows x cols array of message tensors, where the first dimension specifies the spatial direction" + messages::Array{T, 3} +end + +""" +Construct a message tensor on a certain bond of a network, +with bond space specified by `pspaces`. +In the 2-layer case, the message tensor will be +``` + ┌--- pspaces[1] (ket layer) + m + └--- pspaces[2] (bra layer) +``` +Returning axis order is `bra ← ket`. + +When `posdef` is true, the message will be made semi-posdef definite +(when interpreted as a `bra ← ket` TensorMap). +""" +function _message_tensor(f, ::Type{T}, pspaces::P; posdef::Bool = true) where {T, P <: ProductSpace} + N = length(pspaces) + @assert N == 2 "BPEnv is currently only defined for 2-layer InfiniteSquareNetwork." + m = f(T, pspaces[2] ← pspaces[1]') + return posdef ? (m' * m) : m +end + +""" + BPEnv( + [f=randn, T=ComplexF64], Ds_north::A, Ds_east::A; posdef::Bool = true + ) where {A <: AbstractMatrix{<:ProductSpace}} + +Construct a BP environment by specifying matrices of north and east virtual spaces of the +corresponding [`InfiniteSquareNetwork`](@ref). Each matrix entry corresponds to a site in the unit cell. + +When `posdef` is true, all messages will be made semi-positive definite +(when interpreted as a `bra ← ket` TensorMap). + +Each entry of the `Ds_north` and `Ds_east` matrices corresponds to an effective local space +of the network, and can be represented as a `ProductSpace` (e.g. +for the case of a network representing overlaps of PEPSs). +""" +function BPEnv( + Ds_north::A, Ds_east::A; posdef::Bool = true + ) where {A <: AbstractMatrix{<:ProductSpace}} + return BPEnv(randn, ComplexF64, N, Ds_north, Ds_east; posdef) +end +function BPEnv( + f, T, Ds_north::A, Ds_east::A; posdef::Bool = true + ) where {A <: AbstractMatrix{<:ProductSpace}} + # no recursive broadcasting? + Ds_south = _elementwise_dual.(circshift(Ds_north, (-1, 0))) + Ds_west = _elementwise_dual.(circshift(Ds_east, (0, 1))) + messages = map(Iterators.product(1:4, axes(Ds_north, 1), axes(Ds_north, 2))) do (dir, r, c) + msg = if dir == NORTH + _message_tensor(f, T, Ds_north[_next(r, end), c]; posdef) + elseif dir == EAST + _message_tensor(f, T, Ds_east[r, _prev(c, end)]; posdef) + elseif dir == SOUTH + _message_tensor(f, T, Ds_south[_prev(r, end), c]; posdef) + else # WEST + _message_tensor(f, T, Ds_west[r, _next(c, end)]; posdef) + end + if dir == NORTH || dir == EAST + msg = permute(msg, ((2,), (1,))) + end + return msg + end + normalize!.(messages) + return BPEnv(messages) +end + +""" + BPEnv( + [f=randn, T=ComplexF64], D_north::P, D_east::P; + unitcell::Tuple{Int, Int} = (1, 1), posdef::Bool = true + ) where {P <: ProductSpace} + +Construct a BP environment by specifying the north and east virtual spaces of the +corresponding [`InfiniteSquareNetwork`](@ref). The network unit cell can be specified +by the `unitcell` keyword argument. +""" +function BPEnv( + D_north::P, D_east::P; + unitcell::Tuple{Int, Int} = (1, 1), posdef::Bool = true + ) where {P <: ProductSpace} + return BPEnv(randn, ComplexF64, D_north, D_east; unitcell, posdef) +end +function BPEnv( + f, T, D_north::P, D_east::P; + unitcell::Tuple{Int, Int} = (1, 1), posdef::Bool = true + ) where {P <: ProductSpace} + return BPEnv(f, T, N, fill(D_north, unitcell), fill(D_east, unitcell); posdef) +end + +""" + BPEnv([f=randn, T=ComplexF64], network::InfiniteSquareNetwork; posdef::Bool = true) + +Construct a BP environment by specifying a corresponding [`InfiniteSquareNetwork`](@ref). +""" +function BPEnv(f, T, network::InfiniteSquareNetwork; posdef::Bool = true) + Ds_north = _north_edge_physical_spaces(network) + Ds_east = _east_edge_physical_spaces(network) + return BPEnv(f, T, Ds_north, Ds_east; posdef) +end +function BPEnv(network::InfiniteSquareNetwork; posdef::Bool = true) + return BPEnv(randn, scalartype(network), network; posdef) +end + +function BPEnv(state::InfinitePartitionFunction, args...; kwargs...) + return BPEnv(InfiniteSquareNetwork(state), args...; kwargs...) +end +function BPEnv(state::InfinitePEPS, args...; kwargs...) + bp_env = BPEnv(InfiniteSquareNetwork(state), args...; kwargs...) + TensorKit.id!.(bp_env.messages) + return bp_env +end +function BPEnv(f, T, state::Union{InfinitePartitionFunction, InfinitePEPS}, args...; kwargs...) + return BPEnv(f, T, InfiniteSquareNetwork(state), args...; kwargs...) +end + +Base.eltype(::Type{BPEnv{T}}) where {T} = T +Base.size(env::BPEnv, args...) = size(env.messages, args...) +Base.getindex(env::BPEnv, args...) = Base.getindex(env.messages, args...) +Base.axes(env::BPEnv, args...) = Base.axes(env.messages, args...) +Base.eachindex(env::BPEnv) = eachindex(IndexCartesian(), env.messages) +VectorInterface.scalartype(::Type{BPEnv{T}}) where {T} = scalartype(T) +TensorKit.spacetype(::Type{BPEnv{T}}) where {T} = spacetype(T) + +function eachcoordinate(x::BPEnv) + return collect(Iterators.product(axes(x, 2), axes(x, 3))) +end +function eachcoordinate(x::BPEnv, dirs) + return collect(Iterators.product(dirs, axes(x, 2), axes(x, 3))) +end + +# conversion to CTMRGEnv +""" + CTMRGEnv(bp_env::BPEnv) + +Construct a CTMRG environment with bond dimension χ = 1 +from the belief propagation environment `bp_env`. +""" +function CTMRGEnv(bp_env::BPEnv) + edges = map(CartesianIndices(bp_env.messages)) do idx + M = bp_env.messages[idx] + M = if idx[1] == SOUTH || idx[1] == WEST + permute(M, ((2, 1), ())) + else + repartition(M, numind(M), 0; copy = true) + end + return insertleftunit(insertleftunit(M), 1) + end + corners = map(CartesianIndices(edges)) do _ + return TensorKit.id(scalartype(bp_env), oneunit(spacetype(bp_env))) + end + return CTMRGEnv(corners, edges) +end + +# rotation (the behavior is the same as CTMRGEnv edges) + +function Base.rotl90(env::BPEnv{T}) where {T} + messages′ = Array{T, 3}(undef, 4, size(env.messages, 3), size(env.messages, 2)) + for dir in 1:4 + dir2 = _prev(dir, 4) + messages′[dir2, :, :] = rotl90(env.messages[dir, :, :]) + end + env2 = BPEnv(copy(messages′)) + for idx in eachcoordinate(env2, (EAST, WEST)) + cidx = CartesianIndex(idx) + env2.messages[cidx] = permute(env2.messages[cidx], ((2,), (1,))) + end + return env2 +end + +function Base.rotr90(env::BPEnv{T}) where {T} + messages′ = Array{T, 3}(undef, 4, size(env.messages, 3), size(env.messages, 2)) + for dir in 1:4 + dir2 = _next(dir, 4) + messages′[dir2, :, :] = rotr90(env.messages[dir, :, :]) + end + env2 = BPEnv(copy(messages′)) + for idx in eachcoordinate(env2, (NORTH, SOUTH)) + cidx = CartesianIndex(idx) + env2.messages[cidx] = permute(env2.messages[cidx], ((2,), (1,))) + end + return env2 +end + +function Base.rot180(env::BPEnv{T}) where {T} + messages′ = Array{T, 3}(undef, 4, size(env.messages, 2), size(env.messages, 3)) + for dir in 1:4 + dir2 = _next(_next(dir, 4), 4) + messages′[dir2, :, :] = rot180(env.messages[dir, :, :]) + end + env2 = BPEnv(copy(messages′)) + for cidx in CartesianIndices(env2.messages) + env2.messages[cidx] = permute(env2.messages[cidx], ((2,), (1,))) + end + return env2 +end diff --git a/src/environments/ctmrg_environments.jl b/src/environments/ctmrg_environments.jl index 5d4895727..7adfb99e0 100644 --- a/src/environments/ctmrg_environments.jl +++ b/src/environments/ctmrg_environments.jl @@ -341,6 +341,17 @@ end Base.:*(α::Number, e::CTMRGEnv) = CTMRGEnv(α * e.corners, α * e.edges) Base.similar(e::CTMRGEnv) = CTMRGEnv(similar(e.corners), similar(e.edges)) +# (approximate) equality +function Base.isapprox(env1::CTMRGEnv, env2::CTMRGEnv; kwargs...) + for (c1, c2) in zip(env1.corners, env2.corners) + !isapprox(c1, c2; kwargs...) && return false + end + for (e1, e2) in zip(env1.edges, env2.edges) + !isapprox(e1, e2; kwargs...) && return false + end + return true +end + # VectorInterface # --------------- diff --git a/src/networks/local_sandwich.jl b/src/networks/local_sandwich.jl index 955c89a7f..a980bbc6f 100644 --- a/src/networks/local_sandwich.jl +++ b/src/networks/local_sandwich.jl @@ -5,7 +5,7 @@ ## Space utils _elementwise_dual(S::ElementarySpace) = dual(S) -_elementwise_dual(P::ProductSpace) = ProductSpace(dual.(P)...) +_elementwise_dual(P::ProductSpace{S, N}) where {S, N} = ProductSpace{S, N}(dual.(P)...) # route all virtualspace getters through a single method for convenience north_virtualspace(O, args...) = virtualspace(O, args..., NORTH) diff --git a/src/utility/util.jl b/src/utility/util.jl index a63a2a97f..3fa7e0cc1 100644 --- a/src/utility/util.jl +++ b/src/utility/util.jl @@ -48,20 +48,20 @@ end _safe_pow(a::Number, pow::Real, tol::Real) = (pow < 0 && abs(a) < tol) ? zero(a) : a^pow """ - sdiag_pow(s, pow::Real; tol::Real=eps(scalartype(s))^(3 / 4)) + sdiag_pow(s, pow::Real; tol::Real=eps(real(scalartype(s)))^(3 / 4)) Compute `s^pow` for a diagonal matrix `s`. """ -function sdiag_pow(s::DiagonalTensorMap, pow::Real; tol::Real = eps(scalartype(s))^(3 / 4)) - # Relative tol w.r.t. largest singular value (use norm(∘, Inf) to make differentiable) +function sdiag_pow(s::DiagonalTensorMap, pow::Real; tol::Real = eps(real(scalartype(s)))^(3 / 4)) + # Relative tol w.r.t. largest abs value of `s` (use norm(∘, Inf) to make differentiable) tol *= norm(s, Inf) spow = DiagonalTensorMap(_safe_pow.(s.data, pow, tol), space(s, 1)) return spow end function sdiag_pow( - s::AbstractTensorMap{T, S, 1, 1}, pow::Real; tol::Real = eps(scalartype(s))^(3 / 4) + s::AbstractTensorMap{T, S, 1, 1}, pow::Real; tol::Real = eps(real(scalartype(s)))^(3 / 4) ) where {T, S} - # Relative tol w.r.t. largest singular value (use norm(∘, Inf) to make differentiable) + # Relative tol w.r.t. largest abs value of `s` (use norm(∘, Inf) to make differentiable) tol *= norm(s, Inf) spow = similar(s) for (k, b) in blocks(s) @@ -74,7 +74,7 @@ end function ChainRulesCore.rrule( ::typeof(sdiag_pow), s::AbstractTensorMap, pow::Real; - tol::Real = eps(scalartype(s))^(3 / 4), + tol::Real = eps(real(scalartype(s)))^(3 / 4), ) tol *= norm(s, Inf) spow = sdiag_pow(s, pow; tol) @@ -205,3 +205,13 @@ macro showtypeofgrad(x) end ) end + +""" +Randomly take the dual of `ElementarySpace`s in `Vs` with propability `p` +""" +function random_dual!(Vs::AbstractMatrix{E}; p = 0.7) where {E <: ElementarySpace} + for (i, V) in enumerate(Vs) + (rand() < p) && (Vs[i] = V') + end + return Vs +end diff --git a/test/bp/expvals.jl b/test/bp/expvals.jl new file mode 100644 index 000000000..d67cc3ba1 --- /dev/null +++ b/test/bp/expvals.jl @@ -0,0 +1,54 @@ +using Test +using Random +using TensorKit +using PEPSKit +using PEPSKit: random_dual!, SUGauge, gauge_fix + +ds = Dict( + U1Irrep => U1Space(i => d for (i, d) in zip(-1:1, (1, 1, 2))), + FermionParity => Vect[FermionParity](0 => 2, 1 => 1) +) +Ds = Dict( + U1Irrep => U1Space(i => D for (i, D) in zip(-1:1, (1, 3, 2))), + FermionParity => Vect[FermionParity](0 => 3, 1 => 2) +) +Random.seed!(41973582) + +@testset "Expectation values of BPEnv ($S)" for S in keys(ds) + d, D, uc = ds[S], Ds[S], (2, 3) + ψds = fill(d, uc) + ψDNs = random_dual!(fill(D, uc)) + ψDEs = random_dual!(fill(D, uc)) + ψ0 = InfinitePEPS(ψds, ψDNs, ψDEs) + + ψ, wts, _ = gauge_fix(ψ0, SUGauge(; maxiter = 100, tol = 1.0e-10)) + for (a0, a) in zip(ψ0.A, ψ.A) + @test space(a0) == space(a) + end + bp_env = BPEnv(wts) + ctm_env = CTMRGEnv(wts) + @test ctm_env ≈ CTMRGEnv(bp_env) + + # SU fixed point wts should already be a BP fixed point of ψ + bp_alg = BeliefPropagation(; miniter = 1, maxiter = 1, tol = 1.0e-7) + _, err = leading_boundary(bp_env, ψ, bp_alg) + @test err < 1.0e-9 + + op = randn(d → d) + for site in CartesianIndices(size(ψ)) + lo = LocalOperator(ψds, (site,) => op) + val1 = expectation_value(ψ, lo, bp_env) + val2 = expectation_value(ψ, lo, ctm_env) + @test val1 ≈ val2 + end + + op = randn(d ⊗ d → d ⊗ d) + vs = [CartesianIndex(1, 0), CartesianIndex(0, 1)] + for site1 in CartesianIndices(size(ψ)), v in vs + site2 = site1 + v + lo = LocalOperator(ψds, (site1, site2) => op) + val1 = expectation_value(ψ, lo, bp_env) + val2 = expectation_value(ψ, lo, ctm_env) + @test val1 ≈ val2 + end +end diff --git a/test/bp/gaugefix.jl b/test/bp/gaugefix.jl new file mode 100644 index 000000000..e53bc6a8f --- /dev/null +++ b/test/bp/gaugefix.jl @@ -0,0 +1,64 @@ +using Test, TestExtras +using Random +using TensorKit +using PEPSKit +using PEPSKit: SUGauge, gauge_fix, compare_weights, random_dual!, twistdual + +@testset "Compare BP and SU ($S, posdef msgs = $h)" for (S, h) in + Iterators.product([U1Irrep, FermionParity], [true, false]) + unitcell = (2, 3) + elt = ComplexF64 + maxiter, tol = 100, 1.0e-9 + Random.seed!(52840679) + Pspaces, Nspaces, Espaces = if S == U1Irrep + map(zip(rand(1:2, unitcell), rand(1:2, unitcell), rand(1:2, unitcell))) do (d0, d1, d2) + Vect[S](0 => d0, 1 => d1, -1 => d2) + end, + map(zip(rand(2:4, unitcell), rand(2:4, unitcell), rand(2:4, unitcell))) do (d0, d1, d2) + Vect[S](0 => d0, 1 => d1, -1 => d2) + end, + map(zip(rand(2:4, unitcell), rand(2:4, unitcell), rand(2:4, unitcell))) do (d0, d1, d2) + Vect[S](0 => d0, 1 => d1, -1 => d2) + end + else + map(zip(rand(2:3, unitcell), rand(2:3, unitcell))) do (d0, d1) + Vect[S](0 => d0, 1 => d1) + end, + map(zip(rand(2:4, unitcell), rand(2:4, unitcell))) do (d0, d1) + Vect[S](0 => d0, 1 => d1) + end, + map(zip(rand(2:4, unitcell), rand(2:4, unitcell))) do (d0, d1) + Vect[S](0 => d0, 1 => d1) + end + end + Nspaces, Espaces = random_dual!(Nspaces), random_dual!(Espaces) + peps0 = InfinitePEPS(randn, elt, Pspaces, Nspaces, Espaces) + + # 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 + normalize!.(wts1.data) + + # find BP fixed point and SUWeight + bp_alg = BeliefPropagation(; maxiter, tol, project_hermitian = h) + env = BPEnv(randn, elt, peps1; posdef = h) + env, err = leading_boundary(env, peps1, bp_alg) + 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) + for (a1, a2) in zip(peps1.A, peps2.A) + @test space(a1) == space(a2) + end + for (X, Xinv) in XXinv + # X, Xinv should contract to identity + @tensor tmp[-1; -2] := X[-1; 1] * Xinv[1; -2] + @test tmp ≈ twistdual(TensorKit.id(space(X, 1)), 1) + # BP should differ from SU only by a unitary gauge transformation + @test inv(X) ≈ adjoint(X) ≈ Xinv + end +end diff --git a/test/bp/rotation.jl b/test/bp/rotation.jl new file mode 100644 index 000000000..8b6fc7428 --- /dev/null +++ b/test/bp/rotation.jl @@ -0,0 +1,45 @@ +using Test +using Random +using TensorKit +using PEPSKit +using PEPSKit: random_dual! + +ds = Dict( + Trivial => ℂ^2, + U1Irrep => U1Space(i => d for (i, d) in zip(-1:1, (1, 1, 2))), + FermionParity => Vect[FermionParity](0 => 2, 1 => 1) +) +Ds = Dict( + Trivial => ℂ^3, + U1Irrep => U1Space(i => D for (i, D) in zip(-1:1, (1, 3, 2))), + FermionParity => Vect[FermionParity](0 => 3, 1 => 2) +) +Random.seed!(41973582) + +function meas_sites( + op::O, ψ::InfinitePEPS, env::Union{BPEnv, CTMRGEnv} + ) where {O <: AbstractTensorMap{<:Any, <:Any, 1, 1}} + lattice = physicalspace(ψ) + return map(CartesianIndices(ψ.A)) do site1 + lo = LocalOperator(lattice, (site1,) => op) + return expectation_value(ψ, lo, env) + end +end + +@testset "Rotation of BPEnv ($S)" for S in keys(ds) + d, D, unitcell = ds[S], Ds[S], (2, 3) + ψds = fill(d, unitcell) + ψDNs = random_dual!(fill(D, unitcell)) + ψDEs = random_dual!(fill(D, unitcell)) + ψ = InfinitePEPS(ψds, ψDNs, ψDEs) + env = BPEnv(ψ) + + op = randn(d → d) + meas1 = meas_sites(op, ψ, env) + # rotated peps and env + for f in (rotl90, rotr90, rot180) + ψ′, env′ = f(ψ), f(env) + meas1′ = meas_sites(op, ψ′, env′) + @test meas1′ ≈ f(meas1) + end +end diff --git a/test/bp/unitcell.jl b/test/bp/unitcell.jl new file mode 100644 index 000000000..81e2896be --- /dev/null +++ b/test/bp/unitcell.jl @@ -0,0 +1,98 @@ +using Test +using Random +using PEPSKit +using PEPSKit: bp_iteration +using TensorKit + +# settings +Random.seed!(91283219347) +elt = ComplexF64 + +function test_unitcell(unitcell, Pspaces, Nspaces, Espaces) + peps = InfinitePEPS(randn, elt, Pspaces, Nspaces, Espaces) + env0 = BPEnv(ones, elt, peps) + alg = BeliefPropagation() + + # apply one BP iteration + network = InfiniteSquareNetwork(peps) + env1 = bp_iteration(network, env0, alg) + # another iteration to detect bond mismatches + env1 = bp_iteration(network, env1, alg) + + # compute random expecation value to test matching bonds + random_op = LocalOperator( + Pspaces, ( + (c,) => randn(elt, Pspaces[c], Pspaces[c]) + for c in CartesianIndices(unitcell) + )..., + ) + @test expectation_value(peps, random_op, env0) isa Number + @test expectation_value(peps, random_op, env1) isa Number + return +end + +@testset "Random Cartesian spaces with BP" begin + unitcell = (3, 3) + + Pspaces = ComplexSpace.(rand(2:3, unitcell...)) + Nspaces = ComplexSpace.(rand(2:4, unitcell...)) + Espaces = ComplexSpace.(rand(2:4, unitcell...)) + + test_unitcell(unitcell, Pspaces, Nspaces, Espaces) +end + +@testset "Specific U1 spaces with BP" begin + unitcell = (2, 2) + + PA = U1Space(-1 => 1, 0 => 1) + PB = U1Space(0 => 1, 1 => 1) + Vpeps = U1Space(-1 => 2, 0 => 1, 1 => 2) + Venv = U1Space(-2 => 2, -1 => 3, 0 => 4, 1 => 3, 2 => 2) + + Pspaces = [PA PB; PB PA] + Nspaces = [Vpeps Vpeps'; Vpeps' Vpeps] + + test_unitcell(unitcell, Pspaces, Nspaces, Nspaces) + + # 4x4 unit cell with all 32 inequivalent bonds + # + # 10 4 7 32 + # | | | | + # 3--A--1--B--5--C--8--D--3 + # | | | | + # 2 6 9 11 + # | | | | + # 14--E-12--F-15--G-17--H-14 + # | | | | + # 13 16 18 19 + # | | | | + # 22--I-20--J-23--K-25--L-22 + # | | | | + # 21 24 26 27 + # | | | | + # 29--M-28--N-30--O-31--P-29 + # | | | | + # 10 4 7 32 + + phys_space = Vect[U1Irrep](1 => 1, -1 => 1) + corner_space = Vect[U1Irrep](0 => 1, 1 => 1, -1 => 1) + vspaces = map(i -> Vect[U1Irrep](0 => 1 + i % 4, 1 => i ÷ 4 % 4, -2 => i ÷ 16), 1:32) + @test length(Set(vspaces)) == 32 + + Espaces = [ + vspaces[1] vspaces[5] vspaces[8] vspaces[3] + vspaces[12] vspaces[15] vspaces[17] vspaces[14] + vspaces[20] vspaces[23] vspaces[25] vspaces[22] + vspaces[28] vspaces[30] vspaces[31] vspaces[29] + ] + + Nspaces = [ + vspaces[10] vspaces[4] vspaces[7] vspaces[32] + vspaces[2] vspaces[6] vspaces[9] vspaces[11] + vspaces[13] vspaces[16] vspaces[18] vspaces[19] + vspaces[21] vspaces[24] vspaces[26] vspaces[27] + ] + Pspaces = fill(phys_space, (4, 4)) + + test_unitcell(unitcell, Pspaces, Nspaces, Nspaces) +end diff --git a/test/runtests.jl b/test/runtests.jl index 532359c72..5fc589315 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -53,6 +53,20 @@ end include("ctmrg/correlation_length.jl") end end + if GROUP == "ALL" || GROUP == "BP" + @time @safetestset "Unit cell bond matching" begin + include("bp/unitcell.jl") + end + @time @safetestset "Expectation values" begin + include("bp/expvals.jl") + end + @time @safetestset "Rotation of BPEnv" begin + include("bp/rotation.jl") + end + @time @safetestset "Gauge-fixing iPEPS" begin + include("bp/gaugefix.jl") + end + end if GROUP == "ALL" || GROUP == "GRADIENTS" @time @safetestset "CTMRG gradients" begin include("gradients/ctmrg_gradients.jl")