More Ising tests

Project.toml

@@ -1,7 +1,7 @@
 name = "ITensorNetworksNext"
 uuid = "302f2e75-49f0-4526-aef7-d8ba550cb06c"
 authors = ["ITensor developers <support@itensor.org> and contributors"]
-version = "0.1.13"
+version = "0.1.14"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"

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])

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