Add simple and full update algorithms

src/PEPSKit.jl

@@ -17,12 +17,14 @@ include("utility/util.jl")
 include("utility/diffable_threads.jl")
 include("utility/svd.jl")
 include("utility/rotations.jl")
+include("utility/mirror.jl")
 include("utility/diffset.jl")
 include("utility/hook_pullback.jl")
 include("utility/autoopt.jl")
 
 include("states/abstractpeps.jl")
 include("states/infinitepeps.jl")
+include("states/suweight.jl")
 
 include("operators/transferpeps.jl")
 include("operators/infinitepepo.jl")
@@ -38,11 +40,16 @@ include("environments/transferpepo_environments.jl")
 
 include("algorithms/contractions/localoperator.jl")
 include("algorithms/contractions/ctmrg_contractions.jl")
+include("algorithms/contractions/ctmrg_rhos.jl")
+include("algorithms/contractions/measure_rhos.jl")
 
 include("algorithms/ctmrg/sparse_environments.jl")
 include("algorithms/ctmrg/ctmrg.jl")
 include("algorithms/ctmrg/gaugefix.jl")
 
+include("algorithms/time_evolution/simpleupdate.jl")
+include("algorithms/time_evolution/fullupdate.jl")
+
 include("algorithms/toolbox.jl")
 
 include("algorithms/peps_opt.jl")
@@ -59,7 +66,7 @@ include("utility/symmetrization.jl")
         const ctmrgscheme = :simultaneous
         const reuse_env = true
         const trscheme = FixedSpaceTruncation()
-        const fwd_alg = TensorKit.SVD()
+        const fwd_alg = TensorKit.SDD()
         const rrule_alg = Arnoldi(; tol=1e-2fpgrad_tol, krylovdim=48, verbosity=-1)
         const svd_alg = SVDAdjoint(; fwd_alg, rrule_alg)
         const optimizer = LBFGS(32; maxiter=100, gradtol=1e-4, verbosity=2)
@@ -167,6 +174,13 @@ export leading_boundary
 export PEPSOptimize, GeomSum, ManualIter, LinSolver
 export fixedpoint
 
+export absorb_wt, absorb_wt!
+export su_iter!, simpleupdate!
+export fu_iter!, fullupdate!
+export meas_site, meas_bond
+export calrho_site, calrho_bondx, calrho_bondy, calrho_all
+
+export SUWeight
 export InfinitePEPS, InfiniteTransferPEPS
 export InfinitePEPO, InfiniteTransferPEPO
 export initializeMPS, initializePEPS

src/algorithms/contractions/ctmrg_contractions.jl

@@ -530,7 +530,7 @@ Apply bottom projector to southwest corner and south edge.
 function renormalize_bottom_corner((row, col), envs::CTMRGEnv, projectors)
     C_southwest = envs.corners[SOUTHWEST, row, _prev(col, end)]
     E_south = envs.edges[SOUTH, row, col]
-    P_bottom = projectors[1][row, col]
+    P_bottom = projectors[1][row]
     return @autoopt @tensor corner[χ_in; χ_out] :=
         E_south[χ_in D1 D2; χ1] * C_southwest[χ1; χ2] * P_bottom[χ2 D1 D2; χ_out]
 end
@@ -549,7 +549,7 @@ Apply top projector to northwest corner and north edge.
 function renormalize_top_corner((row, col), envs::CTMRGEnv, projectors)
     C_northwest = envs.corners[NORTHWEST, row, _prev(col, end)]
     E_north = envs.edges[NORTH, row, col]
-    P_top = projectors[2][_next(row, end), col]
+    P_top = projectors[2][_next(row, end)]
     return @autoopt @tensor corner[χ_in; χ_out] :=
         P_top[χ_in; χ1 D1 D2] * C_northwest[χ1; χ2] * E_north[χ2 D1 D2; χ_out]
 end
@@ -705,8 +705,8 @@ function renormalize_west_edge(  # For sequential CTMRG scheme
 )
     return renormalize_west_edge(
         envs.edges[WEST, row, _prev(col, end)],
-        projectors[1][row, col],
-        projectors[2][_next(row, end), col],
+        projectors[1][row],
+        projectors[2][_next(row, end)],
         ket[row, col],
         bra[row, col],
     )

src/algorithms/contractions/ctmrg_rhos.jl

@@ -0,0 +1,230 @@
+"""
+Calculate 1-site rho at site `(r,c)`
+```
+    C1 - χ4 - T1 - χ6 - C2  r-1
+    |         ‖         |
+    χ2        DN        χ8
+    |         ‖         |
+    T4 = DW =k/b = DE = T2  r
+    |         ‖         |
+    χ1        DS        χ7
+    |         ‖         |
+    C4 - χ3 - T3 - χ5 - C3  r+1
+    c-1       c        c+1
+```
+Indices d0, d1 are physical indices of ket, bra
+"""
+function calrho_site(
+    row::Int, col::Int, envs::CTMRGEnv, ket::InfinitePEPS, bra::InfinitePEPS=ket
+)
+    N1, N2 = size(ket)
+    @assert 1 <= row <= N1 && 1 <= col <= N2
+    rp1, rm1 = _next(row, N1), _prev(row, N1)
+    cp1, cm1 = _next(col, N2), _prev(col, N2)
+    tket, tbra = ket[row, col], bra[row, col]
+    c1 = envs.corners[1, rm1, cm1]
+    t1 = envs.edges[1, rm1, col]
+    c2 = envs.corners[2, rm1, cp1]
+    t2 = envs.edges[2, row, cp1]
+    c3 = envs.corners[3, rp1, cp1]
+    t3 = envs.edges[3, rp1, col]
+    c4 = envs.corners[4, rp1, cm1]
+    t4 = envs.edges[4, row, cm1]
+    PEPSKit.@autoopt @tensor rho1[d1; d0] := (
+        c4[χ3, χ1] *
+        t4[χ1, DW0, DW1, χ2] *
+        c1[χ2, χ4] *
+        t3[χ5, DS0, DS1, χ3] *
+        tket[d0, DN0, DE0, DS0, DW0] *
+        conj(tbra[d1, DN1, DE1, DS1, DW1]) *
+        t1[χ4, DN0, DN1, χ6] *
+        c3[χ7, χ5] *
+        t2[χ8, DE0, DE1, χ7] *
+        c2[χ6, χ8]
+    )
+    return rho1
+end
+
+"""
+Calculate 2-site rho on sites `(r,c)(r,c+1)`
+```
+    C1 - χ4 - T1 - χ6 - T1 - χ8 - C2    r-1
+    |         ‖         ‖         |
+    χ2       DN1       DN2        χ10
+    |         ‖         ‖         |
+    T4 = DW =k/b = DM =k/b = DE = T2    r
+    |         ‖         ‖         |
+    χ1       DS1       DS2        χ9
+    |         ‖         ‖         |
+    C4 - χ3 - T3 - χ5 - T3 - χ7 - C3    r+1
+    c-1       c        c+1       c+2
+```
+Indices d0, d1 are physical indices of ket, bra
+"""
+function calrho_bondx(
+    row::Int, col::Int, envs::CTMRGEnv, ket::InfinitePEPS, bra::InfinitePEPS=ket
+)
+    N1, N2 = size(ket)
+    @assert 1 <= row <= N1 && 1 <= col <= N2
+    rp1, rm1 = _next(row, N1), _prev(row, N1)
+    cp1, cm1 = _next(col, N2), _prev(col, N2)
+    cp2 = _next(cp1, N2)
+    tket1, tbra1 = ket[row, col], bra[row, col]
+    tket2, tbra2 = ket[row, cp1], bra[row, cp1]
+    c1 = envs.corners[1, rm1, cm1]
+    t11, t12 = envs.edges[1, rm1, col], envs.edges[1, rm1, cp1]
+    c2 = envs.corners[2, rm1, cp2]
+    t2 = envs.edges[2, row, cp2]
+    c3 = envs.corners[3, rp1, cp2]
+    t31, t32 = envs.edges[3, rp1, col], envs.edges[3, rp1, cp1]
+    c4 = envs.corners[4, rp1, cm1]
+    t4 = envs.edges[4, row, cm1]
+    PEPSKit.@autoopt @tensor rho2[d11, d21; d10, d20] := (
+        c4[χ3, χ1] *
+        t4[χ1, DW0, DW1, χ2] *
+        c1[χ2, χ4] *
+        t31[χ5, DS10, DS11, χ3] *
+        tket1[d10, DN10, DM0, DS10, DW0] *
+        conj(tbra1[d11, DN11, DM1, DS11, DW1]) *
+        t11[χ4, DN10, DN11, χ6] *
+        t32[χ7, DS20, DS21, χ5] *
+        tket2[d20, DN20, DE0, DS20, DM0] *
+        conj(tbra2[d21, DN21, DE1, DS21, DM1]) *
+        t12[χ6, DN20, DN21, χ8] *
+        c3[χ9, χ7] *
+        t2[χ10, DE0, DE1, χ9] *
+        c2[χ8, χ10]
+    )
+    return rho2
+end
+
+"""
+Calculate 2-site rho on sites `(r,c)(r-1,c)`
+```
+    C1 - χ9 - T1 -χ10 - C2  r-2
+    |         ‖         |
+    χ7        DN        χ8
+    |         ‖         |
+    T4 = DW2=k/b =DE2 = T2  r-1
+    |         ‖         |
+    χ5        DM        χ6
+    |         ‖         |
+    T4 = DW1=k/b =DE1 = T2  r
+    |         ‖         |
+    χ3        DS        χ4
+    |         ‖         |
+    C4 - χ1 - T3 - χ2 - C3  r+1
+    c-1       c        c+1
+```
+Indices d0, d1 are physical indices of ket, bra
+"""
+function calrho_bondy(
+    row::Int, col::Int, envs::CTMRGEnv, ket::InfinitePEPS, bra::InfinitePEPS=ket
+)
+    N1, N2 = size(ket)
+    @assert 1 <= row <= N1 && 1 <= col <= N2
+    rp1, rm1 = _next(row, N1), _prev(row, N1)
+    cp1, cm1 = _next(col, N2), _prev(col, N2)
+    rm2 = _prev(rm1, N1)
+    tket1, tbra1 = ket[row, col], bra[row, col]
+    tket2, tbra2 = ket[rm1, col], bra[rm1, col]
+    c1 = envs.corners[1, rm2, cm1]
+    t1 = envs.edges[1, rm2, col]
+    c2 = envs.corners[2, rm2, cp1]
+    t21, t22 = envs.edges[2, row, cp1], envs.edges[2, rm1, cp1]
+    c3 = envs.corners[3, rp1, cp1]
+    t3 = envs.edges[3, rp1, col]
+    c4 = envs.corners[4, rp1, cm1]
+    t41, t42 = envs.edges[4, row, cm1], envs.edges[4, rm1, cm1]
+    PEPSKit.@autoopt @tensor rho2[d11, d21; d10, d20] := (
+        c4[χ1, χ3] *
+        t3[χ2, DS0, DS1, χ1] *
+        c3[χ4, χ2] *
+        t41[χ3, DW10, DW11, χ5] *
+        tket1[d10, DM0, DE10, DS0, DW10] *
+        conj(tbra1[d11, DM1, DE11, DS1, DW11]) *
+        t21[χ6, DE10, DE11, χ4] *
+        t42[χ5, DW20, DW21, χ7] *
+        tket2[d20, DN0, DE20, DM0, DW20] *
+        conj(tbra2[d21, DN1, DE21, DM1, DW21]) *
+        t22[χ8, DE20, DE21, χ6] *
+        c1[χ7, χ9] *
+        t1[χ9, DN0, DN1, χ10] *
+        c2[χ10, χ8]
+    )
+    return rho2
+end
+
+# TODO: add rhos on next nearest neighbor bonds
+
+"""
+Calculate 2-site rho on 2nd nearest neighbor sites `(r,c)(r-1,c+1)`
+```
+    C1 -χ10 - T1 -χ11 - T1 -χ12 - C2    r-2
+    |         ‖         ‖         |
+    χ8       DN1       DN2        χ9
+    |         ‖         ‖         |
+    T4 =DW2= k/b =DH2= k/b =DE2== T2    r-1
+    |         ‖         ‖         |
+    χ6       DV1       DV2        χ7
+    |         ‖         ‖         |
+    T4 =DW1= k/b =DH1= k/b =DE1== T2    r
+    |         ‖         ‖         |
+    χ4       DS1       DS2        χ5
+    |         ‖         ‖         |
+    C4 - χ1 - T3 - χ2 - T3 - χ3 - C4    r+1
+    c-1       c        c+1       c+2
+```
+Indices d0, d1 are physical indices of ket, bra
+"""
+function calrho_bondd1(
+    row::Int, col::Int, envs::CTMRGEnv, ket::InfinitePEPS, bra::InfinitePEPS=ket
+)
+    N1, N2 = size(ket)
+    @assert 1 <= row <= N1 && 1 <= col <= N2
+    throw("not implemented")
+end
+
+"""
+Calculate 2-site rho on 2nd nearest neighbor sites `(r,c+1)(r-1,c)`
+"""
+function calrho_bondd2(
+    row::Int, col::Int, envs::CTMRGEnv, ket::InfinitePEPS, bra::InfinitePEPS=ket
+)
+    N1, N2 = size(ket)
+    @assert 1 <= row <= N1 && 1 <= col <= N2
+    throw("not implemented")
+end
+
+"""
+Calculate rho for all sites
+"""
+function calrho_allsites(envs::CTMRGEnv, ket::InfinitePEPS, bra::InfinitePEPS=ket)
+    Nr, Nc = size(ket)
+    return collect(
+        calrho_site(r, c, envs, ket, bra) for (r, c) in Iterators.product(1:Nr, 1:Nc)
+    )
+end
+
+"""
+Calculate rho for all nearest-neighbor bonds
+"""
+function calrho_allnbs(envs::CTMRGEnv, ket::InfinitePEPS, bra::InfinitePEPS=ket)
+    Nr, Nc = size(ket)
+    rhoxss = collect(
+        calrho_bondx(r, c, envs, ket, bra) for (r, c) in Iterators.product(1:Nr, 1:Nc)
+    )
+    rhoyss = collect(
+        calrho_bondy(r, c, envs, ket, bra) for (r, c) in Iterators.product(1:Nr, 1:Nc)
+    )
+    return [rhoxss, rhoyss]
+end
+
+"""
+Calculate rho for all sites and nearest-neighbor bonds
+"""
+function calrho_all(envs::CTMRGEnv, ket::InfinitePEPS, bra::InfinitePEPS=ket)
+    rho1ss = calrho_allsites(envs, ket, bra)
+    rho2sss = calrho_allnbs(envs, ket, bra)
+    return rho1ss, rho2sss
+end

src/algorithms/contractions/measure_rhos.jl

@@ -0,0 +1,39 @@
+"""
+Get identity operator on the physical space
+"""
+function _getid_from_rho(rho::AbstractTensorMap)
+    Pspace = codomain(rho)[1]
+    if isdual(Pspace)
+        Pspace = adjoint(Pspace)
+    end
+    return TensorKit.id(Pspace)
+end
+
+"""
+Measure `<op>` using 1-site rho
+"""
+function meas_site(op::AbstractTensorMap, rho1::AbstractTensorMap)
+    Id = _getid_from_rho(rho1)
+    val = ncon((rho1, op), ([1, 2], [1, 2]))
+    nrm = ncon((rho1, Id), ([1, 2], [1, 2]))
+    meas = first(blocks(val / nrm))[2][1]
+    return meas
+end
+
+"""
+Measure `<op1 op2>` using 2-site rho
+"""
+function meas_bond(op1::AbstractTensorMap, op2::AbstractTensorMap, rho2::AbstractTensorMap)
+    return meas_bond(op1 ⊗ op2, rho2)
+end
+
+"""
+Measure `<gate>` using 2-site rho
+"""
+function meas_bond(gate::AbstractTensorMap, rho2::AbstractTensorMap)
+    Id = _getid_from_rho(rho2)
+    val = ncon((rho2, gate), ([1, 2, 3, 4], [1, 2, 3, 4]))
+    nrm = ncon((rho2, Id ⊗ Id), ([1, 2, 3, 4], [1, 2, 3, 4]))
+    meas = first(blocks(val / nrm))[2][1]
+    return meas
+end