Add test for classical Ising model

Project.toml

@@ -16,14 +16,16 @@ TupleTools = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6"
 MPSKit = "0.12"
 MacroTools = "0.5"
 PrecompileTools = "1"
+QuadGK = "2.11.1"
 TensorKit = "0.13,0.14"
 TensorOperations = "5"
 TupleTools = "1"
 julia = "1.10"
 
 [extras]
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestExtras = "5ed8adda-3752-4e41-b88a-e8b09835ee3a"
 
 [targets]
-test = ["Test", "TestExtras"]
+test = ["Test", "TestExtras", "QuadGK"]

src/models/transfermatrices.jl

@@ -8,7 +8,7 @@
 MPO for the partition function of the two-dimensional classical Ising model, defined as
 
 ```math
-\\mathcal{Z}(\\beta) = \\sum_{\\{s\\}} \\exp(-\\beta H(s)) \\text{ with } H(s) = \\sum_{\\langle i, j \\rangle} s_i s_j
+\\mathcal{Z}(\\beta) = \\sum_{\\{s\\}} \\exp(-\\beta H(s)) \\text{ with } H(s) = -\\sum_{\\langle i, j \\rangle} s_i s_j
 
 ```
 where each classical spin can take the values ``s = \\pm 1``.

test/classical_ising.jl

@@ -0,0 +1,41 @@
+using Test
+using TensorKit
+using QuadGK
+using MPSKit
+using MPSKitModels
+
+## Setup
+
+beta = 0.6
+Vspaces = [ComplexSpace(12), Z2Space(0 => 6, 1 => 6)]
+alg = VUMPS(; tol=1e-8, maxiter=25, verbosity=1)
+
+"""
+    classical_ising_free_energy(; beta)
+
+[Exact Onsager solution](https://en.wikipedia.org/wiki/Square_lattice_Ising_model#Exact_solution)
+for the free energy of the 2D classical Ising Model with partition function
+
+```math
+\\mathcal{Z}(\\beta) = \\sum_{\\{s\\}} \\exp(-\\beta H(s)) \\text{ with } H(s) = - \\sum_{\\langle i, j \\rangle} s_i s_j
+```
+"""
+function classical_ising_free_energy(; beta=log(1 + sqrt(2)) / 2)
+    k = 1 / sinh(2 * beta)^2
+    F = quadgk(theta -> log(cosh(2 * beta)^2 +
+                            1 / k * sqrt(1 + k^2 - 2 * k * cos(2 * theta))),
+               0, pi)[1]
+    return -1 / beta * (log(2) / 2 + 1 / (2 * pi) * F)
+end
+
+## Test
+
+@testset "Classical Ising for $(sectortype(V)) symmetry" for V in Vspaces
+    O = classical_ising(sectortype(V); beta)
+    psi0 = InfiniteMPS(physicalspace(O, 1), V)
+    psi, envs, = leading_boundary(psi0, O, alg)
+    λ = expectation_value(psi, O, envs)
+    f = -log(λ) / beta
+    f_exact = classical_ising_free_energy(; beta)
+    @test abs(f - f_exact) < 1e-4
+end

test/quantum_potts.jl

@@ -0,0 +1,40 @@
+using Test
+using TensorKit
+using MPSKit
+using MPSKitModels
+
+## Setup
+
+χ = 6
+qs = [3, 4, 5]
+symmetries = [Trivial, ZNIrrep]
+Vspaces = [ComplexSpace(12), Z2Space(0 => 6, 1 => 6)]
+alg = VUMPS(; maxiter=25, verbosity=0)
+
+# https://iopscience.iop.org/article/10.1088/0305-4470/14/11/020/meta
+function quantum_potts_energy(Q::Int; maxiter::Int=1000)
+    Q == 3 && return -(4 / 3 + 2sqrt(3) / π)
+    Q == 4 && return 2 - 8 * log(2)
+    summation = sum((-1)^n / (sqrt(Q) / 2 - cosh((2 * n + 1) * acosh(sqrt(Q) / 2)))
+                    for n in 1:maxiter)
+    limit = 2 - Q - sqrt(Q) * (Q - 4) * summation
+    return limit
+end
+
+_sectortype(::Type{Trivial}, q::Int) = Trivial
+_sectortype(::Type{ZNIrrep}, q::Int) = ZNIrrep{q}
+_virtualspace(::Type{Trivial}, q::Int, χ::Int) = ComplexSpace(q * χ)
+_virtualspace(::Type{ZNIrrep}, q::Int, χ::Int) = Rep[ℤ{q}](i => χ for i in 0:(q - 1))
+
+## Test
+
+@testset "$q-state Potts with $(_sectortype(symmetry, q))) symmetry" for q in qs,
+                                                                         symmetry in
+                                                                         symmetries
+
+    H = quantum_potts(_sectortype(symmetry, q); q)
+    ψ = InfiniteMPS(physicalspace(H, 1), _virtualspace(symmetry, q, χ))
+    @test imag(expectation_value(ψ, H)) ≈ 0 atol = 1e-10
+    ψ, envs, δ = find_groundstate(ψ, H, alg)
+    @test quantum_potts_energy(q) ≈ expectation_value(ψ, H, envs) atol = 1e-2
+end

test/runtests.jl

@@ -37,6 +37,14 @@ end
     include("heisenberg.jl")
 end
 
+@testset "quantum potts model" begin
+    include("quantum_potts.jl")
+end
+
+@testset "classical ising model" begin
+    include("classical_ising.jl")
+end
+
 @testset "sixvertex model" begin
     include("sixvertex.jl")
 end