Type Stability and general improvements attempt (#118)

* remove PEPSKit dependency

* comment on twistdual

* formatting

* first experiments

* type annotation

* changing types

* type almost all of the ctm methods

* fix bug in in BTRG

* use eltype of central tensor for bond tensor initialization

* k is a float in BTRG

* fix type in HOTRG constructor

* fix constructor slooptnr

* fix isreal bug

* fix order of operations

* attempt fix for ImpurityHOTRG

* same fix for impurityTRG

---------

Co-authored-by: sanderdemeyer <sander.demeyer@ugent.be>
Co-authored-by: Sander De Meyer <74001142+sanderdemeyer@users.noreply.github.com>

Author: VictorVanthilt

Project.toml

@@ -9,20 +9,17 @@ KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LoggingExtras = "e6f89c97-d47a-5376-807f-9c37f3926c36"
 OptimKit = "77e91f04-9b3b-57a6-a776-40b61faaebe0"
-PEPSKit = "52969e89-939e-4361-9b68-9bc7cde4bdeb"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [compat]
 DocStringExtensions = "0.9.4"
-KrylovKit = "0.9.5"
 LoggingExtras = "~1.0"
 OptimKit = "0.4.0"
-PEPSKit = "0.6"
-TensorKit = "0.14"
 Zygote = "0.7.7"
 julia = "1.11"
+TensorKit = "0.14"
 
 [extras]
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"

src/TNRKit.jl

@@ -3,8 +3,6 @@ using TensorKit, LinearAlgebra
 using LoggingExtras, Printf
 using KrylovKit
 using OptimKit, Zygote
-using PEPSKit: InfinitePartitionFunction, CTMRGEnv
-using PEPSKit: network_value, twistdual, twistdual!
 using DocStringExtensions
 
 # stop criteria
@@ -97,4 +95,6 @@ include("utility/projectors.jl")
 
 include("utility/blocking.jl")
 export block_tensors
+
+include("utility/network_value.jl")
 end

src/schemes/atrg.jl

@@ -23,11 +23,12 @@ $(TYPEDFIELDS)
 ### References
 * [Adachi et. al. Phys. Rev. B 102 (2020)](@cite adachiAnisotropicTensorRenormalization2020)
 """
-mutable struct ATRG <: TNRScheme
-    T::TensorMap
+mutable struct ATRG{E, S, TT <: AbstractTensorMap{E, S, 2, 2}} <: TNRScheme{E, S}
+    "Central tensor"
+    T::TT
 
-    function ATRG(T::TensorMap{E, S, 2, 2}) where {E, S}
-        return new(T)
+    function ATRG(T::TT) where {E, S, TT <: AbstractTensorMap{E, S, 2, 2}}
+        return new{E, S, TT}(T)
     end
 end
 

src/schemes/atrg3d.jl

@@ -23,11 +23,12 @@ $(TYPEDFIELDS)
 ### References
 * [Adachi et. al. Phys. Rev. B 102 (2020)](@cite adachiAnisotropicTensorRenormalization2020)
 """
-mutable struct ATRG_3D <: TNRScheme
-    T::TensorMap
+mutable struct ATRG_3D{E, S, TT <: AbstractTensorMap{E, S, 2, 4}} <: TNRScheme{E, S}
+    "Central tensor"
+    T::TT
 
-    function ATRG_3D(T::TensorMap{E, S, 2, 4}) where {E, S}
-        return new(T)
+    function ATRG_3D(T::TT) where {E, S, TT <: AbstractTensorMap{E, S, 2, 4}}
+        return new{E, S, TT}(T)
     end
 end
 

src/schemes/btrg.jl

@@ -23,20 +23,29 @@ $(TYPEDFIELDS)
 ### References
 * [Adachi et. al. Phys. Rev. B 105 (2022)](@cite adachiBondweightedTensorRenormalization2022)
 """
-mutable struct BTRG <: TNRScheme
-    T::TensorMap
-    S1::TensorMap
-    S2::TensorMap
+mutable struct BTRG{E, S, TT <: AbstractTensorMap{E, S, 2, 2}, BT <: AbstractTensorMap{E, S, 1, 1}} <: TNRScheme{E, S}
+    "Central tensor"
+    T::TT
+
+    "Bond tensor on vertical bonds"
+    S1::BT
+
+    "Bond tensor on horizontal bonds"
+    S2::BT
+
+    "Bond weight exponent"
     k::Float64
 
-    function BTRG(T::TensorMap{E, S, 2, 2}, k::Number) where {E, S}
+    function BTRG(T::TT, k::Number) where {E, S, TT <: AbstractTensorMap{E, S, 2, 2}}
         # Construct S1 and S2 as identity matrices.
-        return new(T, id(space(T, 2)), id(space(T, 1)), k)
+        S1 = id(E, space(T, 2))
+        S2 = id(E, space(T, 1))
+        return new{E, S, TT, typeof(S1)}(T, S1, S2, k)
     end
 end
 
 # Default implementation using the optimal value for k
-BTRG(T::TensorMap) = BTRG(T, -0.5)
+BTRG(T) = BTRG(T, -0.5)
 
 function pseudopow(t::DiagonalTensorMap, a::Real; tol = eps(scalartype(t))^(3 / 4))
     t′ = copy(t)

src/schemes/ctm/c4vctm.jl

@@ -23,15 +23,20 @@ The keyword argument symmetrize makes the tensor C4v symmetric when set to true.
 
 $(TYPEDFIELDS)
 """
-mutable struct c4vCTM{A, S}
-    T::TensorMap{A, S, 0, 4}
-    C::TensorMap{A, S, 1, 1}
-    E::TensorMap{A, S, 2, 1}
+mutable struct c4vCTM{A, S, TT <: AbstractTensorMap{A, S, 0, 4}, TC <: AbstractTensorMap{A, S, 1, 1}, TE <: AbstractTensorMap{A, S, 2, 1}} <: TNRScheme{A, S}
+    "Central tensor"
+    T::TT
 
-    function c4vCTM(T::TensorMap{A, S, 0, 4}) where {A, S}
+    "Corner tensor"
+    C::TC
+
+    "Edge tensor"
+    E::TE
+
+    function c4vCTM(T::TT) where {A, S, TT <: AbstractTensorMap{A, S, 0, 4}}
         C, E = c4vCTM_init(T)
 
-        return new{A, S}(T, C, E)
+        return new{A, S, TT, typeof(C), typeof(E)}(T, C, E)
     end
 end
 
@@ -113,7 +118,7 @@ function step!(scheme::c4vCTM, trunc)
     mat, U, S = find_U_sym(scheme, trunc)
 
     @tensor scheme.C[-1; -2] := mat[1 2; 3 4] * U[3 4; -2] * conj(U[1 2; -1])
-    @tensor scheme.E[-1 -2; -3] := scheme.E[1 5; 3] * flip(scheme.T, (3, 4); inv = true)[5 4 -2 2] *
+    @tensor scheme.E[-1 -2; -3] := scheme.E[1 5; 3] * flip(scheme.T, (3, 4); inv = false)[5 4 -2 2] *
         U[3 4; -3] *
         conj(U[1 2; -1])
 
@@ -123,15 +128,21 @@ function step!(scheme::c4vCTM, trunc)
 end
 
 function lnz(scheme::c4vCTM)
-    Z, env = tensor2env(permute(flip(scheme.T, (3, 4); inv = true), ((4, 3), (1, 2))), scheme.C, scheme.E)
-    # should be inv = false ??
-    return real(log(network_value(Z, env)))
+    T_unflipped = permute(flip(scheme.T, (3, 4); inv = false), ((4, 3), (1, 2)))
+
+    corners = fill(scheme.C, 4)
+    edges = fill(scheme.E, 4)
+
+    edges[3] = flip(scheme.E, 2; inv = true)
+    edges[4] = flip(scheme.E, 2; inv = true)
+
+    return real(log(network_value(T_unflipped, corners, edges)))
 end
 
 function build_corner_matrix(scheme)
     @tensor opt = true mat[-1 -2; -3 -4] := scheme.C[1; 2] * scheme.E[-1 3; 1] *
         scheme.E[2 4; -3] *
-        flip(scheme.T, 3; inv = true)[4 -4 -2 3]
+        flip(scheme.T, 3; inv = false)[4 -4 -2 3]
     return mat
 end
 
@@ -152,20 +163,6 @@ function c4vCTM_init(T::TensorMap{A, S, 0, 4}) where {A, S}
     return C, E
 end
 
-function tensor2env(O, C, T)
-    Z = InfinitePartitionFunction(O)
-    env = CTMRGEnv(Z, space(C)[1])
-
-    for i in 1:4
-        env.corners[i] = C
-        env.edges[i] = T
-    end
-
-    env.edges[3] = flip(T, 2)
-    env.edges[4] = flip(T, 2)
-    return Z, env
-end
-
 function Base.show(io::IO, scheme::c4vCTM)
     println(io, "c4vCTM - C4v symmetric Corner Transfer Matrix")
     println(io, "  * T: $(summary(scheme.T))")

src/schemes/ctm/ctm_hotrg.jl

@@ -1,18 +1,18 @@
-mutable struct ctm_HOTRG{A, S} <: TNRScheme
-    T::TensorMap{A, S, 2, 2}
-    C2::TensorMap{A, S, 1, 1}
-    E1::TensorMap{A, S, 2, 1}
-    E2::TensorMap{A, S, 2, 1}
+mutable struct ctm_HOTRG{E, S, TT <: AbstractTensorMap{E, S, 2, 2}, TC <: AbstractTensorMap{E, S, 1, 1}, TE <: AbstractTensorMap{E, S, 2, 1}} <: TNRScheme{E, S}
+    T::TT
+    C2::TC
+    E1::TE
+    E2::TE
     χenv::Int64
     function ctm_HOTRG(
-            T::TensorMap{A, S, 2, 2},
+            T::TT,
             χenv::Int64;
             ctm_iter = maxiter(2.0e4),
             ctm_tol = trivial_convcrit(1.0e-9),
             ctm_obc = false,
             χenv_ini = 2,
-        ) where {A, S}
-        if eltype(T) != Float64
+        ) where {E, S, TT <: AbstractTensorMap{E, S, 2, 2}}
+        if E != Float64
             @error "This scheme only supports tensors with real numbers"
         end
         scheme_init = TNRKit.rCTM(T)
@@ -27,7 +27,7 @@ mutable struct ctm_HOTRG{A, S} <: TNRScheme
         @info "rCTM finished"
         C2, E1, E2 = scheme_init.C2, scheme_init.E1, scheme_init.E2
         @assert BraidingStyle(sectortype(T)) == Bosonic() "$(summary(BraidingStyle(sectortype(T)))) braiding style is not supported for rCTM"
-        return new{A, S}(T, C2, E1, E2, χenv)
+        return new{E, S, TT, typeof(C2), typeof(E1)}(T, C2, E1, E2, χenv)
     end
 end
 

src/schemes/ctm/ctm_trg.jl

@@ -1,15 +1,15 @@
-mutable struct ctm_TRG{A, S} <: TNRScheme
-    T::TensorMap{A, S, 2, 2}
-    C2::TensorMap{A, S, 1, 1}
-    E1::TensorMap{A, S, 2, 1}
-    E2::TensorMap{A, S, 2, 1}
+mutable struct ctm_TRG{E, S, TT <: AbstractTensorMap{E, S, 2, 2}, TC <: AbstractTensorMap{E, S, 1, 1}, TE <: AbstractTensorMap{E, S, 2, 1}} <: TNRScheme{E, S}
+    T::TT
+    C2::TC
+    E1::TE
+    E2::TE
     χenv::Int64
     function ctm_TRG(
-            T::TensorMap{A, S, 2, 2},
+            T::TT,
             χenv::Int64;
             ctm_iter = 2.0e4,
             ctm_tol = 1.0e-9,
-        ) where {A, S}
+        ) where {E, S, TT <: AbstractTensorMap{E, S, 2, 2}}
         if eltype(T) != Float64
             @error "This scheme only support tensors with real numbers"
         end
@@ -24,7 +24,7 @@ mutable struct ctm_TRG{A, S} <: TNRScheme
         @info "rCTM finished"
         C2, E1, E2 = scheme_init.C2, scheme_init.E1, scheme_init.E2
         @assert BraidingStyle(sectortype(T)) == Bosonic() "$(summary(BraidingStyle(sectortype(T)))) braiding style is not supported for rCTM"
-        return new{A, S}(T, C2, E1, E2, χenv)
+        return new{E, S, TT, typeof(C2), typeof(E1)}(T, C2, E1, E2, χenv)
     end
 end
 

src/schemes/ctm/onesite_ctm.jl

@@ -9,18 +9,20 @@
             │Cbl├─┤Eb├─┤Cbr│
             └───┘ └──┘ └───┘
 """
-mutable struct CTM
-    T::AbstractTensorMap
-    Ctl::AbstractTensorMap
-    Ctr::AbstractTensorMap
-    Cbr::AbstractTensorMap
-    Cbl::AbstractTensorMap
-    El::AbstractTensorMap
-    Eb::AbstractTensorMap
-    Et::AbstractTensorMap
-    Er::AbstractTensorMap
+mutable struct CTM{E, S, TT <: AbstractTensorMap{E, S, 2, 2}, TC <: AbstractTensorMap{E, S, 1, 1}, TE <: AbstractTensorMap{E, S, 2, 1}} <: TNRScheme{E, S}
+    T::TT
+    Ctl::TC
+    Ctr::TC
+    Cbr::TC
+    Cbl::TC
+    El::TE
+    Eb::TE
+    Et::TE
+    Er::TE
 end
 
+#TODO: type everything
+
 CTM(T; bc = ones, bc_free = false) = CTM(T, CTM_init(T; bc, bc_free)...)
 
 function lnz(ctm::CTM)
@@ -34,14 +36,17 @@ function lnz(ctm::CTM)
         ctm.Eb[6 10; 7] *
         ctm.Cbl[7; 8] *
         ctm.El[8 9; 1]
+
     B = tr(ρA(ctm))
+
     @tensor opt = true C =
         ctm.Ctl[1; 2] *
         ctm.Et[2 7; 3] *
         ctm.Ctr[3; 4] *
         ctm.Cbr[4; 5] *
         ctm.Eb[5 7; 6] *
         ctm.Cbl[6; 1]
+
     @tensor opt = true D =
         ctm.Ctl[1; 2] *
         ctm.Ctr[2; 3] *

src/schemes/ctm/rctm.jl

@@ -15,19 +15,19 @@ The corner tensors are related by its mirror images.
                │  
                ▼  
 =#
-mutable struct rCTM{A, S}
-    T::TensorMap{A, S, 2, 2}
-    C2::TensorMap{A, S, 1, 1}
-    E1::TensorMap{A, S, 2, 1}
-    E2::TensorMap{A, S, 2, 1}
-
-    function rCTM(T::TensorMap{A, S, 2, 2}) where {A, S}
-        if typeof(T.data[1]) != Float64
+mutable struct rCTM{E, S, TT <: AbstractTensorMap{E, S, 2, 2}, TC <: AbstractTensorMap{E, S, 1, 1}, TE <: AbstractTensorMap{E, S, 2, 1}} <: TNRScheme{E, S}
+    T::TT
+    C2::TC
+    E1::TE
+    E2::TE
+
+    function rCTM(T::TT) where {E, S, TT <: AbstractTensorMap{E, S, 2, 2}}
+        if !(E <: Real)
             @error "This scheme only support tensors with real numbers"
         end
         C, E1, E2 = rCTM_init(T)
         @assert BraidingStyle(sectortype(T)) == Bosonic() "$(summary(BraidingStyle(sectortype(T)))) braiding style is not supported for rCTM"
-        return new{A, S}(T, C, E1, E2)
+        return new{E, S, TT, typeof(C), typeof(E1)}(T, C, E1, E2)
     end
 end
 
@@ -103,22 +103,10 @@ function run!(
 end
 
 function lnz(scheme::rCTM)
-    Z, env = tensor2env(scheme.T, scheme.C2, scheme.E1, scheme.E2)
-    return real(log(network_value(Z, env)))
-end
+    corners = [adjoint(scheme.C2), scheme.C2, adjoint(scheme.C2), scheme.C2]
+    edges = [scheme.E1 scheme.E2 flip(scheme.E1, 2) flip(scheme.E2, 2)]
 
-function tensor2env(T, C2, E1, E2)
-    Z = InfinitePartitionFunction(T)
-    env = CTMRGEnv(Z, space(C2)[1])
-    for i in 1:2
-        env.corners[2 * i] = C2
-        env.edges[2 * i] = E2
-        env.corners[2 * i - 1] = adjoint(C2)
-        env.edges[2 * i - 1] = E1
-    end
-    env.edges[3] = flip(env.edges[3], 2)
-    env.edges[4] = flip(env.edges[4], 2)
-    return Z, env
+    return real(log(network_value(scheme.T, corners, edges)))
 end
 
 function Base.show(io::IO, scheme::rCTM)