diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index 4a39d6bd3..dd7770541 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -77,6 +77,7 @@ 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/fullupdate.jl") include("algorithms/transfermatrix.jl") include("algorithms/toolbox.jl") @@ -105,7 +106,7 @@ export fixedpoint export absorb_weight export ALSTruncation, FullEnvTruncation -export SimpleUpdate +export SimpleUpdate, FullUpdate export TimeEvolver, timestep, time_evolve export InfiniteSquareNetwork diff --git a/src/algorithms/contractions/bondenv/als_solve.jl b/src/algorithms/contractions/bondenv/als_solve.jl index c244c7370..f86bf3112 100644 --- a/src/algorithms/contractions/bondenv/als_solve.jl +++ b/src/algorithms/contractions/bondenv/als_solve.jl @@ -206,3 +206,13 @@ function _solve_ab( x1, info = linsolve(f, Sx, x0, 0, 1) return x1, info end + +function _solve_ab_pinv!( + Rx::AbstractTensorMap{T, S, 2, 2}, Sx::AbstractTensorMap{T, S, 2, 1}; kwargs... + ) where {T <: Number, S <: ElementarySpace} + Rx_inv, ϵ = _pinv!(copy(Rx); kwargs...) + is = filter(i -> isdual(codomain(Rx_inv, i)), 1:numout(Rx_inv)) + x = Rx_inv * Sx + twist!(x, is) + return x, ϵ +end diff --git a/src/algorithms/contractions/bondenv/benv_ctm.jl b/src/algorithms/contractions/bondenv/benv_ctm.jl index a936a74f5..979aa762e 100644 --- a/src/algorithms/contractions/bondenv/benv_ctm.jl +++ b/src/algorithms/contractions/bondenv/benv_ctm.jl @@ -8,7 +8,8 @@ Construct the environment (norm) tensor C4---T3---------T3---C3 r+1 c-1 c c+1 c+2 ``` -where `XX = X' X` and `YY = Y' Y` (stacked together). +where `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). Axis order: `[DX1 DY1; DX0 DY0]`, as in ``` @@ -21,7 +22,9 @@ Axis order: `[DX1 DY1; DX0 DY0]`, as in └---------------------┘ ``` """ -function bondenv_fu(row::Int, col::Int, X::PEPSOrth, Y::PEPSOrth, env::CTMRGEnv) +function bondenv_fu( + row::Int, col::Int, X::PO, Y::PO, env::CTMRGEnv + ) where {PO <: Union{PEPSOrth, PEPOOrth}} Nr, Nc = size(env.corners)[[2, 3]] cm1 = _prev(col, Nc) cp1 = _next(col, Nc) @@ -49,12 +52,23 @@ function bondenv_fu(row::Int, col::Int, X::PEPSOrth, Y::PEPSOrth, env::CTMRGEnv) C4--χ3--T3X---------χ5---------T3Y--χ7---C3 r+1 c-1 c c+1 c+2 =# - @autoopt @tensor benv[DX1 DY1; DX0 DY0] := - c4[χ3 χ1] * t4[χ1 DWX0 DWX1 χ2] * c1[χ2 χ4] * t3X[χ5 DSX0 DSX1 χ3] * - X[DNX0 DX0 DSX0 DWX0] * conj(X[DNX1 DX1 DSX1 DWX1]) * t1X[χ4 DNX0 DNX1 χ6] * - c3[χ9 χ7] * t2[χ10 DEY0 DEY1 χ9] * c2[χ8 χ10] * t3Y[χ7 DSY0 DSY1 χ5] * - Y[DNY0 DEY0 DSY0 DY0] * conj(Y[DNY1 DEY1 DSY1 DY1]) * t1Y[χ6 DNY0 DNY1 χ8] - + benv = nothing + if PO <: PEPSOrth + @autoopt @tensor benv[DX1 DY1; DX0 DY0] := + c4[χ3 χ1] * t4[χ1 DWX0 DWX1 χ2] * c1[χ2 χ4] * t3X[χ5 DSX0 DSX1 χ3] * + X[DNX0 DX0 DSX0 DWX0] * conj(X[DNX1 DX1 DSX1 DWX1]) * t1X[χ4 DNX0 DNX1 χ6] * + c3[χ9 χ7] * t2[χ10 DEY0 DEY1 χ9] * c2[χ8 χ10] * t3Y[χ7 DSY0 DSY1 χ5] * + Y[DNY0 DEY0 DSY0 DY0] * conj(Y[DNY1 DEY1 DSY1 DY1]) * t1Y[χ6 DNY0 DNY1 χ8] + else + # eliminate fermion sign when contracting the remaining physical axis in X, Y + X2 = isdual(space(X, 1)) ? twist(X, 1) : X + Y2 = isdual(space(Y, 1)) ? twist(Y, 1) : Y + @autoopt @tensor benv[DX1 DY1; DX0 DY0] := + c4[χ3 χ1] * t4[χ1 DWX0 DWX1 χ2] * c1[χ2 χ4] * t3X[χ5 DSX0 DSX1 χ3] * + X2[pX DNX0 DX0 DSX0 DWX0] * conj(X[pX DNX1 DX1 DSX1 DWX1]) * t1X[χ4 DNX0 DNX1 χ6] * + c3[χ9 χ7] * t2[χ10 DEY0 DEY1 χ9] * c2[χ8 χ10] * t3Y[χ7 DSY0 DSY1 χ5] * + Y2[pY DNY0 DEY0 DSY0 DY0] * conj(Y[pY DNY1 DEY1 DSY1 DY1]) * t1Y[χ6 DNY0 DNY1 χ8] + end normalize!(benv, Inf) return benv end diff --git a/src/algorithms/contractions/bondenv/gaugefix.jl b/src/algorithms/contractions/bondenv/gaugefix.jl index c49030c13..911fb159e 100644 --- a/src/algorithms/contractions/bondenv/gaugefix.jl +++ b/src/algorithms/contractions/bondenv/gaugefix.jl @@ -92,12 +92,13 @@ function fixgauge_benv( end """ -When the (half) bond environment `Z` consists of two `PEPSOrth` tensors `X`, `Y` as +When the (half) bond environment `Z` consists of +two `PEPSOrth` or `PEPOOrth` tensors `X`, `Y` as ``` - ┌---------------┐ ┌-------------------┐ - | | = | | , - └---Z-- --┘ └--Z0---X-- --Y--┘ - ↓ ↓ + ┌-----------------------┐ + | | + └---Z---(X)-- --(Y)---┘ + ↓ ``` apply the gauge transformation `Linv`, `Rinv` for `Z` to `X`, `Y`: ``` @@ -106,15 +107,25 @@ apply the gauge transformation `Linv`, `Rinv` for `Z` to `X`, `Y`: -4 - X - 1 - Rinv - -2 -4 - Linv - 1 - Y - -2 | | -3 -3 + + -2 -2 + | | + -5 - X - 1 - Rinv - -3 -5 - Linv - 1 - Y - -3 + | ╲ | ╲ + -4 -1 -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} + X::PEPSOrth, Y::PEPSOrth, Linv::MPSBondTensor, Rinv::MPSBondTensor, + ) @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 end +function _fixgauge_benvXY( + X::PEPOOrth, Y::PEPOOrth, Linv::MPSBondTensor, Rinv::MPSBondTensor, + ) + @plansor X[-1 -2 -3 -4 -5] := X[-1 -2 1 -4 -5] * Rinv[1; -3] + @plansor Y[-1 -2 -3 -4 -5] := Y[-1 -2 -3 -4 1] * Linv[1; -5] + return X, Y +end diff --git a/src/algorithms/ctmrg/sequential.jl b/src/algorithms/ctmrg/sequential.jl index 9818f79b2..7a36e0b29 100644 --- a/src/algorithms/ctmrg/sequential.jl +++ b/src/algorithms/ctmrg/sequential.jl @@ -57,6 +57,65 @@ function ctmrg_leftmove(col::Int, network, env::CTMRGEnv, alg::SequentialCTMRG) return env, info end +""" + ctmrg_rightmove(col::Int, network, env::CTMRGEnv, alg::SequentialCTMRG) + +Perform sequential CTMRG right move on the `col`-th column. +""" +function ctmrg_rightmove(col::Int, network, env::CTMRGEnv, alg::SequentialCTMRG) + #= + right move <--- + ←-- T1 ← C2 r-1 + ‖ ↑ + === M' = T2 r + ‖ ↑ + --→ T3 → C3 r+1 + c c+1 + =# + Nc = size(network)[2] + @assert 1 <= col <= Nc + env, info = ctmrg_leftmove(Nc + 1 - col, rot180(network), rot180(env), alg) + return rot180(env), info +end + +""" + ctmrg_upmove(col::Int, network, env::CTMRGEnv, alg::SequentialCTMRG) + +Perform sequential CTMRG up move on the `row`-th row. +""" +function ctmrg_upmove(row::Int, network, env::CTMRGEnv, alg::SequentialCTMRG) + #= + c-1 c c+1 + r-1 C1 ← T1 ← C2 + ↓ ‖ ↑ | up + r T4 = M == T2 ↓ move + ↓ ‖ ↑ + =# + Nr = size(network)[1] + @assert 1 <= row <= Nr + env, info = ctmrg_leftmove(Nr + 1 - row, rotl90(network), rotl90(env), alg) + return rotr90(env), info +end + +""" + ctmrg_downmove(col::Int, network, env::CTMRGEnv, alg::SequentialCTMRG) + +Perform sequential CTMRG down move on the `row`-th row. +""" +function ctmrg_downmove(row::Int, network, env::CTMRGEnv, alg::SequentialCTMRG) + #= + ↓ ‖ ↑ + r T4 = M == T2 ↑ down + ↓ ‖ ↑ | move + r+1 C4 → T3 → C3 + c-1 c c+1 + =# + Nr = size(network)[1] + @assert 1 <= row <= Nr + env, info = ctmrg_leftmove(row, rotr90(network), rotr90(env), alg) + return rotl90(env), info +end + function ctmrg_iteration(network, env::CTMRGEnv, alg::SequentialCTMRG) truncation_error = zero(real(scalartype(network))) condition_number = zero(real(scalartype(network))) diff --git a/src/algorithms/time_evolution/fullupdate.jl b/src/algorithms/time_evolution/fullupdate.jl new file mode 100644 index 000000000..fe959e4a0 --- /dev/null +++ b/src/algorithms/time_evolution/fullupdate.jl @@ -0,0 +1,391 @@ +""" +$(TYPEDEF) + +Algorithm struct for full update (FU) of InfinitePEPS or InfinitePEPO. + +## Fields + +$(TYPEDFIELDS) + +Reference: Physical Review B 92, 035142 (2015) +""" +@kwdef struct FullUpdate <: TimeEvolution + "When true, fix gauge of bond environment" + fixgauge::Bool = true + "When true (or false), the Trotter gate is `exp(-H dt)` (or `exp(-iH dt)`)" + imaginary_time::Bool = true + "Number of iterations without fully reconverging the environment" + reconverge_interval::Int = 5 + "Do column and row moves after updating every single bond, + instaed of waiting until an entire row/column of bonds are updated" + finer_env_update::Bool = false + "Bond truncation algorithm after applying time evolution gate" + opt_alg::Union{ALSTruncation, FullEnvTruncation} = + FullEnvTruncation(; trunc = truncerror(; atol = 1.0e-10)) + "CTMRG algorithm to reconverge environment. + Its `projector_alg` is also used for the fast update of + the environment after each FU iteration" + ctm_alg::SequentialCTMRG = SequentialCTMRG(; + tol = 1.0e-9, maxiter = 20, verbosity = 1, + trunc = truncerror(; atol = 1.0e-10), + projector_alg = :fullinfinite, + ) +end + +""" +Internal state of full update algorithm +""" +struct FUState{S <: InfiniteState, E <: CTMRGEnv, N <: Number} + "number of performed iterations" + iter::Int + "evolved time" + t::N + "PEPS/PEPO" + psi::S + "CTMRG environment" + env::E + "whether the current environment is reconverged" + reconverged::Bool +end + +""" + TimeEvolver( + psi0::InfiniteState, H::LocalOperator, dt::Number, nstep::Int, + alg::FullUpdate, env0::CTMRGEnv; t0::Number = 0.0 + ) + +Initialize a TimeEvolver with Hamiltonian `H` and full update `alg`, +starting from the initial state `psi0` and CTMRG environment `env0`. + +- The initial time is specified by `t0`. +""" +function TimeEvolver( + psi0::InfiniteState, H::LocalOperator, dt::Number, nstep::Int, + alg::FullUpdate, env0::CTMRGEnv; t0::Number = 0.0 + ) + _timeevol_sanity_check(psi0, physicalspace(H), alg) + dt′ = _get_dt(psi0, dt, alg.imaginary_time) + gate = get_expham(H, dt′ / 2) + state = FUState(0, t0, psi0, env0, true) + return TimeEvolver(alg, dt, nstep, gate, state) +end + +""" +Full update for the bond between `[row, col]` and `[row, col+1]`. +""" +function _fu_xbond!( + state::InfiniteState, gate::AbstractTensorMap{T, S, 2, 2}, env::CTMRGEnv, + row::Int, col::Int, alg::FullUpdate + ) where {T <: Number, S <: ElementarySpace} + cp1 = _next(col, size(state, 2)) + A, B = state[row, col], state[row, cp1] + X, a, b, Y = _qr_bond(A, B) + # positive/negative-definite approximant: benv = ± Z Z† + benv = bondenv_fu(row, col, X, Y, env) + Z = positive_approx(benv) + @debug "cond(benv) before gauge fix: $(LinearAlgebra.cond(Z' * Z))" + # fix gauge of bond environment + if alg.fixgauge + Z, a, b, (Linv, Rinv) = fixgauge_benv(Z, a, b) + X, Y = _fixgauge_benvXY(X, Y, Linv, Rinv) + benv = Z' * Z + @debug "cond(L) = $(LinearAlgebra.cond(Linv)); cond(R): $(LinearAlgebra.cond(Rinv))" + @debug "cond(benv) after gauge fix: $(LinearAlgebra.cond(benv))" + else + benv = Z' * Z + end + # apply gate + a, s, b, = _apply_gate(a, b, gate, truncerror(; atol = 1.0e-15)) + # optimize a, b; s is already normalized + a, s, b, info = bond_truncate(a, b, benv, alg.opt_alg) + A, B = _qr_bond_undo(X, a, b, Y) + normalize!(A, Inf) + normalize!(B, Inf) + state.A[row, col], state.A[row, cp1] = A, B + return s, info +end + +""" +Update all horizontal bonds in the c-th column +(i.e. `(r,c) (r,c+1)` for all `r = 1, ..., Nr`). +To update rows, rotate the network clockwise by 90 degrees. +The iPEPS/iPEPO `state` is modified in place. +""" +function _fu_column!( + state::InfiniteState, gate::LocalOperator, + alg::FullUpdate, env::CTMRGEnv, col::Int + ) + Nr, Nc, = size(state) + fid = 1.0 + wts_col = Vector{PEPSWeight}(undef, Nr) + colmove_alg = SequentialCTMRG() + @reset colmove_alg.projector_alg = alg.ctm_alg.projector_alg + env2 = deepcopy(env) + for row′ in 1:Nr + # update row order: 1, Nr, 2, Nr-1, 3, Nr-2, ... + row = isodd(row′) ? ((row′ + 1) ÷ 2) : (Nr - (row′ ÷ 2) + 1) + term = get_gateterm(gate, (CartesianIndex(row, col), CartesianIndex(row, col + 1))) + wts_col[row], info = _fu_xbond!(state, term, env, row, col, alg) + fid = min(fid, info.fid) + if alg.finer_env_update + network = isa(state, InfinitePEPS) ? InfiniteSquareNetwork(state) : + InfiniteSquareNetwork(InfinitePEPS(state)) + # match bonds between updated state tensors and the environment + env2, info = ctmrg_leftmove(col, network, env2, colmove_alg) + env2, info = ctmrg_rightmove(_next(col, Nc), network, env2, colmove_alg) + # update environment around the bond to be updated next + env2, info = ctmrg_upmove(row, network, env2, colmove_alg) + env2, info = ctmrg_downmove(row, network, env2, colmove_alg) + end + end + if !alg.finer_env_update + network = isa(state, InfinitePEPS) ? InfiniteSquareNetwork(state) : + InfiniteSquareNetwork(InfinitePEPS(state)) + env2, info = ctmrg_leftmove(col, network, env, colmove_alg) + env2, info = ctmrg_rightmove(_next(col, Nc), network, env2, colmove_alg) + end + for c in [col, _next(col, Nc)] + env.corners[:, :, c] = env2.corners[:, :, c] + env.edges[:, :, c] = env2.edges[:, :, c] + end + return wts_col, fid +end + +""" +One round of fast full update on the input InfinitePEPS or InfinitePEPO `state` +and its 2-layer CTMRGEnv `env`, without fully reconverging `env`. +""" +function fu_iter( + state::InfiniteState, gate::LocalOperator, alg::FullUpdate, env::CTMRGEnv + ) + Nr, Nc, = size(state) + fidmin = 1.0 + state2, env2 = deepcopy(state), deepcopy(env) + wts = Array{PEPSWeight}(undef, 2, Nr, Nc) + for i in 1:4 + N = size(state2, 2) + for col in 1:N + wts_col, fid_col = _fu_column!(state2, gate, alg, env2, col) + fidmin = min(fidmin, fid_col) + # assign the weights to the un-rotated `wts` + if i == 1 + wts[1, :, col] = wts_col + elseif i == 2 + wts[2, _next(col, N), :] = reverse(wts_col) + elseif i == 3 + wts[1, :, mod1(N - col, N)] = reverse(wts_col) + else + wts[2, N + 1 - col, :] = wts_col + end + end + gate, state2, env2 = rotl90(gate), rotl90(state2), rotl90(env2) + end + return state2, env2, SUWeight(collect(wt for wt in wts)), fidmin +end + +function Base.iterate(it::TimeEvolver{<:FullUpdate}, state = it.state) + iter, t = state.iter, state.t + (iter == it.nstep) && return nothing + psi, env, wts, fid = fu_iter(state.psi, it.gate, it.alg, state.env) + iter, t = iter + 1, t + it.dt + # reconverge environment for the last step and every `reconverge_interval` steps + reconverged = (iter % it.alg.reconverge_interval == 0) || (iter == it.nstep) + if reconverged + network = isa(psi, InfinitePEPS) ? psi : InfinitePEPS(psi) + env, = leading_boundary(env, network, it.alg.ctm_alg) + end + # update internal state + it.state = FUState(iter, t, psi, env, reconverged) + info = (; t, wts, fid) + return (psi, env, info), it.state +end + +""" + timestep( + it::TimeEvolver{<:FullUpdate}, psi::InfiniteState, env::CTMRGEnv; + iter::Int = it.state.iter, t::Float64 = it.state.t + ) -> (psi, env, info) + +Given the TimeEvolver iterator `it`, perform one step of time evolution +on the input state `psi` and its environment `env`. + +- Using `iter` and `t` to reset the current iteration number and evolved time + respectively of the TimeEvolver `it`. +- Use `reconverge_env` to force reconverging the obtained environment. +""" +function MPSKit.timestep( + it::TimeEvolver{<:FullUpdate}, psi::InfiniteState, env::CTMRGEnv; + iter::Int = it.state.iter, t::Float64 = it.state.t, reconverge_env::Bool = false + ) + _timeevol_sanity_check(psi, physicalspace(it.state.psi), it.alg) + state = FUState(iter, t, psi, env, true) + result = iterate(it, state) + if result === nothing + @warn "TimeEvolver `it` has already reached the end." + return nothing + else + psi, env, info = first(result) + if reconverge_env && !(it.state.reconverged) + network = isa(psi, InfinitePEPS) ? psi : InfinitePEPS(psi) + env, = leading_boundary(env, network, it.alg.ctm_alg) + # update internal state of `it` + state0 = it.state + it.state = (@set state0.env = env) + end + return psi, env, info + end +end + +""" +Imaginary time full update of InfinitePEPS with convergence checking +""" +function _time_evolve_gs( + it::TimeEvolver{<:FullUpdate}, tol::Float64, H::LocalOperator + ) + time_start = time() + @assert (it.state.psi isa InfinitePEPS) && it.alg.imaginary_time "Only imaginary time evolution of InfinitePEPS allows convergence checking." + time0 = time() + # backup variables + iter0, t0 = it.state.iter, it.state.t + psi0, env0, info0 = it.state.psi, it.state.env, nothing + energy0 = expectation_value(psi0, H, psi0, env0) / prod(size(psi0)) + @info "FU: initial state energy = $(energy0)." + for (psi, env, info) in it + iter = it.state.iter + if iter == 1 + # reconverge for the 1st step + network = isa(psi, InfinitePEPS) ? psi : InfinitePEPS(psi) + env, = leading_boundary(env, network, it.alg.ctm_alg) + # update internal state of `it` + # TODO: more elegant to use `Accessors.@set` + it.state = FUState(iter, it.state.t, it.state.psi, env, true) + end + !(it.state.reconverged) && continue + # do the following only when env has been reconverged + energy = expectation_value(psi, H, psi, env) / prod(size(psi)) + diff = energy - energy0 + stop = (iter == it.nstep) || (diff < 0 && abs(diff) < tol) || (diff > 0) + showinfo = (iter == 1) || it.state.reconverged || stop + time1 = time() + if showinfo + 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)." + Δλ = (info0 === nothing) ? NaN : compare_weights(info.wts, info0.wts) + @info @sprintf( + "FU iter %-6d: E = %.5f, ΔE = %.3e, |Δλ| = %.3e. Time: %.2f s", + it.state.iter, energy, diff, Δλ, time1 - time0 + ) + end + if (diff < 0 && abs(diff) < tol) + @info "FU: energy has converged." + end + if diff > 0 + @warn "FU: energy has increased from last check. Abort evolution and return results from last check." + psi, env, info, energy = psi0, env0, info0, energy0 + # also reset internal state of `it` to last check + it.state = FUState(iter0, t0, psi0, env0, true) + end + if iter == it.nstep + @info "FU: energy has not converged." + end + if stop + @assert it.state.reconverged + time_end = time() + @info @sprintf("Full update finished. Total time elasped: %.2f s", time_end - time_start) + return psi, env, info + else + # update backup variables + iter0, t0 = it.state.iter, it.state.t + psi0, env0, info0, energy0 = psi, env, info, energy + end + time0 = time() + end + return +end + +""" +Full update without convergence checking +""" +function _time_evolve(it::TimeEvolver{<:FullUpdate}) + time_start = time0 = time() + info0 = nothing + for (psi, env, info) in it + iter = it.state.iter + !(it.state.reconverged) && continue + # do the following only when env has been reconverged + stop = (iter == it.nstep) + showinfo = (iter == 1) || it.state.reconverged || stop + time1 = time() + if showinfo + 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)." + Δλ = (info0 === nothing) ? NaN : compare_weights(info.wts, info0.wts) + @info @sprintf( + "FU iter %d: t = %.2e, |Δλ| = %.3e. Time: %.2f s", + it.state.iter, it.state.t, Δλ, time1 - time0 + ) + end + if stop + @assert it.state.reconverged + time_end = time() + @info @sprintf("Full update finished. Total time elasped: %.2f s", time_end - time_start) + return psi, env, info + else + info0 = info + end + time0 = time() + end + return +end + +""" + time_evolve( + it::TimeEvolver{<:FullUpdate}; + tol::Float64 = 0.0, H::Union{Nothing, LocalOperator} = nothing + ) + +Perform time evolution to the end of FullUpdate TimeEvolver `it`, +or until convergence of energy set by a positive `tol`. + +- Setting `tol > 0` enables convergence check (for imaginary time evolution of InfinitePEPS only). The Hamiltonian `H` should also be provided to measure the energy. + For other usages they should not be changed. +""" +function MPSKit.time_evolve( + it::TimeEvolver{<:FullUpdate}; + tol::Float64 = 0.0, H::Union{Nothing, LocalOperator} = nothing + ) + if tol > 0 + return _time_evolve_gs(it, tol, H) + else + @assert tol == 0 + return _time_evolve(it) + end +end + +""" + MPSKit.time_evolve( + psi0::InfiniteState, H::LocalOperator, dt::Number, nstep::Int, + alg::FullUpdate, env0::CTMRGEnv; + tol::Float64 = 0.0, t0::Number = 0.0 + ) -> (psi, env, info) + +Perform time evolution on the initial state `psi0` and initial environment `env0` +with Hamiltonian `H`, using FullUpdate algorithm `alg`, time step `dt` for +`nstep` number of steps. + +- Setting `tol > 0` enables convergence check (for imaginary time evolution of InfinitePEPS only). + For other usages it should not be changed. +- Use `t0` to specify the initial time of `psi0`. +- `info` is a NamedTuple containing information of the evolution, + including the time `info.t` evolved since `psi0`. +""" +function MPSKit.time_evolve( + psi0::InfiniteState, H::LocalOperator, dt::Number, nstep::Int, + alg::FullUpdate, env0::CTMRGEnv; + tol::Float64 = 0.0, t0::Number = 0.0 + ) + it = TimeEvolver(psi0, H, dt, nstep, alg, env0; t0) + return time_evolve(it; tol, H) +end diff --git a/src/algorithms/truncation/bond_truncation.jl b/src/algorithms/truncation/bond_truncation.jl index b2f54ac1c..f4e3099fe 100644 --- a/src/algorithms/truncation/bond_truncation.jl +++ b/src/algorithms/truncation/bond_truncation.jl @@ -1,7 +1,7 @@ """ $(TYPEDEF) -Algorithm struct for the alternating least square (ALS) optimization of a bond. +Algorithm struct for the alternating least square (ALS) optimization of a bond between two tensors. ## Fields @@ -15,22 +15,25 @@ The truncation algorithm can be constructed from the following keyword arguments * `trunc::TruncationStrategy`: SVD truncation strategy when initilizing the truncated tensors connected by the bond. * `maxiter::Int=50` : Maximal number of ALS iterations. -* `tol::Float64=1e-15` : ALS converges when fidelity change between two iterations is smaller than `tol`. +* `tol::Float64=1e-9` : ALS converges when bond SVD spectrum change between two iterations is smaller than `tol`. +* `use_pinv::Bool=true`: Use pseudo-inverse (instead of `KrylovKit.linsolve`) to solve linear equations in ALS itertions. * `check_interval::Int=0` : Set number of iterations to print information. Output is suppressed when `check_interval <= 0`. """ @kwdef struct ALSTruncation trunc::TruncationStrategy maxiter::Int = 50 - tol::Float64 = 1.0e-15 + tol::Float64 = 1.0e-9 + use_pinv::Bool = true check_interval::Int = 0 end function _als_message( - iter::Int, cost::Float64, fid::Float64, Δcost::Float64, Δfid::Float64, time_elapsed::Float64, + iter::Int, cost::Float64, fid::Float64, Δcost::Float64, + Δfid::Float64, Δs::Float64, time_elapsed::Float64, ) return @sprintf( "%5d, fid = %.8e, Δfid = %.8e, time = %.4f s\n", iter, fid, Δfid, time_elapsed - ) * @sprintf(" cost = %.3e, Δcost/cost0 = %.3e", cost, Δcost) + ) * @sprintf(" cost = %.3e, Δcost/cost0 = %.3e, |Δs| = %.4e.", cost, Δcost, Δs) end """ @@ -65,14 +68,15 @@ function bond_truncate( @assert !isdual(space(a, 2)) @assert !isdual(space(b, 2)) @assert codomain(benv) == domain(benv) + 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, s, b = svd_trunc(permute(a2b2, perm_ab); trunc = alg.trunc) - s /= norm(s, Inf) - a, b = absorb_s(a, s, b) + a, s0, b = svd_trunc(permute(a2b2, perm_ab); trunc = alg.trunc) + normalize!(s0, Inf) + a, b = absorb_s(a, s0, b) #= temporarily reorder axes of a and b to 1 -a/b- 2 ↓ @@ -84,8 +88,8 @@ function bond_truncate( # cost function will be normalized by initial value cost00 = cost_function_als(benv, ab, a2b2) fid = fidelity(benv, ab, a2b2) - cost0, fid0, Δfid = cost00, fid, 0.0 - verbose && @info "ALS init" * _als_message(0, cost0, fid, NaN, NaN, 0.0) + 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() #= @@ -98,25 +102,39 @@ function bond_truncate( =# Ra = _tensor_Ra(benv, b) Sa = _tensor_Sa(benv, b, a2b2) - a, info_a = _solve_ab(Ra, Sa, a) + a, info_a = if alg.use_pinv + _solve_ab_pinv!(Ra, Sa; atol = 1.0e-10) + else + _solve_ab(Ra, Sa, a) + end # Fixing `a`, solve for `b` from `Rb b = Sb` Rb = _tensor_Rb(benv, a) Sb = _tensor_Sb(benv, a, a2b2) - b, info_b = _solve_ab(Rb, Sb, b) + b, info_b = if alg.use_pinv + _solve_ab_pinv!(Rb, Sb; atol = 1.0e-10) + else + _solve_ab(Rb, Sb, b) + end + @debug "Bond truncation info" info_a info_b ab = _combine_ab(a, b) cost = cost_function_als(benv, ab, a2b2) fid = fidelity(benv, ab, a2b2) + _, s, _ = svd_trunc!(permute(ab, perm_ab); trunc = alg.trunc) + normalize!(s, Inf) + # fidelity, cost and bond-s change + Δs = (space(s) == space(s0)) ? _singular_value_distance((s, s0)) : NaN Δcost = abs(cost - cost0) / cost00 Δfid = abs(fid - fid0) - cost0, fid0 = cost, fid + cost0, fid0, s0 = cost, fid, s time1 = time() - converge = (Δfid < alg.tol) + converge = (Δs < alg.tol) cancel = (iter == alg.maxiter) showinfo = cancel || (verbose && (converge || iter == 1 || iter % alg.check_interval == 0)) if showinfo message = _als_message( - iter, cost, fid, Δcost, Δfid, time1 - ((cancel || converge) ? time00 : time0), + iter, cost, fid, Δcost, Δfid, Δs, + time1 - ((cancel || converge) ? time00 : time0), ) if converge @info "ALS conv" * message @@ -131,7 +149,11 @@ function bond_truncate( a, s, b = svd_trunc!(permute(_combine_ab(a, b), perm_ab); trunc = alg.trunc) # normalize singular value spectrum s /= norm(s, Inf) - return a, s, b, (; fid, Δfid) + a, b = absorb_s(a, s, b) + if need_flip + a, s, b = flip_svd(a, s, b) + end + return a, s, b, (; fid, Δfid, Δs) end function bond_truncate( @@ -144,6 +166,7 @@ function bond_truncate( @assert !isdual(space(a, 2)) @assert !isdual(space(b, 2)) @assert codomain(benv) == domain(benv) + need_flip = isdual(space(b, 1)) #= initialize bond matrix using QR as `Ra Lb` --- a == b --- ==> - Qa - Ra == Rb - Qb - @@ -174,8 +197,12 @@ function bond_truncate( ) # 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] + if need_flip + a, s, b = flip_svd(a, s, b) + end return a, s, b, info end diff --git a/src/algorithms/truncation/fullenv_truncation.jl b/src/algorithms/truncation/fullenv_truncation.jl index 6590787de..763d31ec7 100644 --- a/src/algorithms/truncation/fullenv_truncation.jl +++ b/src/algorithms/truncation/fullenv_truncation.jl @@ -15,7 +15,7 @@ The truncation algorithm can be constructed from the following keyword arguments * `trunc::TruncationStrategy` : SVD truncation strategy when optimizing the new bond matrix. * `maxiter::Int=50` : Maximal number of FET iterations. -* `tol::Float64=1e-15` : FET converges when fidelity change between two FET iterations is smaller than `tol`. +* `tol::Float64=1e-9` : FET converges when bond SVD spectrum change between two FET iterations is smaller than `tol`. * `trunc_init::Bool=true` : Controls whether the initialization of the new bond matrix is obtained from truncated SVD of the old bond matrix. * `check_interval::Int=0` : Set number of iterations to print information. Output is suppressed when `check_interval <= 0`. @@ -26,7 +26,7 @@ The truncation algorithm can be constructed from the following keyword arguments @kwdef struct FullEnvTruncation trunc::TruncationStrategy maxiter::Int = 50 - tol::Float64 = 1.0e-15 + tol::Float64 = 1.0e-9 trunc_init::Bool = true check_interval::Int = 0 end @@ -242,7 +242,7 @@ function fullenv_truncate( @tensor B[-1 -2; -3 -4] := conj(u[1; -1]) * benv[1 -2; 3 -4] * u[3; -3] _linearmap_twist!(p) _linearmap_twist!(B) - r, info_r = linsolve(Base.Fix1(*, B), p, r, 0, 1) + r, info_r = linsolve(Base.Fix1(*, B), p, r) @tensor b1[-1; -2] = u[-1; 1] * r[1 -2] u, s, vh = svd_trunc(b1; trunc = alg.trunc) s /= norm(s, Inf) @@ -252,7 +252,8 @@ function fullenv_truncate( @tensor B[-1 -2; -3 -4] := conj(vh[-2; 2]) * benv[-1 2; -3 4] * vh[-4; 4] _linearmap_twist!(p) _linearmap_twist!(B) - l, info_l = linsolve(Base.Fix1(*, B), p, l, 0, 1) + l, info_l = linsolve(Base.Fix1(*, B), p, l) + @debug "Bond truncation info" info_l info_r @tensor b1[-1; -2] = l[-1 1] * vh[1; -2] fid = fidelity(benv, b0, b1) u, s, vh = svd_trunc!(b1; trunc = alg.trunc) @@ -263,7 +264,7 @@ function fullenv_truncate( s0 = deepcopy(s) fid0 = fid time1 = time() - converge = (Δfid < alg.tol) + converge = (Δs < alg.tol) cancel = (iter == alg.maxiter) showinfo = cancel || (verbose && (converge || iter == 1 || iter % alg.check_interval == 0)) diff --git a/src/environments/ctmrg_environments.jl b/src/environments/ctmrg_environments.jl index 08fa91e62..014dd4082 100644 --- a/src/environments/ctmrg_environments.jl +++ b/src/environments/ctmrg_environments.jl @@ -277,6 +277,11 @@ edgetype(::Type{CTMRGEnv{C, E}}) where {C, E} = E TensorKit.spacetype(::Type{E}) where {E <: CTMRGEnv} = spacetype(cornertype(E)) +## (Approximate) equality +function Base.:(==)(env1::CTMRGEnv, env2::CTMRGEnv) + return env1.corners == env2.corners && env1.edges == env2.edges +end + # In-place update of environment function update!(env::CTMRGEnv{C, T}, env´::CTMRGEnv{C, T}) where {C, T} env.corners .= env´.corners diff --git a/src/utility/svd.jl b/src/utility/svd.jl index 27aa8a6c2..f4299e886 100644 --- a/src/utility/svd.jl +++ b/src/utility/svd.jl @@ -628,3 +628,11 @@ function svd_pullback!( end return ΔA end + +# Calculate the pseudo-inverse using SVD +function _pinv!(a::AbstractTensorMap; kwargs...) + # TODO: replace this with actual truncation error once TensorKit is updated + u, s, vh = svd_compact!(a; alg = LAPACK_QRIteration()) + u, s, vh, ϵ = _truncate_compact((u, s, vh), truncerror(; kwargs...)) + return vh' * sdiag_pow(s, -1) * u', ϵ +end diff --git a/test/bondenv/benv_fu.jl b/test/bondenv/benv_fu.jl index 1a531032a..5fc6c90bf 100644 --- a/test/bondenv/benv_fu.jl +++ b/test/bondenv/benv_fu.jl @@ -7,8 +7,12 @@ using Random Random.seed!(100) Nr, Nc = 2, 2 +Envspace = Vect[FermionParity ⊠ U1Irrep]( + (0, 0) => 4, (1, 1 // 2) => 1, (1, -1 // 2) => 1, (0, 1) => 1, (0, -1) => 1 +) +ctm_alg = SequentialCTMRG(; tol = 1.0e-10, verbosity = 2, trunc = truncerror(; atol = 1.0e-10) & truncrank(8)) # create Hubbard iPEPS using simple update -function get_hubbard_state(t::Float64 = 1.0, U::Float64 = 8.0) +function get_hubbard_peps(t::Float64 = 1.0, U::Float64 = 8.0) H = hubbard_model(ComplexF64, Trivial, U1Irrep, InfiniteSquare(Nr, Nc); t, U, mu = U / 2) Vphy = Vect[FermionParity ⊠ U1Irrep]((0, 0) => 2, (1, 1 // 2) => 1, (1, -1 // 2) => 1) peps = InfinitePEPS(rand, ComplexF64, Vphy, Vphy; unitcell = (Nr, Nc)) @@ -20,30 +24,46 @@ function get_hubbard_state(t::Float64 = 1.0, U::Float64 = 8.0) return peps end -peps = get_hubbard_state() -# calculate CTMRG environment -Envspace = Vect[FermionParity ⊠ U1Irrep]( - (0, 0) => 4, (1, 1 // 2) => 1, (1, -1 // 2) => 1, (0, 1) => 1, (0, -1) => 1 -) -ctm_alg = SequentialCTMRG(; tol = 1.0e-10, verbosity = 2, trunc = truncerror(; atol = 1.0e-10) & truncrank(8)) -env, = leading_boundary(CTMRGEnv(rand, ComplexF64, peps, Envspace), peps, ctm_alg) -for row in 1:Nr, col in 1:Nc - cp1 = PEPSKit._next(col, Nc) - A, B = peps.A[row, col], peps.A[row, cp1] - X, a, b, Y = PEPSKit._qr_bond(A, B) - benv = PEPSKit.bondenv_fu(row, col, X, Y, env) - @assert [isdual(space(benv, ax)) for ax in 1:numind(benv)] == [0, 0, 1, 1] - Z = PEPSKit.positive_approx(benv) - # verify that gauge fixing can greatly reduce - # condition number for physical state bond envs - cond1 = cond(Z' * Z) - Z2, a2, b2, (Linv, Rinv) = PEPSKit.fixgauge_benv(Z, a, b) - benv2 = Z2' * Z2 - cond2 = cond(benv2) - @test 1 <= cond2 < cond1 - @info "benv cond number: (gauge-fixed) $(cond2) ≤ $(cond1) (initial)" - # verify gauge fixing is done correctly - @tensor half[:] := Z[-1; 1 3] * a[1; -2 2] * b[2 -3; 3] - @tensor half2[:] := Z2[-1; 1 3] * a2[1; -2 2] * b2[2 -3; 3] - @test half ≈ half2 +function get_hubbard_pepo(t::Float64 = 1.0, U::Float64 = 8.0) + H = hubbard_model(ComplexF64, Trivial, U1Irrep, InfiniteSquare(Nr, Nc); t, U, mu = U / 2) + pepo = PEPSKit.infinite_temperature_density_matrix(H) + wts = SUWeight(pepo) + alg = SimpleUpdate(; + trunc = truncerror(; atol = 1.0e-10) & truncrank(4), bipartite = false + ) + pepo, = time_evolve(pepo, H, 2.0e-3, 500, alg, wts; check_interval = 100) + normalize!.(pepo.A, Inf) + return pepo +end + +function test_benv_fu(state::Union{InfinitePEPS, InfinitePEPO}) + network = isa(state, InfinitePEPS) ? state : InfinitePEPS(state) + env, = leading_boundary(CTMRGEnv(rand, ComplexF64, network, Envspace), network, ctm_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] + X, a, b, Y = PEPSKit._qr_bond(A, B) + benv = PEPSKit.bondenv_fu(row, col, X, Y, env) + @assert [isdual(space(benv, ax)) for ax in 1:numind(benv)] == [0, 0, 1, 1] + Z = PEPSKit.positive_approx(benv) + # verify that gauge fixing can greatly reduce + # condition number for physical state bond envs + cond1 = cond(Z' * Z) + Z2, a2, b2, (Linv, Rinv) = PEPSKit.fixgauge_benv(Z, a, b) + benv2 = Z2' * Z2 + cond2 = cond(benv2) + @test 1 <= cond2 < cond1 + @info "benv cond number: (gauge-fixed) $(cond2) ≤ $(cond1) (initial)" + # verify gauge fixing is done correctly + @tensor half[:] := Z[-1; 1 3] * a[1; -2 2] * b[2 -3; 3] + @tensor half2[:] := Z2[-1; 1 3] * a2[1; -2 2] * b2[2 -3; 3] + @test half ≈ half2 + end + return +end + +peps = get_hubbard_peps() +pepo = get_hubbard_pepo() +for state in (peps, pepo) + test_benv_fu(state) end diff --git a/test/bondenv/bond_truncate.jl b/test/bondenv/bond_truncate.jl index 647e68d6c..5132b8cba 100644 --- a/test/bondenv/bond_truncate.jl +++ b/test/bondenv/bond_truncate.jl @@ -7,7 +7,7 @@ using LinearAlgebra using KrylovKit Random.seed!(0) -maxiter = 500 +maxiter = 600 check_interval = 20 trunc = truncerror(; atol = 1.0e-10) & truncrank(8) Vext = Vect[FermionParity](0 => 100, 1 => 100) @@ -30,15 +30,16 @@ for Vbondl in (Vint, Vint'), Vbondr in (Vint, Vint') @info "Fidelity of simple SVD truncation = $fid0.\n" ss = Dict{String, DiagonalTensorMap}() for (label, alg) in ( - ("ALS", ALSTruncation(; trunc, maxiter, check_interval)), + ("ALS", ALSTruncation(; trunc, maxiter, check_interval, use_pinv = false)), + ("ALS (pinv)", ALSTruncation(; trunc, maxiter, check_interval, use_pinv = true)), ("FET", FullEnvTruncation(; trunc, maxiter, check_interval, trunc_init = false)), ) a1, ss[label], b1, info = PEPSKit.bond_truncate(a2, b2, benv, alg) @info "$label improved fidelity = $(info.fid)." - display(ss[label]) - a1, b1 = PEPSKit.absorb_s(a1, ss[label], b1) + # display(ss[label]) @test info.fid ≈ PEPSKit.fidelity(benv, PEPSKit._combine_ab(a1, b1), a2b2) @test info.fid > fid0 end + @test isapprox(ss["ALS"], ss["ALS (pinv)"], atol = 1.0e-3) @test isapprox(ss["ALS"], ss["FET"], atol = 1.0e-3) end diff --git a/test/runtests.jl b/test/runtests.jl index 38d4fdd4f..9604ffaef 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -70,6 +70,9 @@ end @time @safetestset "Time evolution with site-dependent truncation" begin include("timeevol/sitedep_truncation.jl") end + @time @safetestset "Transverse Field Ising model: real-time full update" 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/tf_ising_realtime.jl b/test/timeevol/tf_ising_realtime.jl new file mode 100644 index 000000000..d59d57f06 --- /dev/null +++ b/test/timeevol/tf_ising_realtime.jl @@ -0,0 +1,102 @@ +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 χ = 4D = 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 +] + +# redefine tfising Hamiltonian with only 2-site gate +function tfising( + T::Type{<:Number}, + S::Union{Type{Trivial}, Type{Z2Irrep}}, + lattice::InfiniteSquare; + J = 1.0, + g = 1.0, + ) + ZZ = rmul!(4 * S_zz(T, S), -J) + X = rmul!(σˣ(T, S), g * -J) + unit = id(space(X, 1)) + gate = ZZ + (1 / 4) * (unit ⊗ X + X ⊗ unit) + spaces = fill(domain(X)[1], (lattice.Nrows, lattice.Ncols)) + return LocalOperator( + spaces, (neighbor => gate for neighbor in PEPSKit.nearest_neighbours(lattice))... + ) +end + +function tfising_fu(g::Float64, Dcut::Int, chi::Int; als = true, use_pinv = true) + # the fully polarized state + peps = InfinitePEPS(zeros, ComplexF64, ℂ^2, ℂ^1; unitcell = (2, 2)) + for t in peps.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 peps.A) + op = LocalOperator(lattice, ((1, 1),) => σˣ()) + ham = tfising(ComplexF64, Trivial, InfiniteSquare(2, 2); J = 1.0, g = g) + + trunc_peps = truncerror(; atol = 1.0e-10) & truncrank(Dcut) + opt_alg = if als + ALSTruncation(; trunc = trunc_peps, tol = 1.0e-10, use_pinv) + else + FullEnvTruncation(; trunc = trunc_peps, tol = 1.0e-10) + end + + 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 + ) + + env = CTMRGEnv(ones, ComplexF64, peps, ℂ^1) + env, = leading_boundary(env, peps, ctm_alg) + + # do one extra step at the beginning to match benchmark data + fu_alg = FullUpdate(; opt_alg, ctm_alg, imaginary_time = false, reconverge_interval = 5) + evolver = TimeEvolver(peps, ham, 0.01, 30, fu_alg, env) + peps, env, info = timestep(evolver, peps, env; reconverge_env = true) + # ensure the recoverged environment is updated to the internal state of `evolver` + @test env == evolver.state.env + magx = expectation_value(peps, op, env) + @info Printf.format(formatter, info.t, real(magx), imag(magx)) + @test isapprox(magx, data[1, 2]; atol = 0.005) + + # reset the number of performed iterations + state0 = evolver.state + evolver.state = (@set state0.iter = 0) + # continue the remaining evolution + count = 2 + for (peps, env, info) in evolver + !evolver.state.reconverged && continue + # 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)." + + 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 + end + return nothing +end + +tfising_fu(hc, 6, 24; als = false, use_pinv = true) diff --git a/test/timeevol/timestep.jl b/test/timeevol/timestep.jl index 39153a898..fd110a526 100644 --- a/test/timeevol/timestep.jl +++ b/test/timeevol/timestep.jl @@ -15,7 +15,7 @@ using PEPSKit evolver = TimeEvolver(ψ0, H, dt, nstep, alg, env0) ψ1, env1, info1 = deepcopy(ψ0), deepcopy(env0), nothing for iter in 0:(nstep - 1) - ψ1, env1, info1 = timestep(evolver, ψ1, env1) + ψ1, env1, info1 = timestep(evolver, ψ1, env1; iter) end # time_evolve ψ2, env2, info2 = time_evolve(ψ0, H, dt, nstep, alg, env0) @@ -31,3 +31,37 @@ using PEPSKit @test env1 == env2 == env3 @test info1 == info2 == info3 end + +@testset "FullUpdate timestep" begin + Nr, Nc = 2, 2 + H = real(heisenberg_XYZ(ComplexF64, Trivial, InfiniteSquare(Nr, Nc); Jx = 1, Jy = 1, Jz = 1)) + ψ0 = PEPSKit.infinite_temperature_density_matrix(H) + env0 = CTMRGEnv(ones, Float64, InfinitePEPS(ψ0), ℂ^1) + opt_alg = FullEnvTruncation(; trunc = truncerror(; atol = 1.0e-10) & truncrank(4)) + ctm_alg = SequentialCTMRG(; + tol = 1.0e-9, maxiter = 20, verbosity = 2, + trunc = truncerror(; atol = 1.0e-10) & truncrank(8), + projector_alg = :fullinfinite, + ) + alg = FullUpdate(; opt_alg, ctm_alg) + dt, nstep = 1.0e-2, 20 + # manual timestep + evolver = TimeEvolver(ψ0, H, dt, nstep, alg, env0) + ψ1, env1, info1 = deepcopy(ψ0), deepcopy(env0), nothing + for iter in 0:(nstep - 1) + ψ1, env1, info1 = timestep(evolver, ψ1, env1) + end + # time_evolve + ψ2, env2, info2 = time_evolve(ψ0, H, dt, nstep, alg, env0) + # for-loop syntax + ## manually reset internal state of evolver + evolver.state = PEPSKit.FUState(0, 0.0, ψ0, env0, true) + ψ3, env3, info3 = nothing, nothing, nothing + for state in evolver + ψ3, env3, info3 = state + end + # results should be *exactly* the same + @test ψ1 == ψ2 == ψ3 + @test env1 == env2 == env3 + @test info1.wts == info2.wts == info3.wts +end