update tests

Author: lkdvos

Project.toml

@@ -26,6 +26,7 @@ Accessors = "0.1"
 Aqua = "0.8.9"
 BlockTensorKit = "0.1.4"
 Compat = "3.47, 4.10"
+Combinatorics = "1"
 DocStringExtensions = "0.9.3"
 HalfIntegers = "1.6.0"
 KrylovKit = "0.8.3, 0.9.2"
@@ -48,11 +49,12 @@ julia = "1.10"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestExtras = "5ed8adda-3752-4e41-b88a-e8b09835ee3a"
 
 [targets]
-test = ["Aqua", "Pkg", "Test", "TestExtras", "Plots"]
+test = ["Aqua", "Pkg", "Test", "TestExtras", "Plots", "Combinatorics"]

src/operators/mpohamiltonian.jl

@@ -340,7 +340,7 @@ function Base.getproperty(H::MPOHamiltonian, sym::Symbol)
     elseif sym === :C
         return map(h -> h[1, 1, 1, 2:(end - 1)], parent(H))
     elseif sym === :D
-        return map(h -> h[1, 1, 1, end], parent(H))
+        return map(h -> h[1:1, 1, 1, end:end], parent(H))
     else
         return getfield(H, sym)
     end

test/algorithms.jl

@@ -7,7 +7,7 @@ module TestAlgorithms
 
 using ..TestSetup
 using Test, TestExtras
-using MPSKit
+using MPSKit: fuse_mul_mpo
 using TensorKit
 using TensorKit: ℙ
 
@@ -763,6 +763,85 @@ end
 end
 
 @testset "TaylorCluster time evolution" begin
+    L = 4
+    dt = 0.05
+    dτ = -im * dt
+
+    @testset "O(1) exact expression" begin
+        alg = TaylorCluster(; N=1, compression=true, extension=true)
+        for H in (transverse_field_ising(), heisenberg_XXX())
+            # Infinite
+            mpo = make_time_mpo(H, dt, alg)
+            O = mpo[1]
+            I, A, B, C, D = H[1][1], H.A[1], H.B[1], H.C[1], H.D[1]
+
+            @test size(O, 1) == size(D, 1)^2 + (size(D, 1) * size(B, 1))
+            @test size(O, 4) == size(D, 4)^2 + (size(C, 4) * size(D, 4))
+
+            O_exact = similar(O)
+            O_exact[1, 1, 1, 1] = I + dτ * D + dτ^2 / 2 * fuse_mul_mpo(D, D)
+            O_exact[1, 1, 1, 2] = C + dτ * symm_mul_mpo(C, D)
+            O_exact[2, 1, 1, 1] = dτ * (B + dτ * symm_mul_mpo(B, D))
+            O_exact[2, 1, 1, 2] = A + dτ * (symm_mul_mpo(A, D) + symm_mul_mpo(C, B))
+
+            @test all(isapprox.(parent(O), parent(O_exact)))
+
+            # Finite
+            H_fin = open_boundary_conditions(H, L)
+            mpo_fin = make_time_mpo(H_fin, dt, alg)
+            mpo_fin2 = open_boundary_conditions(mpo, L)
+            for i in 1:L
+                @test all(isapprox.(parent(mpo_fin[i]), parent(mpo_fin2[i])))
+            end
+        end
+    end
+
+    @testset "O(2) exact expression" begin
+        alg = TaylorCluster(; N=2, compression=true, extension=true)
+        for H in (transverse_field_ising(), heisenberg_XXX())
+            # Infinite
+            mpo = make_time_mpo(H, dt, alg)
+            O = mpo[1]
+            I, A, B, C, D = H[1][1], H.A[1], H.B[1], H.C[1], H.D[1]
+
+            @test size(O, 1) ==
+                  size(D, 1)^3 + (size(D, 1)^2 * size(B, 1)) + (size(D, 1) * size(B, 1)^2)
+            @test size(O, 4) ==
+                  size(D, 4)^3 + (size(C, 4) * size(D, 4)^2) + (size(D, 4) * size(C, 4)^2)
+
+            O_exact = similar(O)
+            O_exact[1, 1, 1, 1] = I + dτ * D + dτ^2 / 2 * fuse_mul_mpo(D, D) +
+                                  dτ^3 / 6 * fuse_mul_mpo(fuse_mul_mpo(D, D), D)
+            O_exact[1, 1, 1, 2] = C + dτ * symm_mul_mpo(C, D) +
+                                  dτ^2 / 2 * symm_mul_mpo(C, D, D)
+            O_exact[1, 1, 1, 3] = fuse_mul_mpo(C, C) + dτ * symm_mul_mpo(C, C, D)
+            O_exact[2, 1, 1, 1] = dτ *
+                                  (B + dτ * symm_mul_mpo(B, D) +
+                                   dτ^2 / 2 * symm_mul_mpo(B, D, D))
+            O_exact[2, 1, 1, 2] = A + dτ * symm_mul_mpo(A, D) +
+                                  dτ^2 / 2 * symm_mul_mpo(A, D, D) +
+                                  dτ * (symm_mul_mpo(C, B) + dτ * symm_mul_mpo(C, B, D))
+            O_exact[2, 1, 1, 3] = 2 * (symm_mul_mpo(A, C) + dτ * symm_mul_mpo(A, C, D)) +
+                                  dτ * symm_mul_mpo(C, C, B)
+            O_exact[3, 1, 1, 1] = dτ^2 / 2 *
+                                  (fuse_mul_mpo(B, B) + dτ * symm_mul_mpo(B, B, D))
+            O_exact[3, 1, 1, 2] = dτ * (symm_mul_mpo(A, B) + dτ * symm_mul_mpo(A, B, D)) +
+                                  dτ^2 / 2 * symm_mul_mpo(B, B, C)
+            O_exact[3, 1, 1, 3] = fuse_mul_mpo(A, A) + dτ * symm_mul_mpo(A, A, D) +
+                                  2 * dτ * symm_mul_mpo(A, C, B)
+
+            @test all(isapprox.(parent(O), parent(O_exact)))
+
+            # Finite
+            H_fin = open_boundary_conditions(H, L)
+            mpo_fin = make_time_mpo(H_fin, dt, alg)
+            mpo_fin2 = open_boundary_conditions(mpo, L)
+            for i in 1:L
+                @test all(isapprox.(parent(mpo_fin[i]), parent(mpo_fin2[i])))
+            end
+        end
+    end
+
     L = 4
     Hs = [transverse_field_ising(; L=L), heisenberg_XXX(; L=L)]
 
@@ -800,28 +879,4 @@ end
     end
 end
 
-@testset "TaylorCluster for non-square Hamiltonians" begin
-    # Make a MPOHamiltonian with very non-square Tensors
-    # Make different orders and kinds of TaylorCluster time-evolution MPOs
-    # Check the virtual dimensions match up between the MPOs
-
-    chain = fill(ℂ^2, 8)
-    O = randn(ComplexF64, ℂ^2 ← ℂ^2)
-
-    ops1 = fill((1, 2) => O ⊗ O, 5)
-    ops2 = [(2, 3) => O ⊗ O ⊗ O for i in 2:(length(chain) - 1)]
-    ops3 = [(3, 5) => O ⊗ O ⊗ O for i in 3:(length(chain) - 2)]
-
-    H = FiniteMPOHamiltonian(chain, ops1..., ops2..., ops3...)
-
-    for N in 1:3
-        for configuration in ((false, false), (false, true), (true, false), (true, true))
-            O = make_time_mpo(H, 0.1, TaylorCluster(N, configuration...))
-            for i in 1:(length(chain) - 1)
-                @test right_virtualspace(O, i) == left_virtualspace(O, i + 1)
-            end
-        end
-    end
-end
-
 end

test/setup.jl

@@ -8,10 +8,12 @@ using TensorKit
 using TensorKit: PlanarTrivial, ℙ, BraidingTensor
 using BlockTensorKit
 using LinearAlgebra: Diagonal
+using Combinatorics: permutations
 
 # exports
 export S_xx, S_yy, S_zz, S_x, S_y, S_z
 export force_planar
+export symm_mul_mpo
 export transverse_field_ising, heisenberg_XXX, bilinear_biquadratic_model
 export classical_ising, finite_classical_ising, sixvertex
 
@@ -58,6 +60,14 @@ end
 force_planar(mpo::MPOHamiltonian) = MPOHamiltonian(map(force_planar, parent(mpo)))
 force_planar(mpo::MPO) = MPO(map(force_planar, parent(mpo)))
 
+# sum of all permutations: {Os...}
+function symm_mul_mpo(Os::MPSKit.MPOTensor...)
+    N! = factorial(length(Os))
+    return sum(permutations(Os)) do os
+        return foldl(MPSKit.fuse_mul_mpo, os)
+    end / N!
+end
+
 # Toy models
 # ----------------------------