More tests

test/test_contract_mpo.jl

@@ -1,13 +1,26 @@
 @eval module $(gensym())
-using ITensors: ITensors, MPO, OpSum, apply, delta, inner, randomMPS, siteinds, truncate!
+using ITensors:
+  ITensors,
+  MPO,
+  OpSum,
+  apply,
+  dag,
+  delta,
+  denseblocks,
+  inner,
+  randomMPS,
+  siteinds,
+  truncate!
 using ITensorTDVP: ITensorTDVP
 using Test: @test, @testset
 @testset "Contract MPO (eltype=$elt)" for elt in (
-  Float32, Float64, Complex{Float32}, Complex{Float64}
-)
+    Float32, Float64, Complex{Float32}, Complex{Float64}
+  ),
+  conserve_qns in [false, true]
+
   N = 20
-  s = siteinds("S=1/2", N)
-  psi = randomMPS(elt, s; linkdims=8)
+  s = siteinds("S=1/2", N; conserve_qns)
+  psi = randomMPS(elt, s, j -> isodd(j) ? "↑" : "↓"; linkdims=8)
   os = OpSum()
   for j in 1:(N - 1)
     os += 0.5, "S+", j, "S-", j + 1
@@ -20,34 +33,26 @@ using Test: @test, @testset
     os += "Sz", j, "Sz", j + 2
   end
   H = MPO(elt, os, s)
-
-  # Test basic usage with default parameters
-  Hpsi = apply(H, psi; alg="fit")
+  Hpsi = apply(H, psi; alg="fit", nsweeps=2)
   @test ITensors.scalartype(Hpsi) == elt
-  @test inner(psi, Hpsi) ≈ inner(psi', H, psi) atol = 1e-5
-
-  #
+  @test inner(psi, Hpsi) ≈ inner(psi', H, psi) rtol = √eps(real(elt))
   # Change "top" indices of MPO to be a different set
-  #
-  t = siteinds("S=1/2", N)
+  t = siteinds("S=1/2", N; conserve_qns)
   psit = deepcopy(psi)
   for j in 1:N
-    H[j] *= delta(elt, s[j]', t[j])
-    psit[j] *= delta(elt, s[j], t[j])
+    H[j] *= delta(elt, dag(s[j])', t[j])
+    psit[j] *= delta(elt, dag(s[j]), t[j])
   end
-
   # Test with nsweeps=2
   Hpsi = apply(H, psi; alg="fit", nsweeps=2)
   @test ITensors.scalartype(Hpsi) == elt
-  @test inner(psit, Hpsi) ≈ inner(psit, H, psi) atol = 1e-5
-
+  @test inner(psit, Hpsi) ≈ inner(psit, H, psi) rtol = √eps(real(elt))
   # Test with less good initial guess MPS not equal to psi
   psi_guess = copy(psi)
   truncate!(psi_guess; maxdim=2)
   Hpsi = apply(H, psi; alg="fit", nsweeps=4, init_mps=psi_guess)
   @test ITensors.scalartype(Hpsi) == elt
   @test inner(psit, Hpsi) ≈ inner(psit, H, psi) atol = 1e-5
-
   # Test with nsite=1
   Hpsi_guess = apply(H, psi; alg="naive", cutoff=1E-4)
   Hpsi = apply(H, psi; alg="fit", init_mps=Hpsi_guess, nsite=1, nsweeps=2)

test/test_dmrg.jl

@@ -5,19 +5,20 @@ using Test: @test, @testset
 @testset "DMRG (eltype=$elt, nsite=$nsite)" for elt in (
     Float32, Float64, Complex{Float32}, Complex{Float64}
   ),
-  nsite in [1, 2]
+  nsite in [1, 2],
+  conserve_qns in [false, true]
 
   N = 10
   cutoff = eps(real(elt)) * 10^4
-  s = siteinds("S=1/2", N)
+  s = siteinds("S=1/2", N; conserve_qns)
   os = OpSum()
   for j in 1:(N - 1)
     os += 0.5, "S+", j, "S-", j + 1
     os += 0.5, "S-", j, "S+", j + 1
     os += "Sz", j, "Sz", j + 1
   end
   H = MPO(elt, os, s)
-  psi = randomMPS(elt, s; linkdims=20)
+  psi = randomMPS(elt, s, j -> isodd(j) ? "↑" : "↓"; linkdims=20)
   nsweeps = 10
   maxdim = [10, 20, 40, 100]
   psi = ITensorTDVP.dmrg(