diff --git a/Project.toml b/Project.toml index e49759e8d..3d607b569 100644 --- a/Project.toml +++ b/Project.toml @@ -22,7 +22,6 @@ RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec" TensorKitManifolds = "11fa318c-39cb-4a83-b1ed-cdc7ba1e3684" Transducers = "28d57a85-8fef-5791-bfe6-a80928e7c999" -TupleTools = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6" VectorInterface = "409d34a3-91d5-4945-b6ec-7529ddf182d8" [compat] @@ -50,7 +49,6 @@ TensorOperations = "5" Test = "1" TestExtras = "0.3" Transducers = "0.4" -TupleTools = "1.6.0" VectorInterface = "0.2, 0.3, 0.4, 0.5" julia = "1.10" diff --git a/src/MPSKit.jl b/src/MPSKit.jl index ecad4dcf3..a2b7555cb 100644 --- a/src/MPSKit.jl +++ b/src/MPSKit.jl @@ -12,7 +12,6 @@ using RecipesBase using VectorInterface using Accessors using HalfIntegers -import TupleTools as TT using DocStringExtensions using LinearAlgebra: diag, Diagonal diff --git a/src/algorithms/timestep/timeevmpo.jl b/src/algorithms/timestep/timeevmpo.jl index b7e2cde65..b3eea50a2 100644 --- a/src/algorithms/timestep/timeevmpo.jl +++ b/src/algorithms/timestep/timeevmpo.jl @@ -38,12 +38,13 @@ $(TYPEDFIELDS) "include higher-order corrections" extension::Bool = true "approximate compression of corrections, accurate up to order `N`" - compression::Bool = true + compression::Bool = false end const WI = TaylorCluster(; N=1, extension=false, compression=false) -function make_time_mpo(H::MPOHamiltonian, dt::Number, alg::TaylorCluster) +function make_time_mpo(H::MPOHamiltonian, dt::Number, alg::TaylorCluster; + tol=eps(real(scalartype(H)))^(3 / 4)) N = alg.N τ = -1im * dt @@ -80,20 +81,23 @@ function make_time_mpo(H::MPOHamiltonian, dt::Number, alg::TaylorCluster) if alg.extension H_next = H_n * H′ linds_next = LinearIndices(ntuple(i -> V, N + 1)) - cinds_next = CartesianIndices(linds_next) for (i, slice) in enumerate(parent(H_n)) - for Inext in nonzero_keys(H_next[i]) - aₑ = cinds_next[Inext[1]] - bₑ = cinds_next[Inext[4]] - for c in findall(==(1), aₑ.I), d in findall(==(V), bₑ.I) - a = TT.deleteat(Tuple(aₑ), c) - b = TT.deleteat(Tuple(bₑ), d) - if all(>(1), b) && !(all(in((1, V), a)) && any(==(V), a)) - n1 = count(==(1), a) + 1 - n3 = count(==(V), b) + 1 - factor = τ * factorial(N) / (factorial(N + 1) * n1 * n3) - slice[linds[a...], 1, 1, linds[b...]] += factor * H_next[i][Inext] - end + for a in cinds, b in cinds + all(>(1), b.I) || continue + all(in((1, V)), a.I) && any(==(V), a.I) && continue + + n1 = count(==(1), a.I) + 1 + n3 = count(==(V), 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) + + # 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ₑ...]] end end end @@ -101,52 +105,41 @@ function make_time_mpo(H::MPOHamiltonian, dt::Number, alg::TaylorCluster) # loopback step: Algorithm 1 # constructing the Nth order time evolution MPO - if H isa FiniteMPOHamiltonian - mpo = FiniteMPO(parent(H_n)) - else - mpo = InfiniteMPO(parent(H_n)) - end + mpo = MPO(parent(H_n)) for slice in parent(mpo) - for (I, v) in nonzero_pairs(slice) - a = cinds[I[1]] - b = cinds[I[4]] - if all(in((1, V)), a.I) && any(!=(1), a.I) - delete!(slice, I) - elseif any(!=(1), b.I) && all(in((1, V)), b.I) - exponent = count(==(V), b.I) - factor = τ^exponent * factorial(N - exponent) / factorial(N) - Idst = CartesianIndex(Base.setindex(I.I, 1, 4)) - slice[Idst] += factor * v - delete!(slice, I) + for b in cinds[2:end] + all(in((1, V)), b.I) || continue + + b_lin = linds[b] + a = count(==(V), 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 slice in parent(mpo) - for I in nonzero_keys(slice) - a = cinds[I[1]] - a_sort = CartesianIndex(sort(collect(a.I); by=(!=(1)))...) - n1_a = count(==(1), a.I) - n3_a = count(==(V), a.I) - if n3_a <= n1_a && a_sort != a - Idst = CartesianIndex(Base.setindex(I.I, linds[a_sort], 1)) - slice[Idst] += slice[I] - delete!(slice, I) - elseif a != a_sort - delete!(slice, I) - end + for c in cinds + c_lin = linds[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) - b = cinds[I[4]] - b_sort = CartesianIndex(sort(collect(b.I); by=(!=(V)))...) - n1_b = count(==(1), b.I) - n3_b = count(==(V), b.I) - if n3_b > n1_b && b_sort != b - Idst = CartesianIndex(Base.setindex(I.I, linds[b_sort], 4)) - slice[Idst] += slice[I] - delete!(slice, I) - elseif b != b_sort - delete!(slice, I) + if n3 <= n1 && s_c != c + slice[linds[s_c], 1, 1, :] += slice[c_lin, 1, 1, :] + for I in nonzero_keys(slice) + (I[1] == c_lin || I[4] == c_lin) && delete!(slice, I) + end + elseif n3 > n1 && s_r != c + slice[:, 1, 1, linds[s_r]] += slice[:, 1, 1, c_lin] + for I in nonzero_keys(slice) + (I[1] == c_lin || I[4] == c_lin) && delete!(slice, I) + end end end end @@ -161,14 +154,10 @@ function make_time_mpo(H::MPOHamiltonian, dt::Number, alg::TaylorCluster) b = CartesianIndex(replace(a.I, V => 1)) b_lin = linds[b] factor = τ^n1 * factorial(N - n1) / factorial(N) + slice[:, 1, 1, b_lin] += factor * slice[:, 1, 1, a_lin] + for I in nonzero_keys(slice) - if I[4] == a_lin - Idst = CartesianIndex(Base.setindex(I.I, b_lin, 4)) - slice[Idst] += factor * slice[I] - delete!(slice, I) - elseif I[1] == a_lin - delete!(slice, I) - end + (I[1] == a_lin || I[4] == a_lin) && delete!(slice, I) end end end @@ -180,7 +169,7 @@ function make_time_mpo(H::MPOHamiltonian, dt::Number, alg::TaylorCluster) mpo[end] = mpo[end][:, :, :, 1] end - return remove_orphans!(mpo) + return remove_orphans!(mpo; tol=tol) end has_prod_elem(slice, t1, t2) = all(map(x -> contains(slice, x...), zip(t1, t2))) @@ -223,65 +212,6 @@ function interweave(a, b) end end -# function make_time_mpo(H::MPOHamiltonian{S,T}, dt, alg::WII) where {S,T} -# WA = PeriodicArray{T,3}(undef, H.period, H.odim - 2, H.odim - 2) -# WB = PeriodicArray{T,2}(undef, H.period, H.odim - 2) -# WC = PeriodicArray{T,2}(undef, H.period, H.odim - 2) -# WD = PeriodicArray{T,1}(undef, H.period) -# -# δ = dt * (-1im) -# -# for i in 1:(H.period), j in 2:(H.odim - 1), k in 2:(H.odim - 1) -# init_1 = isometry(storagetype(H[i][1, H.odim]), codomain(H[i][1, H.odim]), -# domain(H[i][1, H.odim])) -# init = [init_1, zero(H[i][1, k]), zero(H[i][j, H.odim]), zero(H[i][j, k])] -# -# (y, convhist) = exponentiate(1.0, RecursiveVec(init), -# Arnoldi(; tol=alg.tol, maxiter=alg.maxiter)) do x -# out = similar(x.vecs) -# -# @plansor out[1][-1 -2; -3 -4] := δ * x[1][-1 1; -3 -4] * -# H[i][1, H.odim][2 3; 1 4] * τ[-2 4; 2 3] -# -# @plansor out[2][-1 -2; -3 -4] := δ * x[2][-1 1; -3 -4] * -# H[i][1, H.odim][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, 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, H.odim][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, H.odim][-1 -2; 3 -4] * τ[3 4; 1 2] -# -# @plansor out[4][-1 -2; -3 -4] := δ * x[4][-1 1; -3 -4] * -# H[i][1, H.odim][2 3; 1 4] * τ[-2 4; 2 3] -# @plansor out[4][-1 -2; -3 -4] += x[1][1 2; -3 4] * H[i][j, 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, H.odim][-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, k][2 -2; 1 -4] * τ[3 4; 1 2] -# -# return RecursiveVec(out) -# end -# convhist.converged == 0 && -# @warn "exponentiate failed to converge: normres = $(convhist.normres)" -# -# WA[i, j - 1, k - 1] = y[4] -# WB[i, j - 1] = y[3] -# WC[i, k - 1] = y[2] -# WD[i] = y[1] -# end -# -# W2 = PeriodicArray{Union{T,Missing},3}(missing, H.period, H.odim - 1, H.odim - 1) -# W2[:, 2:end, 2:end] = WA -# W2[:, 2:end, 1] = WB -# W2[:, 1, 2:end] = WC -# W2[:, 1, 1] = WD -# -# return SparseMPO(W2) -# end - function make_time_mpo(H::InfiniteMPOHamiltonian{T}, dt, alg::WII) where {T} WA = H.A WB = H.B diff --git a/test/algorithms.jl b/test/algorithms.jl index 97a924475..fa4dc7e2d 100644 --- a/test/algorithms.jl +++ b/test/algorithms.jl @@ -769,4 +769,41 @@ end end end +@testset "TaylorCluster time evolution" begin + L = 4 + Hs = [transverse_field_ising(; L=L), heisenberg_XXX(; L=L)] + + Ns = [1, 2, 3] + dts = [1e-2, 1e-3] + for H in Hs + ψ = FiniteMPS(L, physicalspace(H, 1), ℂ^16) + for N in Ns + εs = zeros(ComplexF64, 2) + for (i, dt) in enumerate(dts) + ψ₀, _ = find_groundstate(ψ, H, DMRG(; verbosity=0)) + E₀ = expectation_value(ψ₀, H) + + O = make_time_mpo(H, dt, TaylorCluster(; N=N)) + + ψ₁, _ = 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 + 1 atol = 0.1 + end + + for N in Ns + εs = zeros(ComplexF64, 2) + for (i, dt) in enumerate(dts) + ψ₀, _ = find_groundstate(ψ, H, DMRG(; verbosity=0)) + E₀ = expectation_value(ψ₀, H) + + O = make_time_mpo(H, dt, TaylorCluster(; N=N, compression=true)) + + ψ₁, _ = 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 + end + end +end end