A method to compute the conformal spin up to 10 (#74)

* Added functions to compress the tensor network in a single direction

* Added function to calculate conformal spin

* format

* Added generic QR_R,L and the function to calculate spin up to 8

* format

* Deleted junction, added conformal dimension

* changed the corresponding dependent function in symmetric loop tnr

* Deleted garbages

* Format

* Trying to add the loop optimization to transfer MPO

* It seems that preserving translational symmetry is very important. So it may be better to only performing loop in a unit cell...

* The planar_opt is still the best

* format

* The 1x4 case is extremely precise

The loop optimization of the transfer MPO can indeed reduce the truncation error by 1/10, however, the final spectrum becomes worse when loop is added. Still need to find the reason

* Reorganized some functions in cft_data

* back truncation of the virtual bond

* Change back

* Changed another configuration

Now the cost decreases to O(D^6) and can calculate spins from -5 to 5 now!

* Deleted garbage

* Deleted some export

Author: Chenqitrg

Project.toml

@@ -16,13 +16,13 @@ Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [compat]
 DocStringExtensions = "0.9.4"
-julia = "1.10"
 KrylovKit = "0.9.5"
+LoggingExtras = "~1.0"
 OptimKit = "0.4.0"
 PEPSKit = "0.5.0"
 TensorKit = "0.14"
 Zygote = "0.7.7"
-LoggingExtras = "~1.0"
+julia = "1.10"
 
 [extras]
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"

src/TNRKit.jl

@@ -76,6 +76,8 @@ export finalize!, finalize_two_by_two!, finalize_cftdata!, finalize_central_char
 include("utility/cdl.jl")
 export cdl_tensor
 
+include("utility/projectors.jl")
+
 include("utility/blocking.jl")
 export block_tensors
 end

src/schemes/looptnr.jl

@@ -43,189 +43,25 @@ function Ψ_A(scheme::LoopTNR)
     return psi
 end
 
-#Utility functions for QR decomp
-
-# A single step of the QR decomposition from the left with 3 in-coming legs
-#       |     |
-#        2   3
-#         v v
-# --L-1-<--T--<-4-----
-# =
-#       |     |
-#        2   3
-#         v v
-# ----1-<--Q--<-4--Rt-
-
-function QR_L(L::TensorMap, T::AbstractTensorMap{E,S,1,3}) where {E,S}
-    @planar LT[-1; -2 -3 -4] := L[-1; 1] * T[1; -2 -3 -4]
-    temp = transpose(LT, (3, 2, 1), (4,); copy=true)
-    _, Rt = leftorth(temp)
-    return Rt / norm(Rt, Inf)
-end
-
-# A single step of the QR decomposition from the right with 3 in-coming legs
-#        |     |
-#         2   3
-#          v v
-# -----1-<--T--<-4-R--
-# =
-#        |     |
-#         2   3
-#          v v
-# -Lt--1-<--Q--<-4----
-
-function QR_R(R::TensorMap, T::AbstractTensorMap{E,S,1,3}) where {E,S}
-    @planar TR[-1; -2 -3 -4] := T[-1; -2 -3 1] * R[1; -4]
-    Lt, _ = rightorth(TR)
-    return Lt / norm(Lt, Inf)
-end
-
-# A single step of the QR decomposition from the left with 2 in-coming legs
-#          |
-#          2
-#          v
-# --L-1-<--T--<-3----
-# =
-#          | 
-#          2
-#          v
-# ----1-<--Q--<-3--Rt-
-
-function QR_L(L::TensorMap, T::AbstractTensorMap{E,S,1,2}) where {E,S}
-    @planar LT[-1; -2 -3] := L[-1; 1] * T[1; -2 -3]
-    temp = transpose(LT, (2, 1), (3,); copy=true)
-    _, Rt = leftorth(temp)
-    return Rt / norm(Rt, Inf)
-end
-
-# A single step of the QR decomposition from the right with 2 in-coming legs
-#           |
-#           2
-#           v
-# -----1-<--T--<-3-R--
-# =
-#           |
-#           2
-#           v
-# -Lt--1-<--Q--<-3----
-function QR_R(R::TensorMap, T::AbstractTensorMap{E,S,1,2}) where {E,S}
-    @planar TR[-1; -2 -3] := T[-1; -2 1] * R[1; -3]
-    Lt, _ = rightorth(TR)
-    return Lt / norm(Lt, Inf)
-end
-
-# Functions to find the left and right projectors
-
-# Function to find the list of left projectors L_list
-function find_L(psi::Array, entanglement_criterion::stopcrit)
-    type = eltype(psi[1])
-    n = length(psi)
-    L_list = map(x -> id(type, codomain(psi[x])[1]), 1:n)
-    crit = true
-    steps = 0
-    error = [Inf]
-    running_pos = 1
-    while crit
-        pos_next = mod(running_pos, n) + 1
-        L_last_time = L_list[pos_next]
-        L_list[pos_next] = QR_L(L_list[running_pos], psi[running_pos])
-
-        if space(L_list[pos_next]) == space(L_last_time)
-            push!(error, abs(norm(L_list[pos_next] - L_last_time)))
-        end
-
-        running_pos = pos_next
-        steps += 1
-        crit = entanglement_criterion(steps, error)
-    end
-
-    return L_list
-end
-
-# Function to find the list of right projectors R_list
-function find_R(psi::Array, entanglement_criterion::stopcrit)
-    type = eltype(psi[1])
-    n = length(psi)
-    R_list = map(x -> id(type, domain(psi[x]).spaces[end]), 1:n)
-    crit = true
-    steps = 0
-    error = [Inf]
-
-    running_pos = n
-    while crit
-        pos_last = mod(running_pos - 2, n) + 1
-        R_last_time = R_list[pos_last]
-        R_list[pos_last] = QR_R(R_list[running_pos], psi[running_pos])
-
-        if space(R_list[pos_last]) == space(R_last_time)
-            push!(error, abs(norm(R_list[pos_last] - R_last_time)))
-        end
-
-        running_pos = pos_last
-        steps += 1
-        crit = entanglement_criterion(steps, error)
-    end
-
-    return R_list
-end
-
-# Function to find the projector P_L and P_R
-function P_decomp(R::TensorMap, L::TensorMap, trunc::TensorKit.TruncationScheme)
-    U, S, V, _ = tsvd(L * R; trunc=trunc, alg=TensorKit.SVD())
-    re_sq = pseudopow(S, -0.5)
-    PR = R * V' * re_sq
-    PL = re_sq * U' * L
-    return PR, PL
-end
-
-# Function to find the list of projectors
-function find_projectors(psi::Array, entanglement_criterion::stopcrit,
-                         trunc::TensorKit.TruncationScheme)
-    PR_list = []
-    PL_list = []
-
-    n = length(psi)
-    L_list = find_L(psi, entanglement_criterion)
-    R_list = find_R(psi, entanglement_criterion)
-    for i in 1:n
-        pr, pl = P_decomp(R_list[mod(i - 2, n) + 1], L_list[i], trunc)
-
-        push!(PR_list, pr)
-        push!(PL_list, pl)
-    end
-    return PR_list, PL_list
-end
-
-function SVD12(T::AbstractTensorMap{E,S,1,3}, trunc::TensorKit.TruncationScheme) where {E,S}
-    T_trans = transpose(T, (2, 1), (3, 4); copy=true)
-    U, s, V, _ = tsvd(T_trans; trunc=trunc, alg=TensorKit.SVD())
-    @planar S1[-1; -2 -3] := U[-2 -1; 1] * sqrt(s)[1; -3]
-    @planar S2[-1; -2 -3] := sqrt(s)[-1; 1] * V[1; -2 -3]
-    return S1, S2
-end
-
 # Function to construct MPS Ψ_B from MPS Ψ_A. Using a large cut-off dimension in SVD but a small cut-off dimension in loop to increase the precision of initialization.
 function Ψ_B(ΨA, trunc::TensorKit.TruncationScheme,
              truncentanglement::TensorKit.TruncationScheme)
+    NA = length(ΨA)
     ΨB = []
-
-    for i in 1:4
+    for i in 1:NA
         s1, s2 = SVD12(ΨA[i], truncdim(trunc.dim * 2))
         push!(ΨB, s1)
         push!(ΨB, s2)
     end
 
     ΨB_function(steps, data) = abs(data[end])
     criterion = maxiter(10) & convcrit(1e-12, ΨB_function)
-    PR_list, PL_list = find_projectors(ΨB, criterion, trunc & truncentanglement)
-
-    ΨB_disentangled = []
-    for i in 1:8
-        @planar B1[-1; -2 -3] := PL_list[i][-1; 1] * ΨB[i][1; -2 2] *
-                                 PR_list[mod(i, 8) + 1][2; -3]
-        push!(ΨB_disentangled, B1)
-    end
-    return ΨB_disentangled
+    in_inds = ones(Int, 2*NA)
+    out_inds = 2*ones(Int, 2*NA)
+    PR_list, PL_list = find_projectors(ΨB, in_inds, out_inds, criterion,
+                                       trunc & truncentanglement)
+    MPO_disentangled!(ΨB, in_inds, out_inds, PR_list, PL_list)
+    return ΨB
 end
 
 # Construct the list of transfer matrices for ΨAΨA
@@ -235,8 +71,9 @@ end
 #       | |
 # ---2'--A--4'---
 function ΨAΨA(psiA)
+    NA = length(psiA)
     ΨAΨA_list = []
-    for i in 1:4
+    for i in 1:NA
         @planar tmp[-1 -2; -3 -4] := psiA[i][-2; 1 2 -4] * psiA[i]'[1 2 -3; -1]
         push!(ΨAΨA_list, tmp)
     end
@@ -251,8 +88,9 @@ end
 # ---2'--B--4'---
 
 function ΨBΨB(psiB)
+    NB = length(psiB)
     ΨBΨB_list = []
-    for i in 1:8
+    for i in 1:NB
         @planar tmp[-1 -2; -3 -4] := psiB[i][-2; 1 -4] * psiB[i]'[1 -3; -1]
         push!(ΨBΨB_list, tmp)
     end
@@ -267,8 +105,9 @@ end
 # ---2'----A----4'---
 
 function ΨBΨA(psiB, psiA)
+    NA = length(psiA)
     ΨBΨA_list = []
-    for i in 1:4
+    for i in 1:NA
         @planar temp[-1 -2; -3 -4] := psiB[2 * i - 1]'[1 3; -1] * psiA[i][-2; 1 2 -4] *
                                       psiB[2 * i]'[2 -3; 3]
         push!(ΨBΨA_list, temp)
@@ -295,7 +134,8 @@ loop_criterion = maxiter(50) & convcrit(1e-8, entanglement_function)
 function entanglement_filtering!(scheme::LoopTNR, entanglement_criterion::stopcrit,
                                  trunc::TensorKit.TruncationScheme)
     ΨA = Ψ_A(scheme)
-    PR_list, PL_list = find_projectors(ΨA, entanglement_criterion, trunc)
+    PR_list, PL_list = find_projectors(ΨA, [1, 1, 1, 1], [3, 3, 3, 3],
+                                       entanglement_criterion, trunc)
 
     TA = copy(scheme.TA)
     TB = copy(scheme.TB)
@@ -387,20 +227,20 @@ function right_cache(tensor_list)
     return cache
 end
 
-# Function to perform the optimization loop for the LoopTNR scheme. Sweeping from left to right, we optimize the tensors in the loop by minimizing the cost function.
+# A general function to optimize the truncation error of an MPS on a ring.
+# Sweeping from left to right, we optimize the tensors in the loop by minimizing the cost function.
 # Here cache of right-half-chain is used to minimize the number of multiplications to accelerate the sweeping. 
 # The transfer matrix on the left is updated after each optimization step.
 # The cache technique is from Chenfeng Bao's thesis, see http://hdl.handle.net/10012/14674.
-function loop_opt!(scheme::LoopTNR, loop_criterion::stopcrit,
-                   trunc::TensorKit.TruncationScheme,
-                   truncentanglement::TensorKit.TruncationScheme, verbosity::Int)
-    psiA = Ψ_A(scheme)
+function loop_opt(psiA::Array, loop_criterion::stopcrit,
+                  trunc::TensorKit.TruncationScheme,
+                  truncentanglement::TensorKit.TruncationScheme, verbosity::Int)
     psiB = Ψ_B(psiA, trunc, truncentanglement)
+    NB = length(psiB) # Number of tensors in the MPS Ψ_B
     psiBpsiB = ΨBΨB(psiB)
     psiBpsiA = ΨBΨA(psiB, psiA)
     psiApsiA = ΨAΨA(psiA)
     C = to_number(psiApsiA) # Since C is not changed during the optimization, we can compute it once and use it in the cost function.
-
     cost = Float64[Inf]
     sweep = 0
     crit = true
@@ -411,7 +251,19 @@ function loop_opt!(scheme::LoopTNR, loop_criterion::stopcrit,
         left_BA = id(codomain(psiBpsiA[1])) # Initialize the left transfer matrix for ΨBΨA
 
         t_start = time()
-        for pos_psiB in 1:8
+
+        if sweep == 0
+            tNt = tr(psiBpsiB[1]*right_cache_BB[1])
+            tdw = tr(psiBpsiA[1]*right_cache_BA[1])
+            wdt = conj(tdw)
+            cost_this = real((C + tNt - wdt - tdw) / C)
+            if verbosity > 1
+                @infov 3 "Initial cost: $cost_this"
+            end
+            push!(cost, cost_this)
+        end
+
+        for pos_psiB in 1:NB
             pos_psiA = (pos_psiB - 1) ÷ 2 + 1 # Position in the MPS Ψ_A
 
             N = tN(left_BB, right_cache_BB[pos_psiB]) # Compute the half of the matrix N for the current position in the loop, right cache is used to minimize the number of multiplications
@@ -434,31 +286,34 @@ function loop_opt!(scheme::LoopTNR, loop_criterion::stopcrit,
             end
         end
         sweep += 1
-        crit = loop_criterion(sweep, cost)
 
         tNt = tr(left_BB)
         tdw = tr(left_BA)
         wdt = conj(tdw)
         cost_this = real((C + tNt - wdt - tdw) / C)
         push!(cost, cost_this)
-
+        crit = loop_criterion(sweep, cost)
         if verbosity > 1
             @infov 3 "Sweep: $sweep, Cost: $(cost[end]), Time: $(time() - t_start)s" # Included the time taken for the sweep
         end
     end
 
+    return psiB
+end
+
+function loop_opt!(scheme::LoopTNR, loop_criterion::stopcrit,
+                   trunc::TensorKit.TruncationScheme,
+                   truncentanglement::TensorKit.TruncationScheme,
+                   verbosity::Int)
+    psiA = Ψ_A(scheme)
+    psiB = loop_opt(psiA, loop_criterion, trunc, truncentanglement, verbosity)
     @planar scheme.TB[-1 -2; -3 -4] := psiB[1][1; 2 -2] * psiB[4][-4; 2 3] *
                                        psiB[5][3; 4 -3] * psiB[8][-1; 4 1]
     @planar scheme.TA[-1 -2; -3 -4] := psiB[6][-2; 1 2] * psiB[7][2; 3 -4] *
                                        psiB[2][-3; 3 4] * psiB[3][4; 1 -1]
     return scheme
 end
 
-function loop_opt!(scheme::LoopTNR, trunc::TensorKit.TruncationScheme,
-                   verbosity::Int)
-    return loop_opt!(scheme, loop_criterion, trunc, verbosity)
-end
-
 function step!(scheme::LoopTNR, trunc::TensorKit.TruncationScheme,
                truncentanglement::TensorKit.TruncationScheme,
                entanglement_criterion::stopcrit,

src/schemes/symmetric_looptnr.jl

@@ -126,11 +126,14 @@ function entanglement_filtering(T; trunc=truncbelow(1e-12))
     psi_center = Ψ_center(T)
     psi_corner = Ψ_corner(T)
 
-    PR_list, PL_list = TNRKit.find_projectors(psi_center, entanglement_criterion, trunc)
+    PR_list, PL_list = TNRKit.find_projectors(psi_center, [1, 1, 1, 1], [3, 3, 3, 3],
+                                              entanglement_criterion, trunc)
     P_bottom = PL_list[1]
     P_right = PL_list[1]
 
-    PR_list, PL_list = TNRKit.find_projectors(psi_corner, entanglement_criterion, trunc)
+    PR_list, PL_list = TNRKit.find_projectors(psi_corner,
+                                              [1, 1, 1, 1], [3, 3, 3, 3],
+                                              entanglement_criterion, trunc)
     P_top = PL_list[3]
     P_left = PL_list[3]
 
@@ -157,7 +160,8 @@ function ef_oneloop(T, trunc::TensorKit.TruncationScheme)
 
     ΨB_function(steps, data) = abs(data[end])
     criterion = maxiter(100) & convcrit(1e-12, ΨB_function)
-    PR_list, _ = find_projectors(ΨB, criterion, trunc)
+    PR_list, _ = find_projectors(ΨB, [1, 1, 1, 1, 1, 1, 1, 1], [2, 2, 2, 2, 2, 2, 2, 2],
+                                 criterion, trunc)
 
     ΨB_disentangled = []
     for i in 1:1