diff --git a/Project.toml b/Project.toml index 5bb5d04..0392ff1 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ITensorNetworksNext" uuid = "302f2e75-49f0-4526-aef7-d8ba550cb06c" authors = ["ITensor developers and contributors"] -version = "0.1.13" +version = "0.1.14" [deps] AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c" diff --git a/src/TensorNetworkGenerators/ising_network.jl b/src/TensorNetworkGenerators/ising_network.jl index edc7ac6..afd8d74 100644 --- a/src/TensorNetworkGenerators/ising_network.jl +++ b/src/TensorNetworkGenerators/ising_network.jl @@ -2,16 +2,17 @@ using DiagonalArrays: DiagonalArray using Graphs: degree, dst, edges, src using LinearAlgebra: Diagonal, eigen using NamedDimsArrays: apply, dename, inds, operator, randname +using NamedGraphs.GraphsExtensions: vertextype -function sqrt_ising_bond(β; h1 = zero(typeof(β)), h2 = zero(typeof(β))) - f11 = exp(β * (1 + h1 + h2)) - f12 = exp(β * (-1 + h1 - h2)) - f21 = exp(β * (-1 - h1 + h2)) - f22 = exp(β * (1 - h1 - h2)) - m² = eltype(β)[f11 f12; f21 f22] - d², v = eigen(m²) - d = sqrt.(d²) - return v * Diagonal(d) * inv(v) +function sqrt_ising_bond(β; J = one(β), h = zero(β), deg1::Integer, deg2::Integer) + h1 = h / deg1 + h2 = h / deg2 + m = [ + exp(β * (J + h1 + h2)) exp(β * (-J + h1 - h2)); + exp(β * (-J - h1 + h2)) exp(β * (J - h1 - h2)); + ] + d, v = eigen(m) + return v * √(Diagonal(d)) * inv(v) end """ @@ -23,7 +24,8 @@ partition function tensors on each vertex. Link dimensions are defined using the edge. """ function ising_network( - f, β::Number, g::AbstractGraph; h::Number = zero(eltype(β)), sz_vertices = [] + f, β::Number, g::AbstractGraph; J::Number = one(β), h::Number = zero(β), + sz_vertices = vertextype(g)[] ) elt = typeof(β) l̃ = Dict(e => randname(f(e)) for e in edges(g)) @@ -38,7 +40,7 @@ function ising_network( v2 = dst(e) deg1 = degree(tn, v1) deg2 = degree(tn, v2) - m = sqrt_ising_bond(β; h1 = h / deg1, h2 = h / deg2) + m = sqrt_ising_bond(β; J, h, deg1, deg2) t = operator(m, (f̃(e),), (f(e),)) tn[v1] = apply(t, tn[v1]) tn[v2] = apply(t, tn[v2]) diff --git a/test/test_tensornetworkgenerators.jl b/test/test_tensornetworkgenerators.jl index 090a964..152e67b 100644 --- a/test/test_tensornetworkgenerators.jl +++ b/test/test_tensornetworkgenerators.jl @@ -11,13 +11,41 @@ using Test: @test, @testset module TestUtils using QuadGK: quadgk # Exact critical inverse temperature for 2D square lattice Ising model. - βc() = 0.5 * log(1 + √2) + βc_2d_ising(elt::Type{<:Number} = Float64) = elt(log(1 + √2) / 2) # Exact infinite volume free energy density for 2D square lattice Ising model. - function ising_free_energy_density(β::Real) - κ = 2sinh(2β) / cosh(2β)^2 - integrand(θ) = log(0.5 * (1 + sqrt(abs(1 - (κ * sin(θ))^2)))) + function f_2d_ising(β::Real; J::Real = one(β)) + κ = 2sinh(2β * J) / cosh(2β * J)^2 + integrand(θ) = log((1 + √(abs(1 - (κ * sin(θ))^2))) / 2) integral, _ = quadgk(integrand, 0, π) - return (-log(2cosh(2β)) - (1 / (2π)) * integral) / β + return (-log(2cosh(2β * J)) - (1 / (2π)) * integral) / β + end + function f_1d_ising(β::Real; J::Real = one(β), h::Real = zero(β)) + λ⁺ = exp(β * J) * (cosh(β * h) + √(sinh(β * h)^2 + exp(-4β * J))) + return -(log(λ⁺) / β) + end + function f_1d_ising(β::Real, N::Integer; periodic::Bool = true, kwargs...) + return if periodic + f_1d_ising_periodic(β, N; kwargs...) + else + f_1d_ising_open(β, N; kwargs...) + end + end + function f_1d_ising_periodic(β::Real, N::Integer; J::Real = one(β), h::Real = zero(β)) + r = √(sinh(β * h)^2 + exp(-4β * J)) + λ⁺ = exp(β * J) * (cosh(β * h) + r) + λ⁻ = exp(β * J) * (cosh(β * h) - r) + Z = λ⁺^N + λ⁻^N + return -(log(Z) / (β * N)) + end + function f_1d_ising_open(β::Real, N::Integer; J::Real = one(β), h::Real = zero(β)) + isone(N) && return 2cosh(β * h) + T = [ + exp(β * (J + h)) exp(-β * J); + exp(-β * J) exp(β * (J - h)); + ] + b = [exp(β * h / 2), exp(-β * h / 2)] + Z = (b' * (T^(N - 1)) * b)[] + return -(log(Z) / (β * N)) end end @@ -38,25 +66,49 @@ end end end @testset "Ising Network" begin - dims = (4, 4) - β = TestUtils.βc() - g = named_grid(dims; periodic = true) - ldict = Dict(e => Index(2) for e in edges(g)) - l(e) = get(() -> ldict[reverse(e)], ldict, e) - tn = ising_network(l, β, g) - @test nv(tn) == 16 - @test ne(tn) == ne(g) - @test issetequal(vertices(tn), vertices(g)) - @test issetequal(arranged_edges(tn), arranged_edges(g)) - for v in vertices(tn) - is = l.(incident_edges(g, v)) - @test issetequal(is, inds(tn[v])) - @test tn[v] ≠ δ(Tuple(is)) + @testset "1D Ising (periodic = $periodic)" for periodic in (false, true) + dims = (4,) + β = 0.4 + g = named_grid(dims; periodic) + ldict = Dict(e => Index(2) for e in edges(g)) + l(e) = get(() -> ldict[reverse(e)], ldict, e) + tn = ising_network(l, β, g) + @test nv(tn) == 4 + @test ne(tn) == ne(g) + @test issetequal(vertices(tn), vertices(g)) + @test issetequal(arranged_edges(tn), arranged_edges(g)) + for v in vertices(tn) + is = l.(incident_edges(g, v)) + @test issetequal(is, inds(tn[v])) + @test tn[v] ≠ δ(Tuple(is)) + end + # TODO: Use eager contraction sequence finding. + z = contract_network(tn; alg = "exact")[] + f = -log(z) / (β * nv(g)) + f_analytic = TestUtils.f_1d_ising(β, 4; periodic) + @test f ≈ f_analytic + end + @testset "2D Ising" begin + dims = (4, 4) + β = TestUtils.βc_2d_ising() + g = named_grid(dims; periodic = true) + ldict = Dict(e => Index(2) for e in edges(g)) + l(e) = get(() -> ldict[reverse(e)], ldict, e) + tn = ising_network(l, β, g) + @test nv(tn) == 16 + @test ne(tn) == ne(g) + @test issetequal(vertices(tn), vertices(g)) + @test issetequal(arranged_edges(tn), arranged_edges(g)) + for v in vertices(tn) + is = l.(incident_edges(g, v)) + @test issetequal(is, inds(tn[v])) + @test tn[v] ≠ δ(Tuple(is)) + end + # TODO: Use eager contraction sequence finding. + z = contract_network(tn; alg = "exact")[] + f = -log(z) / (β * nv(g)) + f_inf = TestUtils.f_2d_ising(β) + @test f ≈ f_inf rtol = 1.0e-1 end - # TODO: Use eager contraction sequence finding. - z = contract_network(tn; alg = "exact")[] - f = -log(z) / (β * nv(g)) - f_inf = TestUtils.ising_free_energy_density(β) - @test f ≈ f_inf rtol = 1.0e-1 end end