Fix 2D Heisenberg Example (#95)

examples/vumps/vumps_2d_heisenberg.jl

@@ -31,25 +31,6 @@ initstate(n) = isodd(n) ? "↑" : "↓"
 s = infsiteinds("S=1/2", N; conserve_qns, initstate)
 ψ = InfMPS(s, initstate)
 
-function ITensorInfiniteMPS.unit_cell_terms(::Model"heisenberg2D"; width, yperiodic)
-  opsum = OpSum()
-  for i in 1:width
-    # Vertical
-    opsum -= 0.5, "S+", i, "S-", i + 1
-    opsum -= 0.5, "S-", i, "S+", i + 1
-    opsum += "Sz", i, "Sz", i + 1
-    # Horizontal
-    opsum -= 0.5, "S+", i, "S-", i + width
-    opsum -= 0.5, "S-", i, "S+", i + width
-    opsum += "Sz", i, "Sz", i + width
-  end
-  if yperiodic
-    opsum -= 0.5, "S+", 1, "S-", width
-    opsum -= 0.5, "S-", 1, "S+", width
-    opsum += "Sz", 1, "Sz", width
-  end
-  return opsum
-end
 model = Model("heisenberg2D")
 
 # Form the Hamiltonian
@@ -75,6 +56,11 @@ subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim)
 
 energy_infinite = expect(ψ, H)
 @show energy_infinite
+@show sum(energy_infinite) / width
 
+energy_approx_exact = reference(model, Observable("energy"); width, yperiodic)
+@show energy_approx_exact
+
+@show isapprox(sum(energy_infinite) / width, energy_approx_exact, atol=1e-4)
 ## using JLD2
 ## jldsave("infmps.jld2"; ψ)

src/models/heisenberg.jl

@@ -25,3 +25,49 @@ function reference(::Model"heisenberg", ::Observable"energy"; N=∞)
   correction = 1 + 0.375 / log(N)^3
   return (E∞ - Eᶠⁱⁿⁱᵗᵉ * correction) / (2N)
 end
+
+function unit_cell_terms(::Model"heisenberg2D"; width, yperiodic)
+  opsum = OpSum()
+  for i in 1:width
+    # Vertical
+    if (i < width) || (yperiodic && width > 2)
+      opsum -= 0.5, "S+", i, "S-", mod(i, width) + 1
+      opsum -= 0.5, "S-", i, "S+", mod(i, width) + 1
+      opsum += "Sz", i, "Sz", mod(i, width) + 1
+    end
+    # Horizontal
+    opsum -= 0.5, "S+", i, "S-", i + width
+    opsum -= 0.5, "S-", i, "S+", i + width
+    opsum += "Sz", i, "Sz", i + width
+  end
+  return opsum
+end
+
+"""
+    reference(::Model"heisenberg2D", ::Observable"energy"; width, yperiodic)
+
+Report the reference isotropic 2D square heisenberg ground energy per site for length `N`.
+Taken from [1,2]. Note that periodic results have an errorbar of ~1e-4 - 1e-3
+
+[1] Ramos and Xavier. "N-leg spin-S Heisenberg ladders:
+A density-matrix renormalization group study"
+Phys. Rev. B 89, 094424 - Published 27 March 2014
+[2] Frischmuth, Ammon, and Troyer. "Susceptibility and low-temperature
+thermodynamics of spin-½ Heisenberg ladders"
+Phys. Rev. B 54, R3714(R) - Published 1 August 1996
+"""
+function reference(::Model"heisenberg2D", ::Observable"energy"; width, yperiodic)
+  if width > 6
+    error("Ladders of width greater than 6 are not in reference data")
+  end
+
+  if yperiodic
+    (width == ∞) && return -0.66931
+    energies = [-0.4432, -0.5780, -0.6006, -0.6187, -0.6278, -0.635]
+    return energies[width]
+  else
+    (width == ∞) && return -0.6768
+    energies = [-0.4431471, -0.578043140180, -0.600537, -0.618566, -0.62776, -0.6346]
+    return energies[width]
+  end
+end

test/test_modelMPOs.jl

@@ -0,0 +1,89 @@
+using ITensors, ITensorMPS
+using ITensorInfiniteMPS
+using Test
+
+@testset verbose = true "Heisenberg Model Test" begin
+  cell_widths = [2, 3, 4, 5]
+  @testset "cell_width=$width" for width in cell_widths
+    model = Model("heisenberg")
+    os = ITensorInfiniteMPS.opsum_finite(model, width;)
+
+    connections = []
+    for term in os
+      push!(connections, sort(ITensors.sites(term)))
+    end
+    for i in 1:(width - 1)
+      @test [i, i + 1] ∈ connections
+    end
+  end
+end
+
+@testset verbose = true "Heisenberg2D Model Test" begin
+  widths = [2, 3, 4, 5, 6]
+  cell_widths = [2, 3, 4, 5]
+  @testset "cell_width=$cell_width width=$width" for cell_width in cell_widths,
+    width in widths
+
+    model = Model("heisenberg2D")
+
+    os = ITensorInfiniteMPS.opsum_finite(model, cell_width * width; width, yperiodic)
+
+    connections = []
+    for term in os
+      push!(connections, ITensors.sites(term))
+    end
+    for col in 1:(cell_width - 1)
+      for row in 1:(width - 1)
+        i = (col - 1) * width + row
+        @test [i, i + 1] ∈ connections
+        @test [i, i + width] ∈ connections
+      end
+      (yperiodic && width > 2) && @test [i, i + width - 1] ∈ connections
+    end
+    # the above forgets the last horizontal bond
+    for col in 1:(cell_width - 1)
+      i = (col - 1) * width + width
+      @test [i, i + width] ∈ connections
+    end
+  end
+end
+
+@testset verbose = true "Ising Model Test" begin
+  cell_widths = [2, 3, 4, 5]
+  @testset "cell_width=$width" for width in cell_widths
+    model = Model("ising")
+    os = ITensorInfiniteMPS.opsum_finite(model, width;)
+
+    connections = []
+    for term in os
+      push!(connections, sort(ITensors.sites(term)))
+    end
+    for i in 1:(width - 1)
+      @test [i, i + 1] ∈ connections
+      @test [i] ∈ connections
+    end
+    @test [width] ∈ connections
+  end
+end
+
+@testset verbose = true "Hubbard Model Test" begin
+  t, U = 1.0, 4.0
+  cell_widths = [2, 3, 4, 5]
+  @testset "cell_width=$width" for width in cell_widths
+    model = Model("hubbard")
+    os = ITensorInfiniteMPS.opsum_finite(model, width; t, U)
+
+    connections = []
+    for term in os
+      push!(connections, ITensors.sites(term))
+    end
+    for i in 1:(width - 1)
+      @test [i, i + 1] ∈ connections
+      @test [i + 1, i] ∈ connections
+      @test [i] ∈ connections
+    end
+    @test [width] ∈ connections
+  end
+end
+
+nothing