From 1c3b05979284f5296b80b7464ad59b1d55fa750c Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 17 Apr 2025 09:25:51 -0400 Subject: [PATCH 01/27] Add missing transfer operations --- src/transfermatrix/transfer.jl | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/src/transfermatrix/transfer.jl b/src/transfermatrix/transfer.jl index 15cc2237a..27c4eb92c 100644 --- a/src/transfermatrix/transfer.jl +++ b/src/transfermatrix/transfer.jl @@ -81,6 +81,19 @@ function transfer_right(v::MPSTensor{S}, A::MPSTensor{S}, Ab::MPSTensor{S}) wher check_unambiguous_braiding(space(v, 2)) return @plansor v[-1 -2; -3] := A[-1 2; 1] * τ[-2 4; 2 3] * conj(Ab[-3 4; 5]) * v[1 3; 5] end +# similar for Matrix Product Density Operators +function transfer_left(v::MPSTensor{S}, A::GenericMPSTensor{S,3}, + Ab::GenericMPSTensor{S,3}) where {S} + check_unambiguous_braiding(space(v, 2)) + @plansor v[-1 -2; -3] ≔ v[1 3; 8] * A[8 7 6; -3] * τ[3 2; 7 5] * τ[5 4; 6 -2] * + conj(Ab[1 2 4; -1]) +end +function transfer_right(v::MPSTensor{S}, A::GenericMPSTensor{S,3}, + Ab::GenericMPSTensor{S,3}) where {S} + check_unambiguous_braiding(space(v, 2)) + @plansor v[-1 -2; -3] ≔ A[-1 4 2; 1] * τ[-2 6; 4 5] * τ[5 7; 2 3] * + conj(Ab[-3 6 7; 8]) * v[1 3; 8] +end # the transfer operation with a utility leg in both the domain and codomain is also ill defined - only due to the codomain utility space function transfer_left(v::MPOTensor{S}, A::MPSTensor{S}, Ab::MPSTensor{S}) where {S} From 453f9b30d3fcfffd32cc6613489cbf5c5864f35f Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 17 Apr 2025 09:32:08 -0400 Subject: [PATCH 02/27] Add fallback `contract_mpo_expval` definition --- src/algorithms/expval.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/algorithms/expval.jl b/src/algorithms/expval.jl index e8ac68738..50e41c689 100644 --- a/src/algorithms/expval.jl +++ b/src/algorithms/expval.jl @@ -93,6 +93,8 @@ function contract_mpo_expval( ) return @plansor GL[1 2; 3] * AC[3 7; 5] * GR[5 8; 6] * O[2 4; 7 8] * conj(ACbar[1 4; 6]) end +# generic fallback +contract_mpo_expval(AC, GL, O, GR, ACbar=AC) = dot(ACbar, ∂AC(AC, GL, O, GR)) function expectation_value( ψ::FiniteMPS, H::FiniteMPOHamiltonian, From 8d899b48e3e42448e19b9e40357d9d189a922a84 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 18 Sep 2025 11:52:39 +0200 Subject: [PATCH 03/27] fix merge --- src/algorithms/expval.jl | 101 ++++++++++++++++++-------------- src/operators/mpohamiltonian.jl | 11 +++- src/states/finitemps.jl | 10 ++-- 3 files changed, 72 insertions(+), 50 deletions(-) diff --git a/src/algorithms/expval.jl b/src/algorithms/expval.jl index 50e41c689..7fe318c54 100644 --- a/src/algorithms/expval.jl +++ b/src/algorithms/expval.jl @@ -37,53 +37,43 @@ function expectation_value end function expectation_value(ψ::AbstractMPS, (inds, O)::Pair) @boundscheck foreach(Base.Fix1(Base.checkbounds, ψ), inds) - sites, local_mpo = instantiate_operator(physicalspace(ψ), inds => O) - @assert _firstspace(first(local_mpo)) == oneunit(_firstspace(first(local_mpo))) == - dual(_lastspace(last(local_mpo))) - for (site, o) in zip(sites, local_mpo) - if o isa MPOTensor - physicalspace(ψ, site) == physicalspace(o) || - throw(SpaceMismatch("physical space does not match at site $site")) - end + # special cases that can be handled more efficiently + if length(inds) == 1 + return local_expectation_value1(ψ, inds[1], O) + elseif length(inds) == 2 && eltype(inds) <: Integer && inds[1] + 1 == inds[2] + return local_expectation_value2(ψ, inds[1], O) end - Ut = fill_data!(similar(local_mpo[1], _firstspace(first(local_mpo))), one) - - # some special cases that avoid using transfer matrices - if length(sites) == 1 - AC = ψ.AC[sites[1]] - E = @plansor conj(AC[4 5; 6]) * conj(Ut[1]) * local_mpo[1][1 5; 3 2] * Ut[2] * - AC[4 3; 6] - elseif length(sites) == 2 && (sites[1] + 1 == sites[2]) - AC = ψ.AC[sites[1]] - AR = ψ.AR[sites[2]] - O1, O2 = local_mpo - E = @plansor conj(AC[4 5; 10]) * conj(Ut[1]) * O1[1 5; 3 8] * AC[4 3; 6] * - conj(AR[10 9; 11]) * Ut[2] * O2[8 9; 7 2] * AR[6 7; 11] - else - # generic case: write as Vl * T^N * Vr - # left side - T = storagetype(site_type(ψ)) - @plansor Vl[-1 -2; -3] := isomorphism( - T, left_virtualspace(ψ, sites[1]), left_virtualspace(ψ, sites[1]) - )[-1; -3] * - conj(Ut[-2]) - - # middle - M = foldl(zip(sites, local_mpo); init = Vl) do v, (site, o) - if o isa Number - return scale!(v * TransferMatrix(ψ.AL[site], ψ.AL[site]), o) - else - return v * TransferMatrix(ψ.AL[site], o, ψ.AL[site]) - end - end + # generic case: convert to MPO and write as Vl * T^N * Vr + sites, local_mpo = instantiate_operator(ψ, inds => O) - # right side - E = @plansor M[1 2; 3] * Ut[2] * ψ.C[sites[end]][3; 4] * - conj(ψ.C[sites[end]][1; 4]) + # left side + Vl = insertrightunit(l_LL(ψ, sites[1]), 1; dual=true) + + # middle + M = foldl(zip(sites, local_mpo); init=Vl) do v, (site, o) + if o isa Number + return scale!(v * TransferMatrix(ψ.AL[site], ψ.AL[site]), o) + else + return v * TransferMatrix(ψ.AL[site], o, ψ.AL[site]) + end end - return E / norm(ψ)^2 + # right side + E = @plansor removeunit(M, 2)[1; 2] * ψ.C[sites[end]][2; 3] * + conj(ψ.C[sites[end]][1; 3]) + return E / dot(ψ, ψ) +end + +function local_expectation_value1(ψ::AbstractMPS, site, O) + E = contract_mpo_expval1(ψ.AC[site], O, ψ.AC[site]) + return E / dot(ψ, ψ) +end +function local_expectation_value2(ψ::AbstractMPS, site, O) + AC = ψ.AC[site] + AR = ψ.AR[site + 1] + E = contract_mpo_expval2(AC, AR, O, AC, AR) + return E / dot(ψ, ψ) end # MPOHamiltonian @@ -94,7 +84,32 @@ function contract_mpo_expval( return @plansor GL[1 2; 3] * AC[3 7; 5] * GR[5 8; 6] * O[2 4; 7 8] * conj(ACbar[1 4; 6]) end # generic fallback -contract_mpo_expval(AC, GL, O, GR, ACbar=AC) = dot(ACbar, ∂AC(AC, GL, O, GR)) +contract_mpo_expval(AC, GL, O, GR, ACbar=AC) = dot(ACbar, ∂AC(AC, O, GL, GR)) + +function contract_mpo_expval1(AC::MPSTensor, O::AbstractTensorMap, ACbar::MPSTensor=AC) + numin(O) == numout(O) == 1 || throw(ArgumentError("O is not a single-site operator")) + return @plansor conj(ACbar[2 3; 4]) * O[3; 1] * AC[2 1; 4] +end +function contract_mpo_expval1(AC::GenericMPSTensor{S,3}, O::AbstractTensorMap{<:Any,S}, + ACbar::GenericMPSTensor{S,3}=AC) where {S} + numin(O) == numout(O) == 1 || throw(ArgumentError("O is not a single-site operator")) + return @plansor conj(ACbar[2 3 4; 5]) * O[3; 1] * AC[2 1 4; 5] +end + +function contract_mpo_expval2(A1::MPSTensor, A2::MPSTensor, O::AbstractTensorMap, + A1bar::MPSTensor=A1, A2bar::MPSTensor=A2) + numin(O) == numout(O) == 2 || throw(ArgumentError("O is not a two-site operator")) + return @plansor conj(A1bar[4 5; 6]) * conj(A2bar[6 7; 8]) * O[5 7; 2 3] * A1[4 2; 1] * + A2[1 3; 8] +end +function contract_mpo_expval2(A1::GenericMPSTensor{S,3}, A2::GenericMPSTensor{S,3}, + O::AbstractTensorMap{<:Any,S}, + A1bar::GenericMPSTensor{S,3}=A1, + A2bar::GenericMPSTensor{S,3}=A2) where {S} + numin(O) == numout(O) == 2 || throw(ArgumentError("O is not a two-site operator")) + return @plansor conj(A1bar[8 3 4; 11]) * conj(A2bar[11 12 13; 14]) * τ[9 6; 1 2] * + τ[3 4; 9 10] * A1[8 1 2; 5] * A2[5 7 13; 14] * O[10 12; 6 7] +end function expectation_value( ψ::FiniteMPS, H::FiniteMPOHamiltonian, diff --git a/src/operators/mpohamiltonian.jl b/src/operators/mpohamiltonian.jl index 0727f4ca3..2f03a4cda 100644 --- a/src/operators/mpohamiltonian.jl +++ b/src/operators/mpohamiltonian.jl @@ -299,11 +299,15 @@ function _split_mpoham_types(W::Matrix)::Type{<:MPOTensor} end """ + instantiate_operator(state, O::Pair) instantiate_operator(lattice::AbstractArray{<:VectorSpace}, O::Pair) -Instantiate a local operator `O` on a lattice `lattice` as a vector of MPO tensors, and a -vector of linear site indices. +Instantiate a local operator `O` for a `state` or `lattice` as a vector of MPO tensors, and +a vector of linear site indices. """ +function instantiate_operator(state::AbstractMPS, O::Pair) + return instantiate_operator(physicalspace(state), O) +end function instantiate_operator(lattice::AbstractArray{<:VectorSpace}, (inds′, O)::Pair) inds = inds′ isa Int ? tuple(inds′) : inds′ mpo = O isa FiniteMPO ? O : FiniteMPO(O) @@ -320,7 +324,8 @@ function instantiate_operator(lattice::AbstractArray{<:VectorSpace}, (inds′, O i = 1 for j in first(indices):last(indices) if j == indices[i] - @assert space(operators[i], 2) == lattice[j] "operator does not fit into the given Hilbert space: $(space(operators[i], 2)) ≠ $(lattice[j])" + physicalspace(operators[i]) == lattice[j] || + throw(SpaceMismatch("physical space does not match at site $j")) push!(local_mpo, operators[i]) i += 1 else diff --git a/src/states/finitemps.jl b/src/states/finitemps.jl index 1300ceaa4..0a32f965e 100644 --- a/src/states/finitemps.jl +++ b/src/states/finitemps.jl @@ -626,9 +626,11 @@ TensorKit.normalize(ψ::FiniteMPS) = normalize!(copy(ψ)) Fixedpoints ===========================================================================================# -function r_RR(ψ::FiniteMPS{T}) where {T} - return isomorphism(storagetype(T), domain(ψ.AR[end]), domain(ψ.AR[end])) +function r_RR(ψ::FiniteMPS, site::Int=length(ψ)) + Vr = right_virtualspace(ψ.AR[site]) + return isomorphism(storagetype(site_type(ψ)), Vr ← Vr) end -function l_LL(ψ::FiniteMPS{T}) where {T} - return isomorphism(storagetype(T), space(ψ.AL[1], 1), space(ψ.AL[1], 1)) +function l_LL(ψ::FiniteMPS, site::Int=1) + Vl = left_virtualspace(ψ.AL[site]) + return isomorphism(storagetype(site_type(ψ)), Vl ← Vl) end From 5d5de540768f3a2c94cb2b6bd924add7d5037736 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 17 Apr 2025 17:53:17 -0400 Subject: [PATCH 04/27] Add support for imaginary time evolution with real tensors --- src/algorithms/timestep/integrators.jl | 8 ++- src/algorithms/timestep/taylorcluster.jl | 9 +-- src/algorithms/timestep/tdvp.jl | 70 +++++++++++++++++------- src/algorithms/timestep/time_evolve.jl | 37 +++++++++++-- src/algorithms/timestep/wii.jl | 10 ++-- 5 files changed, 98 insertions(+), 36 deletions(-) diff --git a/src/algorithms/timestep/integrators.jl b/src/algorithms/timestep/integrators.jl index 7d19591d8..b324cf1a9 100644 --- a/src/algorithms/timestep/integrators.jl +++ b/src/algorithms/timestep/integrators.jl @@ -17,8 +17,12 @@ function integrate end _eval_t(f, t::Number) = Base.Fix2(f, t) _eval_x(f, x) = Base.Fix1(f, x) -function integrate(f, y₀, t::Number, dt::Number, alg::Union{Arnoldi, Lanczos}) - y, convhist = exponentiate(_eval_t(f, t + dt / 2), -1im * dt, y₀, alg) +function integrate( + f, y₀, t::Number, dt::Number, alg::Union{Arnoldi, Lanczos}; + imaginary_evolution::Bool = false + ) + dt′ = imaginary_evolution ? -dt : -1im * dt + y, convhist = exponentiate(_eval_t(f, t + dt′ / 2), dt′, y₀, alg) convhist.converged == 0 && @warn "integrator failed to converge: normres = $(convhist.normres)" return y diff --git a/src/algorithms/timestep/taylorcluster.jl b/src/algorithms/timestep/taylorcluster.jl index cbad1ee13..977189f89 100644 --- a/src/algorithms/timestep/taylorcluster.jl +++ b/src/algorithms/timestep/taylorcluster.jl @@ -28,10 +28,11 @@ 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))) + H::MPOHamiltonian, dt::Number, alg::TaylorCluster; + tol = eps(real(scalartype(H))), imaginary_evolution::Bool = false ) N = alg.N - τ = -1im * dt + τ = imaginary_evolution ? -dt : -1im * dt # start with H^N H_n = H^N @@ -180,7 +181,7 @@ end # Hack to treat FiniteMPOhamiltonians as Infinite function make_time_mpo( H::FiniteMPOHamiltonian, dt::Number, alg::TaylorCluster; - tol = eps(real(scalartype(H))) + tol = eps(real(scalartype(H))); imaginary_evolution::Bool = false ) H′ = copy(parent(H)) @@ -201,7 +202,7 @@ function make_time_mpo( 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) + mpo = make_time_mpo(InfiniteMPOHamiltonian(H′), dt, alg; tol, imaginary_evolution) # Impose boundary conditions mpo_fin = open_boundary_conditions(mpo, length(H)) diff --git a/src/algorithms/timestep/tdvp.jl b/src/algorithms/timestep/tdvp.jl index b711d022c..dc335e2e1 100644 --- a/src/algorithms/timestep/tdvp.jl +++ b/src/algorithms/timestep/tdvp.jl @@ -26,11 +26,15 @@ $(TYPEDFIELDS) end function timestep( - ψ_::InfiniteMPS, H, t::Number, dt::Number, alg::TDVP, + ψ::InfiniteMPS, H, t::Number, dt::Number, alg::TDVP, envs::AbstractMPSEnvironments = environments(ψ_, H); - leftorthflag = true + leftorthflag = true, imaginary_evolution::Bool = false ) - ψ = complex(ψ_) + # convert state to complex if necessary + if scalartype(ψ) <: Real && (!imaginary_evolution || !isreal(dt)) + return timestep(complex(ψ), H, t, dt, alg, envs; leftorthflag, imaginary_evolution) + end + temp_ACs = similar(ψ.AC) temp_Cs = similar(ψ.C) @@ -38,24 +42,30 @@ function timestep( if scheduler isa SerialScheduler temp_ACs = tmap!(temp_ACs, 1:length(ψ); scheduler) do loc Hac = AC_hamiltonian(loc, ψ, H, ψ, envs) - return integrate(Hac, ψ.AC[loc], t, dt, alg.integrator) + return integrate(Hac, ψ.AC[loc], t, dt, alg.integrator; imaginary_evolution) end temp_Cs = tmap!(temp_Cs, 1:length(ψ); scheduler) do loc Hc = C_hamiltonian(loc, ψ, H, ψ, envs) - return integrate(Hc, ψ.C[loc], t, dt, alg.integrator) + return integrate(Hc, ψ.C[loc], t, dt, alg.integrator; imaginary_evolution) end else @sync begin Threads.@spawn begin temp_ACs = tmap!(temp_ACs, 1:length(ψ); scheduler) do loc Hac = AC_hamiltonian(loc, ψ, H, ψ, envs) - return integrate(Hac, ψ.AC[loc], t, dt, alg.integrator) + return integrate( + Hac, ψ.AC[loc], t, dt, alg.integrator; + imaginary_evolution + ) end end Threads.@spawn begin temp_Cs = tmap!(temp_Cs, 1:length(ψ); scheduler) do loc Hc = C_hamiltonian(loc, ψ, H, ψ, envs) - return integrate(Hc, ψ.C[loc], t, dt, alg.integrator) + return integrate( + Hc, ψ.C[loc], t, dt, alg.integrator; + imaginary_evolution + ) end end end @@ -76,34 +86,47 @@ end function timestep!( ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP, - envs::AbstractMPSEnvironments = environments(ψ, H) + envs::AbstractMPSEnvironments = environments(ψ, H); + imaginary_evolution::Bool = false ) # sweep left to right for i in 1:(length(ψ) - 1) Hac = AC_hamiltonian(i, ψ, H, ψ, envs) - ψ.AC[i] = integrate(Hac, ψ.AC[i], t, dt / 2, alg.integrator) + ψ.AC[i] = integrate(Hac, ψ.AC[i], t, dt / 2, alg.integrator; imaginary_evolution) Hc = C_hamiltonian(i, ψ, H, ψ, envs) - ψ.C[i] = integrate(Hc, ψ.C[i], t + dt / 2, -dt / 2, alg.integrator) + ψ.C[i] = integrate( + Hc, ψ.C[i], t + dt / 2, -dt / 2, alg.integrator; + imaginary_evolution + ) end # edge case Hac = AC_hamiltonian(length(ψ), ψ, H, ψ, envs) - ψ.AC[end] = integrate(Hac, ψ.AC[end], t, dt / 2, alg.integrator) + ψ.AC[end] = integrate(Hac, ψ.AC[end], t, dt / 2, alg.integrator; imaginary_evolution) # sweep right to left for i in length(ψ):-1:2 Hac = AC_hamiltonian(i, ψ, H, ψ, envs) - ψ.AC[i] = integrate(Hac, ψ.AC[i], t + dt / 2, dt / 2, alg.integrator) + ψ.AC[i] = integrate( + Hac, ψ.AC[i], t + dt / 2, dt / 2, alg.integrator; + imaginary_evolution + ) Hc = C_hamiltonian(i - 1, ψ, H, ψ, envs) - ψ.C[i - 1] = integrate(Hc, ψ.C[i - 1], t + dt, -dt / 2, alg.integrator) + ψ.C[i - 1] = integrate( + Hc, ψ.C[i - 1], t + dt, -dt / 2, alg.integrator; + imaginary_evolution + ) end # edge case Hac = AC_hamiltonian(1, ψ, H, ψ, envs) - ψ.AC[1] = integrate(Hac, ψ.AC[1], t + dt / 2, dt / 2, alg.integrator) + ψ.AC[1] = integrate( + Hac, ψ.AC[1], t + dt / 2, dt / 2, alg.integrator; + imaginary_evolution + ) return ψ, envs end @@ -143,14 +166,15 @@ end function timestep!( ψ::AbstractFiniteMPS, H, t::Number, dt::Number, alg::TDVP2, - envs::AbstractMPSEnvironments = environments(ψ, H) + envs::AbstractMPSEnvironments = environments(ψ, H); + imaginary_evolution::Bool = false ) # sweep left to right for i in 1:(length(ψ) - 1) ac2 = _transpose_front(ψ.AC[i]) * _transpose_tail(ψ.AR[i + 1]) Hac2 = AC2_hamiltonian(i, ψ, H, ψ, envs) - ac2′ = integrate(Hac2, ac2, t, dt / 2, alg.integrator) + ac2′ = integrate(Hac2, ac2, t, dt / 2, alg.integrator; imaginary_evolution) nal, nc, nar = tsvd!(ac2′; trunc = alg.trscheme, alg = alg.alg_svd) ψ.AC[i] = (nal, complex(nc)) @@ -158,7 +182,10 @@ function timestep!( if i != (length(ψ) - 1) Hac = AC_hamiltonian(i + 1, ψ, H, ψ, envs) - ψ.AC[i + 1] = integrate(Hac, ψ.AC[i + 1], t + dt / 2, -dt / 2, alg.integrator) + ψ.AC[i + 1] = integrate( + Hac, ψ.AC[i + 1], t + dt / 2, -dt / 2, alg.integrator; + imaginary_evolution + ) end end @@ -166,7 +193,7 @@ function timestep!( for i in length(ψ):-1:2 ac2 = _transpose_front(ψ.AL[i - 1]) * _transpose_tail(ψ.AC[i]) Hac2 = AC2_hamiltonian(i - 1, ψ, H, ψ, envs) - ac2′ = integrate(Hac2, ac2, t + dt / 2, dt / 2, alg.integrator) + ac2′ = integrate(Hac2, ac2, t + dt / 2, dt / 2, alg.integrator; imaginary_evolution) nal, nc, nar = tsvd!(ac2′; trunc = alg.trscheme, alg = alg.alg_svd) ψ.AC[i - 1] = (nal, complex(nc)) @@ -174,14 +201,17 @@ function timestep!( if i != 2 Hac = AC_hamiltonian(i - 1, ψ, H, ψ, envs) - ψ.AC[i - 1] = integrate(Hac, ψ.AC[i - 1], t + dt, -dt / 2, alg.integrator) + ψ.AC[i - 1] = integrate( + Hac, ψ.AC[i - 1], t + dt, -dt / 2, alg.integrator; + imaginary_evolution + ) end end return ψ, envs end -#copying version +# copying version function timestep( ψ::AbstractFiniteMPS, H, time::Number, timestep::Number, alg::Union{TDVP, TDVP2}, envs::AbstractMPSEnvironments...; diff --git a/src/algorithms/timestep/time_evolve.jl b/src/algorithms/timestep/time_evolve.jl index a9b06eeda..a0918ad9f 100644 --- a/src/algorithms/timestep/time_evolve.jl +++ b/src/algorithms/timestep/time_evolve.jl @@ -5,19 +5,29 @@ Time-evolve the initial state `ψ₀` with Hamiltonian `H` over a given time span by stepping through each of the time points obtained by iterating t_span. -# Arguments +## Arguments + - `ψ₀::AbstractMPS`: initial state - `H::AbstractMPO`: operator that generates the time evolution (can be time-dependent). - `t_span::AbstractVector{<:Number}`: time points over which the time evolution is stepped - `[alg]`: algorithm to use for the time evolution. Defaults to [`TDVP`](@ref). - `[envs]`: MPS environment manager + +## Keyword Arguments + +- `verbosity::Int=0`: verbosity level for logging +- `imaginary_evolution::Bool=false`: if true, the time evolution is done with an imaginary time step + instead, (i.e. ``\\exp(-Hdt)`` instead of ``\\exp(-iHdt)``). This can be useful for using this + function to compute the ground state of a Hamiltonian, or to compute finite-temperature + properties of a system. """ function time_evolve end, function time_evolve! end for (timestep, time_evolve) in zip((:timestep, :timestep!), (:time_evolve, :time_evolve!)) @eval function $time_evolve( - ψ, H, t_span::AbstractVector{<:Number}, alg, envs = environments(ψ, H); - verbosity::Int = 0 + ψ, H, t_span::AbstractVector{<:Number}, alg, + envs = environments(ψ, H); + verbosity::Int = 0, imaginary_evolution::Bool = false ) log = IterLog("TDVP") LoggingExtras.withlevel(; verbosity) do @@ -26,7 +36,7 @@ for (timestep, time_evolve) in zip((:timestep, :timestep!), (:time_evolve, :time t = t_span[iter] dt = t_span[iter + 1] - t - ψ, envs = $timestep(ψ, H, t, dt, alg, envs) + ψ, envs = $timestep(ψ, H, t, dt, alg, envs; imaginary_evolution) ψ, envs = alg.finalize(t, ψ, H, envs)::Tuple{typeof(ψ), typeof(envs)} @infov 3 logiter!(log, iter, 0, t) @@ -44,18 +54,33 @@ end Time-step the state `ψ₀` with Hamiltonian `H` over a given time step `dt` at time `t`, solving the Schroedinger equation: ``i ∂ψ/∂t = H ψ``. -# Arguments +## Arguments + - `ψ₀::AbstractMPS`: initial state - `H::AbstractMPO`: operator that generates the time evolution (can be time-dependent). - `t::Number`: starting time of time-step - `dt::Number`: time-step magnitude - `[alg]`: algorithm to use for the time evolution. Defaults to [`TDVP`](@ref). - `[envs]`: MPS environment manager + +## Keyword Arguments + +- `imaginary_evolution::Bool=false`: if true, the time evolution is done with an imaginary time step + instead, (i.e. ``\\exp(-Hdt)`` instead of ``\\exp(-iHdt)``). This can be useful for using this + function to compute the ground state of a Hamiltonian, or to compute finite-temperature + properties of a system. """ function timestep end, function timestep! end @doc """ - make_time_mpo(H::MPOHamiltonian, dt::Number, alg) -> O::MPO + make_time_mpo(H::MPOHamiltonian, dt::Number, alg; kwargs...) -> O::MPO Construct an `MPO` that approximates ``\\exp(-iHdt)``. + +## Keyword Arguments + +- `imaginary_evolution::Bool=false`: if true, the time evolution operator is constructed + with an imaginary time step instead, (i.e. ``\\exp(-Hdt)`` instead of ``\\exp(-iHdt)``). + This can be useful for using this function to compute the ground state of a Hamiltonian, + or to compute finite-temperature properties of a system. """ make_time_mpo diff --git a/src/algorithms/timestep/wii.jl b/src/algorithms/timestep/wii.jl index cd02352d5..874a47d18 100644 --- a/src/algorithms/timestep/wii.jl +++ b/src/algorithms/timestep/wii.jl @@ -19,13 +19,14 @@ $(TYPEDFIELDS) maxiter::Int = Defaults.maxiter end -function make_time_mpo(H::InfiniteMPOHamiltonian, dt::Number, alg::WII) +function make_time_mpo(H::InfiniteMPOHamiltonian, dt::Number, alg::WII; + imaginary_evolution::Bool=false) WA = H.A WB = H.B WC = H.C WD = H.D - δ = dt * (-1im) + δ = imaginary_evolution ? -dt : (-1im * dt) exp_alg = Arnoldi(; alg.tol, alg.maxiter) Wnew = map(1:length(H)) do i @@ -83,7 +84,8 @@ function make_time_mpo(H::InfiniteMPOHamiltonian, dt::Number, alg::WII) end # Hack to treat FiniteMPOhamiltonians as Infinite -function make_time_mpo(H::FiniteMPOHamiltonian, dt::Number, alg::WII) +function make_time_mpo(H::FiniteMPOHamiltonian, dt::Number, alg::WII; + imaginary_evolution::Bool=false) H′ = copy(parent(H)) V_left = left_virtualspace(H[1]) @@ -103,7 +105,7 @@ function make_time_mpo(H::FiniteMPOHamiltonian, dt::Number, alg::WII) 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) + mpo = make_time_mpo(InfiniteMPOHamiltonian(H′), dt, alg; tol, imaginary_evolution) # Impose boundary conditions mpo_fin = open_boundary_conditions(mpo, length(H)) From 9b02e721d8b739e6d7d310ac9937b7397ca5f0e7 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 17 Apr 2025 17:53:27 -0400 Subject: [PATCH 05/27] Fix `complex(::MPOHamiltonian)` --- src/operators/jordanmpotensor.jl | 20 ++++++++++++++++++++ src/operators/mpohamiltonian.jl | 15 ++------------- 2 files changed, 22 insertions(+), 13 deletions(-) diff --git a/src/operators/jordanmpotensor.jl b/src/operators/jordanmpotensor.jl index 7fcc7c6e8..bb7932414 100644 --- a/src/operators/jordanmpotensor.jl +++ b/src/operators/jordanmpotensor.jl @@ -194,6 +194,26 @@ function SparseBlockTensorMap(W::JordanMPOTensor) return W′ end +for f in (:real, :complex) + @eval function Base.$f(W::JordanMPOTensor) + E = $f(scalartype(W)) + W′ = JordanMPOTensor{E}(undef, space(W)) + for (I, v) in nonzero_pairs(W.A) + W′.A[I] = $f(v) + end + for (I, v) in nonzero_pairs(W.B) + W′.B[I] = $f(v) + end + for (I, v) in nonzero_pairs(W.C) + W′.C[I] = $f(v) + end + for (I, v) in nonzero_pairs(W.D) + W′.D[I] = $f(v) + end + return W′ + end +end + # Indexing # -------- diff --git a/src/operators/mpohamiltonian.jl b/src/operators/mpohamiltonian.jl index 2f03a4cda..e32df8c16 100644 --- a/src/operators/mpohamiltonian.jl +++ b/src/operators/mpohamiltonian.jl @@ -634,19 +634,8 @@ end # Base.complex(H::MPOHamiltonian) = MPOHamiltonian(map(complex, parent(H))) function Base.complex(H::MPOHamiltonian) scalartype(H) <: Complex && return H - Ws = map(parent(H)) do W - W′ = jordanmpotensortype(spacetype(W), complex(scalartype(W))) - W′[1] = W[1] - W′[end] = W[end] - for (I, v) in nonzero_pairs(W) - if v isa BraidingTensor - W′[I] = BraidingTensor{scalartype(W′)}(space(v), v.adjoint) - else - W′[I] = complex(v) - end - end - end - return MPOHamiltonian(H) + Ws = map(complex, parent(H)) + return MPOHamiltonian(Ws) end function Base.similar(H::MPOHamiltonian, ::Type{O}, L::Int) where {O <: MPOTensor} From a372067fa27c5cf4c8408a645de3e1a771154aca Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 17 Apr 2025 19:20:44 -0400 Subject: [PATCH 06/27] 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 977189f89..d9cd70247 100644 --- a/src/algorithms/timestep/taylorcluster.jl +++ b/src/algorithms/timestep/taylorcluster.jl @@ -181,7 +181,7 @@ end # Hack to treat FiniteMPOhamiltonians as Infinite function make_time_mpo( H::FiniteMPOHamiltonian, dt::Number, alg::TaylorCluster; - tol = eps(real(scalartype(H))); imaginary_evolution::Bool = false + tol = eps(real(scalartype(H))), imaginary_evolution::Bool = false ) H′ = copy(parent(H)) From 8639e084be51a1c5ec31d1faf89f558e40b67daf Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Fri, 2 May 2025 11:03:38 -0400 Subject: [PATCH 07/27] fix typo --- src/algorithms/timestep/tdvp.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/algorithms/timestep/tdvp.jl b/src/algorithms/timestep/tdvp.jl index dc335e2e1..fae5bb3e0 100644 --- a/src/algorithms/timestep/tdvp.jl +++ b/src/algorithms/timestep/tdvp.jl @@ -27,7 +27,7 @@ end function timestep( ψ::InfiniteMPS, H, t::Number, dt::Number, alg::TDVP, - envs::AbstractMPSEnvironments = environments(ψ_, H); + envs::AbstractMPSEnvironments = environments(ψ, H); leftorthflag = true, imaginary_evolution::Bool = false ) # convert state to complex if necessary From 5ca3248edc5ebc3d731db9d0d230412c0207852d Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Mon, 5 May 2025 11:03:46 -0400 Subject: [PATCH 08/27] Refactor IDMRG to be compatible with density matrix mps From c6b577f205d580123aa335db3042a97dc9a2c949 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Tue, 24 Jun 2025 10:36:26 -0400 Subject: [PATCH 09/27] Fix fallback with missing function --- src/algorithms/expval.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/algorithms/expval.jl b/src/algorithms/expval.jl index 7fe318c54..4b6a66ab5 100644 --- a/src/algorithms/expval.jl +++ b/src/algorithms/expval.jl @@ -84,7 +84,7 @@ function contract_mpo_expval( return @plansor GL[1 2; 3] * AC[3 7; 5] * GR[5 8; 6] * O[2 4; 7 8] * conj(ACbar[1 4; 6]) end # generic fallback -contract_mpo_expval(AC, GL, O, GR, ACbar=AC) = dot(ACbar, ∂AC(AC, O, GL, GR)) +contract_mpo_expval(AC, GL, O, GR, ACbar=AC) = dot(ACbar, MPO_AC_Hamiltonian(GL, O, GR) * AC) function contract_mpo_expval1(AC::MPSTensor, O::AbstractTensorMap, ACbar::MPSTensor=AC) numin(O) == numout(O) == 1 || throw(ArgumentError("O is not a single-site operator")) From 656578eb862e183dd3d88d1341f2978b64d6536d Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Tue, 24 Jun 2025 10:37:09 -0400 Subject: [PATCH 10/27] Real scalartype improvements --- src/algorithms/changebonds/svdcut.jl | 8 ++++---- src/algorithms/timestep/tdvp.jl | 6 +++--- src/states/orthoview.jl | 12 ++++++++++-- 3 files changed, 17 insertions(+), 9 deletions(-) diff --git a/src/algorithms/changebonds/svdcut.jl b/src/algorithms/changebonds/svdcut.jl index 7ad26e706..095e260f6 100644 --- a/src/algorithms/changebonds/svdcut.jl +++ b/src/algorithms/changebonds/svdcut.jl @@ -29,9 +29,9 @@ function changebonds!(ψ::AbstractFiniteMPS, alg::SvdCut; normalize::Bool = true for i in (length(ψ) - 1):-1:1 U, S, V, = tsvd(ψ.C[i]; trunc = alg.trscheme, alg = alg.alg_svd) AL′ = ψ.AL[i] * U - ψ.AC[i] = (AL′, complex(S)) + ψ.AC[i] = (AL′, S) AR′ = _transpose_front(V * _transpose_tail(ψ.AR[i + 1])) - ψ.AC[i + 1] = (complex(S), AR′) + ψ.AC[i + 1] = (S, AR′) end return normalize ? normalize!(ψ) : ψ end @@ -87,7 +87,7 @@ function changebonds(ψ::MultilineMPS, alg::SvdCut) return Multiline(map(x -> changebonds(x, alg), ψ.data)) end function changebonds(ψ::InfiniteMPS, alg::SvdCut) - copied = complex.(ψ.AL) + copied = copy.(ψ.AL) ncr = ψ.C[1] for i in 1:length(ψ) @@ -103,7 +103,7 @@ function changebonds(ψ::InfiniteMPS, alg::SvdCut) ψ = if space(ncr, 1) != space(copied[1], 1) InfiniteMPS(copied) else - C₀ = ncr isa TensorMap ? complex(ncr) : TensorMap(complex(ncr)) + C₀ = ncr isa TensorMap ? ncr : TensorMap(ncr) InfiniteMPS(copied, C₀) end return normalize!(ψ) diff --git a/src/algorithms/timestep/tdvp.jl b/src/algorithms/timestep/tdvp.jl index fae5bb3e0..e372c4fbd 100644 --- a/src/algorithms/timestep/tdvp.jl +++ b/src/algorithms/timestep/tdvp.jl @@ -215,9 +215,9 @@ end function timestep( ψ::AbstractFiniteMPS, H, time::Number, timestep::Number, alg::Union{TDVP, TDVP2}, envs::AbstractMPSEnvironments...; - kwargs... + imaginary_evolution::Bool = false, kwargs... ) - isreal = scalartype(ψ) <: Real + isreal = (scalartype(ψ) <: Real && !imaginary_evolution) ψ′ = isreal ? complex(ψ) : copy(ψ) if length(envs) != 0 && isreal @warn "Currently cannot reuse real environments for complex evolution" @@ -228,5 +228,5 @@ function timestep( @assert length(envs) == 0 "Invalid signature" envs′ = environments(ψ′, H) end - return timestep!(ψ′, H, time, timestep, alg, envs′; kwargs...) + return timestep!(ψ′, H, time, timestep, alg, envs′; imaginary_evolution, kwargs...) end diff --git a/src/states/orthoview.jl b/src/states/orthoview.jl index c3abe3325..fd810e95d 100644 --- a/src/states/orthoview.jl +++ b/src/states/orthoview.jl @@ -160,10 +160,18 @@ function Base.setindex!( a, b = vec return if isa(a, MPSBondTensor) #c/ar - setindex!(v.parent.Cs, a, i) + if !(scalartype(parent(v)) <: Real) && (scalartype(a) <: Real) + setindex!(v.parent.Cs, complex(a), i) + else + setindex!(v.parent.Cs, a, i) + end setindex!(v.parent.ARs, b, i) elseif isa(b, MPSBondTensor) #al/c - setindex!(v.parent.Cs, b, i + 1) + if !(scalartype(parent(v)) <: Real) && (scalartype(b) <: Real) + setindex!(v.parent.Cs, complex(b), i + 1) + else + setindex!(v.parent.Cs, b, i + 1) + end setindex!(v.parent.ALs, a, i) else throw(ArgumentError("invalid value types")) From 9a8653d591cbab485c6993ab69f22884cc791a5b Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Fri, 27 Jun 2025 15:40:10 +0200 Subject: [PATCH 11/27] Disable spacecheck for finite T --- src/operators/mpohamiltonian.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/operators/mpohamiltonian.jl b/src/operators/mpohamiltonian.jl index e32df8c16..b760aff4f 100644 --- a/src/operators/mpohamiltonian.jl +++ b/src/operators/mpohamiltonian.jl @@ -324,8 +324,11 @@ function instantiate_operator(lattice::AbstractArray{<:VectorSpace}, (inds′, O i = 1 for j in first(indices):last(indices) if j == indices[i] - physicalspace(operators[i]) == lattice[j] || + # TODO: fix this check for density matrices + if !(eltype(lattice) isa ProductSpace) && + physicalspace(operators[i]) == lattice[j] throw(SpaceMismatch("physical space does not match at site $j")) + end push!(local_mpo, operators[i]) i += 1 else From 1537c89bff4d99f506a1e927b1df2ff0abc08792 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 3 Jul 2025 12:55:06 -0400 Subject: [PATCH 12/27] Fix spacecheck on operator --- src/operators/mpohamiltonian.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/operators/mpohamiltonian.jl b/src/operators/mpohamiltonian.jl index b760aff4f..ab66a75d2 100644 --- a/src/operators/mpohamiltonian.jl +++ b/src/operators/mpohamiltonian.jl @@ -326,7 +326,7 @@ function instantiate_operator(lattice::AbstractArray{<:VectorSpace}, (inds′, O if j == indices[i] # TODO: fix this check for density matrices if !(eltype(lattice) isa ProductSpace) && - physicalspace(operators[i]) == lattice[j] + physicalspace(operators[i]) != lattice[j] throw(SpaceMismatch("physical space does not match at site $j")) end push!(local_mpo, operators[i]) From b33278bd88f600bad3fcefbff85dfd38a68f50cd Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 9 Jul 2025 14:50:34 -0400 Subject: [PATCH 13/27] Add missing return value --- src/algorithms/approximate/approximate.jl | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/src/algorithms/approximate/approximate.jl b/src/algorithms/approximate/approximate.jl index f905c5405..3c564cde2 100644 --- a/src/algorithms/approximate/approximate.jl +++ b/src/algorithms/approximate/approximate.jl @@ -32,14 +32,13 @@ function approximate( envs = environments(ψ, toapprox) ) envs′ = Multiline([envs]) - multi, envs = approximate( + multi, envs, δ = approximate( convert(MultilineMPS, ψ), (convert(MultilineMPO, toapprox[1]), convert(MultilineMPS, toapprox[2])), - algorithm, - envs′ + algorithm, envs′ ) ψ = convert(InfiniteMPS, multi) - return ψ, envs + return ψ, envs, δ end # dispatch to in-place method @@ -56,12 +55,11 @@ function approximate( algorithm::Union{IDMRG, IDMRG2}, envs = environments(ψ, toapprox) ) envs′ = Multiline([envs]) - multi, envs = approximate( + multi, envs, δ = approximate( convert(MultilineMPS, ψ), (convert(MultilineMPO, toapprox[1]), convert(MultilineMPS, toapprox[2])), - algorithm, - envs′ + algorithm, envs′ ) ψ = convert(InfiniteMPS, multi) - return ψ, envs + return ψ, envs, δ end From 95ad216442aa5ad274e925fb3d328deed2158598 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 9 Jul 2025 14:50:48 -0400 Subject: [PATCH 14/27] Fix wrong derivative --- src/algorithms/expval.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/algorithms/expval.jl b/src/algorithms/expval.jl index 4b6a66ab5..d74ad1a5a 100644 --- a/src/algorithms/expval.jl +++ b/src/algorithms/expval.jl @@ -84,7 +84,9 @@ function contract_mpo_expval( return @plansor GL[1 2; 3] * AC[3 7; 5] * GR[5 8; 6] * O[2 4; 7 8] * conj(ACbar[1 4; 6]) end # generic fallback -contract_mpo_expval(AC, GL, O, GR, ACbar=AC) = dot(ACbar, MPO_AC_Hamiltonian(GL, O, GR) * AC) +function contract_mpo_expval(AC, GL, O, GR, ACbar=AC) + return dot(ACbar, MPO_AC_Hamiltonian(GL, O, GR) * AC) +end function contract_mpo_expval1(AC::MPSTensor, O::AbstractTensorMap, ACbar::MPSTensor=AC) numin(O) == numout(O) == 1 || throw(ArgumentError("O is not a single-site operator")) From 982e4dd152773bb6c8e0fa181b5fa6e85917f0a5 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 9 Jul 2025 14:51:53 -0400 Subject: [PATCH 15/27] fix `dot(::InfiniteMPS, ::InfiniteMPO, ::InfiniteMPS)` --- src/operators/mpo.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/operators/mpo.jl b/src/operators/mpo.jl index 8bc220f58..07b84ea03 100644 --- a/src/operators/mpo.jl +++ b/src/operators/mpo.jl @@ -319,8 +319,8 @@ function TensorKit.dot( randomize!(ρ₀) val, = fixedpoint( - TransferMatrix(ket.AL, parent(mpo), bra.AL), ρ₀, :LM; ishermitian, - krylovdim, kwargs... + TransferMatrix(ket.AL, parent(mpo), bra.AL), ρ₀, :LM; + ishermitian, krylovdim, kwargs... ) return val end From 24c9c11e00e0b9eb64aabf5cc194a1496534461f Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 9 Jul 2025 14:52:06 -0400 Subject: [PATCH 16/27] loosen type restriction --- src/operators/mpo.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/operators/mpo.jl b/src/operators/mpo.jl index 07b84ea03..5934aefda 100644 --- a/src/operators/mpo.jl +++ b/src/operators/mpo.jl @@ -325,7 +325,7 @@ function TensorKit.dot( return val end -function TensorKit.dot(mpo₁::FiniteMPO{TO}, mpo₂::FiniteMPO{TO}) where {TO <: MPOTensor} +function TensorKit.dot(mpo₁::FiniteMPO{<:MPOTensor}, mpo₂::FiniteMPO{<:MPOTensor}) length(mpo₁) == length(mpo₂) || throw(ArgumentError("dimension mismatch")) N = length(mpo₁) Nhalf = N ÷ 2 From 0a934f40408d6522577c63ec0e97e11dd05d5ca4 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 9 Jul 2025 14:54:43 -0400 Subject: [PATCH 17/27] fix type check --- src/operators/mpohamiltonian.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/operators/mpohamiltonian.jl b/src/operators/mpohamiltonian.jl index ab66a75d2..fd2972e4f 100644 --- a/src/operators/mpohamiltonian.jl +++ b/src/operators/mpohamiltonian.jl @@ -325,7 +325,7 @@ function instantiate_operator(lattice::AbstractArray{<:VectorSpace}, (inds′, O for j in first(indices):last(indices) if j == indices[i] # TODO: fix this check for density matrices - if !(eltype(lattice) isa ProductSpace) && + if !(eltype(lattice) <: ProductSpace) && physicalspace(operators[i]) != lattice[j] throw(SpaceMismatch("physical space does not match at site $j")) end From e0fda0d2deedc1d9430cc30d4c86adc4d53c05f5 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 18 Sep 2025 11:45:58 +0200 Subject: [PATCH 18/27] fix possibly-in-place-copy --- src/states/infinitemps.jl | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/src/states/infinitemps.jl b/src/states/infinitemps.jl index b92bd9c40..802602f91 100644 --- a/src/states/infinitemps.jl +++ b/src/states/infinitemps.jl @@ -237,14 +237,19 @@ Base.size(ψ::InfiniteMPS, args...) = size(ψ.AL, args...) Base.length(ψ::InfiniteMPS) = length(ψ.AL) Base.eltype(ψ::InfiniteMPS) = eltype(typeof(ψ)) Base.eltype(::Type{<:InfiniteMPS{A}}) where {A} = A + Base.copy(ψ::InfiniteMPS) = InfiniteMPS(copy(ψ.AL), copy(ψ.AR), copy(ψ.C), copy(ψ.AC)) function Base.copy!(ψ::InfiniteMPS, ϕ::InfiniteMPS) - copy!.(ψ.AL, ϕ.AL) - copy!.(ψ.AR, ϕ.AR) - copy!.(ψ.AC, ϕ.AC) - copy!.(ψ.C, ϕ.C) + ψ.AL .= _copy!!.(ψ.AL, ϕ.AL) + ψ.AR .= _copy!!.(ψ.AR, ϕ.AR) + ψ.AC .= _copy!!.(ψ.AC, ϕ.AC) + ψ.C .= _copy!!.(ψ.C, ϕ.C) return ψ end +# possible in-place copy +function _copy!!(dst::AbstractTensorMap, src::AbstractTensorMap) + return space(dst) == space(src) ? copy!(dst, src) : copy(src) +end function Base.complex(ψ::InfiniteMPS) scalartype(ψ) <: Complex && return ψ From cc8635b0a808b675202a3fe50d379c35fa51249c Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 18 Sep 2025 12:17:36 +0200 Subject: [PATCH 19/27] formatter --- src/algorithms/expval.jl | 36 ++++++++++++++++++--------------- src/algorithms/timestep/wii.jl | 12 +++++++---- src/operators/mpohamiltonian.jl | 3 +-- src/states/finitemps.jl | 4 ++-- src/transfermatrix/transfer.jl | 18 +++++++++-------- 5 files changed, 41 insertions(+), 32 deletions(-) diff --git a/src/algorithms/expval.jl b/src/algorithms/expval.jl index d74ad1a5a..1806f5a41 100644 --- a/src/algorithms/expval.jl +++ b/src/algorithms/expval.jl @@ -48,10 +48,10 @@ function expectation_value(ψ::AbstractMPS, (inds, O)::Pair) sites, local_mpo = instantiate_operator(ψ, inds => O) # left side - Vl = insertrightunit(l_LL(ψ, sites[1]), 1; dual=true) + Vl = insertrightunit(l_LL(ψ, sites[1]), 1; dual = true) # middle - M = foldl(zip(sites, local_mpo); init=Vl) do v, (site, o) + M = foldl(zip(sites, local_mpo); init = Vl) do v, (site, o) if o isa Number return scale!(v * TransferMatrix(ψ.AL[site], ψ.AL[site]), o) else @@ -60,8 +60,7 @@ function expectation_value(ψ::AbstractMPS, (inds, O)::Pair) end # right side - E = @plansor removeunit(M, 2)[1; 2] * ψ.C[sites[end]][2; 3] * - conj(ψ.C[sites[end]][1; 3]) + E = @plansor removeunit(M, 2)[1; 2] * ψ.C[sites[end]][2; 3] * conj(ψ.C[sites[end]][1; 3]) return E / dot(ψ, ψ) end @@ -84,33 +83,38 @@ function contract_mpo_expval( return @plansor GL[1 2; 3] * AC[3 7; 5] * GR[5 8; 6] * O[2 4; 7 8] * conj(ACbar[1 4; 6]) end # generic fallback -function contract_mpo_expval(AC, GL, O, GR, ACbar=AC) +function contract_mpo_expval(AC, GL, O, GR, ACbar = AC) return dot(ACbar, MPO_AC_Hamiltonian(GL, O, GR) * AC) end -function contract_mpo_expval1(AC::MPSTensor, O::AbstractTensorMap, ACbar::MPSTensor=AC) +function contract_mpo_expval1(AC::MPSTensor, O::AbstractTensorMap, ACbar::MPSTensor = AC) numin(O) == numout(O) == 1 || throw(ArgumentError("O is not a single-site operator")) return @plansor conj(ACbar[2 3; 4]) * O[3; 1] * AC[2 1; 4] end -function contract_mpo_expval1(AC::GenericMPSTensor{S,3}, O::AbstractTensorMap{<:Any,S}, - ACbar::GenericMPSTensor{S,3}=AC) where {S} +function contract_mpo_expval1( + AC::GenericMPSTensor{S, 3}, O::AbstractTensorMap{<:Any, S}, + ACbar::GenericMPSTensor{S, 3} = AC + ) where {S} numin(O) == numout(O) == 1 || throw(ArgumentError("O is not a single-site operator")) return @plansor conj(ACbar[2 3 4; 5]) * O[3; 1] * AC[2 1 4; 5] end -function contract_mpo_expval2(A1::MPSTensor, A2::MPSTensor, O::AbstractTensorMap, - A1bar::MPSTensor=A1, A2bar::MPSTensor=A2) +function contract_mpo_expval2( + A1::MPSTensor, A2::MPSTensor, O::AbstractTensorMap, + A1bar::MPSTensor = A1, A2bar::MPSTensor = A2 + ) numin(O) == numout(O) == 2 || throw(ArgumentError("O is not a two-site operator")) return @plansor conj(A1bar[4 5; 6]) * conj(A2bar[6 7; 8]) * O[5 7; 2 3] * A1[4 2; 1] * - A2[1 3; 8] + A2[1 3; 8] end -function contract_mpo_expval2(A1::GenericMPSTensor{S,3}, A2::GenericMPSTensor{S,3}, - O::AbstractTensorMap{<:Any,S}, - A1bar::GenericMPSTensor{S,3}=A1, - A2bar::GenericMPSTensor{S,3}=A2) where {S} +function contract_mpo_expval2( + A1::GenericMPSTensor{S, 3}, A2::GenericMPSTensor{S, 3}, + O::AbstractTensorMap{<:Any, S}, + A1bar::GenericMPSTensor{S, 3} = A1, A2bar::GenericMPSTensor{S, 3} = A2 + ) where {S} numin(O) == numout(O) == 2 || throw(ArgumentError("O is not a two-site operator")) return @plansor conj(A1bar[8 3 4; 11]) * conj(A2bar[11 12 13; 14]) * τ[9 6; 1 2] * - τ[3 4; 9 10] * A1[8 1 2; 5] * A2[5 7 13; 14] * O[10 12; 6 7] + τ[3 4; 9 10] * A1[8 1 2; 5] * A2[5 7 13; 14] * O[10 12; 6 7] end function expectation_value( diff --git a/src/algorithms/timestep/wii.jl b/src/algorithms/timestep/wii.jl index 874a47d18..5f5f06e84 100644 --- a/src/algorithms/timestep/wii.jl +++ b/src/algorithms/timestep/wii.jl @@ -19,8 +19,10 @@ $(TYPEDFIELDS) maxiter::Int = Defaults.maxiter end -function make_time_mpo(H::InfiniteMPOHamiltonian, dt::Number, alg::WII; - imaginary_evolution::Bool=false) +function make_time_mpo( + H::InfiniteMPOHamiltonian, dt::Number, alg::WII; + imaginary_evolution::Bool = false + ) WA = H.A WB = H.B WC = H.C @@ -84,8 +86,10 @@ function make_time_mpo(H::InfiniteMPOHamiltonian, dt::Number, alg::WII; end # Hack to treat FiniteMPOhamiltonians as Infinite -function make_time_mpo(H::FiniteMPOHamiltonian, dt::Number, alg::WII; - imaginary_evolution::Bool=false) +function make_time_mpo( + H::FiniteMPOHamiltonian, dt::Number, alg::WII; + imaginary_evolution::Bool = false + ) H′ = copy(parent(H)) V_left = left_virtualspace(H[1]) diff --git a/src/operators/mpohamiltonian.jl b/src/operators/mpohamiltonian.jl index fd2972e4f..476983658 100644 --- a/src/operators/mpohamiltonian.jl +++ b/src/operators/mpohamiltonian.jl @@ -325,8 +325,7 @@ function instantiate_operator(lattice::AbstractArray{<:VectorSpace}, (inds′, O for j in first(indices):last(indices) if j == indices[i] # TODO: fix this check for density matrices - if !(eltype(lattice) <: ProductSpace) && - physicalspace(operators[i]) != lattice[j] + if !(eltype(lattice) <: ProductSpace) && physicalspace(operators[i]) != lattice[j] throw(SpaceMismatch("physical space does not match at site $j")) end push!(local_mpo, operators[i]) diff --git a/src/states/finitemps.jl b/src/states/finitemps.jl index 0a32f965e..cc877603e 100644 --- a/src/states/finitemps.jl +++ b/src/states/finitemps.jl @@ -626,11 +626,11 @@ TensorKit.normalize(ψ::FiniteMPS) = normalize!(copy(ψ)) Fixedpoints ===========================================================================================# -function r_RR(ψ::FiniteMPS, site::Int=length(ψ)) +function r_RR(ψ::FiniteMPS, site::Int = length(ψ)) Vr = right_virtualspace(ψ.AR[site]) return isomorphism(storagetype(site_type(ψ)), Vr ← Vr) end -function l_LL(ψ::FiniteMPS, site::Int=1) +function l_LL(ψ::FiniteMPS, site::Int = 1) Vl = left_virtualspace(ψ.AL[site]) return isomorphism(storagetype(site_type(ψ)), Vl ← Vl) end diff --git a/src/transfermatrix/transfer.jl b/src/transfermatrix/transfer.jl index 27c4eb92c..aa3d10ed7 100644 --- a/src/transfermatrix/transfer.jl +++ b/src/transfermatrix/transfer.jl @@ -82,17 +82,19 @@ function transfer_right(v::MPSTensor{S}, A::MPSTensor{S}, Ab::MPSTensor{S}) wher return @plansor v[-1 -2; -3] := A[-1 2; 1] * τ[-2 4; 2 3] * conj(Ab[-3 4; 5]) * v[1 3; 5] end # similar for Matrix Product Density Operators -function transfer_left(v::MPSTensor{S}, A::GenericMPSTensor{S,3}, - Ab::GenericMPSTensor{S,3}) where {S} +function transfer_left( + v::MPSTensor{S}, A::GenericMPSTensor{S, 3}, Ab::GenericMPSTensor{S, 3} + ) where {S} check_unambiguous_braiding(space(v, 2)) - @plansor v[-1 -2; -3] ≔ v[1 3; 8] * A[8 7 6; -3] * τ[3 2; 7 5] * τ[5 4; 6 -2] * - conj(Ab[1 2 4; -1]) + return @plansor v[-1 -2; -3] ≔ v[1 3; 8] * A[8 7 6; -3] * τ[3 2; 7 5] * τ[5 4; 6 -2] * + conj(Ab[1 2 4; -1]) end -function transfer_right(v::MPSTensor{S}, A::GenericMPSTensor{S,3}, - Ab::GenericMPSTensor{S,3}) where {S} +function transfer_right( + v::MPSTensor{S}, A::GenericMPSTensor{S, 3}, Ab::GenericMPSTensor{S, 3} + ) where {S} check_unambiguous_braiding(space(v, 2)) - @plansor v[-1 -2; -3] ≔ A[-1 4 2; 1] * τ[-2 6; 4 5] * τ[5 7; 2 3] * - conj(Ab[-3 6 7; 8]) * v[1 3; 8] + return @plansor v[-1 -2; -3] ≔ A[-1 4 2; 1] * τ[-2 6; 4 5] * τ[5 7; 2 3] * + conj(Ab[-3 6 7; 8]) * v[1 3; 8] end # the transfer operation with a utility leg in both the domain and codomain is also ill defined - only due to the codomain utility space From b2c4dafb129f3e3a2ba95a56c03935691851a509 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 18 Sep 2025 13:36:17 +0200 Subject: [PATCH 20/27] fix timedependent integrator --- src/algorithms/timestep/integrators.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/algorithms/timestep/integrators.jl b/src/algorithms/timestep/integrators.jl index b324cf1a9..7b2983c34 100644 --- a/src/algorithms/timestep/integrators.jl +++ b/src/algorithms/timestep/integrators.jl @@ -22,7 +22,7 @@ function integrate( imaginary_evolution::Bool = false ) dt′ = imaginary_evolution ? -dt : -1im * dt - y, convhist = exponentiate(_eval_t(f, t + dt′ / 2), dt′, y₀, alg) + y, convhist = exponentiate(_eval_t(f, t + dt / 2), dt′, y₀, alg) convhist.converged == 0 && @warn "integrator failed to converge: normres = $(convhist.normres)" return y From d039b4ebf3aa59098510d8df4d1d289c1ff4351b Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 18 Sep 2025 13:36:37 +0200 Subject: [PATCH 21/27] fix scalartype determination in InfiniteMPS constructor --- src/states/infinitemps.jl | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/states/infinitemps.jl b/src/states/infinitemps.jl index 802602f91..c9a4cf906 100644 --- a/src/states/infinitemps.jl +++ b/src/states/infinitemps.jl @@ -193,8 +193,10 @@ function InfiniteMPS(AL::AbstractVector{<:GenericMPSTensor}, C₀::MPSBondTensor # initialize tensor storage AC = similar.(AL) AR = similar.(AL) - C = similar(AR, typeof(C₀)) - ψ = InfiniteMPS{eltype(AL), eltype(C)}(AL, AR, C, AC) + T = TensorOperations.promote_contract(scalartype(AL), scalartype(C₀)) + TC = TensorOperations.tensoradd_type(T, C₀, ((1,), (2,)), false) + C = similar(AR, TC) + ψ = InfiniteMPS{eltype(AL), TC}(AL, AR, C, AC) # gaugefix the MPS gaugefix!(ψ, AL, C₀; order = :R, kwargs...) @@ -209,8 +211,10 @@ function InfiniteMPS(C₀::MPSBondTensor, AR::AbstractVector{<:GenericMPSTensor} # initialize tensor storage AC = similar.(AR) AL = similar.(AR) - C = similar(AR, typeof(C₀)) - ψ = InfiniteMPS{eltype(AL), eltype(C)}(AL, AR, C, AC) + T = TensorOperations.promote_contract(eltype(AR), eltype(C₀)) + TC = TensorOperations.tensoradd_type(T, C₀, ((1,), (2,)), false) + C = similar(AR, TC) + ψ = InfiniteMPS{eltype(AL), TC}(AL, AR, C, AC) # gaugefix the MPS gaugefix!(ψ, AR, C₀; order = :L, kwargs...) From f8d213f0fb7faeb92754815331d6260e77b0c69b Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sun, 21 Sep 2025 19:15:06 +0200 Subject: [PATCH 22/27] improve test coverage --- src/algorithms/expval.jl | 9 +++++++++ test/algorithms.jl | 12 ++++++++++-- 2 files changed, 19 insertions(+), 2 deletions(-) diff --git a/src/algorithms/expval.jl b/src/algorithms/expval.jl index 1806f5a41..dfce5bc35 100644 --- a/src/algorithms/expval.jl +++ b/src/algorithms/expval.jl @@ -200,3 +200,12 @@ function expectation_value( n = norm(ψ.AC[end])^2 return sum(ens) / (n * length(ψ)) end + +# Density matrices +# ---------------- +function expectation_value(ρ::FiniteMPO, args...) + return expectation_value(convert(FiniteMPS, ρ), args...) +end +function expectation_value(ρ::InfiniteMPO, args...) + return expectation_value(convert(InfiniteMPS, ρ), args...) +end diff --git a/test/algorithms.jl b/test/algorithms.jl index b009c121b..d309ca486 100644 --- a/test/algorithms.jl +++ b/test/algorithms.jl @@ -962,12 +962,15 @@ module TestAlgorithms Z_dense_2 = tr(exp(-2beta * H_dense))^(1 / L) # taylor cluster - rho_taylor_1 = make_time_mpo(H, -im * beta, TaylorCluster(; N = 2)) + rho_taylor_1 = make_time_mpo(H, beta, TaylorCluster(; N = 2); imaginary_evolution = true) Z_taylor_1 = tr(rho_taylor_1)^(1 / L) @test Z_taylor_1 ≈ Z_dense_1 atol = 1.0e-2 Z_taylor_2 = real(dot(rho_taylor_1, rho_taylor_1))^(1 / L) @test Z_taylor_2 ≈ Z_dense_2 atol = 1.0e-2 + E_x_taylor = @constinferred expectation_value(rho_taylor_1, 1 => S_x()) + E_xx_taylor = @constinferred expectation_value(rho_taylor_1, (1, 2) => S_xx()) + # MPO multiplication rho_mps = convert(FiniteMPS, rho_taylor_1) rho_mps, = approximate(rho_mps, (rho_taylor_1, rho_mps), DMRG2(; trscheme, verbosity)) @@ -977,9 +980,14 @@ module TestAlgorithms # TDVP rho_0 = MPSKit.infinite_temperature_density_matrix(H) rho_0_mps = convert(FiniteMPS, rho_0) - rho_mps, = timestep(rho_0_mps, H, 0.0, im * beta, TDVP2(; trscheme)) + rho_mps, = timestep(rho_0_mps, H, 0.0, beta, TDVP2(; trscheme); imaginary_evolution = true) Z_tdvp = real(dot(rho_mps, rho_mps))^(1 / L) @test Z_tdvp ≈ Z_dense_2 atol = 1.0e-2 + + @test expectation_value(rho_0_mps, 1 => S_x()) ≈ 0 + @test expectation_value(rho_0_mps, (1, 2) => S_xx()) ≈ 0 + @test expectation_value(rho_mps, 1 => S_x()) ≈ E_x_taylor atol = 1.0e-2 + @test expectation_value(rho_mps, (1, 2) => S_xx()) ≈ E_xx_taylor atol = 1.0e-2 end @testset "Sector conventions" begin From 5558e498f3da9553f83196e53a475cb119ffe0b9 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sun, 21 Sep 2025 19:29:51 +0200 Subject: [PATCH 23/27] improve WII coverage and fix --- src/algorithms/timestep/wii.jl | 7 ++++--- test/algorithms.jl | 7 +++++++ 2 files changed, 11 insertions(+), 3 deletions(-) diff --git a/src/algorithms/timestep/wii.jl b/src/algorithms/timestep/wii.jl index 5f5f06e84..3cf65cad6 100644 --- a/src/algorithms/timestep/wii.jl +++ b/src/algorithms/timestep/wii.jl @@ -28,7 +28,8 @@ function make_time_mpo( WC = H.C WD = H.D - δ = imaginary_evolution ? -dt : (-1im * dt) + # needs to be complex if negative because of square root + δ = imaginary_evolution ? -complex(dt) : (-1im * dt) exp_alg = Arnoldi(; alg.tol, alg.maxiter) Wnew = map(1:length(H)) do i @@ -109,9 +110,9 @@ function make_time_mpo( 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, imaginary_evolution) + mpo = make_time_mpo(InfiniteMPOHamiltonian(H′), dt, alg; imaginary_evolution) # Impose boundary conditions mpo_fin = open_boundary_conditions(mpo, length(H)) - return remove_orphans!(mpo_fin; tol) + return remove_orphans!(mpo_fin; alg.tol) end diff --git a/test/algorithms.jl b/test/algorithms.jl index d309ca486..fdfdd9e99 100644 --- a/test/algorithms.jl +++ b/test/algorithms.jl @@ -971,6 +971,13 @@ module TestAlgorithms E_x_taylor = @constinferred expectation_value(rho_taylor_1, 1 => S_x()) E_xx_taylor = @constinferred expectation_value(rho_taylor_1, (1, 2) => S_xx()) + # WII + rho_wii = make_time_mpo(H, beta, WII(); imaginary_evolution = true) + Z_wii = tr(rho_wii)^(1 / L) + @test Z_wii ≈ Z_dense_1 atol = 1.0e-2 + @test expectation_value(rho_wii, 1 => S_x()) ≈ E_x_taylor atol = 1.0e-2 + @test expectation_value(rho_wii, (1, 2) => S_xx()) ≈ E_xx_taylor atol = 1.0e-2 + # MPO multiplication rho_mps = convert(FiniteMPS, rho_taylor_1) rho_mps, = approximate(rho_mps, (rho_taylor_1, rho_mps), DMRG2(; trscheme, verbosity)) From 87e2e471d5fbe1a62b0ba5615418e21ac7ccd06a Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sun, 28 Sep 2025 08:52:37 -0400 Subject: [PATCH 24/27] improve `real(::MPOHamiltonian)` coverage --- test/operators.jl | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/test/operators.jl b/test/operators.jl index 3ea1ae003..61e26c67d 100644 --- a/test/operators.jl +++ b/test/operators.jl @@ -124,7 +124,7 @@ module TestOperators @testset "Finite MPOHamiltonian" begin L = 3 T = ComplexF64 - for V in (ℂ^2, U1Space(-1 => 1, 0 => 1, 1 => 1)) + for T in (Float64, ComplexF64), V in (ℂ^2, U1Space(-1 => 1, 0 => 1, 1 => 1)) lattice = fill(V, L) O₁ = randn(T, V, V) O₁ += O₁' @@ -136,6 +136,20 @@ module TestOperators H2 = FiniteMPOHamiltonian(lattice, (i, i + 1) => O₂ for i in 1:(L - 1)) H3 = FiniteMPOHamiltonian(lattice, 1 => O₁, (2, 3) => O₂, (1, 3) => O₂) + @test scalartype(H1) == scalartype(H2) == scalartype(H3) == T + if !(T <: Complex) + for H in (H1, H2, H3) + Hc = @constinferred complex(H) + @test scalartype(Hc) == complex(T) + # cannot define `real(H)`, so only test elementwise + for (Wc, W) in zip(parent(Hc), parent(H)) + Wr = @constinferred real(Wc) + @test scalartype(Wr) == T + @test Wr ≈ W atol = 1.0e-6 + end + end + end + # check if constructor works by converting back to tensormap H1_tm = convert(TensorMap, H1) operators = vcat(fill(E, L - 1), O₁) From 620da02c243cadee0497168d1f7e692bc84cef42 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sun, 28 Sep 2025 19:17:22 -0400 Subject: [PATCH 25/27] small fix in dot(::InfiniteMPS, ::InfiniteMPO, ::InfiniteMPS) --- src/operators/mpo.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/operators/mpo.jl b/src/operators/mpo.jl index 5934aefda..0ad108a1b 100644 --- a/src/operators/mpo.jl +++ b/src/operators/mpo.jl @@ -319,7 +319,7 @@ function TensorKit.dot( randomize!(ρ₀) val, = fixedpoint( - TransferMatrix(ket.AL, parent(mpo), bra.AL), ρ₀, :LM; + flip(TransferMatrix(ket.AL, parent(mpo), bra.AL)), ρ₀, :LM; ishermitian, krylovdim, kwargs... ) return val From e404ef687fb6ec89550198ab47baadc9363c3982 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sun, 28 Sep 2025 19:18:02 -0400 Subject: [PATCH 26/27] Distinguish between finite and infinite traces of MPO --- src/operators/mpo.jl | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/src/operators/mpo.jl b/src/operators/mpo.jl index 0ad108a1b..c152c578a 100644 --- a/src/operators/mpo.jl +++ b/src/operators/mpo.jl @@ -344,7 +344,7 @@ function TensorKit.dot(mpo₁::FiniteMPO{<:MPOTensor}, mpo₂::FiniteMPO{<:MPOTe return @plansor ρ_left[1; 2] * ρ_right[2; 1] end -function LinearAlgebra.tr(mpo::MPO) +function LinearAlgebra.tr(mpo::FiniteMPO) N = length(mpo) # @assert BraidingStyle(sectortype(mpo)) isa SymmetricBraiding Nhalf = N ÷ 2 @@ -367,6 +367,14 @@ function LinearAlgebra.tr(mpo::MPO) return @plansor ρ_left[(); 1] * ρ_right[1; ()] end +function LinearAlgebra.tr(mpo::InfiniteMPO; ishermitian = false, kwargs...) + T = prod(parent(mpo)) do O + @plansor Tᵢ[-1; -2] := O[-1 3; 1 2] * τ[1 2; -2 3] + end + v₀ = randn(scalartype(mpo), domain(T)) + val, = fixedpoint(Base.Fix1(*, T), v₀, :LM; kwargs...) + return val +end function Base.isapprox( mpo₁::MPO, mpo₂::MPO; From 6020a39d930e408b0254900911335b54035ae429 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sun, 28 Sep 2025 19:18:19 -0400 Subject: [PATCH 27/27] improve coverage some more --- test/algorithms.jl | 131 ++++++++++++++++++++++++++++++--------------- 1 file changed, 87 insertions(+), 44 deletions(-) diff --git a/test/algorithms.jl b/test/algorithms.jl index fdfdd9e99..1589b586d 100644 --- a/test/algorithms.jl +++ b/test/algorithms.jl @@ -950,51 +950,94 @@ module TestAlgorithms end @testset "Finite temperature methods" begin - L = 6 - H = transverse_field_ising(; L) - trscheme = truncdim(20) - verbosity = 1 - beta = 0.1 + imaginary_evolution = true + @testset "Finite-size" begin + L = 6 + H = transverse_field_ising(; L) + trscheme = truncdim(20) + verbosity = 1 + beta = 0.1 + + # exact diagonalization + H_dense = convert(TensorMap, H) + Z_dense_1 = tr(exp(-beta * H_dense))^(1 / L) + Z_dense_2 = tr(exp(-2beta * H_dense))^(1 / L) + + # taylor cluster + rho_taylor_1 = make_time_mpo(H, beta, TaylorCluster(; N = 2); imaginary_evolution) + Z_taylor_1 = tr(rho_taylor_1)^(1 / L) + @test Z_taylor_1 ≈ Z_dense_1 atol = 1.0e-2 + Z_taylor_2 = real(dot(rho_taylor_1, rho_taylor_1))^(1 / L) + @test Z_taylor_2 ≈ Z_dense_2 atol = 1.0e-2 + + E_x_taylor = @constinferred expectation_value(rho_taylor_1, 1 => S_x()) + E_xx_taylor = @constinferred expectation_value(rho_taylor_1, (1, 2) => S_xx()) + + # WII + rho_wii = make_time_mpo(H, beta, WII(); imaginary_evolution) + Z_wii = tr(rho_wii)^(1 / L) + @test Z_wii ≈ Z_dense_1 atol = 1.0e-2 + @test expectation_value(rho_wii, 1 => S_x()) ≈ E_x_taylor atol = 1.0e-2 + @test expectation_value(rho_wii, (1, 2) => S_xx()) ≈ E_xx_taylor atol = 1.0e-2 + + # MPO multiplication + rho_mps = convert(FiniteMPS, rho_taylor_1) + rho_mps, = approximate(rho_mps, (rho_taylor_1, rho_mps), DMRG2(; trscheme, verbosity)) + Z_mpomul = tr(convert(FiniteMPO, rho_mps))^(1 / L) + @test Z_mpomul ≈ Z_dense_2 atol = 1.0e-2 + + # TDVP + rho_0 = MPSKit.infinite_temperature_density_matrix(H) + rho_0_mps = convert(FiniteMPS, rho_0) + rho_mps, = timestep(rho_0_mps, H, 0.0, beta, TDVP2(; trscheme); imaginary_evolution) + Z_tdvp = real(dot(rho_mps, rho_mps))^(1 / L) + @test Z_tdvp ≈ Z_dense_2 atol = 1.0e-2 + + @test expectation_value(rho_0_mps, 1 => S_x()) ≈ 0 + @test expectation_value(rho_0_mps, (1, 2) => S_xx()) ≈ 0 + @test expectation_value(rho_mps, 1 => S_x()) ≈ E_x_taylor atol = 1.0e-2 + @test expectation_value(rho_mps, (1, 2) => S_xx()) ≈ E_xx_taylor atol = 1.0e-2 + end - # exact diagonalization - H_dense = convert(TensorMap, H) - Z_dense_1 = tr(exp(-beta * H_dense))^(1 / L) - Z_dense_2 = tr(exp(-2beta * H_dense))^(1 / L) - - # taylor cluster - rho_taylor_1 = make_time_mpo(H, beta, TaylorCluster(; N = 2); imaginary_evolution = true) - Z_taylor_1 = tr(rho_taylor_1)^(1 / L) - @test Z_taylor_1 ≈ Z_dense_1 atol = 1.0e-2 - Z_taylor_2 = real(dot(rho_taylor_1, rho_taylor_1))^(1 / L) - @test Z_taylor_2 ≈ Z_dense_2 atol = 1.0e-2 - - E_x_taylor = @constinferred expectation_value(rho_taylor_1, 1 => S_x()) - E_xx_taylor = @constinferred expectation_value(rho_taylor_1, (1, 2) => S_xx()) - - # WII - rho_wii = make_time_mpo(H, beta, WII(); imaginary_evolution = true) - Z_wii = tr(rho_wii)^(1 / L) - @test Z_wii ≈ Z_dense_1 atol = 1.0e-2 - @test expectation_value(rho_wii, 1 => S_x()) ≈ E_x_taylor atol = 1.0e-2 - @test expectation_value(rho_wii, (1, 2) => S_xx()) ≈ E_xx_taylor atol = 1.0e-2 - - # MPO multiplication - rho_mps = convert(FiniteMPS, rho_taylor_1) - rho_mps, = approximate(rho_mps, (rho_taylor_1, rho_mps), DMRG2(; trscheme, verbosity)) - Z_mpomul = tr(convert(FiniteMPO, rho_mps))^(1 / L) - @test Z_mpomul ≈ Z_dense_2 atol = 1.0e-2 - - # TDVP - rho_0 = MPSKit.infinite_temperature_density_matrix(H) - rho_0_mps = convert(FiniteMPS, rho_0) - rho_mps, = timestep(rho_0_mps, H, 0.0, beta, TDVP2(; trscheme); imaginary_evolution = true) - Z_tdvp = real(dot(rho_mps, rho_mps))^(1 / L) - @test Z_tdvp ≈ Z_dense_2 atol = 1.0e-2 - - @test expectation_value(rho_0_mps, 1 => S_x()) ≈ 0 - @test expectation_value(rho_0_mps, (1, 2) => S_xx()) ≈ 0 - @test expectation_value(rho_mps, 1 => S_x()) ≈ E_x_taylor atol = 1.0e-2 - @test expectation_value(rho_mps, (1, 2) => S_xx()) ≈ E_xx_taylor atol = 1.0e-2 + @testset "Infinite-size" begin + H = transverse_field_ising() + trscheme = truncdim(20) + verbosity = 1 + beta = 0.1 + + # taylor cluster + alg_taylor = TaylorCluster(; N = 2) + rho_taylor = make_time_mpo(H, beta, alg_taylor; imaginary_evolution) + rho_taylor_2 = make_time_mpo(H, 2beta, alg_taylor; imaginary_evolution) + E_taylor = @constinferred expectation_value(rho_taylor, H) + E_taylor2 = @constinferred expectation_value(rho_taylor_2, H) + + E_z_taylor = @constinferred expectation_value(rho_taylor, 1 => S_z()) + @test E_z_taylor ≈ 0 atol = 1.0e-4 # no spontaneous symmetry breaking at finite T + + # WII + rho_wii = make_time_mpo(H, beta, WII(); imaginary_evolution) + @test expectation_value(rho_wii, H) ≈ E_taylor atol = 1.0e-2 + @test expectation_value(rho_wii, 1 => S_z()) ≈ E_z_taylor atol = 1.0e-2 + + # MPO multiplication + rho_mps = convert(InfiniteMPS, rho_taylor) + rho_mps2, = approximate(rho_mps, (rho_taylor, rho_mps), IDMRG(; verbosity)) + E_mpomul = expectation_value(rho_mps2, H) + @test E_mpomul ≈ E_taylor2 atol = 1.0e-2 + + # TDVP + rho_0 = MPSKit.infinite_temperature_density_matrix(H) + rho_0_mps, = changebonds(convert(InfiniteMPS, rho_0), H, OptimalExpand(; trscheme = truncdim(20))) + rho_mps_tdvp, = timestep(rho_0_mps, H, 0.0, beta, TDVP(); imaginary_evolution) + E_tdvp = expectation_value(rho_mps_tdvp, H) + @test E_tdvp ≈ E_taylor atol = 1.0e-2 + + num_vals = 2 + vals_taylor = @constinferred(transfer_spectrum(convert(InfiniteMPS, rho_taylor); num_vals)) + vals_mps = @constinferred(transfer_spectrum(rho_mps; num_vals)) + @test vals_taylor[1:num_vals] ≈ vals_mps[1:num_vals] + end end @testset "Sector conventions" begin