diff --git a/Project.toml b/Project.toml index b0ac7086..2a380006 100644 --- a/Project.toml +++ b/Project.toml @@ -11,6 +11,7 @@ KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LoggingExtras = "e6f89c97-d47a-5376-807f-9c37f3926c36" MatrixAlgebraKit = "6c742aac-3347-4629-af66-fc926824e5e4" +NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" OptimKit = "77e91f04-9b3b-57a6-a776-40b61faaebe0" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" @@ -28,6 +29,7 @@ FastGaussQuadrature = "1.1.0" KrylovKit = "0.10" LoggingExtras = "~1.0" MatrixAlgebraKit = "0.6.1" +NonlinearSolve = "4.12.0" OptimKit = "0.4.0" QuadGK = "2.11.2" Roots = "3.0.0" diff --git a/src/TNRKit.jl b/src/TNRKit.jl index 3c6243d8..53a0ebc6 100644 --- a/src/TNRKit.jl +++ b/src/TNRKit.jl @@ -11,6 +11,7 @@ using SpecialFunctions using FastGaussQuadrature using QuadGK using Roots +using NonlinearSolve using Base.Threads using Combinatorics: permutations diff --git a/src/schemes/looptnr.jl b/src/schemes/looptnr.jl index e1d1140b..8c5db9f1 100644 --- a/src/schemes/looptnr.jl +++ b/src/schemes/looptnr.jl @@ -433,7 +433,9 @@ function loop_opt( end """ -Coarse-grain `ΨB` to renormalized `TA`, `TB` tensors +Coarse-grain `ΨB` to renormalized `TA`, `TB` tensors. +The lattice is rotated by 135 degrees in counter clockwise direction. +The elementary modular parameter `τ₀ ↦ (1 + τ₀) / (1 - τ₀)`. """ function ΨB_to_TATB(psiB::Vector{T}) where {T <: AbstractTensorMap{<:Any, <:Any, 1, 2}} #= diff --git a/src/schemes/trg.jl b/src/schemes/trg.jl index 021379e4..af1c9f43 100644 --- a/src/schemes/trg.jl +++ b/src/schemes/trg.jl @@ -9,7 +9,9 @@ Tensor Renormalization Group # Running the algorithm run!(::TRG, trunc::TruncationStrategy, stop::Stopcrit[, finalizer=default_Finalizer, finalize_beginning=true, verbosity=1]) -Each step rescales the lattice by a (linear) factor of √2 +Each step rescales the lattice by a (linear) factor of √2, +and rotate the lattice by 45 degrees in counter clockwise direction. +The elementary modular parameter `τ₀ ↦ (τ₀ - 1) / (τ₀ + 1)`. !!! info "verbosity levels" - 0: No output diff --git a/src/utility/cft.jl b/src/utility/cft.jl index e7438470..cab36eb9 100644 --- a/src/utility/cft.jl +++ b/src/utility/cft.jl @@ -1,5 +1,5 @@ """ - CFTData{E, I} where {E, I} + struct CFTData{E, K, V, A <: AbstractVector{E}} A struct to hold conformal data extracted from a TNR scheme. @@ -10,13 +10,15 @@ A struct to hold conformal data extracted from a TNR scheme. # Fields - `central_charge::Union{E, Missing}`: The central charge of the CFT. Will be `nothing` if not calculated. + - `modular_parameter::E`: The elementary modular parameter of a square patch of the lattice network. - `scaling_dimensions::StructuredVector{E, K, V, A}`: The scaling dimensions of the CFT, organized in a `StructuredVector` where the sectors correspond to different spin sectors (or other quantum numbers) and the data contains the scaling dimensions within those sectors """ struct CFTData{E, K, V, A <: AbstractVector{E}} - "Central charge of the CFT. Will be `nothing` if not calculated." + "Central charge of the CFT. Will be `missing` if not calculated." central_charge::Union{E, Missing} - + "Elementary modular parameter for one tensor" + modular_parameter::E "Scaling dimensions of the CFT." scaling_dimensions::StructuredVector{E, K, V, A} end @@ -28,47 +30,52 @@ function Base.show(io::IO, data::CFTData) return nothing end -function CFTData(T::TensorMap{E, S, 2, 2}; shape = [sqrt(2), 2 * sqrt(2), 0], kwargs...) where {E, S} - if shape == [1, 1, 0] # trivial implementation - scaling_dimensions = _scaling_dimensions(T) - return CFTData(missing, StructuredVector(scaling_dimensions, Dict([Trivial => collect(eachindex(scaling_dimensions))]))) - else - CFTData(T, T; shape, kwargs...) - end -end CFTData(scheme::TNRScheme; kwargs...) = CFTData(scheme.T; kwargs...) # simple 1-site unitcell schemes -CFTData(scheme::LoopTNR; kwargs...) = CFTData(scheme.TA, scheme.TB; kwargs...) # simple 1-site unitcell schemes +CFTData(scheme::LoopTNR; kwargs...) = CFTData(scheme.TA, scheme.TB; kwargs...) # 2-site unitcell schemes function CFTData(scheme::BTRG; kwargs...) # merge bond tensors into central tensor @tensor T_unit[-1 -2; -3 -4] := scheme.T[1 2; -3 -4] * scheme.S1[-2; 2] * scheme.S2[-1; 1] return CFTData(T_unit; kwargs...) end +# one-site unitcell +function CFTData(T::TensorMap{E, S, 2, 2}; shape = [sqrt(2), 2 * sqrt(2), 0], kwargs...) where {E, S} + if shape == [1, 1, 0] # trivial implementation + τ0 = elementary_modular_parameter(T, T) + Δs = _scaling_dimensions(T, τ0) + return CFTData(missing, τ0, StructuredVector(Δs, Dict([Trivial => collect(eachindex(Δs))]))) + else + CFTData(T, T; shape, kwargs...) + end +end + # Main implementation, two-site unitcell -function CFTData(TA::TensorMap{E, S, 2, 2}, TB::TensorMap{E, S, 2, 2}; shape = [sqrt(2), 2 * sqrt(2), 0], trunc = truncrank(16), truncentanglement = trunctol(; rtol = 1.0e-14)) where {E, S} +function CFTData( + TA::TensorMap{E, S, 2, 2}, TB::TensorMap{E, S, 2, 2}; shape = [sqrt(2), 2 * sqrt(2), 0], + trunc = truncrank(16), truncentanglement = trunctol(; rtol = 1.0e-14) + ) where {E, S} + norm_const = area_term(TA, TB)^(1 / 4) # canonical normalisation constant + TA′, TB′ = TA / norm_const, TB / norm_const + τ0 = elementary_modular_parameter(TA′, TB′) if shape == [1, 1, 0] throw(ArgumentError("The shape [1, 1, 0] is not compatible with a two-site unit cell.")) elseif (shape ≈ [sqrt(2), 2 * sqrt(2), 0]) || (shape == [1, 4, 1]) # these shapes need no truncation - norm_const = area_term(TA, TB)^(1 / 4) # canonical normalisation constant - return spec(TA / norm_const, TB / norm_const, shape) + return spec(TA′, TB′, shape, τ0) elseif (shape == [1, 8, 1]) || (shape ≈ [4 / sqrt(10), 2 * sqrt(10), 2 / sqrt(10)]) - - norm_const = area_term(TA, TB)^(1 / 4) # canonical normalisation constant - - dl, ur, ul, dr = MPO_opt( - TA / norm_const, TB / norm_const, trunc, truncentanglement - ) + dl, ur, ul, dr = MPO_opt(TA′, TB′, trunc, truncentanglement) T = reduced_MPO(dl, ur, ul, dr, trunc) - return spec(T, T, shape) + return spec(T, T, shape, τ0) else throw(ArgumentError("Shape $shape is not implemented.")) end end -# Trivial diagonalisation of the transfer matrix. Currently the v and unitcell are not acessible from the outside. -# The user should really be using the other shapes anyways. -function _scaling_dimensions(T::TensorMap{E, S, 2, 2}; v = 1, unitcell = 1) where {E, S} - # stack unitcell copies of T and trace +""" +Construct the transfer matrix along vertical direction +with `unitcell` copies of `T` concatenated horizontally. +`τ0` is the modular parameter of a single `T`. +""" +function _scaling_dimensions(T::TensorMap{E, S, 2, 2}, τ0::Number; unitcell = 1) where {E, S} indices = [[i, -i, -(i + unitcell), i + 1] for i in 1:unitcell] indices[end][4] = 1 @@ -84,25 +91,23 @@ function _scaling_dimensions(T::TensorMap{E, S, 2, 2}; v = 1, unitcell = 1) wher data = filter(x -> real(x) > 0, data) # filtering out negative real values data = filter(x -> abs(x) > 1.0e-12, data) # filtering out small values - return unitcell * (1 / (2π * v)) .* log.(data[1] ./ data) + # modular parameter of the constructed transfer matrix + Imτ = imag(τ0) / unitcell + return 1 / (2π * Imτ) .* log.(data[1] ./ data) end """ The "canonical" normalization constant for loop-TNR tensors, -which is the eigenvalue with largest real part of the 2 x 2 transfer matrix. +which is the eigenvalue with largest norm of the 2 x 2 transfer matrix. """ -function area_term(A, B; is_real = true) - a_in = domain(A)[1] - b_in = domain(B)[1] - x0 = ones(a_in ⊗ b_in) - +function area_term(TA, TB; is_real = true) function f0(x) - @plansor fx[-1 -2] := A[c -1; 1 m] * x[1 2] * B[m -2; 2 c] - @plansor ffx[-1 -2] := B[c -1; 1 m] * fx[1 2] * A[m -2; 2 c] + @plansor fx[-1 -2] := TA[c -1; 1 m] * x[1 2] * TB[m -2; 2 c] + @plansor ffx[-1 -2] := TB[c -1; 1 m] * fx[1 2] * TA[m -2; 2 c] return ffx end - - spec0, _, info = eigsolve(f0, x0, 1, :LR; verbosity = 0) + x0 = ones(domain(TA, 1) ⊗ domain(TB, 1)) + spec0, _, info = eigsolve(f0, x0, 1, :LM; verbosity = 0) if info.converged == 0 @warn "The area term eigensolver did not converge." end @@ -114,28 +119,40 @@ function area_term(A, B; is_real = true) end # The case with spin is based on https://arxiv.org/pdf/1512.03846 and some private communications with Yingjie Wei and Atsushi Ueda -function spec(TA::TensorMap, TB::TensorMap, shape::Array; Nh = 25) - area = shape[1] * shape[2] - Imτ = shape[1] / shape[2] - relative_shift = shape[3] / shape[1] - +""" +Internal function to construct transfer matrices and extract conformal data. + +# Arguments +- `TA, TB`: Rank-4 network tensors (may be identical for 1-site unit cells). +- `shape`: A triplet `[h, L, x]` — height, circumference, and horizontal shift + of the tube geometry, in units of the original tensor patch. +- `τ0`: Elementary modular parameter of one tensor. +- `Nh`: Number of eigenvalues to solve for per sector (default 25). +""" +function spec(TA::TensorMap, TB::TensorMap, shape::Vector{<:Number}, τ0::Number; Nh = 25) I = sectortype(TA) 𝔽 = field(TA) if BraidingStyle(I) != Bosonic() throw(ArgumentError("Sectors with non-Bosonic charge $I has not been implemented")) end - xspace, f = if shape ≈ [1, 4, 1] + # eigenvalues of the transfer matrix + xspace, f, τ = if shape ≈ [1, 4, 1] domain(TA)[1] ⊗ domain(TB)[1] ⊗ domain(TA)[1] ⊗ domain(TB)[1], - MPO_action_1x4_twist + MPO_action_1x4_twist, (1 + τ0) / 4 + elseif shape ≈ [sqrt(2), 2 * sqrt(2), 0] + domain(TB) ⊗ domain(TB), MPO_action_2gates, (1 + τ0) / 2 / (1 - τ0) + # in the following cases, (TA, TB) are no longer the original tensor in the network elseif shape ≈ [1, 8, 1] + τ0′ = (τ0 - 1) / 2 domain(TA)[1] ⊗ domain(TB)[1] ⊗ domain(TA)[1] ⊗ domain(TB)[1], - MPO_action_1x4 - elseif shape ≈ [sqrt(2), 2 * sqrt(2), 0] || - shape ≈ [4 / sqrt(10), 2 * sqrt(10), 2 / sqrt(10)] - domain(TB) ⊗ domain(TB), MPO_action_2gates + MPO_action_1x4, τ0′ / 4 + elseif shape ≈ [4 / sqrt(10), 2 * sqrt(10), 2 / sqrt(10)] + τ0′ = (τ0 - 1) / 2 + domain(TB) ⊗ domain(TB), MPO_action_2gates, (1 + τ0′) / 2 / (1 - τ0′) + else + error("Unsupported transfer matrix shape.") end - spec_sector = Dict( map(sectors(fuse(xspace))) do charge V = (I == Trivial) ? 𝔽^1 : Vect[I](charge => 1) @@ -156,19 +173,25 @@ function spec(TA::TensorMap, TB::TensorMap, shape::Array; Nh = 25) end ) + # central charge norm_const_0 = spec_sector[one(I)][1] - central_charge = 6 / pi / (Imτ - area / 4) * log(norm_const_0) + area = shape[1] * shape[2] + central_charge = 6 / pi / (imag(τ) - imag(τ0) * area / 4) * log(norm_const_0) # Construct a StructuredVector from the data of the different sectors data = ComplexF64[] structure = Dict{eltype(sectors(fuse(xspace))), Vector{Int}}() last_index = 1 + relative_shift = real(τ) / imag(τ) for charge in sectors(fuse(xspace)) - DeltaS = -1 / (2 * pi * Imτ) * log.(spec_sector[charge] / norm_const_0) - if !(relative_shift ≈ 0) + # DeltaS = Δ - i s Re(τ) / Im(τ) + DeltaS = -1 / (2 * pi * imag(τ)) * log.(spec_sector[charge] / norm_const_0) + if !isapprox(relative_shift, 0; atol = 1.0e-6) + # save `Δ - i s` in `data` push!(data, (real.(DeltaS) + imag.(DeltaS) / relative_shift * im)...) structure[charge] = [last_index:(last_index + length(DeltaS) - 1)...] else + # not enough precision to resolve conformal spin push!(data, real.(DeltaS)...) structure[charge] = [last_index:(last_index + length(DeltaS) - 1)...] end @@ -179,7 +202,7 @@ function spec(TA::TensorMap, TB::TensorMap, shape::Array; Nh = 25) sv = sort(sv; by = real) sv = filter(x -> real(x) ≤ 1.0e16, sv) - return CFTData(central_charge, sv) + return CFTData(central_charge, τ0, sv) end function MPO_opt( @@ -211,6 +234,10 @@ end # Apply functions for diagonalising different shapes of transfer matrices # ======================================================================= # Fig.25 of https://arxiv.org/pdf/2311.18785. Firstly appear in Chenfeng Bao's thesis, see http://hdl.handle.net/10012/14674. +""" +When the elementary modular parameter for TA, TB is `τ`, +the transfer matrix has `τ_TM = (1 + τ) / 2 / (1 - τ)`. +""" function MPO_action_2gates(TA::TensorMap, TB::TensorMap, x::TensorMap) @tensor fx[-1 -2 -3 -4; 5] := TB[-1 -2; 1 2] * x[1 2 3 4; 5] * TB[-3 -4; 3 4] @tensor ffx[-1 -2 -3 -4; 5] := TA[-3 -4; 2 3] * fx[1 2 3 4; 5] * @@ -218,6 +245,10 @@ function MPO_action_2gates(TA::TensorMap, TB::TensorMap, x::TensorMap) return permute(ffx, ((2, 3, 4, 1), (5,))) end +""" +When the elementary modular parameter for TA, TB is `τ`, +the transfer matrix has `τ_TM = τ / 4`. +""" function MPO_action_1x4(TA::TensorMap, TB::TensorMap, x::TensorMap) @tensor TTTTx[-1 -2 -3 -4; -5] := x[1 2 3 4; -5] * TA[41 -1; 1 12] * TB[12 -2; 2 23] * @@ -225,19 +256,162 @@ function MPO_action_1x4(TA::TensorMap, TB::TensorMap, x::TensorMap) return TTTTx end +""" +When the elementary modular parameter for TA, TB is `τ`, +the transfer matrix has `τ_TM = (τ + 1) / 4`. +""" function MPO_action_1x4_twist(TA::TensorMap, TB::TensorMap, x::TensorMap) TTTTx = MPO_action_1x4(TA, TB, x) return permute(TTTTx, ((2, 3, 4, 1), (5,))) end +""" +This renormalization will change the elementary +modular parameter from τ to (τ - 1) / 2. +""" function reduced_MPO( dl::TensorMap, ur::TensorMap, ul::TensorMap, dr::TensorMap, trunc::TruncationStrategy ) - @plansor temp[-1 -2; -3 -4] := ur[-1; 1 4] * - ul[4; 3 -2] * - dr[-3; 2 1] * dl[2; -4 3] + @plansor temp[-1 -2; -3 -4] := + ur[-1; 1 4] * ul[4; 3 -2] * dr[-3; 2 1] * dl[2; -4 3] D, U = SVD12(temp, trunc) @plansor translate[-1 -2; -3 -4] := U[-2; 1 -4] * D[-1 1; -3] return translate end + +# Elementary modular parameter +# ============================ +# TODO: one-tensor version +""" + elementary_modular_parameter(TA::TensorMap{E, S, 2, 2}, TB::TensorMap{E, S, 2, 2}) where {E, S} + +Extract the elementary modular parameter of one tensor from 2x2 transfer matrices. +""" +function elementary_modular_parameter( + TA::TensorMap{E, S, 2, 2}, TB::TensorMap{E, S, 2, 2} + ) where {E, S} + # TODO: fermions + if BraidingStyle(sectortype(TA)) != Bosonic() + return complex(0.0, 1.0) + end + # leading eigenvalue of each transfer matrix + # corresponding to the Δ = s = 0 identity field + λv = real(_eigsolve_2x2_NtoS(TA, TB)) + λh = real(_eigsolve_2x2_EtoW(TA, TB)) + λa = real(_eigsolve_2x2_NEtoSW(TA, TB)) + λb = real(_eigsolve_2x2_NWtoSE(TA, TB)) + # edge case: c = 0 + if all(isapprox.(λv, (λh, λa, λb); rtol = 1.0e-6)) + return complex(0.0, 1.0) + end + # when c ≠ 0 + a1, a2, a3 = log(λh / λv), log(λa / λv), log(λb / λv) + # c here is actually π c / 6 + c, v, θ = solve_cvtheta(a1, a2, a3) + return v * cis(θ) +end + +function _eigsolve_2x2_NtoS(TA, TB) + function f0(x) + @plansor fx[-1 -2] := TA[c -1; 1 m] * x[1 2] * TB[m -2; 2 c] + @plansor fx[-1 -2] := TB[c -1; 1 m] * fx[1 2] * TA[m -2; 2 c] + return fx + end + x0 = ones(domain(TA, 1) ⊗ domain(TB, 1)) + spec0, _, info = eigsolve(f0, x0, 1, :LM; verbosity = 0) + if info.converged == 0 + @warn "The area term eigensolver did not converge." + end + # TODO: for fermions + return first(spec0) +end + +function _eigsolve_2x2_EtoW(TA, TB) + TA′ = permute(TA, ((3, 1), (4, 2))) + TB′ = permute(TB, ((3, 1), (4, 2))) + return _eigsolve_2x2_NtoS(TB′, TA′) +end + +function _eigsolve_2x2_NEtoSW(TA, TB) + #= + 1 2 + ┌---┬---┐ + | | | + 3'--A---B---┤ 3 + | | | + 4'--B---A---┘ 4 + | | + 1' 2' + =# + function f0(x) + @plansor fx[-1 -2 -3 -4] := TB[-2 -3; a b] * x[-1 a b -4] + @plansor fx[-1 -2 -3 -4] := TA[-1 -2; a b] * fx[a b -3 -4] + @plansor fx[-1 -2 -3 -4] := TA[-3 -4; a b] * fx[-1 -2 a b] + @plansor fx[-1 -2 -3 -4] := TB[-4 -1; a b] * fx[-3 a b -2] + return fx + end + x0 = ones(domain(TA, 1) ⊗ domain(TB, 1) ⊗ domain(TB, 2) ⊗ domain(TA, 2)) + spec0, _, info = eigsolve(f0, x0, 1, :LM; verbosity = 0) + if info.converged == 0 + @warn "The area term eigensolver did not converge." + end + # TODO: for fermions + return first(spec0) +end + +function _eigsolve_2x2_NWtoSE(TA, TB) + TA′ = permute(TA, ((2, 4), (1, 3))) + TB′ = permute(TB, ((2, 4), (1, 3))) + return _eigsolve_2x2_NEtoSW(TB′, TA′) +end + +# Logistic sigmoid and its inverse (logit) +sigmoid(x) = 1 / (1 + exp(-x)) +logit(p) = log(p / (1 - p)) + +""" + solve_cvtheta(a1, a2, a3; c0=0.5, v0=1.0, θ0=π/2) + +Solve for positive (c, v) and θ ∈ (0, π) from the three equations: + + a1 = (1/v - v) * c * sin(θ) + a2 = (v / (1 + v² - 2v cos θ) - v) * c * sin(θ) + a3 = (v / (1 + v² + 2v cos θ) - v) * c * sin(θ) + +Returns `(c, v, θ)`. +""" +function solve_cvtheta(a1, a2, a3; c0 = 0.5, v0 = 1.0, θ0 = π / 2) + # Work in unconstrained coords to keep variables in their natural domain. + # c = exp(xc) → c > 0 + # v = exp(xv) → v > 0 + # θ = π * σ(xθ) → θ ∈ (0, π) + function f!(du, u, p) + xc, xv, xθ = u + c = exp(xc) + v = exp(xv) + θ = π * sigmoid(xθ) + + sinθ = sin(θ) + cosθ = cos(θ) + + du[1] = (1 / v - v) * c * sinθ - a1 + du[2] = (v / (1 + v^2 - 2v * cosθ) - v) * c * sinθ - a2 + return du[3] = (v / (1 + v^2 + 2v * cosθ) - v) * c * sinθ - a3 + end + + # Initial guess in unconstrained space + u0 = [log(c0), log(v0), logit(θ0 / π)] + prob = NonlinearProblem(f!, u0) + sol = solve(prob, NewtonRaphson(; autodiff = AutoForwardDiff())) + + if !SciMLBase.successful_retcode(sol) + @warn "Solver did not converge" retcode = sol.retcode residuals = sol.residuals + end + + xc, xv, xθ = sol.u + c = exp(xc) + v = exp(xv) + θ = π * sigmoid(xθ) + return c, v, θ +end diff --git a/test/schemes/schemes.jl b/test/schemes/schemes.jl index f3c172db..1893a544 100644 --- a/test/schemes/schemes.jl +++ b/test/schemes/schemes.jl @@ -11,31 +11,39 @@ T = classical_ising() T_3D = classical_ising_3D() const f_benchmark3D = ising_3D_free_energy_htse() +const Jx_aniso, Jy_aniso = 1.0, 0.6 +const βc_aniso_test = ising_anisotropic_βc(Jx_aniso, Jy_aniso) +const T_aniso = classical_ising(Z2Irrep, βc_aniso_test; Jx = Jx_aniso, Jy = Jy_aniso) +const f_aniso_exact = f_onsager_anisotropic(βc_aniso_test, Jx_aniso, Jy_aniso) + function cft_finalize!(scheme) finalize!(scheme) return CFTData(scheme) end # TRG -@testset "TRG - Ising Model" begin - @info "TRG ising free energy" - scheme = TRG(T) +@testset "TRG - Anisotropic Ising Model (Jx=1.0, Jy=0.6)" begin + @info "TRG anisotropic ising free energy" + scheme = TRG(T_aniso) data = run!(scheme, truncrank(24), maxiter(25)) - @test free_energy(data, ising_βc) ≈ f_onsager rtol = 2.0e-6 + @test free_energy(data, βc_aniso_test) ≈ f_aniso_exact rtol = 2.0e-6 - @info "TRG ising CFT data" - scheme = TRG(T) + @info "TRG anisotropic ising CFT data — shape [1, 1, 0]" + scheme = TRG(T_aniso) run!(scheme, truncrank(24), maxiter(10)) - cft = sort(CFTData(scheme; shape = [1, 1, 0]).scaling_dimensions[2:end]; by = abs) .|> real + cft = CFTData(scheme; shape = [1, 1, 0]) + sd_all = real(cft.scaling_dimensions[Trivial]) + cft_sorted = sort(sd_all[2:end]; by = abs) - @test cft[1] ≈ ising_cft_exact[1] rtol = 2.0e-4 - @test cft[2] ≈ ising_cft_exact[2] rtol = 1.0e-2 + @test cft_sorted[1] ≈ ising_cft_exact[1] rtol = 2.0e-4 + @test cft_sorted[2] ≈ ising_cft_exact[2] rtol = 2.0e-2 - @info "TRG ising ground state degeneracy" + @info "Obtained scaling dimensions: Δ₁ = $(cft_sorted[1]), Δ₂ = $(cft_sorted[2])" - T1 = classical_ising(ising_βc - 0.01) + @info "TRG anisotropic ising ground state degeneracy" + T1 = classical_ising(βc_aniso_test - 0.01; Jx = Jx_aniso, Jy = Jy_aniso) scheme = TRG(T1) run!(scheme, truncrank(16), maxiter(20)) gsd = ground_state_degeneracy(scheme) @@ -44,7 +52,7 @@ end @test X1 ≈ 1.0 rtol = 1.0e-2 @test X2 ≈ 1.0 rtol = 1.0e-2 - T2 = classical_ising(ising_βc + 0.01) + T2 = classical_ising(βc_aniso_test + 0.01; Jx = Jx_aniso, Jy = Jy_aniso) scheme = TRG(T2) run!(scheme, truncrank(16), maxiter(20)) gsd = ground_state_degeneracy(scheme) @@ -52,7 +60,6 @@ end @test gsd ≈ 2 rtol = 1.0e-2 @test X1 ≈ 2.0 rtol = 1.0e-2 @test X2 ≈ 2.0 rtol = 1.0e-2 - end # BTRG @@ -167,9 +174,10 @@ end end # LoopTNR -@testset "LoopTNR - Ising Model - Dense Solver" begin - @info "LoopTNR ising free energy" - scheme = LoopTNR(T) + +@testset "LoopTNR - Anisotropic Ising Model (Jx=1.0, Jy=0.6)" begin + @info "LoopTNR anisotropic ising free energy" + scheme = LoopTNR(T_aniso) loop_condition = LoopParameters( sweeping = maxiter(5) & convcrit(1.0e-9, (steps, cost) -> abs(cost[end])), @@ -180,30 +188,37 @@ end scheme, truncrank(8), maxiter(25), loop_condition ) - @test free_energy(data, ising_βc) ≈ f_onsager rtol = 1.0e-6 + @test free_energy(data, βc_aniso_test) ≈ f_aniso_exact rtol = 1.0e-6 - @info "LoopTNR ising CFT data" - scheme = LoopTNR(T) + @info "LoopTNR anisotropic ising CFT data" + scheme = LoopTNR(T_aniso) run!(scheme, truncrank(12), maxiter(10)) for shape in [[1, 4, 1], [sqrt(2), 2 * sqrt(2), 0]] - cft = sort(CFTData(scheme; shape = shape).scaling_dimensions; by = abs) - d1, d2 = real(cft[Z2Irrep(1)][1]), real(cft[Z2Irrep(0)][2]) - @info "Obtained lowest scaling dimensions:\n$(d1), $(d2)." - @test d1 ≈ ising_cft_exact[1] rtol = 5.0e-4 - @test d2 ≈ ising_cft_exact[2] rtol = 5.0e-4 + cft = CFTData(scheme; shape = shape) + d_σ = real(cft.scaling_dimensions[Z2Irrep(1)][1]) + d_ε = real(cft.scaling_dimensions[Z2Irrep(0)][2]) + @info "Shape $shape: Δ(σ) = $d_σ, Δ(ε) = $d_ε, c = $(cft.central_charge)" + @test d_σ ≈ ising_cft_exact[1] rtol = 5.0e-4 + @test d_ε ≈ ising_cft_exact[2] rtol = 5.0e-4 + @test cft.central_charge ≈ 0.5 rtol = 5.0e-3 end for shape in [[1, 8, 1], [4 / sqrt(10), 2 * sqrt(10), 2 / sqrt(10)]] - cft = sort(CFTData(scheme; shape = shape, trunc = truncrank(16), truncentanglement = trunctol(atol = 1.0e-10)).scaling_dimensions; by = real) - d1, d2 = real(cft[Z2Irrep(1)][1]), real(cft[Z2Irrep(0)][2]) - @info "Obtained lowest scaling dimensions:\n$(d1), $(d2)." - @test d1 ≈ ising_cft_exact[1] rtol = 1.0e-3 - @test d2 ≈ ising_cft_exact[2] rtol = 1.0e-3 + cft = CFTData( + scheme; shape = shape, trunc = truncrank(16), + truncentanglement = trunctol(atol = 1.0e-10) + ) + d_σ = real(cft.scaling_dimensions[Z2Irrep(1)][1]) + d_ε = real(cft.scaling_dimensions[Z2Irrep(0)][2]) + @info "Shape $shape: Δ(σ) = $d_σ, Δ(ε) = $d_ε, c = $(cft.central_charge)" + @test d_σ ≈ ising_cft_exact[1] rtol = 2.0e-3 + @test d_ε ≈ ising_cft_exact[2] rtol = 2.0e-3 + @test cft.central_charge ≈ 0.5 rtol = 1.0e-2 end - @info "LoopTNR ising ground state degeneracy" - T1 = classical_ising(ising_βc - 0.01) + @info "LoopTNR anisotropic ising ground state degeneracy" + T1 = classical_ising(βc_aniso_test - 0.01; Jx = Jx_aniso, Jy = Jy_aniso) scheme = LoopTNR(T1) run!(scheme, truncrank(12), maxiter(20)) gsd = ground_state_degeneracy(scheme) @@ -212,7 +227,7 @@ end @test X1 ≈ 1.0 rtol = 1.0e-2 @test X2 ≈ 1.0 rtol = 1.0e-2 - T2 = classical_ising(ising_βc + 0.01) + T2 = classical_ising(βc_aniso_test + 0.01; Jx = Jx_aniso, Jy = Jy_aniso) scheme = LoopTNR(T2) run!(scheme, truncrank(12), maxiter(20)) gsd = ground_state_degeneracy(scheme)