From e877f707fbafb4858389dfa1202040bc1149df8d Mon Sep 17 00:00:00 2001 From: victor Date: Wed, 19 Feb 2025 11:27:22 +0100 Subject: [PATCH 01/13] remove all square-forcing code --- src/algorithms/timestep/timeevmpo.jl | 27 +++------------------------ 1 file changed, 3 insertions(+), 24 deletions(-) diff --git a/src/algorithms/timestep/timeevmpo.jl b/src/algorithms/timestep/timeevmpo.jl index 0467fdc82..7cec92bc5 100644 --- a/src/algorithms/timestep/timeevmpo.jl +++ b/src/algorithms/timestep/timeevmpo.jl @@ -59,30 +59,9 @@ function make_time_mpo(H::MPOHamiltonian, dt::Number, alg::TaylorCluster; N = alg.N τ = -1im * dt - # Hack to store FiniteMPOhamiltonians in "square" MPO tensors - if H isa FiniteMPOHamiltonian - H′ = copy(H) - H′[1] = similar(H[2]) - H′[end] = similar(H[end - 1]) - - for i in nonzero_keys(H[1]) - H′[1][i] = H[1][i] - end - for i in nonzero_keys(H[end]) - H′[end][:, 1, 1, end] = H[end][:, 1, 1, 1] - end - H′[1][end, 1, 1, end] += add_util_leg(id(space(H[1][end, 1, 1, end], 2))) - H′[end][1, 1, 1, 1] += add_util_leg(id(space(H[end][1, 1, 1, 1], 2))) - else - H′ = H - end - - # Check if mpo has the same size everywhere. This is assumed in the following. - @assert allequal(size.(H′)) "make_time_mpo assumes all mpo tensors to have equal size. A fix for this is yet to be implemented" - # start with H^N - H_n = H′^N - V = size(H′[1], 1) + H_n = H^N + V = size(H[1], 1) linds = LinearIndices(ntuple(i -> V, N)) cinds = CartesianIndices(linds) @@ -90,7 +69,7 @@ function make_time_mpo(H::MPOHamiltonian, dt::Number, alg::TaylorCluster; # incorporate higher order terms # TODO: don't need to fully construct H_next... if alg.extension - H_next = H_n * H′ + H_next = H_n * H linds_next = LinearIndices(ntuple(i -> V, N + 1)) for (i, slice) in enumerate(parent(H_n)) for a in cinds, b in cinds From d05410630d6f72d4e989b4726830bc66e45eb92b Mon Sep 17 00:00:00 2001 From: victor Date: Wed, 19 Feb 2025 12:03:39 +0100 Subject: [PATCH 02/13] add tests for `make_time_mpo` with non-square `FiniteMPOHamiltonians` --- test/algorithms.jl | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/test/algorithms.jl b/test/algorithms.jl index 199fdfb69..06d954d22 100644 --- a/test/algorithms.jl +++ b/test/algorithms.jl @@ -799,4 +799,29 @@ end 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 From 5f4864223a694b27ff6721330b5b560c6c62ae9d Mon Sep 17 00:00:00 2001 From: victor Date: Thu, 20 Feb 2025 09:08:22 +0100 Subject: [PATCH 03/13] Revert "remove all square-forcing code" This reverts commit bc37174d47cc448f4566ae3c688b3c00e15b3689. --- src/algorithms/timestep/timeevmpo.jl | 27 ++++++++++++++++++++++++--- 1 file changed, 24 insertions(+), 3 deletions(-) diff --git a/src/algorithms/timestep/timeevmpo.jl b/src/algorithms/timestep/timeevmpo.jl index 7cec92bc5..0467fdc82 100644 --- a/src/algorithms/timestep/timeevmpo.jl +++ b/src/algorithms/timestep/timeevmpo.jl @@ -59,9 +59,30 @@ function make_time_mpo(H::MPOHamiltonian, dt::Number, alg::TaylorCluster; N = alg.N τ = -1im * dt + # Hack to store FiniteMPOhamiltonians in "square" MPO tensors + if H isa FiniteMPOHamiltonian + H′ = copy(H) + H′[1] = similar(H[2]) + H′[end] = similar(H[end - 1]) + + for i in nonzero_keys(H[1]) + H′[1][i] = H[1][i] + end + for i in nonzero_keys(H[end]) + H′[end][:, 1, 1, end] = H[end][:, 1, 1, 1] + end + H′[1][end, 1, 1, end] += add_util_leg(id(space(H[1][end, 1, 1, end], 2))) + H′[end][1, 1, 1, 1] += add_util_leg(id(space(H[end][1, 1, 1, 1], 2))) + else + H′ = H + end + + # Check if mpo has the same size everywhere. This is assumed in the following. + @assert allequal(size.(H′)) "make_time_mpo assumes all mpo tensors to have equal size. A fix for this is yet to be implemented" + # start with H^N - H_n = H^N - V = size(H[1], 1) + H_n = H′^N + V = size(H′[1], 1) linds = LinearIndices(ntuple(i -> V, N)) cinds = CartesianIndices(linds) @@ -69,7 +90,7 @@ function make_time_mpo(H::MPOHamiltonian, dt::Number, alg::TaylorCluster; # incorporate higher order terms # TODO: don't need to fully construct H_next... if alg.extension - H_next = H_n * H + H_next = H_n * H′ linds_next = LinearIndices(ntuple(i -> V, N + 1)) for (i, slice) in enumerate(parent(H_n)) for a in cinds, b in cinds From 2bad5031589b59cbe1f226cba61e144f60bfa6c3 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 20 Feb 2025 19:00:06 -0500 Subject: [PATCH 04/13] attempt to fix --- src/algorithms/timestep/timeevmpo.jl | 159 +++++++++++++++++---------- 1 file changed, 100 insertions(+), 59 deletions(-) diff --git a/src/algorithms/timestep/timeevmpo.jl b/src/algorithms/timestep/timeevmpo.jl index 0467fdc82..02a853ff0 100644 --- a/src/algorithms/timestep/timeevmpo.jl +++ b/src/algorithms/timestep/timeevmpo.jl @@ -55,60 +55,73 @@ Construct an MPO that approximates ``\\exp(-iHdt)``. """ make_time_mpo function make_time_mpo(H::MPOHamiltonian, dt::Number, alg::TaylorCluster; - tol=eps(real(scalartype(H)))^(3 / 4)) + tol=eps(real(scalartype(H)))) N = alg.N τ = -1im * dt # Hack to store FiniteMPOhamiltonians in "square" MPO tensors if H isa FiniteMPOHamiltonian H′ = copy(H) - H′[1] = similar(H[2]) - H′[end] = similar(H[end - 1]) - for i in nonzero_keys(H[1]) - H′[1][i] = H[1][i] + V_left = left_virtualspace(H[1]) + V_left′ = BlockTensorKit.oplus(V_left, oneunit(V_left), oneunit(V_left)) + H′[1] = similar(H[1], V_left′ ⊗ space(H[1], 2) ← domain(H[1])) + for (I, v) in nonzero_pairs(H[1]) + H′[1][I] = v end - for i in nonzero_keys(H[end]) - H′[end][:, 1, 1, end] = H[end][:, 1, 1, 1] + + V_right = right_virtualspace(H[end]) + V_right′ = BlockTensorKit.oplus(oneunit(V_right), oneunit(V_right), V_right) + H′[end] = similar(H[end], codomain(H[end]) ← space(H[end], 3)' ⊗ V_right′) + for (I, v) in nonzero_pairs(H[end]) + H′[end][I[1], 1, 1, end] = v end - H′[1][end, 1, 1, end] += add_util_leg(id(space(H[1][end, 1, 1, end], 2))) - H′[end][1, 1, 1, 1] += add_util_leg(id(space(H[end][1, 1, 1, 1], 2))) + + H′[1][end, 1, 1, end] = H′[1][1, 1, 1, 1] + H′[end][1, 1, 1, 1] = H′[end][end, 1, 1, end] else H′ = H end # Check if mpo has the same size everywhere. This is assumed in the following. - @assert allequal(size.(H′)) "make_time_mpo assumes all mpo tensors to have equal size. A fix for this is yet to be implemented" + # @assert allequal(size.(H′)) "make_time_mpo assumes all mpo tensors to have equal size. A fix for this is yet to be implemented" # start with H^N H_n = H′^N - V = size(H′[1], 1) - linds = LinearIndices(ntuple(i -> V, N)) - cinds = CartesianIndices(linds) + virtual_sz = map(0:length(H′)) do i + return i == 0 ? size(H′[1], 1) : size(H′[i], 4) + end # extension step: Algorithm 3 # incorporate higher order terms # TODO: don't need to fully construct H_next... if alg.extension H_next = H_n * H′ - linds_next = LinearIndices(ntuple(i -> V, N + 1)) for (i, slice) in enumerate(parent(H_n)) - for a in cinds, b in cinds + V_left = virtual_sz[i] + V_right = virtual_sz[i + 1] + linds_left = LinearIndices(ntuple(Returns(V_left), N)) + linds_right = LinearIndices(ntuple(Returns(V_right), N)) + linds_next_left = LinearIndices(ntuple(Returns(V_left), N + 1)) + linds_next_right = LinearIndices(ntuple(Returns(V_right), N + 1)) + + for a in CartesianIndices(linds_left), b in CartesianIndices(linds_right) all(>(1), b.I) || continue - all(in((1, V)), a.I) && any(==(V), a.I) && continue + all(in((1, V_left)), a.I) && any(==(V_left), a.I) && continue n1 = count(==(1), a.I) + 1 - n3 = count(==(V), b.I) + 1 + n3 = count(==(V_right), b.I) + 1 factor = τ * factorial(N) / (factorial(N + 1) * n1 * n3) for c in 1:(N + 1), d in 1:(N + 1) aₑ = insert!([a.I...], c, 1) - bₑ = insert!([b.I...], d, V) + bₑ = insert!([b.I...], d, V_right) # TODO: use VectorInterface for memory efficiency - slice[linds[a], 1, 1, linds[b]] += factor * - H_next[i][linds_next[aₑ...], 1, 1, - linds_next[bₑ...]] + slice[linds_left[a], 1, 1, linds_right[b]] += factor * + H_next[i][linds_next_left[aₑ...], + 1, 1, + linds_next_right[bₑ...]] end end end @@ -117,12 +130,15 @@ function make_time_mpo(H::MPOHamiltonian, dt::Number, alg::TaylorCluster; # loopback step: Algorithm 1 # constructing the Nth order time evolution MPO mpo = MPO(parent(H_n)) - for slice in parent(mpo) - for b in cinds[2:end] - all(in((1, V)), b.I) || continue - - b_lin = linds[b] - a = count(==(V), b.I) + for (i, slice) in enumerate(parent(mpo)) + V_right = virtual_sz[i + 1] + linds_right = LinearIndices(ntuple(Returns(V_right), N)) + cinds_right = CartesianIndices(linds_right) + for b in cinds_right[2:end] + all(in((1, V_right)), b.I) || continue + + b_lin = linds_right[b] + a = count(==(V_right), b.I) factor = τ^a * factorial(N - a) / factorial(N) slice[:, 1, 1, 1] = slice[:, 1, 1, 1] + factor * slice[:, 1, 1, b_lin] for I in nonzero_keys(slice) @@ -132,31 +148,43 @@ function make_time_mpo(H::MPOHamiltonian, dt::Number, alg::TaylorCluster; end # Remove equivalent rows and columns: Algorithm 2 - for slice in parent(mpo) - for c in cinds - c_lin = linds[c] + for (i, slice) in enumerate(parent(mpo)) + V_left = virtual_sz[i] + linds_left = LinearIndices(ntuple(Returns(V_left), N)) + for c in CartesianIndices(linds_left) + c_lin = linds_left[c] s_c = CartesianIndex(sort(collect(c.I); by=(!=(1)))...) - s_r = CartesianIndex(sort(collect(c.I); by=(!=(V)))...) n1 = count(==(1), c.I) - n3 = count(==(V), c.I) - - if n3 <= n1 && s_c != c || n3 > n1 && s_r != c - if n3 <= n1 && s_c != c - for k in 1:size(slice, 4) - if CartesianIndex(c_lin, 1, 1, k) in nonzero_keys(slice) - slice[linds[s_c], 1, 1, k] += slice[c_lin, 1, 1, k] - end - end - elseif n3 > n1 && s_r != c - for k in 1:size(slice, 1) - if CartesianIndex(k, 1, 1, c_lin) in nonzero_keys(slice) - slice[k, 1, 1, linds[s_r]] += slice[k, 1, 1, c_lin] - end + n3 = count(==(V_left), c.I) + + if n3 <= n1 && s_c != c + for k in 1:size(slice, 4) + I = CartesianIndex(c_lin, 1, 1, k) + if I in nonzero_keys(slice) + slice[linds_left[s_c], 1, 1, k] += slice[I] + delete!(slice, I) end end - for I in nonzero_keys(slice) - (I[1] == c_lin || I[4] == c_lin) && delete!(slice, I) + end + end + + V_right = virtual_sz[i + 1] + linds_right = LinearIndices(ntuple(Returns(V_right), N)) + for c in CartesianIndices(linds_right) + c_lin = linds_right[c] + s_r = CartesianIndex(sort(collect(c.I); by=(!=(V_right)))...) + + n1 = count(==(1), c.I) + n3 = count(==(V_right), c.I) + + if n3 > n1 && s_r != c + for k in 1:size(slice, 1) + I = CartesianIndex(k, 1, 1, c_lin) + if I in nonzero_keys(slice) + slice[k, 1, 1, linds_right[s_r]] += slice[I] + delete!(slice, I) + end end end end @@ -164,22 +192,35 @@ function make_time_mpo(H::MPOHamiltonian, dt::Number, alg::TaylorCluster; # Approximate compression step: Algorithm 4 if alg.compression - for slice in parent(mpo) - for a in cinds + for (i, slice) in enumerate(parent(mpo)) + V_right = virtual_sz[i + 1] + linds_right = LinearIndices(ntuple(Returns(V_right), N)) + for a in CartesianIndices(linds_right) all(>(1), a.I) || continue - a_lin = linds[a] - n1 = count(==(V), a.I) - b = CartesianIndex(replace(a.I, V => 1)) - b_lin = linds[b] + a_lin = linds_right[a] + n1 = count(==(V_right), a.I) + n1 == 0 && continue + b = CartesianIndex(replace(a.I, V_right => 1)) + b_lin = linds_right[b] factor = τ^n1 * factorial(N - n1) / factorial(N) for k in 1:size(slice, 1) - if CartesianIndex(k, 1, 1, a_lin) in nonzero_keys(slice) - slice[k, 1, 1, b_lin] += factor * slice[k, 1, 1, a_lin] + I = CartesianIndex(k, 1, 1, a_lin) + if I in nonzero_keys(slice) + slice[k, 1, 1, b_lin] += factor * slice[I] + delete!(slice, I) end end - - for I in nonzero_keys(slice) - (I[1] == a_lin || I[4] == a_lin) && delete!(slice, I) + end + V_left = virtual_sz[i] + linds_left = LinearIndices(ntuple(Returns(V_left), N)) + for a in CartesianIndices(linds_left) + all(>(1), a.I) || continue + a_lin = linds_left[a] + n1 = count(==(V_left), a.I) + n1 == 0 && continue + for k in 1:size(slice, 4) + I = CartesianIndex(a_lin, 1, 1, k) + delete!(slice, I) end end end @@ -191,7 +232,7 @@ function make_time_mpo(H::MPOHamiltonian, dt::Number, alg::TaylorCluster; mpo[end] = mpo[end][:, :, :, 1] end - return remove_orphans!(mpo; tol=tol) + return remove_orphans!(mpo; tol) end has_prod_elem(slice, t1, t2) = all(map(x -> contains(slice, x...), zip(t1, t2))) From 913c4e85d3ce9f4724a5ef87759fa11e06c76df5 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 26 Feb 2025 07:45:04 -0500 Subject: [PATCH 05/13] fix size issue in MPO --- src/operators/abstractmpo.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/operators/abstractmpo.jl b/src/operators/abstractmpo.jl index 4236e2a02..b0de56ce0 100644 --- a/src/operators/abstractmpo.jl +++ b/src/operators/abstractmpo.jl @@ -72,7 +72,7 @@ function remove_orphans!(mpo::SparseMPO; tol=eps(real(scalartype(mpo)))^(3 / 4)) for i in 1:length(mpo) # slice empty columns on right or empty rows on left mask = filter(1:size(mpo[i], 4)) do j - return j ∈ getindex.(nonzero_keys(mpo[i]), 1) || + return j ∈ getindex.(nonzero_keys(mpo[i]), 1) && j ∈ getindex.(nonzero_keys(mpo[i + 1]), 4) end changed |= length(mask) == size(mpo[i], 4) From fef41ad9add3cee3513ccd3efd86c1a1248e43b8 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 26 Feb 2025 10:03:36 -0500 Subject: [PATCH 06/13] performance improvement --- src/operators/abstractmpo.jl | 21 +++++++++++++++++---- 1 file changed, 17 insertions(+), 4 deletions(-) diff --git a/src/operators/abstractmpo.jl b/src/operators/abstractmpo.jl index b0de56ce0..938b26a5c 100644 --- a/src/operators/abstractmpo.jl +++ b/src/operators/abstractmpo.jl @@ -196,10 +196,23 @@ end # TODO: diagram function _fuse_mpo_mpo(O1::MPOTensor, O2::MPOTensor, Fₗ, Fᵣ) - return @plansor O′[-1 -2; -3 -4] := Fₗ[-1; 1 2] * - O2[1 3; -3 5] * - O1[2 -2; 3 4] * - conj(Fᵣ[-4; 5 4]) + return if O1 isa BraidingTensor && O2 isa BraidingTensor + elseif O1 isa BraidingTensor + @plansor O′[-1 -2; -3 -4] := Fₗ[-1; 1 2] * + O2[1 3; -3 5] * + τ[2 -2; 3 4] * + conj(Fᵣ[-4; 5 4]) + elseif O2 isa BraidingTensor + @plansor O′[-1 -2; -3 -4] := Fₗ[-1; 1 2] * + τ[1 3; -3 5] * + O1[2 -2; 3 4] * + conj(Fᵣ[-4; 5 4]) + else + @plansor O′[-1 -2; -3 -4] := Fₗ[-1; 1 2] * + O2[1 3; -3 5] * + O1[2 -2; 3 4] * + conj(Fᵣ[-4; 5 4]) + end end """ From f3b97b56880f72ae2403785c46b84881c0a090e3 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 26 Feb 2025 10:26:13 -0500 Subject: [PATCH 07/13] update tests --- Project.toml | 4 +- src/operators/mpohamiltonian.jl | 2 +- test/algorithms.jl | 104 ++++++++++++++++++++++++-------- test/setup.jl | 10 +++ 4 files changed, 94 insertions(+), 26 deletions(-) diff --git a/Project.toml b/Project.toml index e6d048674..1e554931c 100644 --- a/Project.toml +++ b/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,6 +49,7 @@ 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" @@ -55,4 +57,4 @@ 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"] diff --git a/src/operators/mpohamiltonian.jl b/src/operators/mpohamiltonian.jl index f5ba931f0..8e1133393 100644 --- a/src/operators/mpohamiltonian.jl +++ b/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 diff --git a/test/algorithms.jl b/test/algorithms.jl index 06d954d22..26ca46c47 100644 --- a/test/algorithms.jl +++ b/test/algorithms.jl @@ -8,6 +8,7 @@ module TestAlgorithms using ..TestSetup using Test, TestExtras using MPSKit +using MPSKit: fuse_mul_mpo using TensorKit using TensorKit: ℙ @@ -763,6 +764,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 +880,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 diff --git a/test/setup.jl b/test/setup.jl index 138f90549..804ca6e3a 100644 --- a/test/setup.jl +++ b/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 # ---------------------------- From 70c644126e69e976a87efd9a47588a072f8326e6 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 26 Feb 2025 11:45:04 -0500 Subject: [PATCH 08/13] Nth order has (N+1)th order error --- test/algorithms.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/algorithms.jl b/test/algorithms.jl index 26ca46c47..05fd89973 100644 --- a/test/algorithms.jl +++ b/test/algorithms.jl @@ -875,7 +875,7 @@ end ψ₁, _ = approximate(ψ₀, (O, ψ₀), DMRG(; verbosity=0)) εs[i] = norm(dot(ψ₀, ψ₁) - exp(-im * E₀ * dt)) end - @test (log(εs[2]) - log(εs[1])) / (log(dts[2]) - log(dts[1])) ≈ N atol = 0.1 + @test (log(εs[2]) - log(εs[1])) / (log(dts[2]) - log(dts[1])) ≈ N + 1 atol = 0.1 end end end From 80cca0f84b5e5c925e52ae0d5126441b6a632872 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Fri, 28 Feb 2025 17:22:06 -0500 Subject: [PATCH 09/13] clean-up --- src/MPSKit.jl | 3 +- src/algorithms/timestep/taylorcluster.jl | 209 +++++++++++++ src/algorithms/timestep/time_evolve.jl | 6 + src/algorithms/timestep/timeevmpo.jl | 357 ----------------------- src/algorithms/timestep/wii.jl | 101 +++++++ 5 files changed, 318 insertions(+), 358 deletions(-) create mode 100644 src/algorithms/timestep/taylorcluster.jl delete mode 100644 src/algorithms/timestep/timeevmpo.jl create mode 100644 src/algorithms/timestep/wii.jl diff --git a/src/MPSKit.jl b/src/MPSKit.jl index 1fe569db2..b4787e239 100644 --- a/src/MPSKit.jl +++ b/src/MPSKit.jl @@ -143,7 +143,8 @@ include("algorithms/changebonds/svdcut.jl") include("algorithms/changebonds/randexpand.jl") include("algorithms/timestep/tdvp.jl") -include("algorithms/timestep/timeevmpo.jl") +include("algorithms/timestep/taylorcluster.jl") +include("algorithms/timestep/wii.jl") include("algorithms/timestep/integrators.jl") include("algorithms/timestep/time_evolve.jl") diff --git a/src/algorithms/timestep/taylorcluster.jl b/src/algorithms/timestep/taylorcluster.jl new file mode 100644 index 000000000..edab759a9 --- /dev/null +++ b/src/algorithms/timestep/taylorcluster.jl @@ -0,0 +1,209 @@ +""" +$(TYPEDEF) + +Algorithm for constructing the `N`th order time evolution MPO using the Taylor cluster expansion. + +## Fields + +$(TYPEDFIELDS) + +## References + +* [Van Damme et al. SciPost Phys. 17 (2024)](@cite vandamme2024) +""" +@kwdef struct TaylorCluster <: Algorithm + "order of the Taylor expansion" + N::Int = 2 + "include higher-order corrections" + extension::Bool = true + "approximate compression of corrections, accurate up to order `N`" + compression::Bool = false +end + +""" + const WI = TaylorCluster(; N=1, extension=false, compression=false) + +First order Taylor expansion for a time-evolution MPO. +""" +const WI = TaylorCluster(; N=1, extension=false, compression=false) + +function make_time_mpo(H::MPOHamiltonian, dt::Number, alg::TaylorCluster; + tol=eps(real(scalartype(H)))) + N = alg.N + τ = -1im * dt + + # start with H^N + H_n = H^N + virtual_sz = map(0:length(H)) do i + return i == 0 ? size(H[1], 1) : size(H[i], 4) + end + linds = map(virtual_sz) do sz + return LinearIndices(ntuple(Returns(sz), N)) + end + + # extension step: Algorithm 3 + # incorporate higher order terms + # TODO: don't need to fully construct H_next... + if alg.extension + H_next = H_n * H + for (i, slice) in enumerate(parent(H_n)) + linds_left = linds[i] + linds_right = linds[i + 1] + V_left = virtual_sz[i] + V_right = virtual_sz[i + 1] + linds_next_left = LinearIndices(ntuple(Returns(V_left), N + 1)) + linds_next_right = LinearIndices(ntuple(Returns(V_right), N + 1)) + + for a in CartesianIndices(linds_left), b in CartesianIndices(linds_right) + all(>(1), b.I) || continue + all(in((1, V_left)), a.I) && any(==(V_left), a.I) && continue + + n1 = count(==(1), a.I) + 1 + n3 = count(==(V_right), b.I) + 1 + factor = τ * factorial(N) / (factorial(N + 1) * n1 * n3) + + for c in 1:(N + 1), d in 1:(N + 1) + aₑ = insert!([a.I...], c, 1) + bₑ = insert!([b.I...], d, V_right) + + # TODO: use VectorInterface for memory efficiency + slice[linds_left[a], 1, 1, linds_right[b]] += factor * + H_next[i][linds_next_left[aₑ...], + 1, 1, + linds_next_right[bₑ...]] + end + end + end + end + + # loopback step: Algorithm 1 + # constructing the Nth order time evolution MPO + mpo = MPO(parent(H_n)) + for (i, slice) in enumerate(parent(mpo)) + V_right = virtual_sz[i + 1] + linds_right = linds[i + 1] + cinds_right = CartesianIndices(linds_right) + for b in cinds_right[2:end] + all(in((1, V_right)), b.I) || continue + + b_lin = linds_right[b] + a = count(==(V_right), b.I) + factor = τ^a * factorial(N - a) / factorial(N) + slice[:, 1, 1, 1] = slice[:, 1, 1, 1] + factor * slice[:, 1, 1, b_lin] + for I in nonzero_keys(slice) + (I[1] == b_lin || I[4] == b_lin) && delete!(slice, I) + end + end + end + + # Remove equivalent rows and columns: Algorithm 2 + for (i, slice) in enumerate(parent(mpo)) + V_left = virtual_sz[i] + linds_left = linds[i] + for c in CartesianIndices(linds_left) + c_lin = linds_left[c] + s_c = CartesianIndex(sort(collect(c.I); by=(!=(1)))...) + + n1 = count(==(1), c.I) + n3 = count(==(V_left), c.I) + + if n3 <= n1 && s_c != c + for k in 1:size(slice, 4) + I = CartesianIndex(c_lin, 1, 1, k) + if I in nonzero_keys(slice) + slice[linds_left[s_c], 1, 1, k] += slice[I] + delete!(slice, I) + end + end + end + end + + V_right = virtual_sz[i + 1] + linds_right = linds[i + 1] + for c in CartesianIndices(linds_right) + c_lin = linds_right[c] + s_r = CartesianIndex(sort(collect(c.I); by=(!=(V_right)))...) + + n1 = count(==(1), c.I) + n3 = count(==(V_right), c.I) + + if n3 > n1 && s_r != c + for k in 1:size(slice, 1) + I = CartesianIndex(k, 1, 1, c_lin) + if I in nonzero_keys(slice) + slice[k, 1, 1, linds_right[s_r]] += slice[I] + delete!(slice, I) + end + end + end + end + end + + # Approximate compression step: Algorithm 4 + if alg.compression + for (i, slice) in enumerate(parent(mpo)) + V_right = virtual_sz[i + 1] + linds_right = linds[i + 1] + for a in CartesianIndices(linds_right) + all(>(1), a.I) || continue + a_lin = linds_right[a] + n1 = count(==(V_right), a.I) + n1 == 0 && continue + b = CartesianIndex(replace(a.I, V_right => 1)) + b_lin = linds_right[b] + factor = τ^n1 * factorial(N - n1) / factorial(N) + for k in 1:size(slice, 1) + I = CartesianIndex(k, 1, 1, a_lin) + if I in nonzero_keys(slice) + slice[k, 1, 1, b_lin] += factor * slice[I] + delete!(slice, I) + end + end + end + V_left = virtual_sz[i] + linds_left = linds[i] + for a in CartesianIndices(linds_left) + all(>(1), a.I) || continue + a_lin = linds_left[a] + n1 = count(==(V_left), a.I) + n1 == 0 && continue + for k in 1:size(slice, 4) + I = CartesianIndex(a_lin, 1, 1, k) + delete!(slice, I) + end + end + end + end + + return remove_orphans!(mpo; tol) +end + +# Hack to treat FiniteMPOhamiltonians as Infinite +function make_time_mpo(H::FiniteMPOHamiltonian, dt::Number, alg::TaylorCluster; + tol=eps(real(scalartype(H)))) + H′ = copy(parent(H)) + + V_left = left_virtualspace(H[1]) + V_left′ = BlockTensorKit.oplus(V_left, oneunit(V_left), oneunit(V_left)) + H′[1] = similar(H[1], V_left′ ⊗ space(H[1], 2) ← domain(H[1])) + for (I, v) in nonzero_pairs(H[1]) + H′[1][I] = v + end + + V_right = right_virtualspace(H[end]) + V_right′ = BlockTensorKit.oplus(oneunit(V_right), oneunit(V_right), V_right) + H′[end] = similar(H[end], codomain(H[end]) ← space(H[end], 3)' ⊗ V_right′) + for (I, v) in nonzero_pairs(H[end]) + H′[end][I[1], 1, 1, end] = v + end + + H′[1][end, 1, 1, end] = H′[1][1, 1, 1, 1] + H′[end][1, 1, 1, 1] = H′[end][end, 1, 1, end] + + mpo = make_time_mpo(InfiniteMPOHamiltonian(H′), dt, alg; tol) + + # Impose boundary conditions + mpo[1] = mpo[1][1, :, :, :] + mpo[end] = mpo[end][:, :, :, 1] + return remove_orphans(mpo; tol) +end diff --git a/src/algorithms/timestep/time_evolve.jl b/src/algorithms/timestep/time_evolve.jl index 5c0753510..a5d22dc2f 100644 --- a/src/algorithms/timestep/time_evolve.jl +++ b/src/algorithms/timestep/time_evolve.jl @@ -52,3 +52,9 @@ solving the Schroedinger equation: ``i ∂ψ/∂t = H ψ``. - `[envs]`: MPS environment manager """ function timestep end, function timestep! end + +@doc """ + make_time_mpo(H::MPOHamiltonian, dt::Number, alg) -> O::MPO + +Construct an `MPO` that approximates ``\\exp(-iHdt)``. +""" make_time_mpo diff --git a/src/algorithms/timestep/timeevmpo.jl b/src/algorithms/timestep/timeevmpo.jl deleted file mode 100644 index 02a853ff0..000000000 --- a/src/algorithms/timestep/timeevmpo.jl +++ /dev/null @@ -1,357 +0,0 @@ -""" -$(TYPEDEF) - -Generalization of the Euler approximation of the operator exponential for MPOs. - -## Fields - -$(TYPEDFIELDS) - -## References - -* [Zaletel et al. Phys. Rev. B 91 (2015)](@cite zaletel2015) -* [Paeckel et al. Ann. of Phys. 411 (2019)](@cite paeckel2019) -""" -@kwdef struct WII <: Algorithm - "tolerance for convergence criterium" - tol::Float64 = Defaults.tol - "maximal number of iterations" - maxiter::Int = Defaults.maxiter -end - -""" -$(TYPEDEF) - -Algorithm for constructing the `N`th order time evolution MPO using the Taylor cluster expansion. - -## Fields - -$(TYPEDFIELDS) - -## References - -* [Van Damme et al. SciPost Phys. 17 (2024)](@cite vandamme2024) -""" -@kwdef struct TaylorCluster <: Algorithm - "order of the Taylor expansion" - N::Int = 2 - "include higher-order corrections" - extension::Bool = true - "approximate compression of corrections, accurate up to order `N`" - compression::Bool = false -end - -""" - const WI = TaylorCluster(; N=1, extension=false, compression=false) - -First order Taylor expansion for a time-evolution MPO. -""" -const WI = TaylorCluster(; N=1, extension=false, compression=false) - -@doc """ - make_time_mpo(H::MPOHamiltonian, dt::Number, alg) - -Construct an MPO that approximates ``\\exp(-iHdt)``. -""" make_time_mpo - -function make_time_mpo(H::MPOHamiltonian, dt::Number, alg::TaylorCluster; - tol=eps(real(scalartype(H)))) - N = alg.N - τ = -1im * dt - - # Hack to store FiniteMPOhamiltonians in "square" MPO tensors - if H isa FiniteMPOHamiltonian - H′ = copy(H) - - V_left = left_virtualspace(H[1]) - V_left′ = BlockTensorKit.oplus(V_left, oneunit(V_left), oneunit(V_left)) - H′[1] = similar(H[1], V_left′ ⊗ space(H[1], 2) ← domain(H[1])) - for (I, v) in nonzero_pairs(H[1]) - H′[1][I] = v - end - - V_right = right_virtualspace(H[end]) - V_right′ = BlockTensorKit.oplus(oneunit(V_right), oneunit(V_right), V_right) - H′[end] = similar(H[end], codomain(H[end]) ← space(H[end], 3)' ⊗ V_right′) - for (I, v) in nonzero_pairs(H[end]) - H′[end][I[1], 1, 1, end] = v - end - - H′[1][end, 1, 1, end] = H′[1][1, 1, 1, 1] - H′[end][1, 1, 1, 1] = H′[end][end, 1, 1, end] - else - H′ = H - end - - # Check if mpo has the same size everywhere. This is assumed in the following. - # @assert allequal(size.(H′)) "make_time_mpo assumes all mpo tensors to have equal size. A fix for this is yet to be implemented" - - # start with H^N - H_n = H′^N - virtual_sz = map(0:length(H′)) do i - return i == 0 ? size(H′[1], 1) : size(H′[i], 4) - end - - # extension step: Algorithm 3 - # incorporate higher order terms - # TODO: don't need to fully construct H_next... - if alg.extension - H_next = H_n * H′ - for (i, slice) in enumerate(parent(H_n)) - V_left = virtual_sz[i] - V_right = virtual_sz[i + 1] - linds_left = LinearIndices(ntuple(Returns(V_left), N)) - linds_right = LinearIndices(ntuple(Returns(V_right), N)) - linds_next_left = LinearIndices(ntuple(Returns(V_left), N + 1)) - linds_next_right = LinearIndices(ntuple(Returns(V_right), N + 1)) - - for a in CartesianIndices(linds_left), b in CartesianIndices(linds_right) - all(>(1), b.I) || continue - all(in((1, V_left)), a.I) && any(==(V_left), a.I) && continue - - n1 = count(==(1), a.I) + 1 - n3 = count(==(V_right), b.I) + 1 - factor = τ * factorial(N) / (factorial(N + 1) * n1 * n3) - - for c in 1:(N + 1), d in 1:(N + 1) - aₑ = insert!([a.I...], c, 1) - bₑ = insert!([b.I...], d, V_right) - - # TODO: use VectorInterface for memory efficiency - slice[linds_left[a], 1, 1, linds_right[b]] += factor * - H_next[i][linds_next_left[aₑ...], - 1, 1, - linds_next_right[bₑ...]] - end - end - end - end - - # loopback step: Algorithm 1 - # constructing the Nth order time evolution MPO - mpo = MPO(parent(H_n)) - for (i, slice) in enumerate(parent(mpo)) - V_right = virtual_sz[i + 1] - linds_right = LinearIndices(ntuple(Returns(V_right), N)) - cinds_right = CartesianIndices(linds_right) - for b in cinds_right[2:end] - all(in((1, V_right)), b.I) || continue - - b_lin = linds_right[b] - a = count(==(V_right), b.I) - factor = τ^a * factorial(N - a) / factorial(N) - slice[:, 1, 1, 1] = slice[:, 1, 1, 1] + factor * slice[:, 1, 1, b_lin] - for I in nonzero_keys(slice) - (I[1] == b_lin || I[4] == b_lin) && delete!(slice, I) - end - end - end - - # Remove equivalent rows and columns: Algorithm 2 - for (i, slice) in enumerate(parent(mpo)) - V_left = virtual_sz[i] - linds_left = LinearIndices(ntuple(Returns(V_left), N)) - for c in CartesianIndices(linds_left) - c_lin = linds_left[c] - s_c = CartesianIndex(sort(collect(c.I); by=(!=(1)))...) - - n1 = count(==(1), c.I) - n3 = count(==(V_left), c.I) - - if n3 <= n1 && s_c != c - for k in 1:size(slice, 4) - I = CartesianIndex(c_lin, 1, 1, k) - if I in nonzero_keys(slice) - slice[linds_left[s_c], 1, 1, k] += slice[I] - delete!(slice, I) - end - end - end - end - - V_right = virtual_sz[i + 1] - linds_right = LinearIndices(ntuple(Returns(V_right), N)) - for c in CartesianIndices(linds_right) - c_lin = linds_right[c] - s_r = CartesianIndex(sort(collect(c.I); by=(!=(V_right)))...) - - n1 = count(==(1), c.I) - n3 = count(==(V_right), c.I) - - if n3 > n1 && s_r != c - for k in 1:size(slice, 1) - I = CartesianIndex(k, 1, 1, c_lin) - if I in nonzero_keys(slice) - slice[k, 1, 1, linds_right[s_r]] += slice[I] - delete!(slice, I) - end - end - end - end - end - - # Approximate compression step: Algorithm 4 - if alg.compression - for (i, slice) in enumerate(parent(mpo)) - V_right = virtual_sz[i + 1] - linds_right = LinearIndices(ntuple(Returns(V_right), N)) - for a in CartesianIndices(linds_right) - all(>(1), a.I) || continue - a_lin = linds_right[a] - n1 = count(==(V_right), a.I) - n1 == 0 && continue - b = CartesianIndex(replace(a.I, V_right => 1)) - b_lin = linds_right[b] - factor = τ^n1 * factorial(N - n1) / factorial(N) - for k in 1:size(slice, 1) - I = CartesianIndex(k, 1, 1, a_lin) - if I in nonzero_keys(slice) - slice[k, 1, 1, b_lin] += factor * slice[I] - delete!(slice, I) - end - end - end - V_left = virtual_sz[i] - linds_left = LinearIndices(ntuple(Returns(V_left), N)) - for a in CartesianIndices(linds_left) - all(>(1), a.I) || continue - a_lin = linds_left[a] - n1 = count(==(V_left), a.I) - n1 == 0 && continue - for k in 1:size(slice, 4) - I = CartesianIndex(a_lin, 1, 1, k) - delete!(slice, I) - end - end - end - end - - # Impose boundary conditions as in arXiv:2302.14181 - if mpo isa FiniteMPO - mpo[1] = mpo[1][1, :, :, :] - mpo[end] = mpo[end][:, :, :, 1] - end - - return remove_orphans!(mpo; tol) -end - -has_prod_elem(slice, t1, t2) = all(map(x -> contains(slice, x...), zip(t1, t2))) -function calc_prod_elem(slice, t1, t2) - return calc_prod_elem(slice[first(t1), first(t2)], slice, t1[2:end], t2[2:end]) -end -function calc_prod_elem(o, slice, t1, t2) - isempty(t1) && return o - - nel = slice[first(t1), first(t2)] - fuse_front = isomorphism(fuse(_firstspace(o) * _firstspace(nel)), - _firstspace(o) * _firstspace(nel)) - fuse_back = isomorphism(fuse(_lastspace(o)' * _lastspace(nel)'), - _lastspace(o)' * _lastspace(nel)') - - @plansor o[-1 -2; -3 -4] := fuse_front[-1; 1 2] * o[1 3; -3 4] * nel[2 -2; 3 5] * - conj(fuse_back[-4; 4 5]) - - return calc_prod_elem(o, slice, t1[2:end], t2[2:end]) -end - -function interweave(a, b) - map(filter(x -> sum(x .== 1) == length(a) && sum(x .== 2) == length(b), - collect(Iterators.product(fill((1, 2), length(a) + length(b))...)))) do key - ia = 1 - ib = 1 - - output = Vector{eltype(a)}(undef, length(a) + length(b)) - for k in key - if k == 1 - el = a[ia] - ia += 1 - else - el = b[ib] - ib += 1 - end - output[ia + ib - 2] = el - end - return output - end -end - -function make_time_mpo(H::InfiniteMPOHamiltonian{T}, dt, alg::WII) where {T} - WA = H.A - WB = H.B - WC = H.C - WD = H.D - - δ = dt * (-1im) - Wnew = map(1:length(H)) do i - for j in 2:(size(H[i], 1) - 1), k in 2:(size(H[i], 4) - 1) - init_1 = isometry(storagetype(WD[i]), codomain(WD[i]), domain(WD[i])) - init = [init_1, - zero(H[i][1, 1, 1, k]), - zero(H[i][j, 1, 1, end]), - zero(H[i][j, 1, 1, k])] - - y, convhist = exponentiate(1.0, init, - Arnoldi(; tol=alg.tol, maxiter=alg.maxiter)) do x - out = similar(x) - - @plansor out[1][-1 -2; -3 -4] := δ * x[1][-1 1; -3 -4] * - H[i][1, 1, 1, end][2 3; 1 4] * - τ[-2 4; 2 3] - - @plansor out[2][-1 -2; -3 -4] := δ * x[2][-1 1; -3 -4] * - H[i][1, 1, 1, end][2 3; 1 4] * - τ[-2 4; 2 3] - @plansor out[2][-1 -2; -3 -4] += sqrt(δ) * - x[1][1 2; -3 4] * - H[i][1, 1, 1, k][-1 -2; 3 -4] * - τ[3 4; 1 2] - - @plansor out[3][-1 -2; -3 -4] := δ * x[3][-1 1; -3 -4] * - H[i][1, 1, 1, end][2 3; 1 4] * - τ[-2 4; 2 3] - @plansor out[3][-1 -2; -3 -4] += sqrt(δ) * - x[1][1 2; -3 4] * - H[i][j, 1, 1, end][-1 -2; 3 -4] * - τ[3 4; 1 2] - - @plansor out[4][-1 -2; -3 -4] := δ * x[4][-1 1; -3 -4] * - H[i][1, 1, 1, end][2 3; 1 4] * - τ[-2 4; 2 3] - @plansor out[4][-1 -2; -3 -4] += x[1][1 2; -3 4] * - H[i][j, 1, 1, k][-1 -2; 3 -4] * - τ[3 4; 1 2] - @plansor out[4][-1 -2; -3 -4] += sqrt(δ) * - x[2][1 2; -3 -4] * - H[i][j, 1, 1, end][-1 -2; 3 4] * - τ[3 4; 1 2] - @plansor out[4][-1 -2; -3 -4] += sqrt(δ) * - x[3][-1 4; -3 3] * - H[i][1, 1, 1, k][2 -2; 1 -4] * - τ[3 4; 1 2] - - return out - end - convhist.converged == 0 && - @warn "failed to exponentiate $(convhist.normres)" - - WA[i][j - 1, 1, 1, k - 1] = y[4] - WB[i][j - 1, 1, 1, 1] = y[3] - WC[i][1, 1, 1, k - 1] = y[2] - WD[i] = y[1] - end - - Vₗ = left_virtualspace(H, i)[1:(end - 1)] - Vᵣ = right_virtualspace(H, i)[1:(end - 1)] - P = physicalspace(H, i) - - h′ = similar(H[i], Vₗ ⊗ P ← P ⊗ Vᵣ) - h′[2:end, 1, 1, 2:end] = WA[i] - h′[2:end, 1, 1, 1] = WB[i] - h′[1, 1, 1, 2:end] = WC[i] - h′[1, 1, 1, 1] = WD[i] - - return h′ - end - - return InfiniteMPO(PeriodicArray(Wnew)) -end diff --git a/src/algorithms/timestep/wii.jl b/src/algorithms/timestep/wii.jl new file mode 100644 index 000000000..8b85f0321 --- /dev/null +++ b/src/algorithms/timestep/wii.jl @@ -0,0 +1,101 @@ +""" +$(TYPEDEF) + +Generalization of the Euler approximation of the operator exponential for MPOs. + +## Fields + +$(TYPEDFIELDS) + +## References + +* [Zaletel et al. Phys. Rev. B 91 (2015)](@cite zaletel2015) +* [Paeckel et al. Ann. of Phys. 411 (2019)](@cite paeckel2019) +""" +@kwdef struct WII <: Algorithm + "tolerance for convergence criterium" + tol::Float64 = Defaults.tol + "maximal number of iterations" + maxiter::Int = Defaults.maxiter +end + +function make_time_mpo(H::InfiniteMPOHamiltonian{T}, dt, alg::WII) where {T} + WA = H.A + WB = H.B + WC = H.C + WD = H.D + + δ = dt * (-1im) + Wnew = map(1:length(H)) do i + for j in 2:(size(H[i], 1) - 1), k in 2:(size(H[i], 4) - 1) + init_1 = isometry(storagetype(WD[i]), codomain(WD[i]), domain(WD[i])) + init = [init_1, + zero(H[i][1, 1, 1, k]), + zero(H[i][j, 1, 1, end]), + zero(H[i][j, 1, 1, k])] + + y, convhist = exponentiate(1.0, init, + Arnoldi(; tol=alg.tol, maxiter=alg.maxiter)) do x + out = similar(x) + + @plansor out[1][-1 -2; -3 -4] := δ * x[1][-1 1; -3 -4] * + H[i][1, 1, 1, end][2 3; 1 4] * + τ[-2 4; 2 3] + + @plansor out[2][-1 -2; -3 -4] := δ * x[2][-1 1; -3 -4] * + H[i][1, 1, 1, end][2 3; 1 4] * + τ[-2 4; 2 3] + @plansor out[2][-1 -2; -3 -4] += sqrt(δ) * + x[1][1 2; -3 4] * + H[i][1, 1, 1, k][-1 -2; 3 -4] * + τ[3 4; 1 2] + + @plansor out[3][-1 -2; -3 -4] := δ * x[3][-1 1; -3 -4] * + H[i][1, 1, 1, end][2 3; 1 4] * + τ[-2 4; 2 3] + @plansor out[3][-1 -2; -3 -4] += sqrt(δ) * + x[1][1 2; -3 4] * + H[i][j, 1, 1, end][-1 -2; 3 -4] * + τ[3 4; 1 2] + + @plansor out[4][-1 -2; -3 -4] := δ * x[4][-1 1; -3 -4] * + H[i][1, 1, 1, end][2 3; 1 4] * + τ[-2 4; 2 3] + @plansor out[4][-1 -2; -3 -4] += x[1][1 2; -3 4] * + H[i][j, 1, 1, k][-1 -2; 3 -4] * + τ[3 4; 1 2] + @plansor out[4][-1 -2; -3 -4] += sqrt(δ) * + x[2][1 2; -3 -4] * + H[i][j, 1, 1, end][-1 -2; 3 4] * + τ[3 4; 1 2] + @plansor out[4][-1 -2; -3 -4] += sqrt(δ) * + x[3][-1 4; -3 3] * + H[i][1, 1, 1, k][2 -2; 1 -4] * + τ[3 4; 1 2] + + return out + end + convhist.converged == 0 && + @warn "failed to exponentiate $(convhist.normres)" + + WA[i][j - 1, 1, 1, k - 1] = y[4] + WB[i][j - 1, 1, 1, 1] = y[3] + WC[i][1, 1, 1, k - 1] = y[2] + WD[i] = y[1] + end + + Vₗ = left_virtualspace(H, i)[1:(end - 1)] + Vᵣ = right_virtualspace(H, i)[1:(end - 1)] + P = physicalspace(H, i) + + h′ = similar(H[i], Vₗ ⊗ P ← P ⊗ Vᵣ) + h′[2:end, 1, 1, 2:end] = WA[i] + h′[2:end, 1, 1, 1] = WB[i] + h′[1, 1, 1, 2:end] = WC[i] + h′[1, 1, 1, 1] = WD[i] + + return h′ + end + + return InfiniteMPO(PeriodicArray(Wnew)) +end From 261ad0e7c01dc7da4b8647c2c22445450668e4a0 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Fri, 28 Feb 2025 17:33:52 -0500 Subject: [PATCH 10/13] WII for `FiniteMPOHamiltonian` --- src/algorithms/timestep/wii.jl | 118 ++++++++++++++++++++------------- 1 file changed, 73 insertions(+), 45 deletions(-) diff --git a/src/algorithms/timestep/wii.jl b/src/algorithms/timestep/wii.jl index 8b85f0321..684a72cb6 100644 --- a/src/algorithms/timestep/wii.jl +++ b/src/algorithms/timestep/wii.jl @@ -19,62 +19,61 @@ $(TYPEDFIELDS) maxiter::Int = Defaults.maxiter end -function make_time_mpo(H::InfiniteMPOHamiltonian{T}, dt, alg::WII) where {T} +function make_time_mpo(H::InfiniteMPOHamiltonian, dt::Number, alg::WII) WA = H.A WB = H.B WC = H.C WD = H.D δ = dt * (-1im) + exp_alg = Arnoldi(; alg.tol, alg.maxiter) + Wnew = map(1:length(H)) do i for j in 2:(size(H[i], 1) - 1), k in 2:(size(H[i], 4) - 1) - init_1 = isometry(storagetype(WD[i]), codomain(WD[i]), domain(WD[i])) - init = [init_1, + init = (isometry(storagetype(WD[i]), codomain(WD[i]), domain(WD[i])), zero(H[i][1, 1, 1, k]), zero(H[i][j, 1, 1, end]), - zero(H[i][j, 1, 1, k])] - - y, convhist = exponentiate(1.0, init, - Arnoldi(; tol=alg.tol, maxiter=alg.maxiter)) do x - out = similar(x) - - @plansor out[1][-1 -2; -3 -4] := δ * x[1][-1 1; -3 -4] * - H[i][1, 1, 1, end][2 3; 1 4] * - τ[-2 4; 2 3] - - @plansor out[2][-1 -2; -3 -4] := δ * x[2][-1 1; -3 -4] * - H[i][1, 1, 1, end][2 3; 1 4] * - τ[-2 4; 2 3] - @plansor out[2][-1 -2; -3 -4] += sqrt(δ) * - x[1][1 2; -3 4] * - H[i][1, 1, 1, k][-1 -2; 3 -4] * - τ[3 4; 1 2] - - @plansor out[3][-1 -2; -3 -4] := δ * x[3][-1 1; -3 -4] * - H[i][1, 1, 1, end][2 3; 1 4] * - τ[-2 4; 2 3] - @plansor out[3][-1 -2; -3 -4] += sqrt(δ) * - x[1][1 2; -3 4] * - H[i][j, 1, 1, end][-1 -2; 3 -4] * - τ[3 4; 1 2] - - @plansor out[4][-1 -2; -3 -4] := δ * x[4][-1 1; -3 -4] * - H[i][1, 1, 1, end][2 3; 1 4] * - τ[-2 4; 2 3] - @plansor out[4][-1 -2; -3 -4] += x[1][1 2; -3 4] * - H[i][j, 1, 1, k][-1 -2; 3 -4] * - τ[3 4; 1 2] - @plansor out[4][-1 -2; -3 -4] += sqrt(δ) * - x[2][1 2; -3 -4] * - H[i][j, 1, 1, end][-1 -2; 3 4] * - τ[3 4; 1 2] - @plansor out[4][-1 -2; -3 -4] += sqrt(δ) * - x[3][-1 4; -3 3] * - H[i][1, 1, 1, k][2 -2; 1 -4] * - τ[3 4; 1 2] - - return out + zero(H[i][j, 1, 1, k])) + + y, convhist = exponentiate(1.0, init, exp_alg) do (x1, x2, x3, x4) + @plansor y1[-1 -2; -3 -4] := δ * x1[-1 1; -3 -4] * + H[i][1, 1, 1, end][2 3; 1 4] * + τ[-2 4; 2 3] + + @plansor y2[-1 -2; -3 -4] := δ * x2[-1 1; -3 -4] * + H[i][1, 1, 1, end][2 3; 1 4] * + τ[-2 4; 2 3] + + sqrt(δ) * + x1[1 2; -3 4] * + H[i][1, 1, 1, k][-1 -2; 3 -4] * + τ[3 4; 1 2] + + @plansor y3[-1 -2; -3 -4] := δ * x3[-1 1; -3 -4] * + H[i][1, 1, 1, end][2 3; 1 4] * + τ[-2 4; 2 3] + + sqrt(δ) * + x1[1 2; -3 4] * + H[i][j, 1, 1, end][-1 -2; 3 -4] * + τ[3 4; 1 2] + + @plansor y4[-1 -2; -3 -4] := (δ * x4[-1 1; -3 -4] * + H[i][1, 1, 1, end][2 3; 1 4] * + τ[-2 4; 2 3] + + x1[1 2; -3 4] * + H[i][j, 1, 1, k][-1 -2; 3 -4] * + τ[3 4; 1 2]) + + (sqrt(δ) * + x2[1 2; -3 -4] * + H[i][j, 1, 1, end][-1 -2; 3 4] * + τ[3 4; 1 2]) + + (sqrt(δ) * + x3[-1 4; -3 3] * + H[i][1, 1, 1, k][2 -2; 1 -4] * + τ[3 4; 1 2]) + + return y1, y2, y3, y4 end + convhist.converged == 0 && @warn "failed to exponentiate $(convhist.normres)" @@ -99,3 +98,32 @@ function make_time_mpo(H::InfiniteMPOHamiltonian{T}, dt, alg::WII) where {T} return InfiniteMPO(PeriodicArray(Wnew)) end + +# Hack to treat FiniteMPOhamiltonians as Infinite +function make_time_mpo(H::FiniteMPOHamiltonian, dt::Number, alg::WII) + H′ = copy(parent(H)) + + V_left = left_virtualspace(H[1]) + V_left′ = BlockTensorKit.oplus(V_left, oneunit(V_left), oneunit(V_left)) + H′[1] = similar(H[1], V_left′ ⊗ space(H[1], 2) ← domain(H[1])) + for (I, v) in nonzero_pairs(H[1]) + H′[1][I] = v + end + + V_right = right_virtualspace(H[end]) + V_right′ = BlockTensorKit.oplus(oneunit(V_right), oneunit(V_right), V_right) + H′[end] = similar(H[end], codomain(H[end]) ← space(H[end], 3)' ⊗ V_right′) + for (I, v) in nonzero_pairs(H[end]) + H′[end][I[1], 1, 1, end] = v + end + + H′[1][end, 1, 1, end] = H′[1][1, 1, 1, 1] + H′[end][1, 1, 1, 1] = H′[end][end, 1, 1, end] + + mpo = make_time_mpo(InfiniteMPOHamiltonian(H′), dt, alg; tol) + + # Impose boundary conditions + mpo[1] = mpo[1][1, :, :, :] + mpo[end] = mpo[end][:, :, :, 1] + return remove_orphans!(mpo; tol) +end From 0822843540c0d36d1445795c707be0b54e0060e3 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sat, 1 Mar 2025 16:02:42 -0500 Subject: [PATCH 11/13] fix typo --- src/algorithms/timestep/taylorcluster.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/algorithms/timestep/taylorcluster.jl b/src/algorithms/timestep/taylorcluster.jl index edab759a9..f802e3758 100644 --- a/src/algorithms/timestep/taylorcluster.jl +++ b/src/algorithms/timestep/taylorcluster.jl @@ -205,5 +205,5 @@ function make_time_mpo(H::FiniteMPOHamiltonian, dt::Number, alg::TaylorCluster; # Impose boundary conditions mpo[1] = mpo[1][1, :, :, :] mpo[end] = mpo[end][:, :, :, 1] - return remove_orphans(mpo; tol) + return remove_orphans!(mpo; tol) end From 739df362752cab7dc315ba67aba8ec9644aeace4 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sun, 2 Mar 2025 07:32:23 -0500 Subject: [PATCH 12/13] reuse boundary code --- src/algorithms/timestep/taylorcluster.jl | 5 ++--- src/algorithms/timestep/wii.jl | 5 ++--- 2 files changed, 4 insertions(+), 6 deletions(-) diff --git a/src/algorithms/timestep/taylorcluster.jl b/src/algorithms/timestep/taylorcluster.jl index f802e3758..9e3b62034 100644 --- a/src/algorithms/timestep/taylorcluster.jl +++ b/src/algorithms/timestep/taylorcluster.jl @@ -203,7 +203,6 @@ function make_time_mpo(H::FiniteMPOHamiltonian, dt::Number, alg::TaylorCluster; mpo = make_time_mpo(InfiniteMPOHamiltonian(H′), dt, alg; tol) # Impose boundary conditions - mpo[1] = mpo[1][1, :, :, :] - mpo[end] = mpo[end][:, :, :, 1] - return remove_orphans!(mpo; tol) + mpo_fin = open_boundary_conditions(mpo, length(H)) + return remove_orphans!(mpo_fin; tol) end diff --git a/src/algorithms/timestep/wii.jl b/src/algorithms/timestep/wii.jl index 684a72cb6..82aec465f 100644 --- a/src/algorithms/timestep/wii.jl +++ b/src/algorithms/timestep/wii.jl @@ -123,7 +123,6 @@ function make_time_mpo(H::FiniteMPOHamiltonian, dt::Number, alg::WII) mpo = make_time_mpo(InfiniteMPOHamiltonian(H′), dt, alg; tol) # Impose boundary conditions - mpo[1] = mpo[1][1, :, :, :] - mpo[end] = mpo[end][:, :, :, 1] - return remove_orphans!(mpo; tol) + mpo_fin = open_boundary_conditions(mpo, length(H)) + return remove_orphans!(mpo_fin; tol) end From 4439d55ac966c075605baea819ea1c22501a481c Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sun, 2 Mar 2025 07:33:51 -0500 Subject: [PATCH 13/13] fix remove_orphans to avoid infinite loop --- src/operators/abstractmpo.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/operators/abstractmpo.jl b/src/operators/abstractmpo.jl index 938b26a5c..07fff6d30 100644 --- a/src/operators/abstractmpo.jl +++ b/src/operators/abstractmpo.jl @@ -72,8 +72,8 @@ function remove_orphans!(mpo::SparseMPO; tol=eps(real(scalartype(mpo)))^(3 / 4)) for i in 1:length(mpo) # slice empty columns on right or empty rows on left mask = filter(1:size(mpo[i], 4)) do j - return j ∈ getindex.(nonzero_keys(mpo[i]), 1) && - j ∈ getindex.(nonzero_keys(mpo[i + 1]), 4) + return j ∈ getindex.(nonzero_keys(mpo[i]), 4) && + j ∈ getindex.(nonzero_keys(mpo[i + 1]), 1) end changed |= length(mask) == size(mpo[i], 4) mpo[i] = mpo[i][:, :, :, mask]