diff --git a/docs/src/assets/tnrkit.bib b/docs/src/assets/tnrkit.bib index 9a509c5c..98e6bb59 100644 --- a/docs/src/assets/tnrkit.bib +++ b/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}, diff --git a/src/TNRKit.jl b/src/TNRKit.jl index 09eb5bb3..a6c4b36f 100644 --- a/src/TNRKit.jl +++ b/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 diff --git a/src/models/ising.jl b/src/models/ising.jl index 9a69bb7e..58905db2 100644 --- a/src/models/ising.jl +++ b/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 `β`. diff --git a/src/models/potts.jl b/src/models/potts.jl index 640ef86d..4447c115 100644 --- a/src/models/potts.jl +++ b/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) diff --git a/src/schemes/hotrg.jl b/src/schemes/hotrg.jl index 8d8aebe7..bc7426c6 100644 --- a/src/schemes/hotrg.jl +++ b/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 diff --git a/src/schemes/impurityhotrg.jl b/src/schemes/impurityhotrg.jl new file mode 100644 index 00000000..e32f8ac2 --- /dev/null +++ b/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 diff --git a/src/utility/finalize.jl b/src/utility/finalize.jl index d8f3529b..c9435758 100644 --- a/src/utility/finalize.jl +++ b/src/utility/finalize.jl @@ -71,6 +71,19 @@ function finalize!(scheme::SLoopTNR) return tr_norm^0.25 end +# finalize! for ImpurityHOTRG +function finalize!(scheme::ImpurityHOTRG) + n = norm(@tensor scheme.T[1 2; 2 1]) + n_11 = norm(@tensor scheme.T_imp_order1_1[1 2; 2 1]) + n_12 = norm(@tensor scheme.T_imp_order1_2[1 2; 2 1]) + n_2 = norm(@tensor scheme.T_imp_order2[1 2; 2 1]) + scheme.T /= n + scheme.T_imp_order1_1 /= n + scheme.T_imp_order1_2 /= n + scheme.T_imp_order2 /= n + return n, n_11, n_12, n_2 +end + # cft data finalize function finalize_cftdata!(scheme::LoopTNR) finalize!(scheme) diff --git a/src/utility/free_energy.jl b/src/utility/free_energy.jl index 316dd523..65604d04 100644 --- a/src/utility/free_energy.jl +++ b/src/utility/free_energy.jl @@ -7,7 +7,7 @@ and computes the free energy. !!! info The `scalefactor` should be set to the rescaling factor of the area of the tensor network after each iteration of the TNR algorithm. - The `initial_size` should be set to the intial size of the physical lattice, which is typically `1.0` (or `2.0` for the 2D non-symmetric classical Ising model). + The `initial_size` should be set to the intial size of the physical lattice, which is typically `1.0`. """ function free_energy(data, β; scalefactor = 2.0, initial_size = 1.0) diff --git a/test/models.jl b/test/models.jl index 08634e43..ae986fc5 100644 --- a/test/models.jl +++ b/test/models.jl @@ -36,12 +36,10 @@ answers = [ 3 / 2 * log(3 / 4), ] -lattice_sizes = vcat([2.0], fill(1.0, length(models_2D) - 1)) - @testset "2D Models" begin - for (model, temp, answer, unitcell) in zip(models_2D, temperatures, answers, lattice_sizes) + for (model, temp, answer) in zip(models_2D, temperatures, answers) scheme = TRG(model) data = run!(scheme, truncdim(16), maxiter(25)) - @test free_energy(data, temp; initial_size = unitcell) ≈ answer rtol = 1.0e-3 + @test free_energy(data, temp) ≈ answer rtol = 1.0e-3 end end diff --git a/test/schemes.jl b/test/schemes.jl index 7f5aa78c..366e62f9 100644 --- a/test/schemes.jl +++ b/test/schemes.jl @@ -194,3 +194,46 @@ end @info "Calculated f = $(fs)." @test fs ≈ f_benchmark3D rtol = 1.0e-3 end + +# ImpurityHOTRG +@testset "ImpurityHOTRG - Ising Model" begin + + T = classical_ising() + T_imp1 = classical_ising_impurity() + + scheme = ImpurityHOTRG(T, T_imp1, T_imp1, T) + + data = run!(scheme, truncdim(16), maxiter(25)) + + @test free_energy(getindex.(data, 1), ising_βc; scalefactor = 4.0) ≈ f_onsager rtol = 6.0e-7 +end + +@testset "Impurity HOTRG - Magnetisation" begin + # High temperature limit ( -> 0) + β = 0.2 + + T = classical_ising(β) + T_imp_order1_1 = classical_ising_impurity(β) + T_imp_order2 = classical_ising(β) + + scheme = ImpurityHOTRG(T, T_imp_order1_1, T_imp_order1_1, T_imp_order2) + + data = run!(scheme, truncdim(8), maxiter(25)) + + m2_highT = data[end][4] / data[end][1] + @test m2_highT ≈ 0.0 atol = 1.0e-14 + + # Low temperature limit ( -> 1) + β = 1.0 + + T = classical_ising(β) + T_imp_order1_1 = classical_ising_impurity(β) + T_imp_order2 = classical_ising(β) + + scheme = ImpurityHOTRG(T, T_imp_order1_1, T_imp_order1_1, T_imp_order2) + + data = run!(scheme, truncdim(8), maxiter(25)) + + m2_lowT = data[end][4] / data[end][1] + @test m2_lowT ≈ 1 rtol = 1.0e-2 +end