Impurity HOTRG

docs/src/assets/tnrkit.bib

@@ -1,7 +1,7 @@
 @article{adachiAnisotropicTensorRenormalization2020,
   title     = {Anisotropic Tensor Renormalization Group},
   author    = {Adachi, Daiki and Okubo, Tsuyoshi and Todo, Synge},
-  year      = {2020},
+  year      = 2020,
   month     = aug,
   journal   = {Physical Review B},
   volume    = {102},
@@ -13,7 +13,7 @@ @article{adachiAnisotropicTensorRenormalization2020
 @article{adachiBondweightedTensorRenormalization2022,
   title     = {Bond-Weighted Tensor Renormalization Group},
   author    = {Adachi, Daiki and Okubo, Tsuyoshi and Todo, Synge},
-  year      = {2022},
+  year      = 2022,
   month     = feb,
   journal   = {Physical Review B},
   volume    = {105},
@@ -26,7 +26,7 @@ @article{adachiBondweightedTensorRenormalization2022
 @article{akiyamaTensorRenormalizationGroup2024a,
   title     = {Tensor Renormalization Group for Fermions},
   author    = {Akiyama, Shinichiro and Meurice, Yannick and Sakai, Ryo},
-  year      = {2024},
+  year      = 2024,
   month     = may,
   journal   = {Journal of Physics: Condensed Matter},
   volume    = {36},
@@ -40,7 +40,7 @@ @article{akiyamaTensorRenormalizationGroup2024a
 @article{levinTensorRenormalizationGroup2007,
   title     = {Tensor {{Renormalization Group Approach}} to {{Two-Dimensional Classical Lattice Models}}},
   author    = {Levin, Michael and Nave, Cody P.},
-  year      = {2007},
+  year      = 2007,
   month     = sep,
   journal   = {Physical Review Letters},
   volume    = {99},
@@ -50,10 +50,25 @@ @article{levinTensorRenormalizationGroup2007
   keywords  = {TNR Algorithm}
 }
 
+@article{moritaCalculationHigherorderMoments2019,
+  title         = {Calculation of Higher-Order Moments by Higher-Order Tensor Renormalization Group},
+  author        = {Morita, Satoshi and Kawashima, Naoki},
+  year          = 2019,
+  month         = mar,
+  journal       = {Computer Physics Communications},
+  volume        = {236},
+  eprint        = {1806.10275},
+  primaryclass  = {cond-mat},
+  issn          = {00104655},
+  doi           = {10.1016/j.cpc.2018.10.014},
+  archiveprefix = {arXiv},
+  keywords      = {Condensed Matter - Statistical Mechanics}
+}
+
 @article{moritaGlobalOptimizationTensor2021,
   title     = {Global Optimization of Tensor Renormalization Group Using the Corner Transfer Matrix},
   author    = {Morita, Satoshi and Kawashima, Naoki},
-  year      = {2021},
+  year      = 2021,
   month     = jan,
   journal   = {Physical Review B},
   volume    = {103},
@@ -65,7 +80,7 @@ @article{moritaGlobalOptimizationTensor2021
 @article{xieCoarsegrainingRenormalizationHigherorder2012,
   title     = {Coarse-Graining Renormalization by Higher-Order Singular Value Decomposition},
   author    = {Xie, Z. Y. and Chen, J. and Qin, M. P. and Zhu, J. W. and Yang, L. P. and Xiang, T.},
-  year      = {2012},
+  year      = 2012,
   month     = jul,
   journal   = {Physical Review B},
   volume    = {86},
@@ -77,7 +92,7 @@ @article{xieCoarsegrainingRenormalizationHigherorder2012
 @article{yangLoopOptimizationTensor2017,
   title     = {Loop {{Optimization}} for {{Tensor Network Renormalization}}},
   author    = {Yang, Shuo and Gu, Zheng-Cheng and Wen, Xiao-Gang},
-  year      = {2017},
+  year      = 2017,
   month     = mar,
   journal   = {Physical Review Letters},
   volume    = {118},

src/TNRKit.jl

@@ -20,7 +20,7 @@ include("schemes/hotrg.jl")
 include("schemes/hotrg3d.jl")
 include("schemes/atrg.jl")
 include("schemes/atrg3d.jl")
-
+include("schemes/impurityhotrg.jl")
 # CTM methods
 include("schemes/ctm/utility.jl")
 include("schemes/ctm/c4ctm.jl")
@@ -43,6 +43,7 @@ export HOTRG
 export HOTRG_3D
 export ATRG
 export ATRG_3D
+export ImpurityHOTRG
 
 export CTM
 export Sublattice_CTM
@@ -60,7 +61,7 @@ export run!
 # models
 include("models/ising.jl")
 export classical_ising, classical_ising_symmetric, ising_βc, f_onsager, ising_cft_exact,
-    ising_βc_3D, classical_ising_symmetric_3D, classical_ising_3D
+    ising_βc_3D, classical_ising_symmetric_3D, classical_ising_3D, classical_ising_impurity
 
 include("models/gross-neveu.jl")
 export gross_neveu_start
@@ -69,7 +70,7 @@ include("models/sixvertex.jl")
 export sixvertex
 
 include("models/potts.jl")
-export classical_potts, classical_potts_symmetric, potts_βc
+export classical_potts, classical_potts_symmetric, potts_βc, classical_potts_impurity
 
 include("models/clock.jl")
 export classical_clock

src/models/ising.jl

@@ -17,27 +17,22 @@ for the classical Ising model with a given inverse temperature `β` and external
     classical_ising() # Default inverse temperature is `ising_βc`
     classical_ising(0.5; h = 1.0) # Custom inverse temperature and magnetic field.
 ```
-!!! info
-    When calculating the free energy with `free_energy()`, set the `initial_size` keyword argument to `2.0`.
-    The initial lattice holds 2 spins.
 
 See also: [`classical_ising_symmetric`](@ref), [`classical_ising_symmetric_3D`](@ref), [`classical_ising_3D`](@ref).
 """
 function classical_ising(β::Number; h = 0)
-    function σ(i::Int64)
-        return 2i - 3
+    init = zeros(Float64, 2, 2, 2, 2)
+    for (i, j, k, l) in Iterators.product([1:2 for _ in 1:4]...)
+        init[i, j, k, l] = mod(i + j + k + l, 2) == 0 ? cosh(h * β) : sinh(h * β)
     end
+    init = TensorMap(init, ℂ^2 ⊗ ℂ^2 ← ℂ^2 ⊗ ℂ^2)
 
-    T_array = Float64[
-        exp(
-                β * (σ(i)σ(j) + σ(j)σ(l) + σ(l)σ(k) + σ(k)σ(i)) +
-                h / 2 * β * (σ(i) + σ(j) + σ(k) + σ(l))
-            )
-            for i in 1:2, j in 1:2, k in 1:2, l in 1:2
-    ]
-
-    T = TensorMap(T_array, ℂ^2 ⊗ ℂ^2 ← ℂ^2 ⊗ ℂ^2)
+    bond_tensor = zeros(Float64, 2, 2)
+    bond_tensor[1, 1] = sqrt(cosh(β))
+    bond_tensor[2, 2] = sqrt(sinh(β))
+    bond_tensor = TensorMap(bond_tensor, ℂ^2 ← ℂ^2)
 
+    @tensor T[-1 -2; -3 -4] := 2 * init[1 2; 3 4] * bond_tensor[-1; 1] * bond_tensor[-2; 2] * bond_tensor[3; -3] * bond_tensor[4; -4]
     return T
 end
 classical_ising() = classical_ising(ising_βc)
@@ -76,6 +71,40 @@ const f_onsager::BigFloat = -2.1096511446082074596677792835110847808254932754354
 """
 $(SIGNATURES)
 
+Constructs the partition function tensor for a 2D square lattice
+for the classical Ising model with a given inverse temperature `β` and external magnetic field `h` with a magnetisation impurity
+
+### Examples
+```julia
+    classical_ising_impurity() # Default inverse temperature is `ising_βc`
+    classical_ising_impurity(0.5; h = 1.0) # Custom inverse temperature and magnetic field
+```
+!!! info
+    When calculating the free energy with `free_energy()`, set the `initial_size` keyword argument to `2.0`.
+    The initial lattice holds 2 spins.
+
+See also: [`classical_ising_symmetric`](@ref), [`classical_ising_symmetric_3D`](@ref), [`classical_ising_3D`](@ref).
+"""
+function classical_ising_impurity(β::Number; h = 0)
+    init = zeros(Float64, 2, 2, 2, 2)
+    for (i, j, k, l) in Iterators.product([1:2 for _ in 1:4]...)
+        init[i, j, k, l] = mod(i + j + k + l, 2) == 0 ? sinh(h * β) : cosh(h * β)
+    end
+    init = TensorMap(init, ℂ^2 ⊗ ℂ^2 ← ℂ^2 ⊗ ℂ^2)
+
+    bond_tensor = zeros(Float64, 2, 2)
+    bond_tensor[1, 1] = sqrt(cosh(β))
+    bond_tensor[2, 2] = sqrt(sinh(β))
+    bond_tensor = TensorMap(bond_tensor, ℂ^2 ← ℂ^2)
+
+    @tensor T[-1 -2; -3 -4] := 2 * init[1 2; 3 4] * bond_tensor[-1; 1] * bond_tensor[-2; 2] * bond_tensor[3; -3] * bond_tensor[4; -4]
+    return T
+end
+classical_ising_impurity() = classical_ising_impurity(ising_βc)
+
+"""
+$(SIGNATURES)
+
 Constructs the partition function tensor for a symmetric 3D cubic lattice
 for the classical Ising model with a given inverse temperature `β`.
 

src/models/potts.jl

@@ -94,3 +94,43 @@ function classical_potts_symmetric(q::Int64, β::Float64)
     return A_potts
 end
 classical_potts_symmetric(q::Int) = classical_potts_symmetric(q, potts_βc(q))
+
+
+"""
+$(SIGNATURES)
+
+Constructs the partition function tensor for a Potts model with `q` states
+and a given inverse temperature `β` with impurities in sectors `k1` and `k2`.
+
+### Examples
+```julia
+    classical_potts_impurity(3) # Default inverse temperature is `potts_βc(3)`
+    classical_potts_impurity(3, 1, 2) # Custom inverse temperature and impurity sectors.
+    classical_potts_impurity(3, 0.5, 2, 3) # Custom inverse temperature and impurity sectors.
+```
+
+See also: [`classical_potts`](@ref), [`potts_βc`](@ref).
+"""
+function classical_potts_impurity(q::Int64, β::Float64, k1::Int64 = 1, k2::Int64 = 1)
+    bond_tensor = zeros(Float64, q, q)
+    for i in 0:(q - 1)
+        bond_tensor[i + 1, i + 1] = sqrt(exp(β) - 1 + q * (i == 0))
+    end
+    Vp = ℂ^q
+    bond_tensor = TensorMap(bond_tensor, Vp ← Vp)
+
+    core_tensor = zeros(Float64, q, q, q, q)
+    for (i, j, k, l) in Iterators.product(0:(q - 1), 0:(q - 1), 0:(q - 1), 0:(q - 1))
+        core_tensor[i + 1, j + 1, k + 1, l + 1] =
+            mod(i + j - k - l + k1 - k2, q) == 0 ? 1 : 0
+    end
+
+    core_tensor = TensorMap(core_tensor, Vp ⊗ Vp ← Vp ⊗ Vp)
+
+    @tensor T[-1 -2; -3 -4] := core_tensor[1 2; 3 4] * bond_tensor[-1; 1] * bond_tensor[-2; 2] * bond_tensor[3; -3] * bond_tensor[4; -4] * (1 / q)
+
+    return T
+end
+
+classical_potts_impurity(q::Int64, k1::Int64, k2::Int64) = classical_potts_impurity(q, potts_βc(q), k1, k2)
+classical_potts_impurity(q::Int64) = classical_potts_impurity(q, potts_βc(q), 1, 1)

src/schemes/hotrg.jl

@@ -39,6 +39,49 @@ in order to reuse the y-compression code for x-compression.
 Hence both are written explicitly.
 =#
 
+
+function _step_hotrg_x(
+        A1::TensorMap{E, S, 2, 2}, A2::TensorMap{E, S, 2, 2},
+        U::TensorMap{E, S, 2, 1}
+    ) where {E, S}
+    #= compression along the x-direction
+                -3
+                |
+            ┌3--U--4┐
+            |       |
+        -1--A1--5---A2-- -4
+            |       |
+            └1--U†-2┘
+                |
+                -2
+    =#
+
+    @tensor T[-1 -2; -3 -4] :=
+        A1[-1 1; 3 5] * A2[5 2; 4 -4] * conj(U[1 2; -2]) * U[3 4; -3]
+    return T
+end
+
+function _step_hotrg_y(
+        A1::TensorMap{E, S, 2, 2}, A2::TensorMap{E, S, 2, 2},
+        U::TensorMap{E, S, 2, 1}
+    ) where {E, S}
+    #= compression along the y-direction
+                    -3
+                    |
+            ┌---1---A2---3--┐
+            |       |       |
+        -1--U†      5       U-- -4
+            |       |       |
+            └---2---A1---4--┘
+                    |
+                    -2
+    =#
+
+    @tensor T[-1 -2; -3 -4] :=
+        conj(U[1 2; -1]) * U[3 4; -4] * A2[1 5; -3 3] * A1[2 -2; 5 4]
+    return T
+end
+
 function _get_hotrg_xproj(
         A1::AbstractTensorMap{E, S, 2, 2}, A2::AbstractTensorMap{E, S, 2, 2},
         trunc::TensorKit.TruncationScheme

src/schemes/impurityhotrg.jl

@@ -0,0 +1,96 @@
+"""
+$(TYPEDEF)
+
+Single impurity method for Higher-Order Tensor Renormalization Group (for 2nd order)
+
+### Constructors
+    $(FUNCTIONNAME)(T, T_imp_order1_1, T_imp_order1_2, T_imp_order2 [, finalize=finalize!])
+
+### Running the algorithm
+    run!(::ImpurityHOTRG, trunc::TensorKit.TruncationSheme, stop::Stopcrit[, finalize_beginning=true, verbosity=1])
+
+Each step rescales the lattice by a (linear) factor of 2
+
+!!! info "verbosity levels"
+    - 0: No output
+    - 1: Print information at start and end of the algorithm
+    - 2: Print information at each step
+
+### Fields
+
+$(TYPEDFIELDS)
+
+### References
+* [Morita et al 10.1016/j.cpc.2018.10.014 (2018)](@cite moritaCalculationHigherorderMoments2019)
+
+"""
+
+mutable struct ImpurityHOTRG <: TNRScheme
+    T::TensorMap
+    T_imp_order1_1::TensorMap
+    T_imp_order1_2::TensorMap
+    T_imp_order2::TensorMap
+    finalize!::Function
+    function ImpurityHOTRG(
+            T::TensorMap{E, S, 2, 2},
+            T_imp_order1_1::TensorMap{E, S, 2, 2},
+            T_imp_order1_2::TensorMap{E, S, 2, 2},
+            T_imp_order2::TensorMap{E, S, 2, 2},
+            ;
+            finalize = (finalize!),
+        ) where {E, S}
+
+        @assert space(T, 1) == space(T_imp_order1_1, 1) == space(T_imp_order1_2, 1) "First space of T, T_imp_order1_1 and T_imp_order1_2 must be the same"
+        @assert space(T, 2) == space(T_imp_order1_1, 2) == space(T_imp_order1_2, 2) "Second space of T, T_imp_order1_1 and T_imp_order1_2 must be the same"
+        @assert space(T, 3) == space(T_imp_order1_1, 3) == space(T_imp_order1_2, 3) "Third space of T, T_imp_order1_1 and T_imp_order1_2 must be the same"
+        @assert space(T, 4) == space(T_imp_order1_1, 4) == space(T_imp_order1_2, 4) "Fourth space of T, T_imp_order1_1 and T_imp_order1_2 must be the same"
+        return new(T, T_imp_order1_1, T_imp_order1_2, T_imp_order2, finalize)
+    end
+end
+
+function step!(scheme::ImpurityHOTRG, trunc::TensorKit.TruncationScheme)
+    Uy, _ = _get_hotrg_yproj(scheme.T, scheme.T, trunc)
+
+    T = _step_hotrg_x(scheme.T, scheme.T, Uy)
+    T_imp_order1_1 = 0.5 * (_step_hotrg_x(scheme.T_imp_order1_1, scheme.T, Uy) + _step_hotrg_x(scheme.T, scheme.T_imp_order1_1, Uy))
+    T_imp_order1_2 = 0.5 * (_step_hotrg_x(scheme.T_imp_order1_2, scheme.T, Uy) + _step_hotrg_x(scheme.T, scheme.T_imp_order1_2, Uy))
+    T_imp_order2 = 0.25 * (
+        _step_hotrg_x(scheme.T_imp_order2, scheme.T, Uy) +
+            _step_hotrg_x(scheme.T, scheme.T_imp_order2, Uy) +
+            _step_hotrg_x(scheme.T_imp_order1_1, scheme.T_imp_order1_2, Uy) +
+            _step_hotrg_x(scheme.T_imp_order1_2, scheme.T_imp_order1_1, Uy)
+    )
+
+    scheme.T = T
+    scheme.T_imp_order1_1 = T_imp_order1_1
+    scheme.T_imp_order1_2 = T_imp_order1_2
+    scheme.T_imp_order2 = T_imp_order2
+
+    Ux, _ = _get_hotrg_xproj(scheme.T, scheme.T, trunc)
+
+    T = _step_hotrg_y(scheme.T, scheme.T, Ux)
+    T_imp_order1_1 = 0.5 * (_step_hotrg_y(scheme.T_imp_order1_1, scheme.T, Ux) + _step_hotrg_y(scheme.T, scheme.T_imp_order1_1, Ux))
+    T_imp_order1_2 = 0.5 * (_step_hotrg_y(scheme.T_imp_order1_2, scheme.T, Ux) + _step_hotrg_y(scheme.T, scheme.T_imp_order1_2, Ux))
+    T_imp_order2 = 0.25 * (
+        _step_hotrg_y(scheme.T_imp_order2, scheme.T, Ux) +
+            _step_hotrg_y(scheme.T, scheme.T_imp_order2, Ux) +
+            _step_hotrg_y(scheme.T_imp_order1_1, scheme.T_imp_order1_2, Ux) +
+            _step_hotrg_y(scheme.T_imp_order1_2, scheme.T_imp_order1_1, Ux)
+    )
+
+    scheme.T = T
+    scheme.T_imp_order1_1 = T_imp_order1_1
+    scheme.T_imp_order1_2 = T_imp_order1_2
+    scheme.T_imp_order2 = T_imp_order2
+
+    return scheme
+end
+
+function Base.show(io::IO, scheme::ImpurityHOTRG)
+    println(io, "ImpurityHOTRG - Impurity Higher Order TRG")
+    println(io, "  * T: $(summary(scheme.T))")
+    println(io, "  * T_imp_order1_1: $(summary(scheme.T_imp_order1_1))")
+    println(io, "  * T_imp_order1_2: $(summary(scheme.T_imp_order1_2))")
+    println(io, "  * T_imp_order2: $(summary(scheme.T_imp_order2))")
+    return nothing
+end