add Triangular CTMRG

src/PEPSKit.jl

@@ -52,8 +52,12 @@ include("utility/retractions.jl")
 include("networks/tensors.jl")
 include("networks/local_sandwich.jl")
 include("networks/infinitesquarenetwork.jl")
+include("networks/tensors_triangular.jl")
+include("networks/local_triangular_sandwich.jl")
+include("networks/infinitetriangularnetwork.jl")
 
 include("states/infinitepeps.jl")
+include("states/infinitepepstriangular.jl")
 include("states/infinitepartitionfunction.jl")
 
 include("utility/symmetrization.jl")
@@ -65,6 +69,7 @@ include("operators/lattices/squarelattice.jl")
 include("operators/models.jl")
 
 include("environments/ctmrg_environments.jl")
+include("environments/ctmrg_environments_triangular.jl")
 include("environments/vumps_environments.jl")
 include("environments/suweight.jl")
 include("environments/bp_environments.jl")
@@ -100,6 +105,11 @@ include("algorithms/ctmrg/sequential.jl")
 include("algorithms/ctmrg/gaugefix.jl")
 include("algorithms/ctmrg/c4v.jl")
 
+include("algorithms/contractions/ctmrg_contractions_triangular.jl")
+include("algorithms/ctmrg_triangular/ctmrg.jl")
+include("algorithms/ctmrg_triangular/utility.jl")
+include("algorithms/ctmrg_triangular/simultaneous.jl")
+
 include("algorithms/truncation/truncationschemes.jl")
 include("algorithms/truncation/fullenv_truncation.jl")
 include("algorithms/truncation/bond_truncation.jl")
@@ -115,6 +125,7 @@ include("algorithms/bp/gaugefix.jl")
 
 include("algorithms/transfermatrix.jl")
 include("algorithms/toolbox.jl")
+include("algorithms/toolbox_triangular.jl")
 include("algorithms/correlators.jl")
 
 include("algorithms/optimization/fixed_point_differentiation.jl")
@@ -126,6 +137,7 @@ using .Defaults: set_scheduler!
 export set_scheduler!
 export EighAdjoint, IterEigh, SVDAdjoint, IterSVD, QRAdjoint
 export CTMRGEnv, SequentialCTMRG, SimultaneousCTMRG
+export CTMRGEnvTriangular, SimultaneousCTMRGTriangular # CTMRGTriaEnv, SimultaneousCTMRGTria
 export FixedSpaceTruncation, SiteDependentTruncation
 export HalfInfiniteProjector, FullInfiniteProjector
 export C4vCTMRG, C4vEighProjector, C4vQRProjector
@@ -143,9 +155,10 @@ export ALSTruncation, FullEnvTruncation
 export SimpleUpdate
 export TimeEvolver, timestep, time_evolve
 
-export InfiniteSquareNetwork
+export InfiniteSquareNetwork, InfiniteTriangularNetwork
 export InfinitePartitionFunction
 export InfinitePEPS, InfiniteTransferPEPS
+export InfinitePEPSTriangular
 export SUWeight
 export InfinitePEPO, InfiniteTransferPEPO
 

src/algorithms/contractions/ctmrg_contractions_triangular.jl

@@ -0,0 +1,58 @@
+function build_double_corner_matrix_triangular(network::InfiniteTriangularNetwork{P}, env::CTMRGEnvTriangular, dir::Int, r::Int, c::Int) where {P <: PFTensorTriangular}
+    @tensor opt = true mat[-1 -2; -3 -4] := env.C[_coordinates(dir, 0, r, c, size(network))...][6 5; 1] * env.C[_coordinates(dir + 1, 0, r, c, size(network))...][1 3; 2] *
+        env.Ea[_coordinates(dir - 1, 0, r, c, size(network))...][-1 7; 6] * env.Eb[_coordinates(dir + 1, 0, r, c, size(network))...][2 4; -3] * rotl60(network[r, c], mod(dir - 1, 6))[7 -2 -4; 5 3 4]
+    return mat
+end
+
+function build_double_corner_matrix_triangular(network::InfiniteTriangularNetwork{P}, env::CTMRGEnvTriangular, dir::Int, r::Int, c::Int) where {P <: PEPSSandwichTriangular}
+    @tensor opt = true mat[-1 -2 -3; -4 -5 -6] := env.C[_coordinates(dir, 0, r, c, size(network))...][6 5 8; 1] * env.C[_coordinates(dir + 1, 0, r, c, size(network))...][1 3 9; 2] *
+        env.Ea[_coordinates(dir - 1, 0, r, c, size(network))...][-1 7 10; 6] * env.Eb[_coordinates(dir + 1, 0, r, c, size(network))...][2 4 11; -4] *
+        rotl60(ket(network[r, c]), mod(dir - 1, 6))[12; 5 3 4 -5 -2 7] * conj(rotl60(bra(network[r, c]), mod(dir - 1, 6))[12; 8 9 11 -6 -3 10])
+    return mat
+end
+
+function semi_renormalize_edge(network::InfiniteTriangularNetwork{P}, env::CTMRGEnvTriangular, Pas, Pbs, dir, r, c) where {P <: PEPSSandwichTriangular}
+    @tensor opt = true mat[-1 -2 -3; -4 -5 -6] := Pas[dir, r, c][-1; 1 2 8] * Pbs[dir, r, c][6 7 9; -4] *
+        env.Ea[_coordinates(dir, 0, r, c, size(network))...][1 3 10; 4] * env.Eb[_coordinates(dir, 1, r, c, size(network))...][4 5 11; 6] *
+        rotl60(ket(network[r, c]), mod(dir - 1, 6))[12; 3 5 7 -5 -2 2] * conj(rotl60(bra(network[r, c]), mod(dir - 1, 6))[12; 10 11 9 -6 -3 8])
+    return mat / norm(mat)
+end
+
+function semi_renormalize_edge(network::InfiniteTriangularNetwork{P}, env::CTMRGEnvTriangular, Pas, Pbs, dir, r, c) where {P <: PFTensorTriangular}
+    @tensor opt = true mat[-1 -2; -3 -4] := Pas[dir, r, c][-1; 1 2] * Pbs[dir, r, c][6 7; -3] *
+        env.Ea[_coordinates(dir, 0, r, c, size(network))...][1 3; 4] * env.Eb[_coordinates(dir, 1, r, c, size(network))...][4 5; 6] * rotl60(network[r, c], mod(dir - 1, 6))[2 -2 -4; 3 5 7]
+    return mat / norm(mat)
+end
+
+function build_halfinfinite_projectors(network::InfiniteTriangularNetwork{P}, env::CTMRGEnvTriangular, Ẽas, Ẽbs, Ẽastr, Ẽbstr, dir, r, c) where {P <: PEPSSandwichTriangular}
+    @tensor opt = true σL[-1 -2 -3; -4] := env.C[_coordinates(dir, 0, r, c, size(network))...][1 2 10; 8] * env.C[_coordinates(dir - 1, 0, r, c, size(network))...][4 3 11; 1] * env.C[_coordinates(dir - 2, 0, r, c, size(network))...][6 5 12; 4] *
+        Ẽbs[_coordinates(dir, 1, r, c, size(network))...][8 9 13; -4] * Ẽastr[_coordinates(dir - 3, 0, r, c, size(network))...][-1 7 14; 6] *
+        rotl60(ket(network[r, c]), mod(dir - 1, 6))[15 2 9 -2 7 5 3] * conj(rotl60(bra(network[r, c]), mod(dir - 1, 6))[15; 10 13 -3 14 12 11])
+    @tensor opt = true σR[-1; -2 -3 -4] := env.C[_coordinates(dir + 1, 0, r, c, size(network))...][8 2 10; 1] * env.C[_coordinates(dir + 2, 0, r, c, size(network))...][1 3 11; 4] * env.C[_coordinates(dir + 3, 0, r, c, size(network))...][4 5 12; 6] *
+        Ẽbstr[_coordinates(dir + 3, -1, r, c, size(network))...][6 7 13; -2] * Ẽas[_coordinates(dir, 0, r, c, size(network))...][-1 9 14; 8] *
+        rotl60(ket(network[_coordinates(dir + 2, 0, r, c, size(network))[2:3]...]), mod(dir - 1, 6))[15; 9 2 3 5 7 -3] * conj(rotl60(bra(network[_coordinates(dir + 2, 0, r, c, size(network))[2:3]...]), mod(dir - 1, 6))[15; 14 10 11 12 13 -4])
+    return σL, σR
+end
+
+function build_halfinfinite_projectors(network::InfiniteTriangularNetwork{P}, env::CTMRGEnvTriangular, Ẽas, Ẽbs, Ẽastr, Ẽbstr, dir, r, c) where {P <: PFTensorTriangular}
+    @tensor opt = true σL[-1 -2; -3] := env.C[_coordinates(dir, 0, r, c, size(network))...][1 2; 8] * env.C[_coordinates(dir - 1, 0, r, c, size(network))...][4 3; 1] * env.C[_coordinates(dir - 2, 0, r, c, size(network))...][6 5; 4] *
+        Ẽbs[_coordinates(dir, 1, r, c, size(network))...][8 9; -3] * Ẽastr[_coordinates(dir - 3, 0, r, c, size(network))...][-1 7; 6] *
+        rotl60(network[r, c], mod(dir - 1, 6))[3 5 7; 2 9 -2]
+    @tensor opt = true σR[-1; -2 -3] := env.C[_coordinates(dir + 1, 0, r, c, size(network))...][8 2; 1] * env.C[_coordinates(dir + 2, 0, r, c, size(network))...][1 3; 4] * env.C[_coordinates(dir + 3, 0, r, c, size(network))...][4 5; 6] *
+        Ẽbstr[_coordinates(dir + 3, -1, r, c, size(network))...][6 7; -2] * Ẽas[_coordinates(dir, 0, r, c, size(network))...][-1 9; 8] *
+        rotl60(network[_coordinates(dir + 2, 0, r, c, size(network))[2:3]...], mod(dir - 1, 6))[-3 7 5; 9 2 3]
+    return σL, σR
+end
+
+function renormalize_corner_triangular((dir, r, c), network::InfiniteTriangularNetwork{P}, env, Pas, Pbs) where {P <: PFTensorTriangular}
+    @tensor opt = true env_C_new[-1 -2; -3] := env.C[_coordinates(dir, 0, r, c, size(network))...][1 3; 6] * env.Ea[_coordinates(dir - 1, 0, r, c, size(network))...][4 2; 1] * env.Eb[_coordinates(dir, 1, r, c, size(network))...][6 7; 8] *
+        rotl60(network[r, c], mod(dir - 1, 6))[2 5 -2; 3 7 9] * Pas[mod1(dir - 1, 6), r, c][-1; 4 5] * Pbs[dir, r, c][8 9; -3]
+    return env_C_new
+end
+
+function renormalize_corner_triangular((dir, r, c), network::InfiniteTriangularNetwork{P}, env, Pas, Pbs) where {P <: PEPSSandwichTriangular}
+    @tensor opt = true env_C_new[-1 -2 -3; -4] := env.C[_coordinates(dir, 0, r, c, size(network))...][1 3 10; 6] * env.Ea[_coordinates(dir - 1, 0, r, c, size(network))...][4 2 11; 1] * env.Eb[_coordinates(dir, 1, r, c, size(network))...][6 7 12; 8] *
+        rotl60(ket(network[r, c]), mod(dir - 1, 6))[15; 3 7 9 -2 5 2] * conj(rotl60(bra(network[r, c]), mod(dir - 1, 6))[15; 10 12 14 -3 13 11]) *
+        Pas[mod1(dir - 1, 6), r, c][-1; 4 5 13] * Pbs[dir, r, c][8 9 14; -4]
+    return env_C_new
+end

src/algorithms/ctmrg_triangular/ctmrg.jl

@@ -0,0 +1,50 @@
+abstract type CTMRGAlgorithmTriangular end
+
+function leading_boundary(env₀::CTMRGEnvTriangular, network::InfiniteTriangularNetwork; kwargs...)
+    alg = select_algorithm(leading_boundary, env₀; kwargs...)
+    return leading_boundary(env₀, network, alg)
+end
+function leading_boundary(
+        env₀::CTMRGEnvTriangular, network::InfiniteTriangularNetwork, alg::CTMRGAlgorithmTriangular
+    )
+    log = ignore_derivatives(() -> MPSKit.IterLog("CTMRG"))
+    return LoggingExtras.withlevel(; alg.verbosity) do
+        env = deepcopy(env₀)
+        S_old = DiagonalTensorMap.(id.(domain.(env₀.C)))
+        η = one(real(scalartype(network)))
+        ctmrg_loginit!(log, η, network, env₀)
+        local info
+        for iter in 1:(alg.maxiter)
+            env, S = ctmrg_iteration(network, env, alg)  # Grow and renormalize in all 4 directions
+            η = calc_convergence(S, S_old)
+            S_old = copy(S)
+
+            if η ≤ alg.tol && iter ≥ alg.miniter
+                ctmrg_logfinish!(log, iter, η, network, env)
+                break
+            end
+            if iter == alg.maxiter
+                ctmrg_logcancel!(log, iter, η, network, env)
+            else
+                ctmrg_logiter!(log, iter, η, network, env)
+            end
+            info = (; convergence_metric = η)
+        end
+        return env, info
+    end
+end
+function leading_boundary(env₀::CTMRGEnvTriangular, state, args...; kwargs...)
+    return leading_boundary(env₀, InfiniteTriangularNetwork(state), args...; kwargs...)
+end
+
+function calc_convergence(S_new::Array{T, 3}, S_old::Array{T, 3}) where {E, S, T <: DiagonalTensorMap{E, S}}
+    ε = Inf
+    function dist(S1, S2)
+        if space(S1) == space(S2)
+            return norm(S1^4 - S2^4)
+        else
+            return Inf
+        end
+    end
+    return maximum(dist.(S_new, S_old))
+end