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/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..9e3b62034 --- /dev/null +++ b/src/algorithms/timestep/taylorcluster.jl @@ -0,0 +1,208 @@ +""" +$(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_fin = open_boundary_conditions(mpo, length(H)) + return remove_orphans!(mpo_fin; 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 0467fdc82..000000000 --- a/src/algorithms/timestep/timeevmpo.jl +++ /dev/null @@ -1,316 +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)))^(3 / 4)) - 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) - linds = LinearIndices(ntuple(i -> V, N)) - cinds = CartesianIndices(linds) - - # 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 - 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 - end - - # 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) - 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 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) - - 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 - end - end - for I in nonzero_keys(slice) - (I[1] == c_lin || I[4] == c_lin) && delete!(slice, I) - end - end - end - end - - # Approximate compression step: Algorithm 4 - if alg.compression - for slice in parent(mpo) - for a in cinds - 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] - 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] - end - end - - for I in nonzero_keys(slice) - (I[1] == a_lin || I[4] == a_lin) && 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=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..82aec465f --- /dev/null +++ b/src/algorithms/timestep/wii.jl @@ -0,0 +1,128 @@ +""" +$(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, 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 = (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, 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)" + + 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 + +# 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_fin = open_boundary_conditions(mpo, length(H)) + return remove_orphans!(mpo_fin; tol) +end diff --git a/src/operators/abstractmpo.jl b/src/operators/abstractmpo.jl index 4236e2a02..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] @@ -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 """ 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 199fdfb69..05fd89973 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)] @@ -795,8 +875,9 @@ 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 + 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 # ----------------------------