diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index abb2fdc2a..a39e84e81 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -104,7 +104,10 @@ include("algorithms/truncation/truncationschemes.jl") include("algorithms/truncation/fullenv_truncation.jl") include("algorithms/truncation/bond_truncation.jl") -include("algorithms/time_evolution/evoltools.jl") +include("algorithms/time_evolution/trotter_gate.jl") +include("algorithms/time_evolution/trotter_mpo.jl") +include("algorithms/time_evolution/apply_gate.jl") +include("algorithms/time_evolution/apply_mpo.jl") include("algorithms/time_evolution/time_evolve.jl") include("algorithms/time_evolution/simpleupdate.jl") include("algorithms/time_evolution/simpleupdate3site.jl") diff --git a/src/algorithms/contractions/bondenv/gaugefix.jl b/src/algorithms/contractions/bondenv/gaugefix.jl index d30c1b2ce..bfd0f7f01 100644 --- a/src/algorithms/contractions/bondenv/gaugefix.jl +++ b/src/algorithms/contractions/bondenv/gaugefix.jl @@ -61,8 +61,8 @@ function fixgauge_benv( ↓ 2 =# - QL, L = left_orth(permute(Z, ((2, 1), (3,))); positive = true) - QR, R = left_orth(permute(Z, ((3, 1), (2,))); positive = true) + QL, L = left_orth!(permute(Z, ((2, 1), (3,)); copy = true); positive = true) + QR, R = left_orth!(permute(Z, ((3, 1), (2,)); copy = true); positive = true) @debug "cond(L) = $(LinearAlgebra.cond(L)); cond(R) = $(LinearAlgebra.cond(R))" # TODO: find a better way to fix gauge that avoids `inv` Linv, Rinv = inv(L), inv(R) diff --git a/src/algorithms/time_evolution/apply_gate.jl b/src/algorithms/time_evolution/apply_gate.jl new file mode 100644 index 000000000..08981666a --- /dev/null +++ b/src/algorithms/time_evolution/apply_gate.jl @@ -0,0 +1,128 @@ +""" +$(SIGNATURES) + +Use QR decomposition on two tensors `A`, `B` connected by a bond to get the reduced tensors. +When `A`, `B` are PEPSTensors, +``` + 2 1 1 + | | | + 5 -A/B- 3 ====> 4 - X ← 2 1 ← a - 3 1 - b → 3 4 → Y - 2 + | ↘ | ↘ ↘ | + 4 1 3 2 2 3 +``` +When `A`, `B` are PEPOTensors, +- If `gate_ax = 1` +``` + 2 3 1 2 1 2 + ↘ | ↘ | ↘ | + 6 -A/B- 4 ====> 5 - X ← 3 1 ← a - 3 1 - b → 3 5 → Y - 3 + | ↘ | ↘ ↘ | + 5 1 4 2 2 4 +``` +- If `gate_ax = 2` +``` + 2 3 2 2 2 2 + ↘ | | ↘ ↘ | + 6 -A/B- 4 ====> 5 - X ← 3 1 ← a - 3 1 - b → 3 5 → Y - 3 + | ↘ | ↘ | ↘ + 5 1 4 1 4 1 +``` +""" +function _qr_bond(A::PT, B::PT; gate_ax::Int = 1) where {PT <: Union{PEPSTensor, PEPOTensor}} + @assert 1 <= gate_ax <= numout(A) + permA, permB, permX, permY = if A isa PEPSTensor + ((2, 4, 5), (1, 3)), ((2, 3, 4), (1, 5)), (1, 4, 2, 3), Tuple(1:4) + else + if gate_ax == 1 + ((2, 3, 5, 6), (1, 4)), ((2, 3, 4, 5), (1, 6)), (1, 2, 5, 3, 4), Tuple(1:5) + else + ((1, 3, 5, 6), (2, 4)), ((1, 3, 4, 5), (2, 6)), (1, 2, 5, 3, 4), Tuple(1:5) + end + end + X, a = left_orth!(permute(A, permA; copy = true); positive = true) + Y, b = left_orth!(permute(B, permB; copy = true); positive = true) + X, Y = permute(X, permX), permute(Y, permY) + b = permute(b, ((3, 2), (1,))) + return X, a, b, Y +end + +""" +$(SIGNATURES) + +Reconstruct the tensors connected by a bond from their `_qr_bond` results. +For PEPSTensors, +``` + -2 -2 + | | + -5- X - 1 - a - -3 -5 - b - 1 - Y - -3 + | ↘ ↘ | + -4 -1 -1 -4 +``` +For PEPOTensors +``` + -2 -3 -2 -3 + ↘ | ↘ | + -6- X - 1 - a - -4 -6 - b - 1 - Y - -4 + | ↘ ↘ | + -5 -1 -1 -5 + + -3 -2 -2 -3 + | ↘ ↘ | + -6- X - 1 - a - -4 -6 - b - 1 - Y - -4 + | ↘ | ↘ + -5 -1 -5 -1 +``` +""" +function _qr_bond_undo(X::PEPSOrth, a::AbstractTensorMap, b::AbstractTensorMap, Y::PEPSOrth) + @tensor A[-1; -2 -3 -4 -5] := X[-2 1 -4 -5] * a[1 -1 -3] + @tensor B[-1; -2 -3 -4 -5] := b[-5 -1 1] * Y[-2 -3 -4 1] + return A, B +end +function _qr_bond_undo(X::PEPOOrth, a::AbstractTensorMap, b::AbstractTensorMap, Y::PEPOOrth) + if !isdual(space(a, 2)) + @tensor A[-1 -2; -3 -4 -5 -6] := X[-2 -3 1 -5 -6] * a[1 -1 -4] + @tensor B[-1 -2; -3 -4 -5 -6] := b[-6 -1 1] * Y[-2 -3 -4 -5 1] + else + @tensor A[-1 -2; -3 -4 -5 -6] := X[-1 -3 1 -5 -6] * a[1 -2 -4] + @tensor B[-1 -2; -3 -4 -5 -6] := b[-6 -2 1] * Y[-1 -3 -4 -5 1] + end + return A, B +end + +""" +$(SIGNATURES) + +Apply 2-site `gate` on the reduced matrices `a`, `b` +``` + -1← a --- 3 --- b ← -4 -2 -3 + ↓ ↓ ↓ ↓ + 1 2 |----gate---| + ↓ ↓ or ↓ ↓ + |----gate---| 1 2 + ↓ ↓ ↓ ↓ + -2 -3 -1← a --- 3 --- b ← -4 +``` +""" +function _apply_gate( + a::AbstractTensorMap, b::AbstractTensorMap, + gate::AbstractTensorMap{T, S, 2, 2}, trunc::TruncationStrategy + ) where {T <: Number, S <: ElementarySpace} + V = space(b, 1) + need_flip = isdual(V) + if isdual(space(a, 2)) + @tensor a2b2[-1 -2; -3 -4] := gate[1 2; -2 -3] * a[-1 1 3] * b[3 2 -4] + else + @tensor a2b2[-1 -2; -3 -4] := gate[-2 -3; 1 2] * a[-1 1 3] * b[3 2 -4] + end + trunc = if trunc isa FixedSpaceTruncation + need_flip ? truncspace(flip(V)) : truncspace(V) + else + trunc + end + a, s, b, ϵ = svd_trunc!(a2b2; trunc, alg = LAPACK_QRIteration()) + a, b = absorb_s(a, s, b) + if need_flip + a, s, b = flip(a, numind(a)), _fliptwist_s(s), flip(b, 1) + end + return a, s, b, ϵ +end diff --git a/src/algorithms/time_evolution/apply_mpo.jl b/src/algorithms/time_evolution/apply_mpo.jl new file mode 100644 index 000000000..963540eb3 --- /dev/null +++ b/src/algorithms/time_evolution/apply_mpo.jl @@ -0,0 +1,395 @@ +#= +# Mixed canonical form of an open boundary MPS +``` + |ψ⟩ = M[1]-←-...-←-M[N] + ↓ ↓ +``` +For convenience, assume all virtual arrows are ←. + +We perform QR and LQ decompositions: starting from +``` + M[1]-←- = Qa[1]-←-R[1]-←- + ↓ ↓ + + -←-M[N] = -←-L[N-1]-←-Qb[N] + ↓ ↓ +``` +we successively calculate +``` + -←-R[n-1]-←-M[n]-←- = -←-Qa[n]-←-R[n]--←-- (n = 2, ..., N - 1) + ↓ ↓ + + -←-M[n+1]-←-L[n+1]-←- = -←-L[n]-←-Qb[n+1]-←- (n = N - 2, ..., 1) + ↓ ↓ +``` + +For each bond (n = 1, ..., N - 1), we perform SVD +``` + R[n] L[n] = U[n] s[n] V†[n] (n = 1, ..., N - 1) +``` +Then we define the projectors together with the Schmidt weight +``` + -←-Pa[n]-←- = L[n] V[n]-←-(1/√s[n])-←- + -←-Pb[n]-←- = -←-(1/√s[n])-←-U†[n] R[n] +``` +The product `Pa Pb` is the identity operator: +``` + Pa[n]-←-Pb[n] = L[n] (R[n] L[n])⁻¹ R[n] = 1 +``` + +The canonical form is then defined by +``` + -←-M̃[n]-←- = -←-Pb[n-1]-←-M[n]-←-Pa[n]-←- + ↓ ↓ +``` + +Note that +``` + M̃[n] + = 1/√s[n-1] U†[n-1] (R[n-1] M[n]) L[n] V[n] 1/√s[n] + = 1/√s[n-1] U†[n-1] Qa[n] (R[n] L[n]) V[n] 1/√s[n] + = 1/√s[n-1] U†[n-1] Qa[n] U[n] s[n] (V†[n] V[n]) 1/√s[n] + = 1/√s[n-1] U†[n-1] Qa[n] U[n] √s[n] +``` +Then `M̃[n]` (n = 1, ..., N - 1) satisfies the (generalized) left-orthogonal condition +``` + ┌---←--M̃[n]--←- ┌-←- 2 + | | | + s[n-1] ↓ = s[n] (s[0] = 1) + | | | + └---→--M̃†[n]-→- └-→- 1 +``` +Similarly, we can express M̃ using Qb +``` + M̃[n] + = 1/√s[n-1] U†[n-1] R[n-1] (M[n] L[n]) V[n] 1/√s[n] + = 1/√s[n-1] U†[n-1] (R[n-1] L[n-1]) Qb[n] V[n] 1/√s[n] + = 1/√s[n-1] U†[n-1] (R[n-1] L[n-1]) Qb[n] V[n] 1/√s[n] + = 1/√s[n-1] (U†[n-1] U[n-1]) s[n-1] V†[n-1] Qb[n] V[n] 1/√s[n] + = √s[n-1] V†[n-1] Qb[n] V[n] 1/√s[n] +``` +Then `M̃[n]` (n = 2, ..., N) satisfies the (generalized) right-orthogonal condition +``` + -←-M̃[n]-←┐ 1 -←-┐ + ↓ | | + * s[n] = s[n-1] (s[N] = 1) + ↓ | | + -→M̃†[n]-→┘ 2 -→-┘ +``` +Here `-*-` is the twist on the physical axis. + +# Truncation of a bond on OBC-MPS + +Suppose we want to truncate the bond between +the n-th and the (n+1)-th sites such that the truncated state +``` + |ψ̃⟩ = M[1]-←-...-←-M̃[n]-←-M̃[n+1]-←-...-←-M[N] + ↓ ↓ ↓ ↓ +``` +maximizes the fidelity +``` + ⟨ψ|ψ̃⟩ ⟨ψ̃|ψ⟩ + F(ψ̃) = ------------- + ⟨ψ̃|ψ̃⟩ ⟨ψ|ψ⟩ +``` +This is simply done by using a truncated `Pa[n], Pb[n]`; then +``` + M[1]--...-M̃[n]--M̃[n+1]--...-M[N] + ⟨ψ|ψ̃⟩ = | | | | + M†[1]-...-M†[n]-M†[n+1]-...-M†[N] + + ┌----M̃[n]--M̃[n+1]--┐ + = s[n-1] | | s[n+1] + └----M†[n]-M†[n+1]-┘ + + ┌--s̃[n]--┐ + = | | = s[n].s̃[n] + └--s[n]--┘ +``` +Then the fidelity is just +``` + F(ψ̃) = (norm(s̃[n], 2) / norm(s[n], 2))^2 +``` +=# +""" +Perform QR decomposition through a PEPS tensor +``` + ╱ ╱ + -←-R0-←-M-←- => ---Q-←-R1-←- + ╱ | ╱ | +``` +""" +function qr_through( + R0::MPSBondTensor, M::GenericMPSTensor{S, 4}; normalize::Bool = true + ) where {S <: ElementarySpace} + @assert !isdual(codomain(R0, 1)) + @assert !isdual(domain(M, 1)) && !isdual(codomain(M, 1)) + @tensor A[-1 -2 -3 -4; -5] := R0[-1; 1] * M[1 -2 -3 -4; -5] + _, r = left_orth!(A; positive = true) + normalize && normalize!(r, Inf) + return r +end +# for `M` at the left end of the MPS +function qr_through( + ::Nothing, M::GenericMPSTensor{S, 4}; normalize::Bool = true + ) where {S <: ElementarySpace} + @assert !isdual(domain(M, 1)) + _, r = left_orth(M; positive = true) + normalize && normalize!(r, Inf) + return r +end + +""" +Perform LQ decomposition through a tensor +``` + ╱ ╱ + -←-L0-←-Q-←- <= -←-M-←-L1-←- + ╱ | ╱ | +``` +""" +function lq_through( + M::GenericMPSTensor{S, 4}, L1::MPSBondTensor; normalize::Bool = true + ) where {S <: ElementarySpace} + @assert !isdual(domain(L1, 1)) + @assert !isdual(codomain(M, 1)) && !isdual(domain(M, 1)) + @tensor A[-1; -2 -3 -4 -5] := M[-1 -2 -3 -4; 1] * L1[1; -5] + l, _ = right_orth!(A; positive = true) + normalize && normalize!(l, Inf) + return l +end +# for `M` at the right end of the MPS +function lq_through( + M::GenericMPSTensor{S, 4}, ::Nothing; normalize::Bool = true + ) where {S <: ElementarySpace} + @assert !isdual(codomain(M, 1)) + A = permute(M, ((1,), (2, 3, 4, 5))) + l, _ = right_orth!(A; positive = true) + normalize && normalize!(l, Inf) + return l +end + +""" +Given a cluster `Ms`, find all `R`, `L` matrices on each internal bond +""" +function _get_allRLs(Ms::Vector{T}) where {T <: GenericMPSTensor{<:ElementarySpace, 4}} + # M1 -- (R1,L1) -- M2 -- (R2,L2) -- M3 + N = length(Ms) + # get the first R and the last L + R_first = qr_through(nothing, Ms[1]; normalize = true) + L_last = lq_through(Ms[N], nothing; normalize = true) + Rs = Vector{typeof(R_first)}(undef, N - 1) + Ls = Vector{typeof(L_last)}(undef, N - 1) + Rs[1], Ls[end] = R_first, L_last + # get remaining R, L matrices + for n in 2:(N - 1) + m = N - n + 1 + Rs[n] = qr_through(Rs[n - 1], Ms[n]; normalize = true) + Ls[m - 1] = lq_through(Ms[m], Ls[m]; normalize = true) + end + return Rs, Ls +end + +""" +Given the tensors `R`, `L` on a bond, construct +the projectors `Pa`, `Pb` and the new bond weight `s` +such that the contraction of `Pa`, `s`, `Pb` is identity when `trunc = notrunc`, + +The arrows between `Pa`, `s`, `Pb` are +``` + - Pa --←-- Pb - + 1 ← s ← 2 +``` +""" +function _proj_from_RL( + r::MPSBondTensor, l::MPSBondTensor; + trunc::TruncationStrategy = notrunc() + ) + @assert isdual(domain(r, 1)) == isdual(codomain(r, 1)) == false + @assert isdual(domain(l, 1)) == isdual(codomain(l, 1)) == false + rl = r * l + u, s, vh, ϵ = svd_trunc!(rl; trunc) + sinv = sdiag_pow(s, -1 / 2) + Pa, Pb = l * vh' * sinv, sinv * u' * r + return Pa, s, Pb, ϵ +end + +""" +Given a cluster `Ms` and the pre-calculated `R`, `L` bond matrices, +find all projectors `Pa`, `Pb` and Schmidt weights `wts` on internal bonds. +""" +function _get_allprojs(Ms, Rs, Ls, truncs::Vector{E}) where {E <: TruncationStrategy} + N = length(Ms) + @assert length(truncs) == N - 1 + projs_errs = map(1:(N - 1)) do i + trunc = if isa(truncs[i], FixedSpaceTruncation) + tspace = space(Ms[i + 1], 1) + isdual(tspace) ? truncspace(flip(tspace)) : truncspace(tspace) + else + truncs[i] + end + return _proj_from_RL(Rs[i], Ls[i]; trunc) + end + Pas = map(Base.Fix2(getindex, 1), projs_errs) + wts = map(Base.Fix2(getindex, 2), projs_errs) + Pbs = map(Base.Fix2(getindex, 3), projs_errs) + # local truncation error on each bond + ϵs = map(Base.Fix2(getindex, 4), projs_errs) + return Pas, Pbs, wts, ϵs +end + +""" +Flip the virtual arrows in the MPS `Ms` +""" +function _flip_virtuals!( + Ms::Vector{T}, flips::Vector{Bool}; inv::Bool = false + ) where {T <: GenericMPSTensor} + @assert length(flips) == length(Ms) - 1 + for (n, flip) in enumerate(flips) + !flip && continue + M1, M2 = Ms[n], Ms[n + 1] + Ms[n] = TensorKit.flip(M1, numind(M1); inv) + Ms[n + 1] = TensorKit.flip(M2, 1; inv) + end + return Ms +end + +""" +Find projectors to truncate internal bonds of the cluster `Ms`. +""" +function _cluster_truncate!( + Ms::Vector{T}, truncs::Vector{E} + ) where {T <: GenericMPSTensor{<:ElementarySpace, 4}, E <: TruncationStrategy} + Rs, Ls = _get_allRLs(Ms) + Pas, Pbs, wts, ϵs = _get_allprojs(Ms, Rs, Ls, truncs) + # apply projectors + # M1 -- (Pa1,wt1,Pb1) -- M2 -- (Pa2,wt2,Pb2) -- M3 + for (i, (Pa, Pb)) in enumerate(zip(Pas, Pbs)) + @tensor (Ms[i])[-1 -2 -3 -4; -5] := (Ms[i])[-1 -2 -3 -4; 1] * Pa[1; -5] + @tensor (Ms[i + 1])[-1 -2 -3 -4; -5] := Pb[-1; 1] * (Ms[i + 1])[1 -2 -3 -4; -5] + end + return wts, ϵs, Pas, Pbs +end + +""" +Apply the gate MPO `gs` on the cluster `Ms`. +When `gate_ax` is 1 or 2, the gate acts from the physical codomain or domain side. + +e.g. Cluster in PEPS with `gate_ax = 1`: +``` + ╱ ╱ ╱ + --- M1 -←-- M2 -←-- M3 --- + ╱ | ╱ | ╱ | + ↓ ↓ ↓ + g1 -←-- g2 -←-- g3 + ↓ ↓ ↓ +``` + +In the cluster, the axes of each tensor use the MPS order +``` + PEPS: PEPO: + 3 3 4 + ╱ | ╱ + 1 -- M -- 5 1 -- M -- 6 + ╱ | ╱ | + 4 2 5 2 + M[1 2 3 4; 5] M[1 2 3 4 5; 6] +``` +""" +function _apply_gatempo!( + Ms::Vector{T1}, gs::Vector{T2}; gate_ax::Int = 1 + ) where {T1 <: GenericMPSTensor{<:ElementarySpace, 4}, T2 <: AbstractTensorMap} + @assert length(Ms) == length(gs) + @assert gate_ax == 1 + @assert all(!isdual(space(g, 1)) for g in gs[2:end]) + @assert all(!isdual(space(M, 1)) for M in Ms[2:end]) + # fusers to merge axes on bonds in the gate-cluster product + # M1 == f1† -- f1 == M2 == f2† -- f2 == M3 + fusers = map(@view(Ms[2:end]), @view(gs[2:end])) do M, g + V1, V2 = space(M, 1), space(g, 1) + return isomorphism(fuse(V1, V2) ← V1 ⊗ V2) + end + #= gate on codomain of PEPS + -3 -3 -3 + ╱ ┌-┐ ┌-┐ ╱ ┌-┐ ┌-┐ ╱ + -1 --M--2--┤ | | ├--2--M--4--┤ | | ├--2--M-- -5 + ╱ | | ├- -5 -1 -┤ | ╱ | | ├- -5 -1 -┤ | ╱ | + -4 1 | | | | -4 1 | | | | -4 1 + ├--3--┤ | | ├--3--┼--5--┤ | | ├--3--┤ + -2 └-┘ └-┘ -2 └-┘ └-┘ -2 + =# + for (i, (g, M)) in enumerate(zip(gs, Ms)) + @assert !isdual(space(M, 2)) + if i == 1 + fr = fusers[i] + @tensor (Ms[i])[-1 -2 -3 -4; -5] := M[-1 1 -3 -4; 2] * g[-2 1 3] * fr'[2 3; -5] + elseif i == length(Ms) + fl = fusers[i - 1] + @tensor (Ms[i])[-1 -2 -3 -4; -5] := fl[-1; 2 3] * M[2 1 -3 -4; -5] * g[3 -2 1] + else + fl, fr = fusers[i - 1], fusers[i] + @tensor (Ms[i])[-1 -2 -3 -4; -5] := fl[-1; 2 3] * M[2 1 -3 -4; 4] * g[3 -2 1 5] * fr'[4 5; -5] + end + end + return Ms +end + +function _apply_gatempo!( + Ms::Vector{T1}, gs::Vector{T2}; gate_ax::Int = 1 + ) where {T1 <: GenericMPSTensor{<:ElementarySpace, 5}, T2 <: AbstractTensorMap} + @assert length(Ms) == length(gs) + @assert gate_ax == 1 || gate_ax == 2 + @assert all(!isdual(space(g, 1)) for g in gs[2:end]) + @assert all(!isdual(space(M, 1)) for M in Ms[2:end]) + # fusers to merge axes on bonds in the gate-cluster product + # M1 == f1† -- f1 == M2 == f2† -- f2 == M3 + fusers = map(@view(Ms[2:end]), @view(gs[2:end])) do M, g + V1, V2 = space(M, 1), space(g, 1) + return isomorphism(fuse(V1, V2) ← V1 ⊗ V2) + end + #= gate on codomain of PEPO (gate_ax = 1) + + -3 -4 -3 -4 -3 -4 + | ╱ ┌-┐ ┌-┐ | ╱ ┌-┐ ┌-┐ | ╱ + -1 --M--2--┤ | | ├--2--M--4--┤ | | ├--2--M-- -6 + ╱ | | ├- -6 -1 -┤ | ╱ | | ├- -6 -1 -┤ | ╱ | + -5 1 | | | | -5 1 | | | | -5 1 + ├--3--┤ | | ├--3--┼--5--┤ | | ├--3--┤ + -2 └-┘ └-┘ -2 └-┘ └-┘ -2 + + gate on domain of PEPO (gate_ax = 2) + + -3 ┌-┐ ┌-┐ -3 ┌-┐ ┌-┐ -3 + ├--3--┤ | | ├--3--┼--5--┤ | | ├--3--┤ + 1 -4 | ├- -6 -1 -┤ | 1 -4 | ├- -6 -1 -┤ | 1 -4 + | ╱ | | | | | ╱ | | | | | ╱ + -1 --M--2--┤ | | ├--2--M--4--┤ | | ├--2--M-- -6 + ╱ | └-┘ └-┘ ╱ | └-┘ └-┘ ╱ | + -5 -2 -5 -2 -5 -2 + =# + for (i, (g, M)) in enumerate(zip(gs, Ms)) + @assert !isdual(space(M, 2)) + if i == 1 + fr = fusers[i] + if gate_ax == 1 + @tensor (Ms[i])[-1 -2 -3 -4 -5; -6] := M[-1 1 -3 -4 -5; 2] * g[-2 1 3] * fr'[2 3; -6] + else + @tensor (Ms[i])[-1 -2 -3 -4 -5; -6] := M[-1 -2 1 -4 -5; 2] * g[1 -3 3] * fr'[2 3; -6] + end + elseif i == length(Ms) + fl = fusers[i - 1] + if gate_ax == 1 + @tensor (Ms[i])[-1 -2 -3 -4 -5; -6] := fl[-1; 2 3] * M[2 1 -3 -4 -5; -6] * g[3 -2 1] + else + @tensor (Ms[i])[-1 -2 -3 -4 -5; -6] := fl[-1; 2 3] * M[2 -2 1 -4 -5; -6] * g[3 1 -3] + end + else + fl, fr = fusers[i - 1], fusers[i] + if gate_ax == 1 + @tensor (Ms[i])[-1 -2 -3 -4 -5; -6] := fl[-1; 2 3] * M[2 1 -3 -4 -5; 4] * g[3 -2 1 5] * fr'[4 5; -6] + else + @tensor (Ms[i])[-1 -2 -3 -4 -5; -6] := fl[-1; 2 3] * M[2 -2 1 -4 -5; 4] * g[3 1 -3 5] * fr'[4 5; -6] + end + end + end + return Ms +end diff --git a/src/algorithms/time_evolution/evoltools.jl b/src/algorithms/time_evolution/evoltools.jl deleted file mode 100644 index 8c9d34945..000000000 --- a/src/algorithms/time_evolution/evoltools.jl +++ /dev/null @@ -1,302 +0,0 @@ -""" -Process the Trotter time step `dt` according to the intended usage. -""" -function _get_dt( - state::InfiniteState, dt::Number, imaginary_time::Bool - ) - # PEPS update: exp(-H dt)|ψ⟩ - # PEPO update (purified): exp(-H dt/2)|ρ⟩ - # PEPO update (not purified): exp(-H dt/2) ρ exp(-H dt/2) - dt′ = (state isa InfinitePEPS) ? dt : (dt / 2) - if (state isa InfinitePEPO) - @assert size(state)[3] == 1 - end - if !imaginary_time - @assert (state isa InfinitePEPS) "Real time evolution of InfinitePEPO (Heisenberg picture) is not implemented." - dt′ = 1.0im * dt′ - end - return dt′ -end - -function MPSKit.infinite_temperature_density_matrix(H::LocalOperator) - T = scalartype(H) - A = map(physicalspace(H)) do Vp - ψ = permute(TensorKit.id(T, Vp), (1, 2)) - Vv = oneunit(Vp) # trivial (1D) virtual space - virt = ones(T, domain(ψ) ← Vv ⊗ Vv ⊗ Vv' ⊗ Vv') - return ψ * virt - end - return InfinitePEPO(cat(A; dims = 3)) -end - -""" - get_expham(H::LocalOperator, dt::Number) - -Compute `exp(-dt * op)` for each term `op` in `H`, -and combine them into a new LocalOperator. -Each `op` in `H` must be a single `TensorMap`. -""" -function get_expham(H::LocalOperator, dt::Number) - return LocalOperator( - physicalspace(H), (sites => exp(-dt * op) for (sites, op) in H.terms)... - ) -end - -""" - is_nearest_neighbour(H::LocalOperator) - -Check if an operator `H` contains only nearest neighbor terms. -""" -function is_nearest_neighbour(H::LocalOperator) - return all(H.terms) do (sites, op) - return numin(op) == 2 && sum(abs, Tuple(sites[2] - sites[1])) == 1 - end -end - -""" - is_equivalent_bond(bond1::NTuple{2,CartesianIndex{2}}, bond2::NTuple{2,CartesianIndex{2}}, (Nrow, Ncol)::NTuple{2,Int}) - -Check if two 2-site bonds are related by a (periodic) lattice translation. -""" -function is_equivalent_bond( - bond1::NTuple{2, CartesianIndex{2}}, bond2::NTuple{2, CartesianIndex{2}}, - (Nrow, Ncol)::NTuple{2, Int}, - ) - r1 = bond1[1] - bond1[2] - r2 = bond2[1] - bond2[2] - shift_row = bond1[1][1] - bond2[1][1] - shift_col = bond1[1][2] - bond2[1][2] - return r1 == r2 && mod(shift_row, Nrow) == 0 && mod(shift_col, Ncol) == 0 -end - -""" - get_gateterm(gate::LocalOperator, bond::NTuple{2,CartesianIndex{2}}) - -Get the term of a 2-site gate acting on a certain bond. -Input `gate` should only include one term for each nearest neighbor bond. -""" -function get_gateterm(gate::LocalOperator, bond::NTuple{2, CartesianIndex{2}}) - bonds = findall(p -> is_equivalent_bond(p.first, bond, size(gate.lattice)), gate.terms) - if length(bonds) == 0 - # try reversed site order - bonds = findall( - p -> is_equivalent_bond(p.first, reverse(bond), size(gate.lattice)), gate.terms - ) - if length(bonds) == 1 - return permute(gate.terms[bonds[1]].second, ((2, 1), (4, 3))) - elseif length(bonds) == 0 - # if term not found, return the zero operator on this bond - dtype = scalartype(gate) - r1, c1 = (mod1(bond[1][i], n) for (i, n) in zip(1:2, size(gate))) - r2, c2 = (mod1(bond[2][i], n) for (i, n) in zip(1:2, size(gate))) - V1 = physicalspace(gate)[r1, c1] - V2 = physicalspace(gate)[r2, c2] - return zeros(dtype, V1 ⊗ V2 ← V1 ⊗ V2) - else - error("There are multiple terms in `gate` corresponding to the bond $(bond).") - end - else - (length(bonds) == 1) || - error("There are multiple terms in `gate` corresponding to the bond $(bond).") - return gate.terms[bonds[1]].second - end -end - -""" -$(SIGNATURES) - -Use QR decomposition on two tensors `A`, `B` connected by a bond to get the reduced tensors. -When `A`, `B` are PEPSTensors, -``` - 2 1 1 - | | | - 5 -A/B- 3 ====> 4 - X ← 2 1 ← a - 3 1 - b → 3 4 → Y - 2 - | ↘ | ↘ ↘ | - 4 1 3 2 2 3 -``` -When `A`, `B` are PEPOTensors, -- If `gate_ax = 1` -``` - 2 3 1 2 1 2 - ↘ | ↘ | ↘ | - 6 -A/B- 4 ====> 5 - X ← 3 1 ← a - 3 1 - b → 3 5 → Y - 3 - | ↘ | ↘ ↘ | - 5 1 4 2 2 4 -``` -- If `gate_ax = 2` -``` - 2 3 2 2 2 2 - ↘ | | ↘ ↘ | - 6 -A/B- 4 ====> 5 - X ← 3 1 ← a - 3 1 - b → 3 5 → Y - 3 - | ↘ | ↘ | ↘ - 5 1 4 1 4 1 -``` -""" -function _qr_bond(A::PT, B::PT; gate_ax::Int = 1) where {PT <: Union{PEPSTensor, PEPOTensor}} - @assert 1 <= gate_ax <= numout(A) - permA, permB, permX, permY = if A isa PEPSTensor - ((2, 4, 5), (1, 3)), ((2, 3, 4), (1, 5)), (1, 4, 2, 3), Tuple(1:4) - else - if gate_ax == 1 - ((2, 3, 5, 6), (1, 4)), ((2, 3, 4, 5), (1, 6)), (1, 2, 5, 3, 4), Tuple(1:5) - else - ((1, 3, 5, 6), (2, 4)), ((1, 3, 4, 5), (2, 6)), (1, 2, 5, 3, 4), Tuple(1:5) - end - end - X, a = left_orth(permute(A, permA); positive = true) - Y, b = left_orth(permute(B, permB); positive = true) - # no longer needed after TensorKit 0.15 - # @assert !isdual(space(a, 1)) - # @assert !isdual(space(b, 1)) - X, Y = permute(X, permX), permute(Y, permY) - b = permute(b, ((3, 2), (1,))) - return X, a, b, Y -end - -""" -$(SIGNATURES) - -Reconstruct the tensors connected by a bond from their `_qr_bond` results. -For PEPSTensors, -``` - -2 -2 - | | - -5- X - 1 - a - -3 -5 - b - 1 - Y - -3 - | ↘ ↘ | - -4 -1 -1 -4 -``` -For PEPOTensors -``` - -2 -3 -2 -3 - ↘ | ↘ | - -6- X - 1 - a - -4 -6 - b - 1 - Y - -4 - | ↘ ↘ | - -5 -1 -1 -5 - - -3 -2 -2 -3 - | ↘ ↘ | - -6- X - 1 - a - -4 -6 - b - 1 - Y - -4 - | ↘ | ↘ - -5 -1 -5 -1 -``` -""" -function _qr_bond_undo(X::PEPSOrth, a::AbstractTensorMap, b::AbstractTensorMap, Y::PEPSOrth) - @tensor A[-1; -2 -3 -4 -5] := X[-2 1 -4 -5] * a[1 -1 -3] - @tensor B[-1; -2 -3 -4 -5] := b[-5 -1 1] * Y[-2 -3 -4 1] - return A, B -end -function _qr_bond_undo(X::PEPOOrth, a::AbstractTensorMap, b::AbstractTensorMap, Y::PEPOOrth) - if !isdual(space(a, 2)) - @tensor A[-1 -2; -3 -4 -5 -6] := X[-2 -3 1 -5 -6] * a[1 -1 -4] - @tensor B[-1 -2; -3 -4 -5 -6] := b[-6 -1 1] * Y[-2 -3 -4 -5 1] - else - @tensor A[-1 -2; -3 -4 -5 -6] := X[-1 -3 1 -5 -6] * a[1 -2 -4] - @tensor B[-1 -2; -3 -4 -5 -6] := b[-6 -2 1] * Y[-1 -3 -4 -5 1] - end - return A, B -end - -""" -$(SIGNATURES) - -Apply 2-site `gate` on the reduced matrices `a`, `b` -``` - -1← a --- 3 --- b ← -4 -2 -3 - ↓ ↓ ↓ ↓ - 1 2 |----gate---| - ↓ ↓ or ↓ ↓ - |----gate---| 1 2 - ↓ ↓ ↓ ↓ - -2 -3 -1← a --- 3 --- b ← -4 -``` -""" -function _apply_gate( - a::AbstractTensorMap, b::AbstractTensorMap, - gate::AbstractTensorMap{T, S, 2, 2}, trunc::TruncationStrategy - ) where {T <: Number, S <: ElementarySpace} - V = space(b, 1) - need_flip = isdual(V) - if isdual(space(a, 2)) - @tensor a2b2[-1 -2; -3 -4] := gate[1 2; -2 -3] * a[-1 1 3] * b[3 2 -4] - else - @tensor a2b2[-1 -2; -3 -4] := gate[-2 -3; 1 2] * a[-1 1 3] * b[3 2 -4] - end - trunc = if trunc isa FixedSpaceTruncation - need_flip ? truncspace(flip(V)) : truncspace(V) - else - trunc - end - a, s, b, ϵ = svd_trunc!(a2b2; trunc, alg = LAPACK_QRIteration()) - a, b = absorb_s(a, s, b) - if need_flip - a, s, b = flip(a, numind(a)), _fliptwist_s(s), flip(b, 1) - end - return a, s, b, ϵ -end - -""" -Convert a 3-site gate to MPO form by SVD, -in which the axes are ordered as -``` - 2 3 3 - ↓ ↓ ↓ - g1 ←- 3 1 ←- g2 ←- 4 1 ←- g3 - ↓ ↓ ↓ - 1 2 2 -``` -""" -function gate_to_mpo3( - gate::AbstractTensorMap{T, S, 3, 3}, trunc = trunctol(; atol = MPSKit.Defaults.tol) - ) where {T <: Number, S <: ElementarySpace} - Os = MPSKit.decompose_localmpo(MPSKit.add_util_leg(gate), trunc) - g1 = removeunit(Os[1], 1) - g2 = Os[2] - g3 = removeunit(Os[3], 4) - return [g1, g2, g3] -end - -""" -Obtain the 3-site gate MPO on the southeast cluster at position `[row, col]` -``` - r-1 g3 - | - ↓ - r g1 -←- g2 - c c+1 -``` -""" -function _get_gatempo_se(ham::LocalOperator, dt::Number, row::Int, col::Int) - Nr, Nc = size(ham) - @assert 1 <= row <= Nr && 1 <= col <= Nc - sites = [ - CartesianIndex(row, col), - CartesianIndex(row, col + 1), - CartesianIndex(row - 1, col + 1), - ] - nb1x = get_gateterm(ham, (sites[1], sites[2])) - nb1y = get_gateterm(ham, (sites[2], sites[3])) - nb2 = get_gateterm(ham, (sites[1], sites[3])) - # identity operator at each site - units = map(sites) do site - site_ = CartesianIndex(mod1(site[1], Nr), mod1(site[2], Nc)) - return id(physicalspace(ham)[site_]) - end - # when iterating through ┘, └, ┌, ┐ clusters in the unit cell, - # NN / NNN bonds are counted 4 / 2 times, respectively. - @tensor Odt[i' j' k'; i j k] := - -dt * ( - (nb1x[i' j'; i j] * units[3][k' k] + units[1][i'; i] * nb1y[j' k'; j k]) / 4 + - (nb2[i' k'; i k] * units[2][j'; j]) / 2 - ) - op = exp(Odt) - return gate_to_mpo3(op) -end - -""" -Construct the 3-site gate MPOs on the southeast cluster -for 3-site simple update on square lattice. -""" -function _get_gatempos_se(ham::LocalOperator, dt::Number) - Nr, Nc = size(ham.lattice) - return collect(_get_gatempo_se(ham, dt, r, c) for r in 1:Nr, c in 1:Nc) -end diff --git a/src/algorithms/time_evolution/simpleupdate3site.jl b/src/algorithms/time_evolution/simpleupdate3site.jl index bbe432807..ea3d0e2e5 100644 --- a/src/algorithms/time_evolution/simpleupdate3site.jl +++ b/src/algorithms/time_evolution/simpleupdate3site.jl @@ -1,399 +1,3 @@ -#= -# Mixed canonical form of an open boundary MPS -``` - |ψ⟩ = M[1]-←-...-←-M[N] - ↓ ↓ -``` -For convenience, assume all virtual arrows are ←. - -We perform QR and LQ decompositions: starting from -``` - M[1]-←- = Qa[1]-←-R[1]-←- - ↓ ↓ - - -←-M[N] = -←-L[N-1]-←-Qb[N] - ↓ ↓ -``` -we successively calculate -``` - -←-R[n-1]-←-M[n]-←- = -←-Qa[n]-←-R[n]--←-- (n = 2, ..., N - 1) - ↓ ↓ - - -←-M[n+1]-←-L[n+1]-←- = -←-L[n]-←-Qb[n+1]-←- (n = N - 2, ..., 1) - ↓ ↓ -``` - -For each bond (n = 1, ..., N - 1), we perform SVD -``` - R[n] L[n] = U[n] s[n] V†[n] (n = 1, ..., N - 1) -``` -Then we define the projectors together with the Schmidt weight -``` - -←-Pa[n]-←- = L[n] V[n]-←-(1/√s[n])-←- - -←-Pb[n]-←- = -←-(1/√s[n])-←-U†[n] R[n] -``` -The product `Pa Pb` is the identity operator: -``` - Pa[n]-←-Pb[n] = L[n] (R[n] L[n])⁻¹ R[n] = 1 -``` - -The canonical form is then defined by -``` - -←-M̃[n]-←- = -←-Pb[n-1]-←-M[n]-←-Pa[n]-←- - ↓ ↓ -``` - -Note that -``` - M̃[n] - = 1/√s[n-1] U†[n-1] (R[n-1] M[n]) L[n] V[n] 1/√s[n] - = 1/√s[n-1] U†[n-1] Qa[n] (R[n] L[n]) V[n] 1/√s[n] - = 1/√s[n-1] U†[n-1] Qa[n] U[n] s[n] (V†[n] V[n]) 1/√s[n] - = 1/√s[n-1] U†[n-1] Qa[n] U[n] √s[n] -``` -Then `M̃[n]` (n = 1, ..., N - 1) satisfies the (generalized) left-orthogonal condition -``` - ┌---←--M̃[n]--←- ┌-←- 2 - | | | - s[n-1] ↓ = s[n] (s[0] = 1) - | | | - └---→--M̃†[n]-→- └-→- 1 -``` -Similarly, we can express M̃ using Qb -``` - M̃[n] - = 1/√s[n-1] U†[n-1] R[n-1] (M[n] L[n]) V[n] 1/√s[n] - = 1/√s[n-1] U†[n-1] (R[n-1] L[n-1]) Qb[n] V[n] 1/√s[n] - = 1/√s[n-1] U†[n-1] (R[n-1] L[n-1]) Qb[n] V[n] 1/√s[n] - = 1/√s[n-1] (U†[n-1] U[n-1]) s[n-1] V†[n-1] Qb[n] V[n] 1/√s[n] - = √s[n-1] V†[n-1] Qb[n] V[n] 1/√s[n] -``` -Then `M̃[n]` (n = 2, ..., N) satisfies the (generalized) right-orthogonal condition -``` - -←-M̃[n]-←┐ 1 -←-┐ - ↓ | | - * s[n] = s[n-1] (s[N] = 1) - ↓ | | - -→M̃†[n]-→┘ 2 -→-┘ -``` -Here `-*-` is the twist on the physical axis. - -# Truncation of a bond on OBC-MPS - -Suppose we want to truncate the bond between -the n-th and the (n+1)-th sites such that the truncated state -``` - |ψ̃⟩ = M[1]-←-...-←-M̃[n]-←-M̃[n+1]-←-...-←-M[N] - ↓ ↓ ↓ ↓ -``` -maximizes the fidelity -``` - ⟨ψ|ψ̃⟩ ⟨ψ̃|ψ⟩ - F(ψ̃) = ------------- - ⟨ψ̃|ψ̃⟩ ⟨ψ|ψ⟩ -``` -This is simply done by using a truncated `Pa[n], Pb[n]`; then -``` - M[1]--...-M̃[n]--M̃[n+1]--...-M[N] - ⟨ψ|ψ̃⟩ = | | | | - M†[1]-...-M†[n]-M†[n+1]-...-M†[N] - - ┌----M̃[n]--M̃[n+1]--┐ - = s[n-1] | | s[n+1] - └----M†[n]-M†[n+1]-┘ - - ┌--s̃[n]--┐ - = | | = s[n].s̃[n] - └--s[n]--┘ -``` -Then the fidelity is just -``` - F(ψ̃) = (norm(s̃[n], 2) / norm(s[n], 2))^2 -``` -=# -""" -Perform QR decomposition through a PEPS tensor -``` - ╱ ╱ - -←-R0-←-M-←- => ---Q-←-R1-←- - ╱ | ╱ | -``` -""" -function qr_through( - R0::MPSBondTensor, M::GenericMPSTensor{S, 4}; normalize::Bool = true - ) where {S <: ElementarySpace} - @assert !isdual(codomain(R0, 1)) - @assert !isdual(domain(M, 1)) && !isdual(codomain(M, 1)) - @tensor A[-1 -2 -3 -4; -5] := R0[-1; 1] * M[1 -2 -3 -4; -5] - _, r = left_orth!(A; positive = true) - normalize && normalize!(r, Inf) - return r -end -# for `M` at the left end of the MPS -function qr_through( - ::Nothing, M::GenericMPSTensor{S, 4}; normalize::Bool = true - ) where {S <: ElementarySpace} - @assert !isdual(domain(M, 1)) - _, r = left_orth(M; positive = true) - normalize && normalize!(r, Inf) - return r -end - -""" -Perform LQ decomposition through a tensor -``` - ╱ ╱ - -←-L0-←-Q-←- <= -←-M-←-L1-←- - ╱ | ╱ | -``` -""" -function lq_through( - M::GenericMPSTensor{S, 4}, L1::MPSBondTensor; normalize::Bool = true - ) where {S <: ElementarySpace} - @assert !isdual(domain(L1, 1)) - @assert !isdual(codomain(M, 1)) && !isdual(domain(M, 1)) - @tensor A[-1; -2 -3 -4 -5] := M[-1 -2 -3 -4; 1] * L1[1; -5] - l, _ = right_orth!(A; positive = true) - normalize && normalize!(l, Inf) - return l -end -# for `M` at the right end of the MPS -function lq_through( - M::GenericMPSTensor{S, 4}, ::Nothing; normalize::Bool = true - ) where {S <: ElementarySpace} - @assert !isdual(codomain(M, 1)) - A = permute(M, ((1,), (2, 3, 4, 5))) - l, _ = right_orth!(A; positive = true) - normalize && normalize!(l, Inf) - return l -end - -""" -Given a cluster `Ms`, find all `R`, `L` matrices on each internal bond -""" -function _get_allRLs(Ms::Vector{T}) where {T <: GenericMPSTensor{<:ElementarySpace, 4}} - # M1 -- (R1,L1) -- M2 -- (R2,L2) -- M3 - N = length(Ms) - # get the first R and the last L - R_first = qr_through(nothing, Ms[1]; normalize = true) - L_last = lq_through(Ms[N], nothing; normalize = true) - Rs = Vector{typeof(R_first)}(undef, N - 1) - Ls = Vector{typeof(L_last)}(undef, N - 1) - Rs[1], Ls[end] = R_first, L_last - # get remaining R, L matrices - for n in 2:(N - 1) - m = N - n + 1 - Rs[n] = qr_through(Rs[n - 1], Ms[n]; normalize = true) - Ls[m - 1] = lq_through(Ms[m], Ls[m]; normalize = true) - end - return Rs, Ls -end - -""" -Given the tensors `R`, `L` on a bond, construct -the projectors `Pa`, `Pb` and the new bond weight `s` -such that the contraction of `Pa`, `s`, `Pb` is identity when `trunc = notrunc`, - -The arrows between `Pa`, `s`, `Pb` are -``` - - Pa --←-- Pb - - 1 ← s ← 2 -``` -""" -function _proj_from_RL( - r::MPSBondTensor, l::MPSBondTensor; - trunc::TruncationStrategy = notrunc() - ) - @assert isdual(domain(r, 1)) == isdual(codomain(r, 1)) == false - @assert isdual(domain(l, 1)) == isdual(codomain(l, 1)) == false - rl = r * l - u, s, vh, ϵ = svd_trunc!(rl; trunc) - sinv = sdiag_pow(s, -1 / 2) - Pa, Pb = l * vh' * sinv, sinv * u' * r - return Pa, s, Pb, ϵ -end - -""" -Given a cluster `Ms` and the pre-calculated `R`, `L` bond matrices, -find all projectors `Pa`, `Pb` and Schmidt weights `wts` on internal bonds. -""" -function _get_allprojs(Ms, Rs, Ls, truncs::Vector{E}) where {E <: TruncationStrategy} - N = length(Ms) - @assert length(truncs) == N - 1 - projs_errs = map(1:(N - 1)) do i - trunc = if isa(truncs[i], FixedSpaceTruncation) - tspace = space(Ms[i + 1], 1) - isdual(tspace) ? truncspace(flip(tspace)) : truncspace(tspace) - else - truncs[i] - end - return _proj_from_RL(Rs[i], Ls[i]; trunc) - end - Pas = map(Base.Fix2(getindex, 1), projs_errs) - wts = map(Base.Fix2(getindex, 2), projs_errs) - Pbs = map(Base.Fix2(getindex, 3), projs_errs) - # local truncation error on each bond - ϵs = map(Base.Fix2(getindex, 4), projs_errs) - return Pas, Pbs, wts, ϵs -end - -""" -Flip the virtual arrows in the MPS `Ms` -""" -function _flip_virtuals!( - Ms::Vector{T}, flips::Vector{Bool}; inv::Bool = false - ) where {T <: GenericMPSTensor} - @assert length(flips) == length(Ms) - 1 - for (n, flip) in enumerate(flips) - !flip && continue - M1, M2 = Ms[n], Ms[n + 1] - Ms[n] = TensorKit.flip(M1, numind(M1); inv) - Ms[n + 1] = TensorKit.flip(M2, 1; inv) - end - return Ms -end - -""" -Find projectors to truncate internal bonds of the cluster `Ms`. -""" -function _cluster_truncate!( - Ms::Vector{T}, truncs::Vector{E} - ) where {T <: GenericMPSTensor{<:ElementarySpace, 4}, E <: TruncationStrategy} - Rs, Ls = _get_allRLs(Ms) - Pas, Pbs, wts, ϵs = _get_allprojs(Ms, Rs, Ls, truncs) - # apply projectors - # M1 -- (Pa1,wt1,Pb1) -- M2 -- (Pa2,wt2,Pb2) -- M3 - for (i, (Pa, Pb)) in enumerate(zip(Pas, Pbs)) - @tensor (Ms[i])[-1 -2 -3 -4; -5] := (Ms[i])[-1 -2 -3 -4; 1] * Pa[1; -5] - @tensor (Ms[i + 1])[-1 -2 -3 -4; -5] := Pb[-1; 1] * (Ms[i + 1])[1 -2 -3 -4; -5] - end - return wts, ϵs, Pas, Pbs -end - -""" -Apply the gate MPO `gs` on the cluster `Ms`. -When `gate_ax` is 1 or 2, the gate acts from the physical codomain or domain side. - -e.g. Cluster in PEPS with `gate_ax = 1`: -``` - ╱ ╱ ╱ - --- M1 -←-- M2 -←-- M3 --- - ╱ | ╱ | ╱ | - ↓ ↓ ↓ - g1 -←-- g2 -←-- g3 - ↓ ↓ ↓ -``` - -In the cluster, the axes of each tensor use the MPS order -``` - PEPS: PEPO: - 3 3 4 - ╱ | ╱ - 1 -- M -- 5 1 -- M -- 6 - ╱ | ╱ | - 4 2 5 2 - M[1 2 3 4; 5] M[1 2 3 4 5; 6] -``` -""" -function _apply_gatempo!( - Ms::Vector{T1}, gs::Vector{T2}; gate_ax::Int = 1 - ) where {T1 <: GenericMPSTensor{<:ElementarySpace, 4}, T2 <: AbstractTensorMap} - @assert length(Ms) == length(gs) - @assert gate_ax == 1 - @assert all(!isdual(space(g, 1)) for g in gs[2:end]) - @assert all(!isdual(space(M, 1)) for M in Ms[2:end]) - # fusers to merge axes on bonds in the gate-cluster product - # M1 == f1† -- f1 == M2 == f2† -- f2 == M3 - fusers = map(Ms[2:end], gs[2:end]) do M, g - V1, V2 = space(M, 1), space(g, 1) - return isomorphism(fuse(V1, V2) ← V1 ⊗ V2) - end - #= gate on codomain of PEPS - -3 -3 -3 - ╱ ┌-┐ ┌-┐ ╱ ┌-┐ ┌-┐ ╱ - -1 --M--2--┤ | | ├--2--M--4--┤ | | ├--2--M-- -5 - ╱ | | ├- -5 -1 -┤ | ╱ | | ├- -5 -1 -┤ | ╱ | - -4 1 | | | | -4 1 | | | | -4 1 - ├--3--┤ | | ├--3--┼--5--┤ | | ├--3--┤ - -2 └-┘ └-┘ -2 └-┘ └-┘ -2 - =# - for (i, (g, M)) in enumerate(zip(gs, Ms)) - @assert !isdual(space(M, 2)) - if i == 1 - fr = fusers[i] - @tensor (Ms[i])[-1 -2 -3 -4; -5] := M[-1 1 -3 -4; 2] * g[-2 1 3] * fr'[2 3; -5] - elseif i == length(Ms) - fl = fusers[i - 1] - @tensor (Ms[i])[-1 -2 -3 -4; -5] := fl[-1; 2 3] * M[2 1 -3 -4; -5] * g[3 -2 1] - else - fl, fr = fusers[i - 1], fusers[i] - @tensor (Ms[i])[-1 -2 -3 -4; -5] := fl[-1; 2 3] * M[2 1 -3 -4; 4] * g[3 -2 1 5] * fr'[4 5; -5] - end - end - return Ms -end - -function _apply_gatempo!( - Ms::Vector{T1}, gs::Vector{T2}; gate_ax::Int = 1 - ) where {T1 <: GenericMPSTensor{<:ElementarySpace, 5}, T2 <: AbstractTensorMap} - @assert length(Ms) == length(gs) - @assert gate_ax == 1 || gate_ax == 2 - @assert all(!isdual(space(g, 1)) for g in gs[2:end]) - @assert all(!isdual(space(M, 1)) for M in Ms[2:end]) - # fusers to merge axes on bonds in the gate-cluster product - # M1 == f1† -- f1 == M2 == f2† -- f2 == M3 - fusers = map(Ms[2:end], gs[2:end]) do M, g - V1, V2 = space(M, 1), space(g, 1) - return isomorphism(fuse(V1, V2) ← V1 ⊗ V2) - end - #= gate on codomain of PEPO (gate_ax = 1) - - -3 -4 -3 -4 -3 -4 - | ╱ ┌-┐ ┌-┐ | ╱ ┌-┐ ┌-┐ | ╱ - -1 --M--2--┤ | | ├--2--M--4--┤ | | ├--2--M-- -6 - ╱ | | ├- -6 -1 -┤ | ╱ | | ├- -6 -1 -┤ | ╱ | - -5 1 | | | | -5 1 | | | | -5 1 - ├--3--┤ | | ├--3--┼--5--┤ | | ├--3--┤ - -2 └-┘ └-┘ -2 └-┘ └-┘ -2 - - gate on domain of PEPO (gate_ax = 2) - - -3 ┌-┐ ┌-┐ -3 ┌-┐ ┌-┐ -3 - ├--3--┤ | | ├--3--┼--5--┤ | | ├--3--┤ - 1 -4 | ├- -6 -1 -┤ | 1 -4 | ├- -6 -1 -┤ | 1 -4 - | ╱ | | | | | ╱ | | | | | ╱ - -1 --M--2--┤ | | ├--2--M--4--┤ | | ├--2--M-- -6 - ╱ | └-┘ └-┘ ╱ | └-┘ └-┘ ╱ | - -5 -2 -5 -2 -5 -2 - =# - for (i, (g, M)) in enumerate(zip(gs, Ms)) - @assert !isdual(space(M, 2)) - if i == 1 - fr = fusers[i] - if gate_ax == 1 - @tensor (Ms[i])[-1 -2 -3 -4 -5; -6] := M[-1 1 -3 -4 -5; 2] * g[-2 1 3] * fr'[2 3; -6] - else - @tensor (Ms[i])[-1 -2 -3 -4 -5; -6] := M[-1 -2 1 -4 -5; 2] * g[1 -3 3] * fr'[2 3; -6] - end - elseif i == length(Ms) - fl = fusers[i - 1] - if gate_ax == 1 - @tensor (Ms[i])[-1 -2 -3 -4 -5; -6] := fl[-1; 2 3] * M[2 1 -3 -4 -5; -6] * g[3 -2 1] - else - @tensor (Ms[i])[-1 -2 -3 -4 -5; -6] := fl[-1; 2 3] * M[2 -2 1 -4 -5; -6] * g[3 1 -3] - end - else - fl, fr = fusers[i - 1], fusers[i] - if gate_ax == 1 - @tensor (Ms[i])[-1 -2 -3 -4 -5; -6] := fl[-1; 2 3] * M[2 1 -3 -4 -5; 4] * g[3 -2 1 5] * fr'[4 5; -6] - else - @tensor (Ms[i])[-1 -2 -3 -4 -5; -6] := fl[-1; 2 3] * M[2 -2 1 -4 -5; 4] * g[3 1 -3 5] * fr'[4 5; -6] - end - end - end - return Ms -end - function _fuse_physicalspaces(O::GenericMPSTensor{S, 5}) where {S <: ElementarySpace} V1, V2 = codomain(O, 2), codomain(O, 3) F = isomorphism(Int, fuse(V1, V2), V1 ⊗ V2) diff --git a/src/algorithms/time_evolution/time_evolve.jl b/src/algorithms/time_evolution/time_evolve.jl index 38d65eb07..d2b208542 100644 --- a/src/algorithms/time_evolution/time_evolve.jl +++ b/src/algorithms/time_evolution/time_evolve.jl @@ -46,6 +46,26 @@ function _state_bipartite_check(psi::InfiniteState) return true end +""" +Process the Trotter time step `dt` according to the intended usage. +""" +function _get_dt( + state::InfiniteState, dt::Number, imaginary_time::Bool + ) + # PEPS update: exp(-H dt)|ψ⟩ + # PEPO update (purified): exp(-H dt/2)|ρ⟩ + # PEPO update (not purified): exp(-H dt/2) ρ exp(-H dt/2) + dt′ = (state isa InfinitePEPS) ? dt : (dt / 2) + if (state isa InfinitePEPO) + @assert size(state)[3] == 1 + end + if !imaginary_time + @assert (state isa InfinitePEPS) "Real time evolution of InfinitePEPO (Heisenberg picture) is not implemented." + dt′ = complex(zero(dt′), dt′) + end + return dt′ +end + function _timeevol_sanity_check( ψ₀::InfiniteState, Pspaces::M, alg::A ) where {A <: TimeEvolution, M <: AbstractMatrix{<:ElementarySpace}} @@ -60,3 +80,14 @@ function _timeevol_sanity_check( end return nothing end + +function MPSKit.infinite_temperature_density_matrix(H::LocalOperator) + T = scalartype(H) + A = map(physicalspace(H)) do Vp + ψ = permute(TensorKit.id(T, Vp), (1, 2)) + Vv = oneunit(Vp) # trivial (1D) virtual space + virt = ones(T, domain(ψ) ← Vv ⊗ Vv ⊗ Vv' ⊗ Vv') + return ψ * virt + end + return InfinitePEPO(cat(A; dims = 3)) +end diff --git a/src/algorithms/time_evolution/trotter_gate.jl b/src/algorithms/time_evolution/trotter_gate.jl new file mode 100644 index 000000000..dfe151e42 --- /dev/null +++ b/src/algorithms/time_evolution/trotter_gate.jl @@ -0,0 +1,72 @@ +""" + get_expham(H::LocalOperator, dt::Number) + +Compute `exp(-dt * op)` for each term `op` in `H`, +and combine them into a new `LocalOperator`. +Each `op` in `H` must be a single `TensorMap`. +""" +function get_expham(H::LocalOperator, dt::Number) + return LocalOperator( + physicalspace(H), (sites => exp(-dt * op) for (sites, op) in H.terms)... + ) +end + +""" + is_nearest_neighbour(H::LocalOperator) + +Check if an operator `H` contains only nearest neighbor terms. +""" +function is_nearest_neighbour(H::LocalOperator) + return all(H.terms) do (sites, op) + return numin(op) == 2 && sum(abs, Tuple(sites[2] - sites[1])) == 1 + end +end + +""" + is_equivalent_bond(bond1::NTuple{2,CartesianIndex{2}}, bond2::NTuple{2,CartesianIndex{2}}, (Nrow, Ncol)::NTuple{2,Int}) + +Check if two 2-site bonds are related by a (periodic) lattice translation. +""" +function is_equivalent_bond( + bond1::NTuple{2, CartesianIndex{2}}, bond2::NTuple{2, CartesianIndex{2}}, + (Nrow, Ncol)::NTuple{2, Int}, + ) + r1 = bond1[1] - bond1[2] + r2 = bond2[1] - bond2[2] + shift_row = bond1[1][1] - bond2[1][1] + shift_col = bond1[1][2] - bond2[1][2] + return r1 == r2 && mod(shift_row, Nrow) == 0 && mod(shift_col, Ncol) == 0 +end + +""" + get_gateterm(gate::LocalOperator, bond::NTuple{2,CartesianIndex{2}}) + +Get the term of a 2-site gate acting on a certain bond. +Input `gate` should only include one term for each nearest neighbor bond. +""" +function get_gateterm(gate::LocalOperator, bond::NTuple{2, CartesianIndex{2}}) + bonds = findall(p -> is_equivalent_bond(p.first, bond, size(gate.lattice)), gate.terms) + if length(bonds) == 0 + # try reversed site order + bonds = findall( + p -> is_equivalent_bond(p.first, reverse(bond), size(gate.lattice)), gate.terms + ) + if length(bonds) == 1 + return permute(gate.terms[bonds[1]].second, ((2, 1), (4, 3))) + elseif length(bonds) == 0 + # if term not found, return the zero operator on this bond + dtype = scalartype(gate) + r1, c1 = (mod1(bond[1][i], n) for (i, n) in zip(1:2, size(gate))) + r2, c2 = (mod1(bond[2][i], n) for (i, n) in zip(1:2, size(gate))) + V1 = physicalspace(gate)[r1, c1] + V2 = physicalspace(gate)[r2, c2] + return zeros(dtype, V1 ⊗ V2 ← V1 ⊗ V2) + else + error("There are multiple terms in `gate` corresponding to the bond $(bond).") + end + else + (length(bonds) == 1) || + error("There are multiple terms in `gate` corresponding to the bond $(bond).") + return gate.terms[bonds[1]].second + end +end diff --git a/src/algorithms/time_evolution/trotter_mpo.jl b/src/algorithms/time_evolution/trotter_mpo.jl new file mode 100644 index 000000000..d8f1d3e87 --- /dev/null +++ b/src/algorithms/time_evolution/trotter_mpo.jl @@ -0,0 +1,66 @@ +""" +Convert a 3-site gate to MPO form by SVD, +in which the axes are ordered as +``` + 2 3 3 + ↓ ↓ ↓ + g1 ←- 3 1 ←- g2 ←- 4 1 ←- g3 + ↓ ↓ ↓ + 1 2 2 +``` +""" +function gate_to_mpo3( + gate::AbstractTensorMap{T, S, 3, 3}, trunc = trunctol(; atol = MPSKit.Defaults.tol) + ) where {T <: Number, S <: ElementarySpace} + Os = MPSKit.decompose_localmpo(MPSKit.add_util_leg(gate), trunc) + g1 = removeunit(Os[1], 1) + g2 = Os[2] + g3 = removeunit(Os[3], 4) + return [g1, g2, g3] +end + +""" +Obtain the 3-site gate MPO on the southeast cluster at position `[row, col]` +``` + r-1 g3 + | + ↓ + r g1 -←- g2 + c c+1 +``` +""" +function _get_gatempo_se(ham::LocalOperator, dt::Number, row::Int, col::Int) + Nr, Nc = size(ham) + @assert 1 <= row <= Nr && 1 <= col <= Nc + sites = [ + CartesianIndex(row, col), + CartesianIndex(row, col + 1), + CartesianIndex(row - 1, col + 1), + ] + nb1x = get_gateterm(ham, (sites[1], sites[2])) + nb1y = get_gateterm(ham, (sites[2], sites[3])) + nb2 = get_gateterm(ham, (sites[1], sites[3])) + # identity operator at each site + units = map(sites) do site + site_ = CartesianIndex(mod1(site[1], Nr), mod1(site[2], Nc)) + return id(physicalspace(ham)[site_]) + end + # when iterating through ┘, └, ┌, ┐ clusters in the unit cell, + # NN / NNN bonds are counted 4 / 2 times, respectively. + @tensor Odt[i' j' k'; i j k] := + -dt * ( + (nb1x[i' j'; i j] * units[3][k' k] + units[1][i'; i] * nb1y[j' k'; j k]) / 4 + + (nb2[i' k'; i k] * units[2][j'; j]) / 2 + ) + op = exp(Odt) + return gate_to_mpo3(op) +end + +""" +Construct the 3-site gate MPOs on the southeast cluster +for 3-site simple update on square lattice. +""" +function _get_gatempos_se(ham::LocalOperator, dt::Number) + Nr, Nc = size(ham.lattice) + return collect(_get_gatempo_se(ham, dt, r, c) for r in 1:Nr, c in 1:Nc) +end