diff --git a/Project.toml b/Project.toml
index 4bd03cf..0ac17ae 100644
--- a/Project.toml
+++ b/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"]
diff --git a/src/models/transfermatrices.jl b/src/models/transfermatrices.jl
index c467e81..c9ff600 100644
--- a/src/models/transfermatrices.jl
+++ b/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``.
diff --git a/test/classical_ising.jl b/test/classical_ising.jl
new file mode 100644
index 0000000..9daaaa6
--- /dev/null
+++ b/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
diff --git a/test/potts.jl b/test/potts.jl
deleted file mode 100644
index 1526563..0000000
--- a/test/potts.jl
+++ /dev/null
@@ -1,66 +0,0 @@
-using MPSKit
-using TensorKit
-
-alg = VUMPS(; maxiter=25, verbosity=0)
-
-# https://iopscience.iop.org/article/10.1088/0305-4470/14/11/020/meta
-function E₀(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
-
-# Q = 3
-@testset "Q=3 Trivial" begin
-    H = quantum_potts(; q=3)
-    ψ = InfiniteMPS(3, 32)
-    @test imag(expectation_value(ψ, H)) ≈ 0 atol = 1e-10
-    ψ, envs, δ = find_groundstate(ψ, H, alg)
-    @test E₀(3) ≈ expectation_value(ψ, H, envs) atol = 1e-2
-end
-
-@testset "Z3Irrep" begin
-    H = quantum_potts(Z3Irrep; q=3)
-    ψ = InfiniteMPS(Rep[ℤ₃](i => 1 for i in 0:2), Rep[ℤ₃](i => 12 for i in 0:2))
-    @test imag(expectation_value(ψ, H)) ≈ 0 atol = 1e-10
-    ψ, envs, δ = find_groundstate(ψ, H, alg)
-    @test E₀(3) ≈ expectation_value(ψ, H, envs) atol = 1e-2
-end
-
-# Q = 4
-@testset "Q=4 Trivial" begin
-    H = quantum_potts(; q=4)
-    ψ = InfiniteMPS(4, 45)
-    @test imag(expectation_value(ψ, H)) ≈ 0 atol = 1e-10
-    ψ, envs, δ = find_groundstate(ψ, H, alg)
-    @test E₀(4) ≈ expectation_value(ψ, H, envs) atol = 1e-2
-end
-
-@testset "Z4Irrep" begin
-    H = quantum_potts(Z4Irrep; q=4)
-    ψ = InfiniteMPS(Vect[Z4Irrep](i => 1 for i in 0:3), Vect[Z4Irrep](i => 12 for i in 0:3))
-    @test imag(expectation_value(ψ, H)) ≈ 0 atol = 1e-10
-    ψ, envs, δ = find_groundstate(ψ, H, alg)
-    @test E₀(4) ≈ expectation_value(ψ, H, envs) atol = 1e-2
-end
-
-# Q = 5
-@testset "Q=5 Trivial" begin
-    H = quantum_potts(; q=5)
-    ψ = InfiniteMPS(5, 60)
-    @test imag(expectation_value(ψ, H)) ≈ 0 atol = 1e-10
-    ψ, envs, δ = find_groundstate(ψ, H, alg)
-    @test E₀(5) ≈ expectation_value(ψ, H, envs) atol = 1e-2
-end
-
-@testset "ZNIrrep{5}" begin
-    H = quantum_potts(ZNIrrep{5}; q=5)
-    ψ = InfiniteMPS(Vect[ZNIrrep{5}](i => 1 for i in 0:4),
-                    Vect[ZNIrrep{5}](i => 12 for i in 0:4))
-    @test imag(expectation_value(ψ, H)) ≈ 0 atol = 1e-10
-    ψ, envs, δ = find_groundstate(ψ, H, alg)
-    @test E₀(5) ≈ expectation_value(ψ, H, envs) atol = 1e-2
-end
diff --git a/test/quantum_potts.jl b/test/quantum_potts.jl
new file mode 100644
index 0000000..07bef37
--- /dev/null
+++ b/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
diff --git a/test/runtests.jl b/test/runtests.jl
index 08e9b52..853d174 100644
--- a/test/runtests.jl
+++ b/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