From 0943a3e95488ec60e90aa1206fe651fd6daad91c Mon Sep 17 00:00:00 2001 From: Jan Schneider Date: Wed, 28 May 2025 11:43:20 +0200 Subject: [PATCH 01/15] add vscode to gitignore --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index ace59e2..b73f6ce 100644 --- a/.gitignore +++ b/.gitignore @@ -2,3 +2,4 @@ /examples/vumps/Manifest.toml /test/Manifest.toml *.swp +.vscode/ \ No newline at end of file From b9cf4783fb76147fc57f9ee42c75d9e56adb4436 Mon Sep 17 00:00:00 2001 From: Jan Schneider Date: Fri, 5 Sep 2025 14:21:05 +0200 Subject: [PATCH 02/15] bump ITensor version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index dce4615..2357c4b 100644 --- a/Project.toml +++ b/Project.toml @@ -20,7 +20,7 @@ SplitApplyCombine = "03a91e81-4c3e-53e1-a0a4-9c0c8f19dd66" Compat = "3, 4" HDF5 = "0.15, 0.16, 0.17" ITensorMPS = "0.3" -ITensors = "0.7, 0.8" +ITensors = "0.7, 0.8, 0.9" Infinities = "0.1" IterTools = "1" KrylovKit = "0.5, 0.6, 0.7, 0.8, 0.9" From f0ddcc68c95ff51785fef579895c6df11acd6096 Mon Sep 17 00:00:00 2001 From: Jan Schneider Date: Wed, 10 Sep 2025 12:22:48 +0200 Subject: [PATCH 03/15] working code for VUMPS of an exponentially interacting Hamiltonian with InfiniteBlockMPO --- src/vumps_mpo.jl | 98 ++++++++++++++++++++++++++++++++++++------------ 1 file changed, 73 insertions(+), 25 deletions(-) diff --git a/src/vumps_mpo.jl b/src/vumps_mpo.jl index 4a7c561..067493f 100644 --- a/src/vumps_mpo.jl +++ b/src/vumps_mpo.jl @@ -103,8 +103,12 @@ function left_environment(H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS; tol=1e- s = siteinds(only, ψ) δʳ(n) = δ(dag(r[n]), prime(r[n])) δˡ(n) = δ(l[n], l′[n]) + # op("Id", s[1]) is permuted w.r.t. below line δˢ(n) = δ(dag(s[n]), prime(s[n])) + # this code is limited to homogeneous physical spaces on each site within the unit cell! + phys_dim = dim(s[1]) + eₗ = [0.0] dₕ = size(H[1])[1] #Ls = [Vector{ITensor}(undef, dₕ) for j in 1:N] @@ -112,11 +116,11 @@ function left_environment(H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS; tol=1e- #Building the L vector for n_1 = 1 # TM is 2 3 ... N 1 localR = ψ.C[1] * δʳ(1) * ψ′.C[1] #to revise - for n in reverse(1:(dₕ - 1)) + for a in reverse(1:(dₕ - 1)) temp_Ls = apply_left_transfer_matrix( - translatecell(translator(ψ), Ls[1][n + 1], -1), n + 1, H, ψ, 2 - N + translatecell(translator(ψ), Ls[1][a + 1], -1), a + 1, H, ψ, 2 - N ) - for j in 1:n + for j in 1:a if isassigned(temp_Ls, j) if isassigned(Ls[1], j) Ls[1][j] += temp_Ls[j] @@ -125,25 +129,45 @@ function left_environment(H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS; tol=1e- end end end - if !isempty(H[1][n, n]) - λ = H[1][n, n][1, 1] - δiag = δˢ(1) + # starting and ending identity OPs are order 2, + # as they have NO link dimension, + # because they terminate and do not transport QN numbers. + # Intermediate OPs proportional to the identity DO transport QNs and must be order 4! + if (!isempty(H[1][a, a])) && (order(H[1][a, a]) == 2 ) + λ = H[1][a, a][1, 1] + # δiag = δˢ(1) # not used #@assert norm(H[1][n, n] - λ * δiag) == 0 "Non identity diagonal not implemented in MPO" @assert abs(λ) <= 1 "Diverging term" - eₗ[1] = (Ls[1][n] * localR)[] - Ls[1][n] += -(eₗ[1] * denseblocks(δˡ(1))) + eₗ[1] = (Ls[1][a] * localR)[] + Ls[1][a] += -(eₗ[1] * denseblocks(δˡ(1))) if λ == 1 - A = AOᴸ(ψ, H, n) - Ls[1][n], info = linsolve(A, Ls[1][n], 1, -1; tol=tol) + A = AOᴸ(ψ, H, a) + Ls[1][a], info = linsolve(A, Ls[1][a], 1, -1; tol=tol) else - println("Not implemented") + # println("Not implemented") + error("This case cannot exist! The last and first term of the diagonal of the MPO must be the identity.") flush(stdout) flush(stderr) end + elseif (!isempty(H[1][a, a])) && (order(H[1][a, a]) == 4 ) + # assert diagonal operator! + λ = H[1][a, a][1, 1, 1, 1] + test_Tensor = reshape(diagm(λ*ones(phys_dim)), (1,phys_dim,phys_dim,1)) + # verify that this OP = λ⋅Id + @assert norm(test_Tensor - array(H[1][a, a])) ≈ 0.0 "Only operators proportional to the identity are allowed on the diagonal!" + + @assert (λ <= 1.0) && (λ >= 0.0) "The proportionality factor of the operator to the identity must be 0 ≤ λ ≤ 1" + + # solve equations with linsolve + eₗ[1] = (Ls[1][a] * localR)[] + solo_link = inds(Ls[1][a])[1] + Ls[1][a] += -(eₗ[1] * denseblocks(δˡ(1))*setelt(solo_link[1])) + A = AOᴸ(ψ, H, a) + Ls[1][a], info = linsolve(A, Ls[1][a], 1, -1; tol=tol) end end - for n in 2:N - Ls[n] = apply_local_left_transfer_matrix(Ls[n - 1], H, ψ, n) + for a in 2:N + Ls[a] = apply_local_left_transfer_matrix(Ls[a - 1], H, ψ, a) end return CelledVector(Ls), eₗ[1] end @@ -192,12 +216,14 @@ function initialize_right_environment( Rs[1] = ITensor(Float64, link, dag(prime(link))) end for j in 2:(dₕ - 1) - mpo_link = only(uniqueinds(H[n - 1][1, j], sit)) + mpo_link = only(uniqueinds(H[n - 1][dₕ,j], sit)) Rs[j] = ITensor(Float64, dag(mpo_link), link, dag(prime(link))) end return Rs end +Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}} + function apply_local_right_transfer_matrix( Lstart::Vector{ITensor}, H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS, n_1::Int64 ) @@ -256,17 +282,20 @@ function right_environment(H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS; tol=1e δˡ(n) = δ(l[n], dag(prime(l[n]))) δˢ(n) = δ(dag(s[n]), prime(s[n])) + # this code is limited to homogeneous physical spaces on each site within the unit cell! + phys_dim = dim(s[1]) + eᵣ = [0.0] dₕ = size(H[1])[1] Rs = [initialize_right_environment(H, ψ, j; init_first=true) for j in 1:N] #Building the R vector for n_1 = 1 # TM is 2-N 3-N ... 0 localL = ψ.C[0] * δˡ(0) * dag(prime(ψ.C[0])) - for n in 2:dₕ + for a in 2:dₕ temp_Rs = apply_right_transfer_matrix( - translatecell(translator(ψ), Rs[1][n - 1], 1), n - 1, H, ψ, N + translatecell(translator(ψ), Rs[1][a - 1], 1), a - 1, H, ψ, N ) - for j in n:dₕ + for j in a:dₕ if isassigned(temp_Rs, j) if isassigned(Rs[1], j) Rs[1][j] += temp_Rs[j] @@ -275,21 +304,40 @@ function right_environment(H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS; tol=1e end end end - if !isempty(H[1][n, n]) - λ = H[1][n, n][1, 1] - δiag = δˢ(1) + # starting and ending identity OPs are order 2, + # as they have NO link dimension, + # because they terminate and do not transport QN numbers. + # Intermediate OPs proportional to the identity DO transport QNs and must be order 4! + if (!isempty(H[1][a, a])) && (order(H[1][a, a]) == 2 ) + λ = H[1][a, a][1, 1] + # δiag = δˢ(1) #@assert norm(H[1][n, n] - λ * δiag) == 0 "Non identity diagonal not implemented in MPO" @assert abs(λ) <= 1 "Diverging term" - eᵣ[1] = (localL * Rs[1][n])[] - Rs[1][n] += -(eᵣ[1] * denseblocks(δʳ(0))) + eᵣ[1] = (localL * Rs[1][a])[] + Rs[1][a] += -(eᵣ[1] * denseblocks(δʳ(0))) if λ == 1 - A = AOᴿ(ψ, H, n) - Rs[1][n], info = linsolve(A, Rs[1][n], 1, -1; tol=tol) + A = AOᴿ(ψ, H, a) + Rs[1][a], info = linsolve(A, Rs[1][a], 1, -1; tol=tol) else - println("Not yet implemented") + error("This case cannot exist! The last and first term of the diagonal of the MPO must be the identity.") flush(stdout) flush(stderr) end + elseif (!isempty(H[1][a, a])) && (order(H[1][a, a]) == 4 ) + # assert diagonal operator! + λ = H[1][a, a][1, 1, 1, 1] + test_Tensor = reshape(diagm(λ*ones(phys_dim)), (1,phys_dim,phys_dim,1)) + # verify that this OP = λ⋅Id + @assert norm(test_Tensor - array(H[1][a, a])) ≈ 0.0 "Only operators proportional to the identity are allowed on the diagonal!" + + @assert (λ <= 1.0) && (λ >= 0.0) "The proportionality factor of the operator to the identity must be 0 ≤ λ ≤ 1" + + # solve euqations with linsolve + eᵣ[1] = (localL * Rs[1][a])[] + solo_link = inds(Rs[1][a])[1] + Rs[1][a] += -(eᵣ[1] * denseblocks(δʳ(0)) * setelt(solo_link[1])) + A = AOᴿ(ψ, H, a) + Rs[1][a], info = linsolve(A, Rs[1][a], 1, -1; tol=tol) end end if N > 1 From 9babf6239f097120648aa88ea244f8f7b99e9c0d Mon Sep 17 00:00:00 2001 From: Jan Schneider Date: Thu, 11 Sep 2025 11:31:29 +0200 Subject: [PATCH 04/15] update element type of infMPS to be also Complex for real-time evo --- src/infinitecanonicalmps.jl | 20 +++++++++++++------- src/vumps_mpo.jl | 27 +++++++++++++++++++++------ 2 files changed, 34 insertions(+), 13 deletions(-) diff --git a/src/infinitecanonicalmps.jl b/src/infinitecanonicalmps.jl index 6d1dd6a..a995132 100644 --- a/src/infinitecanonicalmps.jl +++ b/src/infinitecanonicalmps.jl @@ -171,8 +171,12 @@ function ITensorMPS.linkinds(ψ::InfiniteMPS) return CelledVector([linkinds(ψ, (n, n + 1)) for n in 1:N], translator(ψ)) end +function InfMPS(eltype::Type{<:Number}, s::Vector, f::Function, translator::Function=translatecelltags, ) + return InfMPS(eltype, infsiteinds(s, translator), f) +end + function InfMPS(s::Vector, f::Function, translator::Function=translatecelltags) - return InfMPS(infsiteinds(s, translator), f) + return InfMPS(Float64, infsiteinds(s, translator), f) end function indval(iv::Pair) @@ -219,9 +223,9 @@ function insert_linkinds!(A; left_dir=ITensors.Out) return A end -function UniformMPS(s::CelledVector, f::Function; left_dir=ITensors.Out) +function UniformMPS(eltype::Type{<:Number}, s::CelledVector, f::Function; left_dir=ITensors.Out) sᶜ¹ = s[Cell(1)] - A = InfiniteMPS([ITensor(sⁿ) for sⁿ in sᶜ¹], translator(s)) + A = InfiniteMPS([ITensor(eltype,sⁿ) for sⁿ in sᶜ¹], translator(s)) #A.data.translator = translator(s) N = length(sᶜ¹) for n in 1:N @@ -233,16 +237,18 @@ function UniformMPS(s::CelledVector, f::Function; left_dir=ITensors.Out) return A end -function InfMPS(s::CelledVector, f::Function) +InfMPS(s::CelledVector, f::Function) = InfMPS(Float64, s::CelledVector, f::Function) + +function InfMPS(eltype::Type{<:Number}, s::CelledVector, f::Function) # TODO: rename cell_length N = length(s) - ψL = UniformMPS(s, f; left_dir=ITensors.Out) - ψR = UniformMPS(s, f; left_dir=ITensors.In) + ψL = UniformMPS(eltype,s, f; left_dir=ITensors.Out) + ψR = UniformMPS(eltype,s, f; left_dir=ITensors.In) ψC = InfiniteMPS(N, translator(s)) l = linkinds(ψL) r = linkinds(ψR) for n in 1:N - ψCₙ = ITensor(dag(l[n])..., r[n]...) + ψCₙ = ITensor(eltype,dag(l[n])..., r[n]...) ψCₙ[l[n]... => 1, r[n]... => 1] = 1.0 ψC[n] = ψCₙ end diff --git a/src/vumps_mpo.jl b/src/vumps_mpo.jl index 067493f..96a6f1a 100644 --- a/src/vumps_mpo.jl +++ b/src/vumps_mpo.jl @@ -109,7 +109,9 @@ function left_environment(H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS; tol=1e- # this code is limited to homogeneous physical spaces on each site within the unit cell! phys_dim = dim(s[1]) - eₗ = [0.0] + eltype = ITensorMPS.promote_itensor_eltype(ψ) + + eₗ = [zero(eltype)] dₕ = size(H[1])[1] #Ls = [Vector{ITensor}(undef, dₕ) for j in 1:N] Ls = [initialize_left_environment(H, ψ, j; init_last=true) for j in 1:N] @@ -285,7 +287,9 @@ function right_environment(H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS; tol=1e # this code is limited to homogeneous physical spaces on each site within the unit cell! phys_dim = dim(s[1]) - eᵣ = [0.0] + eltype = ITensorMPS.promote_itensor_eltype(ψ) + + eᵣ = [zero(eltype)] dₕ = size(H[1])[1] Rs = [initialize_right_environment(H, ψ, j; init_first=true) for j in 1:N] #Building the R vector for n_1 = 1 @@ -413,8 +417,14 @@ function tdvp_iteration_sequential( Ãᴸ = InfiniteMPS(Vector{ITensor}(undef, N)) Ãᴿ = InfiniteMPS(Vector{ITensor}(undef, N)) - eL = zeros(N) - eR = zeros(N) + eltype_ψ = ITensorMPS.promote_itensor_eltype(ψ) + + eltype_t = typeof(time_step) + + eltype = typeof(one(eltype_ψ) * one(eltype_t)) + + eL = zeros(eltype,N) + eR = zeros(eltype,N) for n in 1:N L, eL[n] = left_environment(H, ψ; tol=_solver_tol) #TODO currently computing two many of them R, eR[n] = right_environment(H, ψ; tol=_solver_tol) #TODO currently computing two many of them @@ -475,8 +485,13 @@ function tdvp_iteration_parallel( Ãᴸ = InfiniteMPS(Vector{ITensor}(undef, N)) Ãᴿ = InfiniteMPS(Vector{ITensor}(undef, N)) - eL = zeros(1) - eR = zeros(1) + eltype_ψ = ITensorMPS.promote_itensor_eltype(ψ) + eltype_t = typeof(time_step) + eltype = typeof(one(eltype_ψ) * one(eltype_t)) + + eL = zeros(eltype,1) + eR = zeros(eltype,1) + L, eL[1] = left_environment(H, ψ; tol=_solver_tol) #TODO currently computing two many of them R, eR[1] = right_environment(H, ψ; tol=_solver_tol) #TODO currently computing two many of them for n in 1:N From f99532343be8fec537703224b412b3567e898af3 Mon Sep 17 00:00:00 2001 From: Jan Schneider Date: Thu, 11 Sep 2025 11:46:57 +0200 Subject: [PATCH 05/15] add base tests for general blockmpo --- test/test_vumps_GeneralMPO.jl | 258 ++++++++++++++++++++++++++++++++++ 1 file changed, 258 insertions(+) create mode 100644 test/test_vumps_GeneralMPO.jl diff --git a/test/test_vumps_GeneralMPO.jl b/test/test_vumps_GeneralMPO.jl new file mode 100644 index 0000000..4f1393a --- /dev/null +++ b/test/test_vumps_GeneralMPO.jl @@ -0,0 +1,258 @@ +using ITensors, ITensorMPS +using ITensorInfiniteMPS + +base_path = joinpath(pkgdir(ITensorInfiniteMPS), "examples", "vumps", "src") +src_files = ["vumps_subspace_expansion.jl", "entropy.jl"] +for f in src_files + include(joinpath(base_path, f)) +end + + +""" +build the transverse field Ising model with exponential interactions in the thermodynamic limit 𝑁 → ∞ +use convention H = -J ∑_{n,m=−∞}^∞ σˣₙ σˣₘ ⋅ λ^(-|n-m-1|) - hz ∑_{n=−∞}^∞ σᶻₙ +with +""" +function InfiniteExpHTFI( + sites::CelledVector{<:Index}, + λ, # exponential decay base + J, # interaction kinetic (tunneling) coupling + hz; # interaction strength of the transverse field + kwargs... +) + if abs(λ) > 1.0 + throw(ArgumentError("cannot implement exponential decay with base larger than 1, λ = $(λ)!")) + end + + link_dimension = 3 + + N = length(sites) + + EType = eltype(union(λ, J, hz)) + + linkindices = CelledVector( + hasqns(sites.data) ? + [Index([QN("SzParity",1,2) => 1], "Link,c=1,n=$n") for n in 1:N] : [Index(1, "Link,c=1,n=$(n)") for n in 1:N] + ) + + mpos = [Matrix{ITensor}(undef, link_dimension, link_dimension) for i in 1:N] + for n in 1:N + # define local matrix Hmat with empty tensors as local operators + Hmat = fill( + ITensor(EType, dag(sites[n]), prime(sites[n])), link_dimension, link_dimension + ) + # left link index ll with daggered QN conserving direction (if applicable) + ll = dag(linkindices[n-1]) + # right link index rl + rl = linkindices[n] + + # add both Identities as netral elements in the MPO + # replace all known tensors from empty to known interactions + Hmat[1,1] = op("Id", sites[n]) + Hmat[3,3] = op("Id", sites[n]) + # local nearest neighbour and exp. decaying interaction terms + Hmat[2,1] = op("X", sites[n]) + if !iszero(λ) + Hmat[2,2] = op("Id", sites[n]) * λ # λ Id, on the diagonal + end + Hmat[3,2] = op("X", sites[n]) * -J # Jxx σˣ + if !iszero(hz) + Hmat[3,1] = op("Z", sites[n]) * -hz # hz σᶻ + end + + # add all missing links that connect the interaction + # operators in the unit cell + Hmat[2,1] = setelt(ll[1]) * Hmat[2,1] + Hmat[1,2] = Hmat[1,2] * setelt(rl[1]) + Hmat[2,2] = setelt(ll[1]) * Hmat[2,2] * setelt(rl[1]) + Hmat[3,2] = setelt(rl[1]) * Hmat[3,2] + Hmat[2,3] = setelt(ll[1]) * Hmat[2,3] + mpos[n] = Hmat + end + + return InfiniteBlockMPO(mpos, sites.translator) +end + + +# sanity check if I can construct the MPO for the normal TFI correctly +""" +build the transverse field Ising model with nearest neighbour interactions in the thermodynamic limit 𝑁 → ∞ +use convention H = -J ∑_{n=−∞}^∞ σˣₙ σˣₙ₊₁ - hz ∑_{n=−∞}^∞ σᶻₙ +with +""" +function InfiniteHTFI( + sites::CelledVector{<:Index}, + kinetic_coupling, + hz; # interaction kinetic (tunneling) coupling + kwargs... +) + if abs(λ) > 1.0 + throw(ArgumentError("cannot implement exponential decay with base larger than 1, λ = $(λ)!")) + end + link_dimension = 3 + + N = length(sites) + + EType = eltype(union(J, hz)) + + linkindices = CelledVector( + hasqns(sites.data) ? + [Index([QN("SzParity",1,2) => 1], "Link,c=1,n=$n") for n in 1:N] : [Index(1, "Link,c=1,n=$(n)") for n in 1:N] + ) + + mpos = [Matrix{ITensor}(undef, link_dimension, link_dimension) for i in 1:N] + for n in 1:N + # define local matrix Hmat with empty tensors as local operators + Hmat = fill( + ITensor(EType, dag(sites[n]), prime(sites[n])), link_dimension, link_dimension + ) + # left link index ll with daggered QN conserving direction (if applicable) + ll = dag(linkindices[n-1]) + # right link index rl + rl = linkindices[n] + + # add both Identities as netral elements in the MPO + # replace all known tensors from empty to known interactions + Hmat[1,1] = op("Id", sites[n]) + Hmat[3,3] = op("Id", sites[n]) + # local nearest neighbour and exp. decaying interaction terms + Hmat[2,1] = op("X", sites[n]) + Hmat[3,2] = op("X", sites[n]) * -J # Jxx σˣ + if !iszero(hz) + Hmat[3,1] = op("Z", sites[n]) * -hz # hz σᶻ + end + + # add all missing links that connect the + # interaction operators in the unit cell + Hmat[2,1] = setelt(ll[1]) * Hmat[2,1] + # Hmat[1,2] = Hmat[1,2] * setelt(rl[1]) + Hmat[3,2] = setelt(rl[1]) * Hmat[3,2] + # Hmat[2,3] = setelt(ll[1]) * Hmat[2,3] + mpos[n] = Hmat + end + return InfiniteBlockMPO(mpos, sites.translator) +end + +function expect_two_site(ψ::InfiniteCanonicalMPS, h::ITensor, n1n2) + n1, n2 = n1n2 + ϕ = ψ.AL[n1] * ψ.AL[n2] * ψ.C[n2] + return inner(ϕ, apply(h, ϕ)) +end + + +function expect_two_site(ψ::InfiniteCanonicalMPS, h::MPO, n1n2) + return expect_two_site(ψ, contract(h), n1n2) +end + + +function energy_local(ψ1, ψ2, h::ITensor) + ϕ = ψ1 * ψ2 + return (noprime(ϕ * h) * dag(ϕ))[] +end + +energy_local(ψ1, ψ2, h::MPO) = energy_local(ψ1, ψ2, prod(h)) +# Check translational invariance +function ITensorMPS.expect(ψ, o) + return (noprime(ψ * op(o, filterinds(ψ, "Site")...)) * dag(ψ))[] +end + + +maxdim = 16 # Maximum bond dimension +cutoff = 1e-10 # Singular value cutoff when increasing the bond dimension +max_vumps_iters = 100 # Maximum number of iterations of the VUMPS/TDVP algorithm at a fixed bond dimension +tol = 1e-8 # Precision error tolerance for outer loop of VUMPS or TDVP +outer_iters = 4 # Number of times to increase the bond dimension +time_step = -Inf # -Inf corresponds to VUMPS, finite time_step corresponds to TDVP +solver_tol = (x -> x / 100) # Tolerance for the local solver (eigsolve in VUMPS and exponentiate in TDVP) +multisite_update_alg = "parallel" # Choose between ["sequential", "parallel"]. Only parallel works with TDVP. +conserve_qns = true # Whether or not to conserve spin parity + + +nsite = 2 # Number of sites in the unit cell +initstate(n) = "↑" +s = infsiteinds("S=1/2", nsite; initstate, conserve_szparity=conserve_qns) + +ψ = InfMPS(s, initstate) + +@show norm(contract(ψ.AL[1:nsite]..., ψ.C[nsite]) - contract(ψ.C[0], ψ.AR[1:nsite]...)) + +# J = -0.25 +# hz = 2.0 +# λ = 0.4 + +J = 1 +hz = 1.0 +λ = 0.4 + +# H_test = InfiniteHTFI(s,J,hz) +H_test = InfiniteHTFI(s,J,hz) +H_test0 = InfiniteExpHTFI(s,0.0,J,hz) +H_test_exp = InfiniteExpHTFI(s,λ,J,hz) + +vumps_kwargs = ( + tol=tol, + maxiter=max_vumps_iters, + solver_tol=solver_tol, + multisite_update_alg=multisite_update_alg, +) + +subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim) + +H_ref = InfiniteSum{MPO}(Model("ising"), s; J=J, h=hz) + + +ψ_ref = vumps_subspace_expansion(H_ref, ψ; outer_iters, subspace_expansion_kwargs, vumps_kwargs) + +ψ_test_NN = vumps_subspace_expansion(H_test, ψ; outer_iters, subspace_expansion_kwargs, vumps_kwargs) +ψ_test0_NN = vumps_subspace_expansion(H_test0, ψ; outer_iters, subspace_expansion_kwargs, vumps_kwargs) + + +E_ref = energy_local(ψ_ref.AL[1], ψ_ref.AL[2] * ψ_ref.C[2], H_ref[(1, 2)]) + +L_test, energy_test = ITensorInfiniteMPS.left_environment(H_test, ψ_test_NN; tol=1e-10); +L_test0, energy_test0 = ITensorInfiniteMPS.left_environment(H_test0, ψ_test0_NN; tol=1e-10); + +energy_test /= length(ψ_test_NN) +energy_test0 /= length(ψ_test0_NN) + +@test energy_test ≈ E_ref +@test energy_test0 ≈ E_ref + +@show E_ref + 4/π +@show energy_test + 4/π +@show energy_test0 + 4/π + + +# energy_local(ψ_NN.AL[1], ψ_NN.AL[2] * ψ.C[2], H_test[(1, 2)]) +# energy_local(ψ_test_NN,ψ_test_NN,H_test) + +function ITensorMPS.expect(ψ, o) + return (noprime(ψ * op(o, filterinds(ψ, "Site")...)) * dag(ψ))[] +end + + +@test isapprox(expect(ψ_ref.AL[1] * ψ_ref.C[1], "Z"), expect(ψ_test_NN.AL[1] * ψ_test_NN.C[1], "Z"); rtol=1e-6) +@test isapprox(expect(ψ_ref.AL[1] * ψ_ref.C[1], "X"), expect(ψ_test_NN.AL[1] * ψ_test_NN.C[1], "X"); atol=1e-6) + +@test isapprox(expect(ψ_ref.AL[1] * ψ_ref.C[1], "Z"), expect(ψ_test0_NN.AL[1] * ψ_test0_NN.C[1], "Z"); rtol=1e-6) +@test isapprox(expect(ψ_ref.AL[1] * ψ_ref.C[1], "X"), expect(ψ_test0_NN.AL[1] * ψ_test0_NN.C[1], "X"); rtol=1e-6) + +entropy(ψ_ref,1) +entropy(ψ_test_NN,1) +entropy(ψ_test0_NN,1) + +@test isapprox(entropy(ψ_ref,1), entropy(ψ_test_NN,1); rtol=1e-6) +@test isapprox(entropy(ψ_ref,1), entropy(ψ_test0_NN,1); rtol=1e-6) + +################## +#### actually exponential interactions +################# +ψ_exp = vumps_subspace_expansion(H_test_exp, ψ; outer_iters, subspace_expansion_kwargs, vumps_kwargs) + +L_test, energy_test_exp = ITensorInfiniteMPS.left_environment(H_test_exp, ψ_exp; tol=1e-10); + + +E_gs_04_MPSKit = -1.8177300191088515 # computed at χ=16, λ=0.4, Jₓₓ=1.0, hz=1.0 +E_gs_04_ITensor = energy_test_exp / length(ψ_exp) +@test E_gs_04_MPSKit ≈ E_gs_04_ITensor + From f6b1ae196d83124960cbf8da726f1ef55316cee8 Mon Sep 17 00:00:00 2001 From: Jan Schneider Date: Fri, 12 Sep 2025 18:33:21 +0200 Subject: [PATCH 06/15] update element type for InfMPS --- src/vumps_mpo.jl | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/src/vumps_mpo.jl b/src/vumps_mpo.jl index 96a6f1a..090d1b1 100644 --- a/src/vumps_mpo.jl +++ b/src/vumps_mpo.jl @@ -109,9 +109,9 @@ function left_environment(H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS; tol=1e- # this code is limited to homogeneous physical spaces on each site within the unit cell! phys_dim = dim(s[1]) - eltype = ITensorMPS.promote_itensor_eltype(ψ) + EType = ITensorMPS.promote_itensor_eltype(ψ) - eₗ = [zero(eltype)] + eₗ = zeros(EType,1) dₕ = size(H[1])[1] #Ls = [Vector{ITensor}(undef, dₕ) for j in 1:N] Ls = [initialize_left_environment(H, ψ, j; init_last=true) for j in 1:N] @@ -287,9 +287,9 @@ function right_environment(H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS; tol=1e # this code is limited to homogeneous physical spaces on each site within the unit cell! phys_dim = dim(s[1]) - eltype = ITensorMPS.promote_itensor_eltype(ψ) + EType = ITensorMPS.promote_itensor_eltype(ψ) - eᵣ = [zero(eltype)] + eᵣ = zeros(EType,1) dₕ = size(H[1])[1] Rs = [initialize_right_environment(H, ψ, j; init_first=true) for j in 1:N] #Building the R vector for n_1 = 1 @@ -417,14 +417,14 @@ function tdvp_iteration_sequential( Ãᴸ = InfiniteMPS(Vector{ITensor}(undef, N)) Ãᴿ = InfiniteMPS(Vector{ITensor}(undef, N)) - eltype_ψ = ITensorMPS.promote_itensor_eltype(ψ) + EType_ψ = ITensorMPS.promote_itensor_eltype(ψ) - eltype_t = typeof(time_step) + EType_t = typeof(time_step) - eltype = typeof(one(eltype_ψ) * one(eltype_t)) + EType = typeof(one(EType_ψ) * one(EType_t)) - eL = zeros(eltype,N) - eR = zeros(eltype,N) + eL = zeros(EType,N) + eR = zeros(EType,N) for n in 1:N L, eL[n] = left_environment(H, ψ; tol=_solver_tol) #TODO currently computing two many of them R, eR[n] = right_environment(H, ψ; tol=_solver_tol) #TODO currently computing two many of them @@ -485,12 +485,12 @@ function tdvp_iteration_parallel( Ãᴸ = InfiniteMPS(Vector{ITensor}(undef, N)) Ãᴿ = InfiniteMPS(Vector{ITensor}(undef, N)) - eltype_ψ = ITensorMPS.promote_itensor_eltype(ψ) - eltype_t = typeof(time_step) - eltype = typeof(one(eltype_ψ) * one(eltype_t)) + EType_ψ = ITensorMPS.promote_itensor_eltype(ψ) + EType_t = typeof(time_step) + EType = typeof(one(EType_ψ) * one(EType_t)) - eL = zeros(eltype,1) - eR = zeros(eltype,1) + eL = zeros(EType,1) + eR = zeros(EType,1) L, eL[1] = left_environment(H, ψ; tol=_solver_tol) #TODO currently computing two many of them R, eR[1] = right_environment(H, ψ; tol=_solver_tol) #TODO currently computing two many of them From 21446374b8b806571aa46667a272c24595aa8d75 Mon Sep 17 00:00:00 2001 From: Jan Schneider <9451057+jtschneider@users.noreply.github.com> Date: Sat, 13 Sep 2025 11:16:59 +0200 Subject: [PATCH 07/15] delete orphan line --- src/vumps_mpo.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/vumps_mpo.jl b/src/vumps_mpo.jl index 090d1b1..f9f3ada 100644 --- a/src/vumps_mpo.jl +++ b/src/vumps_mpo.jl @@ -224,7 +224,6 @@ function initialize_right_environment( return Rs end -Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}} function apply_local_right_transfer_matrix( Lstart::Vector{ITensor}, H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS, n_1::Int64 From 4689e87e0fb74799efcebe16ec5a8a698fdfcd5b Mon Sep 17 00:00:00 2001 From: Jan Schneider Date: Sat, 13 Sep 2025 11:22:39 +0200 Subject: [PATCH 08/15] JuliaFormatter run --- src/infinitecanonicalmps.jl | 16 +-- src/vumps_mpo.jl | 42 ++++---- test/test_vumps_GeneralMPO.jl | 183 +++++++++++++++++++++------------- 3 files changed, 148 insertions(+), 93 deletions(-) diff --git a/src/infinitecanonicalmps.jl b/src/infinitecanonicalmps.jl index a995132..af974f8 100644 --- a/src/infinitecanonicalmps.jl +++ b/src/infinitecanonicalmps.jl @@ -171,7 +171,9 @@ function ITensorMPS.linkinds(ψ::InfiniteMPS) return CelledVector([linkinds(ψ, (n, n + 1)) for n in 1:N], translator(ψ)) end -function InfMPS(eltype::Type{<:Number}, s::Vector, f::Function, translator::Function=translatecelltags, ) +function InfMPS( + eltype::Type{<:Number}, s::Vector, f::Function, translator::Function=translatecelltags +) return InfMPS(eltype, infsiteinds(s, translator), f) end @@ -223,9 +225,11 @@ function insert_linkinds!(A; left_dir=ITensors.Out) return A end -function UniformMPS(eltype::Type{<:Number}, s::CelledVector, f::Function; left_dir=ITensors.Out) +function UniformMPS( + eltype::Type{<:Number}, s::CelledVector, f::Function; left_dir=ITensors.Out +) sᶜ¹ = s[Cell(1)] - A = InfiniteMPS([ITensor(eltype,sⁿ) for sⁿ in sᶜ¹], translator(s)) + A = InfiniteMPS([ITensor(eltype, sⁿ) for sⁿ in sᶜ¹], translator(s)) #A.data.translator = translator(s) N = length(sᶜ¹) for n in 1:N @@ -242,13 +246,13 @@ InfMPS(s::CelledVector, f::Function) = InfMPS(Float64, s::CelledVector, f::Funct function InfMPS(eltype::Type{<:Number}, s::CelledVector, f::Function) # TODO: rename cell_length N = length(s) - ψL = UniformMPS(eltype,s, f; left_dir=ITensors.Out) - ψR = UniformMPS(eltype,s, f; left_dir=ITensors.In) + ψL = UniformMPS(eltype, s, f; left_dir=ITensors.Out) + ψR = UniformMPS(eltype, s, f; left_dir=ITensors.In) ψC = InfiniteMPS(N, translator(s)) l = linkinds(ψL) r = linkinds(ψR) for n in 1:N - ψCₙ = ITensor(eltype,dag(l[n])..., r[n]...) + ψCₙ = ITensor(eltype, dag(l[n])..., r[n]...) ψCₙ[l[n]... => 1, r[n]... => 1] = 1.0 ψC[n] = ψCₙ end diff --git a/src/vumps_mpo.jl b/src/vumps_mpo.jl index 090d1b1..f4cb2d6 100644 --- a/src/vumps_mpo.jl +++ b/src/vumps_mpo.jl @@ -111,7 +111,7 @@ function left_environment(H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS; tol=1e- EType = ITensorMPS.promote_itensor_eltype(ψ) - eₗ = zeros(EType,1) + eₗ = zeros(EType, 1) dₕ = size(H[1])[1] #Ls = [Vector{ITensor}(undef, dₕ) for j in 1:N] Ls = [initialize_left_environment(H, ψ, j; init_last=true) for j in 1:N] @@ -135,7 +135,7 @@ function left_environment(H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS; tol=1e- # as they have NO link dimension, # because they terminate and do not transport QN numbers. # Intermediate OPs proportional to the identity DO transport QNs and must be order 4! - if (!isempty(H[1][a, a])) && (order(H[1][a, a]) == 2 ) + if (!isempty(H[1][a, a])) && (order(H[1][a, a]) == 2) λ = H[1][a, a][1, 1] # δiag = δˢ(1) # not used #@assert norm(H[1][n, n] - λ * δiag) == 0 "Non identity diagonal not implemented in MPO" @@ -147,23 +147,25 @@ function left_environment(H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS; tol=1e- Ls[1][a], info = linsolve(A, Ls[1][a], 1, -1; tol=tol) else # println("Not implemented") - error("This case cannot exist! The last and first term of the diagonal of the MPO must be the identity.") + error( + "This case cannot exist! The last and first term of the diagonal of the MPO must be the identity.", + ) flush(stdout) flush(stderr) end - elseif (!isempty(H[1][a, a])) && (order(H[1][a, a]) == 4 ) + elseif (!isempty(H[1][a, a])) && (order(H[1][a, a]) == 4) # assert diagonal operator! λ = H[1][a, a][1, 1, 1, 1] - test_Tensor = reshape(diagm(λ*ones(phys_dim)), (1,phys_dim,phys_dim,1)) + test_Tensor = reshape(diagm(λ * ones(phys_dim)), (1, phys_dim, phys_dim, 1)) # verify that this OP = λ⋅Id @assert norm(test_Tensor - array(H[1][a, a])) ≈ 0.0 "Only operators proportional to the identity are allowed on the diagonal!" @assert (λ <= 1.0) && (λ >= 0.0) "The proportionality factor of the operator to the identity must be 0 ≤ λ ≤ 1" - + # solve equations with linsolve eₗ[1] = (Ls[1][a] * localR)[] solo_link = inds(Ls[1][a])[1] - Ls[1][a] += -(eₗ[1] * denseblocks(δˡ(1))*setelt(solo_link[1])) + Ls[1][a] += -(eₗ[1] * denseblocks(δˡ(1)) * setelt(solo_link[1])) A = AOᴸ(ψ, H, a) Ls[1][a], info = linsolve(A, Ls[1][a], 1, -1; tol=tol) end @@ -218,13 +220,13 @@ function initialize_right_environment( Rs[1] = ITensor(Float64, link, dag(prime(link))) end for j in 2:(dₕ - 1) - mpo_link = only(uniqueinds(H[n - 1][dₕ,j], sit)) + mpo_link = only(uniqueinds(H[n - 1][dₕ, j], sit)) Rs[j] = ITensor(Float64, dag(mpo_link), link, dag(prime(link))) end return Rs end -Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}} +Index{Vector{Pair{QN,Int64}}}, Index{Vector{Pair{QN,Int64}}} function apply_local_right_transfer_matrix( Lstart::Vector{ITensor}, H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS, n_1::Int64 @@ -289,7 +291,7 @@ function right_environment(H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS; tol=1e EType = ITensorMPS.promote_itensor_eltype(ψ) - eᵣ = zeros(EType,1) + eᵣ = zeros(EType, 1) dₕ = size(H[1])[1] Rs = [initialize_right_environment(H, ψ, j; init_first=true) for j in 1:N] #Building the R vector for n_1 = 1 @@ -312,7 +314,7 @@ function right_environment(H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS; tol=1e # as they have NO link dimension, # because they terminate and do not transport QN numbers. # Intermediate OPs proportional to the identity DO transport QNs and must be order 4! - if (!isempty(H[1][a, a])) && (order(H[1][a, a]) == 2 ) + if (!isempty(H[1][a, a])) && (order(H[1][a, a]) == 2) λ = H[1][a, a][1, 1] # δiag = δˢ(1) #@assert norm(H[1][n, n] - λ * δiag) == 0 "Non identity diagonal not implemented in MPO" @@ -323,19 +325,21 @@ function right_environment(H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS; tol=1e A = AOᴿ(ψ, H, a) Rs[1][a], info = linsolve(A, Rs[1][a], 1, -1; tol=tol) else - error("This case cannot exist! The last and first term of the diagonal of the MPO must be the identity.") + error( + "This case cannot exist! The last and first term of the diagonal of the MPO must be the identity.", + ) flush(stdout) flush(stderr) end - elseif (!isempty(H[1][a, a])) && (order(H[1][a, a]) == 4 ) + elseif (!isempty(H[1][a, a])) && (order(H[1][a, a]) == 4) # assert diagonal operator! λ = H[1][a, a][1, 1, 1, 1] - test_Tensor = reshape(diagm(λ*ones(phys_dim)), (1,phys_dim,phys_dim,1)) + test_Tensor = reshape(diagm(λ * ones(phys_dim)), (1, phys_dim, phys_dim, 1)) # verify that this OP = λ⋅Id @assert norm(test_Tensor - array(H[1][a, a])) ≈ 0.0 "Only operators proportional to the identity are allowed on the diagonal!" @assert (λ <= 1.0) && (λ >= 0.0) "The proportionality factor of the operator to the identity must be 0 ≤ λ ≤ 1" - + # solve euqations with linsolve eᵣ[1] = (localL * Rs[1][a])[] solo_link = inds(Rs[1][a])[1] @@ -423,8 +427,8 @@ function tdvp_iteration_sequential( EType = typeof(one(EType_ψ) * one(EType_t)) - eL = zeros(EType,N) - eR = zeros(EType,N) + eL = zeros(EType, N) + eR = zeros(EType, N) for n in 1:N L, eL[n] = left_environment(H, ψ; tol=_solver_tol) #TODO currently computing two many of them R, eR[n] = right_environment(H, ψ; tol=_solver_tol) #TODO currently computing two many of them @@ -489,8 +493,8 @@ function tdvp_iteration_parallel( EType_t = typeof(time_step) EType = typeof(one(EType_ψ) * one(EType_t)) - eL = zeros(EType,1) - eR = zeros(EType,1) + eL = zeros(EType, 1) + eR = zeros(EType, 1) L, eL[1] = left_environment(H, ψ; tol=_solver_tol) #TODO currently computing two many of them R, eR[1] = right_environment(H, ψ; tol=_solver_tol) #TODO currently computing two many of them diff --git a/test/test_vumps_GeneralMPO.jl b/test/test_vumps_GeneralMPO.jl index 4f1393a..34d5629 100644 --- a/test/test_vumps_GeneralMPO.jl +++ b/test/test_vumps_GeneralMPO.jl @@ -7,21 +7,22 @@ for f in src_files include(joinpath(base_path, f)) end - """ build the transverse field Ising model with exponential interactions in the thermodynamic limit 𝑁 → ∞ use convention H = -J ∑_{n,m=−∞}^∞ σˣₙ σˣₘ ⋅ λ^(-|n-m-1|) - hz ∑_{n=−∞}^∞ σᶻₙ with """ -function InfiniteExpHTFI( +function InfiniteExpHTFI( sites::CelledVector{<:Index}, λ, # exponential decay base J, # interaction kinetic (tunneling) coupling hz; # interaction strength of the transverse field - kwargs... + kwargs..., ) if abs(λ) > 1.0 - throw(ArgumentError("cannot implement exponential decay with base larger than 1, λ = $(λ)!")) + throw( + ArgumentError("cannot implement exponential decay with base larger than 1, λ = $(λ)!") + ) end link_dimension = 3 @@ -31,10 +32,13 @@ function InfiniteExpHTFI( EType = eltype(union(λ, J, hz)) linkindices = CelledVector( - hasqns(sites.data) ? - [Index([QN("SzParity",1,2) => 1], "Link,c=1,n=$n") for n in 1:N] : [Index(1, "Link,c=1,n=$(n)") for n in 1:N] + if hasqns(sites.data) + [Index([QN("SzParity", 1, 2) => 1], "Link,c=1,n=$n") for n in 1:N] + else + [Index(1, "Link,c=1,n=$(n)") for n in 1:N] + end, ) - + mpos = [Matrix{ITensor}(undef, link_dimension, link_dimension) for i in 1:N] for n in 1:N # define local matrix Hmat with empty tensors as local operators @@ -42,52 +46,53 @@ function InfiniteExpHTFI( ITensor(EType, dag(sites[n]), prime(sites[n])), link_dimension, link_dimension ) # left link index ll with daggered QN conserving direction (if applicable) - ll = dag(linkindices[n-1]) + ll = dag(linkindices[n - 1]) # right link index rl rl = linkindices[n] # add both Identities as netral elements in the MPO # replace all known tensors from empty to known interactions - Hmat[1,1] = op("Id", sites[n]) - Hmat[3,3] = op("Id", sites[n]) + Hmat[1, 1] = op("Id", sites[n]) + Hmat[3, 3] = op("Id", sites[n]) # local nearest neighbour and exp. decaying interaction terms - Hmat[2,1] = op("X", sites[n]) + Hmat[2, 1] = op("X", sites[n]) if !iszero(λ) - Hmat[2,2] = op("Id", sites[n]) * λ # λ Id, on the diagonal + Hmat[2, 2] = op("Id", sites[n]) * λ # λ Id, on the diagonal end - Hmat[3,2] = op("X", sites[n]) * -J # Jxx σˣ + Hmat[3, 2] = op("X", sites[n]) * -J # Jxx σˣ if !iszero(hz) - Hmat[3,1] = op("Z", sites[n]) * -hz # hz σᶻ + Hmat[3, 1] = op("Z", sites[n]) * -hz # hz σᶻ end # add all missing links that connect the interaction # operators in the unit cell - Hmat[2,1] = setelt(ll[1]) * Hmat[2,1] - Hmat[1,2] = Hmat[1,2] * setelt(rl[1]) - Hmat[2,2] = setelt(ll[1]) * Hmat[2,2] * setelt(rl[1]) - Hmat[3,2] = setelt(rl[1]) * Hmat[3,2] - Hmat[2,3] = setelt(ll[1]) * Hmat[2,3] + Hmat[2, 1] = setelt(ll[1]) * Hmat[2, 1] + Hmat[1, 2] = Hmat[1, 2] * setelt(rl[1]) + Hmat[2, 2] = setelt(ll[1]) * Hmat[2, 2] * setelt(rl[1]) + Hmat[3, 2] = setelt(rl[1]) * Hmat[3, 2] + Hmat[2, 3] = setelt(ll[1]) * Hmat[2, 3] mpos[n] = Hmat end return InfiniteBlockMPO(mpos, sites.translator) end - # sanity check if I can construct the MPO for the normal TFI correctly """ build the transverse field Ising model with nearest neighbour interactions in the thermodynamic limit 𝑁 → ∞ use convention H = -J ∑_{n=−∞}^∞ σˣₙ σˣₙ₊₁ - hz ∑_{n=−∞}^∞ σᶻₙ with """ -function InfiniteHTFI( +function InfiniteHTFI( sites::CelledVector{<:Index}, kinetic_coupling, hz; # interaction kinetic (tunneling) coupling - kwargs... + kwargs..., ) if abs(λ) > 1.0 - throw(ArgumentError("cannot implement exponential decay with base larger than 1, λ = $(λ)!")) + throw( + ArgumentError("cannot implement exponential decay with base larger than 1, λ = $(λ)!") + ) end link_dimension = 3 @@ -96,10 +101,13 @@ function InfiniteHTFI( EType = eltype(union(J, hz)) linkindices = CelledVector( - hasqns(sites.data) ? - [Index([QN("SzParity",1,2) => 1], "Link,c=1,n=$n") for n in 1:N] : [Index(1, "Link,c=1,n=$(n)") for n in 1:N] + if hasqns(sites.data) + [Index([QN("SzParity", 1, 2) => 1], "Link,c=1,n=$n") for n in 1:N] + else + [Index(1, "Link,c=1,n=$(n)") for n in 1:N] + end, ) - + mpos = [Matrix{ITensor}(undef, link_dimension, link_dimension) for i in 1:N] for n in 1:N # define local matrix Hmat with empty tensors as local operators @@ -107,26 +115,26 @@ function InfiniteHTFI( ITensor(EType, dag(sites[n]), prime(sites[n])), link_dimension, link_dimension ) # left link index ll with daggered QN conserving direction (if applicable) - ll = dag(linkindices[n-1]) + ll = dag(linkindices[n - 1]) # right link index rl rl = linkindices[n] # add both Identities as netral elements in the MPO # replace all known tensors from empty to known interactions - Hmat[1,1] = op("Id", sites[n]) - Hmat[3,3] = op("Id", sites[n]) + Hmat[1, 1] = op("Id", sites[n]) + Hmat[3, 3] = op("Id", sites[n]) # local nearest neighbour and exp. decaying interaction terms - Hmat[2,1] = op("X", sites[n]) - Hmat[3,2] = op("X", sites[n]) * -J # Jxx σˣ + Hmat[2, 1] = op("X", sites[n]) + Hmat[3, 2] = op("X", sites[n]) * -J # Jxx σˣ if !iszero(hz) - Hmat[3,1] = op("Z", sites[n]) * -hz # hz σᶻ + Hmat[3, 1] = op("Z", sites[n]) * -hz # hz σᶻ end # add all missing links that connect the # interaction operators in the unit cell - Hmat[2,1] = setelt(ll[1]) * Hmat[2,1] + Hmat[2, 1] = setelt(ll[1]) * Hmat[2, 1] # Hmat[1,2] = Hmat[1,2] * setelt(rl[1]) - Hmat[3,2] = setelt(rl[1]) * Hmat[3,2] + Hmat[3, 2] = setelt(rl[1]) * Hmat[3, 2] # Hmat[2,3] = setelt(ll[1]) * Hmat[2,3] mpos[n] = Hmat end @@ -139,12 +147,10 @@ function expect_two_site(ψ::InfiniteCanonicalMPS, h::ITensor, n1n2) return inner(ϕ, apply(h, ϕ)) end - function expect_two_site(ψ::InfiniteCanonicalMPS, h::MPO, n1n2) return expect_two_site(ψ, contract(h), n1n2) end - function energy_local(ψ1, ψ2, h::ITensor) ϕ = ψ1 * ψ2 return (noprime(ϕ * h) * dag(ϕ))[] @@ -156,7 +162,6 @@ function ITensorMPS.expect(ψ, o) return (noprime(ψ * op(o, filterinds(ψ, "Site")...)) * dag(ψ))[] end - maxdim = 16 # Maximum bond dimension cutoff = 1e-10 # Singular value cutoff when increasing the bond dimension max_vumps_iters = 100 # Maximum number of iterations of the VUMPS/TDVP algorithm at a fixed bond dimension @@ -165,8 +170,11 @@ outer_iters = 4 # Number of times to increase the bond dimension time_step = -Inf # -Inf corresponds to VUMPS, finite time_step corresponds to TDVP solver_tol = (x -> x / 100) # Tolerance for the local solver (eigsolve in VUMPS and exponentiate in TDVP) multisite_update_alg = "parallel" # Choose between ["sequential", "parallel"]. Only parallel works with TDVP. -conserve_qns = true # Whether or not to conserve spin parity +# conserve_qns = true # Whether or not to conserve spin parity + +# for conserve_qns in [true, false] +println("Test with Sx conservation = $(conserve_qns)") nsite = 2 # Number of sites in the unit cell initstate(n) = "↑" @@ -180,14 +188,12 @@ s = infsiteinds("S=1/2", nsite; initstate, conserve_szparity=conserve_qns) # hz = 2.0 # λ = 0.4 -J = 1 -hz = 1.0 -λ = 0.4 +J = 1 +hz = 1.1 # H_test = InfiniteHTFI(s,J,hz) -H_test = InfiniteHTFI(s,J,hz) -H_test0 = InfiniteExpHTFI(s,0.0,J,hz) -H_test_exp = InfiniteExpHTFI(s,λ,J,hz) +H_test = InfiniteHTFI(s, J, hz) +H_test0 = InfiniteExpHTFI(s, 0.0, J, hz) vumps_kwargs = ( tol=tol, @@ -200,28 +206,33 @@ subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim) H_ref = InfiniteSum{MPO}(Model("ising"), s; J=J, h=hz) +ψ_ref = vumps_subspace_expansion( + H_ref, ψ; outer_iters, subspace_expansion_kwargs, vumps_kwargs +) -ψ_ref = vumps_subspace_expansion(H_ref, ψ; outer_iters, subspace_expansion_kwargs, vumps_kwargs) - -ψ_test_NN = vumps_subspace_expansion(H_test, ψ; outer_iters, subspace_expansion_kwargs, vumps_kwargs) -ψ_test0_NN = vumps_subspace_expansion(H_test0, ψ; outer_iters, subspace_expansion_kwargs, vumps_kwargs) - +ψ_test_NN = vumps_subspace_expansion( + H_test, ψ; outer_iters, subspace_expansion_kwargs, vumps_kwargs +) +ψ_test0_NN = vumps_subspace_expansion( + H_test0, ψ; outer_iters, subspace_expansion_kwargs, vumps_kwargs +) E_ref = energy_local(ψ_ref.AL[1], ψ_ref.AL[2] * ψ_ref.C[2], H_ref[(1, 2)]) L_test, energy_test = ITensorInfiniteMPS.left_environment(H_test, ψ_test_NN; tol=1e-10); L_test0, energy_test0 = ITensorInfiniteMPS.left_environment(H_test0, ψ_test0_NN; tol=1e-10); -energy_test /= length(ψ_test_NN) +energy_test /= length(ψ_test_NN) energy_test0 /= length(ψ_test0_NN) @test energy_test ≈ E_ref @test energy_test0 ≈ E_ref -@show E_ref + 4/π -@show energy_test + 4/π -@show energy_test0 + 4/π +E_exact = reference(Model("ising"), Observable("energy"); h=hz, J) +@test isapprox(E_ref, E_exact; rtol=1e-10) +@test isapprox(energy_test, E_exact; rtol=1e-10) +@test isapprox(energy_test0, E_exact; rtol=1e-10) # energy_local(ψ_NN.AL[1], ψ_NN.AL[2] * ψ.C[2], H_test[(1, 2)]) # energy_local(ψ_test_NN,ψ_test_NN,H_test) @@ -230,29 +241,65 @@ function ITensorMPS.expect(ψ, o) return (noprime(ψ * op(o, filterinds(ψ, "Site")...)) * dag(ψ))[] end +@test isapprox( + expect(ψ_ref.AL[1] * ψ_ref.C[1], "Z"), + expect(ψ_test_NN.AL[1] * ψ_test_NN.C[1], "Z"); + atol=1e-7, +) +@test isapprox( + abs(expect(ψ_ref.AL[1] * ψ_ref.C[1], "X")), + abs(expect(ψ_test_NN.AL[1] * ψ_test_NN.C[1], "X")); + atol=1e-7, +) -@test isapprox(expect(ψ_ref.AL[1] * ψ_ref.C[1], "Z"), expect(ψ_test_NN.AL[1] * ψ_test_NN.C[1], "Z"); rtol=1e-6) -@test isapprox(expect(ψ_ref.AL[1] * ψ_ref.C[1], "X"), expect(ψ_test_NN.AL[1] * ψ_test_NN.C[1], "X"); atol=1e-6) - -@test isapprox(expect(ψ_ref.AL[1] * ψ_ref.C[1], "Z"), expect(ψ_test0_NN.AL[1] * ψ_test0_NN.C[1], "Z"); rtol=1e-6) -@test isapprox(expect(ψ_ref.AL[1] * ψ_ref.C[1], "X"), expect(ψ_test0_NN.AL[1] * ψ_test0_NN.C[1], "X"); rtol=1e-6) +@test isapprox( + expect(ψ_ref.AL[1] * ψ_ref.C[1], "Z"), + expect(ψ_test0_NN.AL[1] * ψ_test0_NN.C[1], "Z"); + atol=1e-7, +) +@test isapprox( + abs(expect(ψ_ref.AL[1] * ψ_ref.C[1], "X")), + abs(expect(ψ_test0_NN.AL[1] * ψ_test0_NN.C[1], "X")); + atol=1e-7, +) -entropy(ψ_ref,1) -entropy(ψ_test_NN,1) -entropy(ψ_test0_NN,1) +entropy(ψ_ref, 1) +entropy(ψ_test_NN, 1) +entropy(ψ_test0_NN, 1) -@test isapprox(entropy(ψ_ref,1), entropy(ψ_test_NN,1); rtol=1e-6) -@test isapprox(entropy(ψ_ref,1), entropy(ψ_test0_NN,1); rtol=1e-6) +@test isapprox(entropy(ψ_ref, 1), entropy(ψ_test_NN, 1); rtol=1e-7) +@test isapprox(entropy(ψ_ref, 1), entropy(ψ_test0_NN, 1); rtol=1e-7) ################## #### actually exponential interactions ################# -ψ_exp = vumps_subspace_expansion(H_test_exp, ψ; outer_iters, subspace_expansion_kwargs, vumps_kwargs) +λ = 0.6 +s = infsiteinds("S=1/2", nsite; initstate, conserve_szparity=false) + +vumps_kwargs = ( + tol=1e-6, + maxiter=max_vumps_iters, + solver_tol=solver_tol, + multisite_update_alg=multisite_update_alg, +) +ψ = InfMPS(s, initstate) + +H_test_exp = InfiniteExpHTFI(s, λ, J, hz) + +ψ_exp = vumps_subspace_expansion( + H_test_exp, ψ; outer_iters, subspace_expansion_kwargs, vumps_kwargs +) -L_test, energy_test_exp = ITensorInfiniteMPS.left_environment(H_test_exp, ψ_exp; tol=1e-10); +L_test, energy_test_exp_L = ITensorInfiniteMPS.left_environment(H_test_exp, ψ_exp; tol=1e-6); +R_test, energy_test_exp_R = ITensorInfiniteMPS.right_environment( + H_test_exp, ψ_exp; tol=1e-6 +); +expect(ψ_exp) -E_gs_04_MPSKit = -1.8177300191088515 # computed at χ=16, λ=0.4, Jₓₓ=1.0, hz=1.0 -E_gs_04_ITensor = energy_test_exp / length(ψ_exp) -@test E_gs_04_MPSKit ≈ E_gs_04_ITensor +E_gs_λ04_J1_h1_MPSKit = -1.8177300191088515 # computed at χ=16, λ=0.4, Jₓₓ=1.0, hz=1.0 +E_gs_λ04_J1_h11_MPSKit = -1.8497463013720623 # computed at χ=16, λ=0.4, Jₓₓ=1.0, hz=1.1 +E_gs_04_ITensor = energy_test_exp_L / length(ψ_exp) +@test E_gs_λ04_J1_h11_MPSKit ≈ E_gs_04_ITensor +# end From 0b81ed3fddff213ceebada4e7617a17c327d9814 Mon Sep 17 00:00:00 2001 From: Jan Schneider Date: Sat, 13 Sep 2025 11:23:54 +0200 Subject: [PATCH 09/15] delete orphan line --- src/vumps_mpo.jl | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/vumps_mpo.jl b/src/vumps_mpo.jl index 61abcf7..dc1d207 100644 --- a/src/vumps_mpo.jl +++ b/src/vumps_mpo.jl @@ -226,10 +226,6 @@ function initialize_right_environment( return Rs end -<<<<<<< HEAD -Index{Vector{Pair{QN,Int64}}}, Index{Vector{Pair{QN,Int64}}} -======= ->>>>>>> 21446374b8b806571aa46667a272c24595aa8d75 function apply_local_right_transfer_matrix( Lstart::Vector{ITensor}, H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS, n_1::Int64 From 22c82ce9e7a00039913421ceb233b8e5bfaaafa4 Mon Sep 17 00:00:00 2001 From: Jan Schneider Date: Sat, 13 Sep 2025 11:26:51 +0200 Subject: [PATCH 10/15] delete unwanted files --- src/vumps_mpo.jl | 1 - test/MPO-time-evo.jl | 536 ------------------------------- test/test-contraction-MPO-MPS.jl | 185 ----------- 3 files changed, 722 deletions(-) delete mode 100644 test/MPO-time-evo.jl delete mode 100644 test/test-contraction-MPO-MPS.jl diff --git a/src/vumps_mpo.jl b/src/vumps_mpo.jl index dc1d207..abc2ff8 100644 --- a/src/vumps_mpo.jl +++ b/src/vumps_mpo.jl @@ -226,7 +226,6 @@ function initialize_right_environment( return Rs end - function apply_local_right_transfer_matrix( Lstart::Vector{ITensor}, H::InfiniteBlockMPO, ψ::InfiniteCanonicalMPS, n_1::Int64 ) diff --git a/test/MPO-time-evo.jl b/test/MPO-time-evo.jl deleted file mode 100644 index 853d1a7..0000000 --- a/test/MPO-time-evo.jl +++ /dev/null @@ -1,536 +0,0 @@ -# Reference: Efficient higher-order matrix product operators for time evolution -# Maarten Van Damme, Jutho Haegeman, Ian McCulloch, Laurens Vanderstraeten -# SciPost Phys. 17, 135 (2024) · published 15 November 2024 -# https://scipost.org/SciPostPhys.17.5.135 -# -# -# Given an MPO W of a Hamiltonian H in form -# | Id | C | D | -# W = | 0 | A | B | -# | 0 | 0 | Id | - -# then we find that the time evolution opertator U(t) -# U(t) = exp(-ı t H) -# | Id + τD + τ^2/2 D^2 + τ^3/6 D^3 | C + τ^2/2 {CD} + τ^3/6 {CDD} | CC + τ/3 {CCD} -# ≈ | τB + τ^2/2 {B D} + τ^3/6 {BDD} | A + τ/2({BC} + {AD}) + τ^2/6( {CBD} + {ADD}) | {AC} + τ/3 ({ACD} + {CCB}) -# | τ^2/2 BB + τ^3/6 {BBD} | τ/2 {AB} + τ^2/6 ({ABD} + {BBC}) | AA + τ/3 ({ABC} +{AAD}) - -# NB: -# {AB} = AB + BA -# {ABB} = {BBA} = ABB + BAB + BBA -# {ABC} = ABC + ACB + BAC + BCA + CAB + CBA -# measure τ in imaginary units for actualy time evolution - -braket(A::ITensor, B::ITensor) = replaceprime(prime(A) * B + prime(B) * A, 2 => 1) -function braket(A::ITensor, B::ITensor, C::ITensor) - return replaceprime( - A'' * B' * C + A'' * C' * B + B'' * A' * C + B'' * C' * A + C'' * A' * B + C'' * B' * A, - 3 => 1, - ) -end -function braket2(A::ITensor, B::ITensor) - return replaceprime(A'' * B' * B + B'' * A' * B + B'' * B' * A, 3 => 1) -end - -function timeEvo_MPO_2ndOrder( - s::Union{CelledVector{Index{Int64}},CelledVector{Index{Vector{Pair{QN,Int64}}}}}, - AStrings::Vector{String}, - JAs::Vector{<:Number}, - BStrings::Vector{String}, - JBs::Vector{<:Number}, - CStrings::Vector{String}, - JCs::Vector{<:Number}, - DString::String, - JD::Number, - t::Number; -) - return MPO( - timeEvo_ITensors_2ndOrder( - s.data, AStrings, JAs, BStrings, JBs, CStrings, JCs, DString, JD, t; - ), - ) -end - -function timeEvo_ITensor_2ndOrder( - s::Union{CelledVector{Index{Int64}},CelledVector{Index{Vector{Pair{QN,Int64}}}}}, - # linkindices, - AStrings::Vector{String}, - JAs::Vector{<:Number}, - BStrings::Vector{String}, - JBs::Vector{<:Number}, - CStrings::Vector{String}, - JCs::Vector{<:Number}, - DString::String, - JD::Number, - dt::Number; -) - linkindices = get_linkindices_timeEvo_MPO(s, BStrings, CStrings) - - # note that - τ = -1.0im * dt - - return CelledVector( - map( - x -> bulk_timeEvo_ITensor_2ndOrder( - x[1], s, linkindices, AStrings, JAs, BStrings, JBs, CStrings, JCs, DString, JD, τ - ), - enumerate(s), - ), - translator(s), - ) -end - -function bulk_timeEvo_ITensor_2ndOrder( - n, - sites, - linkindices, - AStrings::Vector{String}, - JAs::Vector{<:Number}, - BStrings::Vector{String}, - JBs::Vector{<:Number}, - CStrings::Vector{String}, - JCs::Vector{<:Number}, - DString::String, - JD::Number, - τ::Number; -) - s_n = sites[n] - # left link index ll with daggered QN conserving direction (if applicable) - ll = (linkindices[n - 1]) - # right link index rl - rl = dag(linkindices[n]) - # Id = op(sites, "Id", n) - - NrOfTerms = length(AStrings) - - # A is possible exponential decay so test for "0" - As = map(x -> x[1] * op(s_n, x[2]), zip(JAs, AStrings)) - Bs = map(x -> x[1] * op(s_n, x[2]), zip(JBs, BStrings)) # JB * op(sites, BString, n) - Cs = map(x -> x[1] * op(s_n, x[2]), zip(JCs, CStrings)) # JC * op(sites, CString, n) - D = JD * op(s_n, DString) - # D = JD * Op(DString, n) - - # Init ITensor inside MPO - # U_t[n] = ITensor(ComplexF64, ll, dag(s), s', rl) - - # first element - firstelement = if iszero(D) - setelt(ll[1]) * (setelt(rl[1])) * op(s_n, "Id") - else - setelt(ll[1]) * - setelt(rl[1]) * - ( - op(s_n, "Id") + - τ * D + - (τ^2 / 2) * replaceprime(D' * D, 2, 1) + - (τ^3 / 6) * replaceprime(D'' * D' * D, 3, 1) - ) - end - - # first row - Cterm = mapreduce( - x -> - setelt(ll[1]) * - setelt(rl[1 + x[1]]) * - (x[2] + (τ / 2) * braket(x[2], D) + (τ^2 / 6) * braket2(x[2], D)), - +, - enumerate(Cs), - ) - - # CHECK FOR NILL-POTENT OPERATORS or operators proportional to zero - # avoid setting entries explicitly zero because that counters the purpose of sparse matrices and - # reduces the runtime efficiency - # (this applies only to Sparse Matrices in the case of QN conservation) - # Note that the empty ITensor is not equal to the ITensor with only zero entries! - Cterm_checked = iszero(Cterm) ? emptyITensor(ll, rl, s_n', dag(s_n)) : Cterm - - Csquared_term = mapreduce( - x -> - setelt(ll[1]) * - setelt(rl[1 + x[1] + NrOfTerms]) * - (replaceprime(x[2]' * x[2], 2, 1) + (τ / 3) * braket2(D, x[2])), - +, - enumerate(Cs), - ) - Csquared_term_checked = - iszero(Csquared_term) ? emptyITensor(ll, rl, s_n', dag(s_n)) : Csquared_term - - # first column (exept first row) - Bterm = mapreduce( - x -> - setelt(ll[1 + x[1]]) * - setelt(rl[1]) * - (τ * x[2] + (τ^2 / 2) * braket(x[2], D) + (τ^3 / 6) * braket2(x[2], D)), - +, - enumerate(Bs), - ) - Bterm_checked = iszero(Bterm) ? emptyITensor(ll, rl, s_n', dag(s_n)) : Bterm - - Bsquared_term = mapreduce( - x -> - setelt(ll[1 + x[1] + NrOfTerms]) * - setelt(rl[1]) * - ((τ^2 / 2) * replaceprime(x[2]' * x[2], 2, 1) + (τ^3 / 6) * braket2(D, x[2])), - +, - enumerate(Bs), - ) - Bsquared_term_checked = - iszero(Bsquared_term) ? emptyITensor(ll, rl, s_n', dag(s_n)) : Bsquared_term - - # diagonal - diagterm = mapreduce( - x -> - setelt(ll[1 + x[1]]) * - setelt(rl[1 + x[1]]) * - ( - x[2][1] + - (τ / 2) * (braket(x[2][2], x[2][3]) + braket(x[2][1], D)) + - (τ^2 / 6) * (braket(x[2][3], x[2][2], D) + braket2(x[2][1], D)) - ), - +, - enumerate(zip(As, Bs, Cs)), - ) - diagterm_checked = iszero(diagterm) ? emptyITensor(ll, rl, s_n', dag(s_n)) : diagterm - - diagsquaredterm = mapreduce( - x -> - setelt(ll[1 + x[1] + NrOfTerms]) * - setelt(rl[1 + x[1] + NrOfTerms]) * - ( - replaceprime(x[2][1]' * x[2][1], 2, 1) + - (τ / 3) * (braket(x[2][1], x[2][2], x[2][3]) + braket2(D, x[2][1])) - ), - +, - enumerate(zip(As, Bs, Cs)), - ) - diagsquaredterm_checked = - iszero(diagsquaredterm) ? emptyITensor(ll, rl, s_n', dag(s_n)) : diagsquaredterm - - # the "rest" of mixed terms - mixedterm = mapreduce( - x -> - setelt(ll[1 + x[1]]) * - setelt(rl[1 + x[1] + NrOfTerms]) * - ( - braket(x[2][1], x[2][3]) + - (τ / 3) * (braket(x[2][1], x[2][3], D) + braket2(x[2][2], x[2][3])) - ), - +, - enumerate(zip(As, Bs, Cs)), - ) - mixedterm_checked = iszero(mixedterm) ? emptyITensor(ll, rl, s_n', dag(s_n)) : mixedterm - - mixedsquareterm = mapreduce( - x -> - setelt(ll[1 + x[1] + NrOfTerms]) * - (setelt(rl[1 + x[1]])) * - ( - (τ / 2) * braket(x[2][1], x[2][2]) + - (τ^2 / 6) * (braket(x[2][1], x[2][2], D) + braket2(x[2][3], x[2][2])) - ), - +, - enumerate(zip(As, Bs, Cs)), - ) - mixedsquareterm_checked = - iszero(mixedsquareterm) ? emptyITensor(ll, rl, s_n', dag(s_n)) : mixedsquareterm - - return sum(( - firstelement, - Cterm_checked, - Csquared_term_checked, - Bterm_checked, - Bsquared_term_checked, - diagterm_checked, - diagsquaredterm_checked, - mixedterm_checked, - mixedsquareterm_checked, - )) -end - -function Left_timeEvo_ITensor_2ndOrder( - sites, - linkindices, - CStrings::Vector{String}, - JCs::Vector{<:Number}, - DString::String, - JD::Number, - τ::Number; -) - NrOfTerms = length(CStrings) - # right link index rl - rl = linkindices[1] - n = 1 - - Cs = map(x -> x[1] * op(sites, x[2], 1), zip(JCs, CStrings)) # JC * op(sites, CString, n) - D = JD * op(sites, DString, 1) - - # first element - firstelement = if iszero(D) - setelt(rl[1]) * op(sites, "Id", n) - else - setelt(rl[1]) * ( - op(sites, "Id", n) + - τ * D + - (τ^2 / 2) * replaceprime(D' * D, 2, 1) + - (τ^3 / 6) * replaceprime(D'' * D' * D, 3, 1) - ) - end - - # first row - Cterm = mapreduce( - x -> - setelt(rl[1 + x[1]]) * - (x[2] + (τ / 2) * braket(x[2], D) + (τ^2 / 6) * braket2(x[2], D)), - +, - enumerate(Cs), - ) - Csquared_term = mapreduce( - x -> - setelt(rl[1 + x[1] + NrOfTerms]) * - (replaceprime(x[2]' * x[2], 2, 1) + (τ / 3) * braket2(D, x[2])), - +, - enumerate(Cs), - ) - - return sum( - # filter( - # !(iszero), - [firstelement, Cterm, Csquared_term] - # ) -) -end - -function Right_timeEvo_ITensor_2ndOrder( - sites, - linkindices, - BStrings::Vector{String}, - JBs::Vector{<:Number}, - DString::String, - JD::Number, - τ::Number; -) - L = length(sites) - n = L - # left link index ll with daggered QN conserving direction (if applicable) - ll = dag(linkindices[L - 1]) - - NrOfTerms = length(BStrings) - - # A is possible exponential decay so test for "0" - # As = map(x -> x[1] * op(sites, x[2], n), zip(JAs, AStrings)) - Bs = map(x -> x[1] * op(sites, x[2], n), zip(JBs, BStrings)) # JB * op(sites, BString, n) - # Cs = map(x -> x[1] * op(sites, x[2], n), zip(JCs, CStrings)) # JC * op(sites, CString, n) - D = JD * op(sites, DString, n) - - # first element - firstelement = if iszero(D) - setelt(ll[1]) * op(sites, "Id", n) - else - setelt(ll[1]) * ( - op(sites, "Id", n) + - τ * D + - (τ^2 / 2) * replaceprime(D' * D, 2, 1) + - (τ^3 / 6) * replaceprime(D'' * D' * D, 3, 1) - ) - end - - # first column (exept first row) - Bterm = mapreduce( - x -> - setelt(ll[1 + x[1]]) * - (τ * x[2] + (τ^2 / 2) * braket(x[2], D) + (τ^3 / 6) * braket2(x[2], D)), - +, - enumerate(Bs), - ) - Bsquared_term = mapreduce( - x -> - setelt(ll[1 + x[1] + NrOfTerms]) * - ((τ^2 / 2) * replaceprime(x[2]' * x[2], 2, 1) + (τ^3 / 6) * braket2(D, x[2])), - +, - enumerate(Bs), - ) - - return sum( - # filter( - # !(iszero), - [firstelement, Bterm, Bsquared_term] - # ) -) -end - -function get_linkindices_timeEvo_MPO( - sites::Vector{Index{Vector{Pair{QN,Int64}}}}, - BStrings::Vector{String}, - CStrings::Vector{String}, -) - N = length(sites) - NOps = length(BStrings) - if NOps != length(CStrings) - throw(ArgumentError("Input BStrings and CString do not share same length!")) - end - - # save QN flux of each operator, note that multiple operators may have same flux, - # thus beloning to the same block and increasing the local dimension of this block - # reuse the vector to save memory - QNFlux_vector = Vector{QN}(undef, length(BStrings)) - - linkindices = Vector{Index{Vector{Pair{QN,Int64}}}}(undef, N) - - for n in 1:N - # save the flux and the corresponding dimension in a dynamically sized vector as it is ordered in historical order - # same as the corresponding operators - QN_local_index_dim = Vector{Pair{QN,Int64}}(undef, 1) - nameQN = String(qn(sites[n][1]).data[1].name) - QNmodulus = qn(sites[n][1]).data[1].modulus - # loop over all interaction operators - for (indexOP, (BString, CString)) in enumerate(zip(BStrings, CStrings)) - # get the flux of interaction - QNFlux_vector[indexOP] = flux(op(sites[n], BString)) - checkflux = flux(op(sites[n], CString)) - if !(checkflux == -QNFlux_vector[indexOP]) - error( - "Operators B and C are not conserving the total QN in the system as their flux is not opposite", - ) - end - # local dimension of flux is at least 1 - if indexOP == 1 - QN_local_index_dim[indexOP] = Pair(QNFlux_vector[indexOP], 1) - elseif indexOP > 1 && QNFlux_vector[indexOP - 1] == QNFlux_vector[indexOP] - # increase tuple of number and dimension by one in dimension - QN_local_index_dim[end] = Pair( - QN_local_index_dim[end][1], QN_local_index_dim[end][2] + 1 - ) - elseif indexOP > 1 - # if new flux is encountered, save in vector - push!(QN_local_index_dim, Pair(QNFlux_vector[indexOP], 1)) - end - end # loop over interaction operators - - # second part of the entries are generated by the double of the fluxes from before - # generated by the square of the interaction operators... - # Note this is because we are using the second order approximation! - - QN_local_index_dim_final = Vector{Pair{QN,Int64}}(undef, 2 * length(QN_local_index_dim)) - - for (index, element) in enumerate(QN_local_index_dim) - dim_of_flux = element[2] - flux_val = val(element[1], nameQN) - QN_local_index_dim_final[index] = element - QN_local_index_dim_final[index + length(QN_local_index_dim)] = Pair( - QN(nameQN, 2 * flux_val, QNmodulus), dim_of_flux - ) - end - - # construct the local link index from all flux blocks and their dimension - linkindices[n] = Index( - [QN(nameQN, 0, QNmodulus) => 1, QN_local_index_dim_final...], "Link,c=1,l=$(n)" - ) - end - - return linkindices -end - -function get_linkindices_timeEvo_MPO( - sites::Vector{Index{Int64}}, BStrings::Vector{String}, CStrings::Vector{String} -) - - # L = - NOps = length(BStrings) - if NOps != length(CStrings) - throw(ArgumentError("Input BStrings and CString do not share same length!")) - end - - return [Index(1 + 2 * length(BStrings), "Link,l=$(n)") for n in 1:(length(sites) - 1)] -end - -function get_linkindices_timeEvo_MPO( - sites::CelledVector{Index{Int64}}, BStrings::Vector{String}, CStrings::Vector{String} -) - - # L = - NOps = length(BStrings) - if NOps != length(CStrings) - throw(ArgumentError("Input BStrings and CString do not share same length!")) - end - - return CelledVector( - [Index(1 + 2 * length(BStrings), "Link,c=1,l=$(n)") for n in 1:length(sites)], - translator(s), - ) -end - -function get_linkindices_timeEvo_MPO( - sites::CelledVector{Index{Vector{Pair{QN,Int64}}}}, - BStrings::Vector{String}, - CStrings::Vector{String}, -) - linkindices = get_linkindices_timeEvo_MPO(sites.data, BStrings, CStrings) - return CelledVector(linkindices, translator(sites)) -end - -function timeEvo_ITensors_2ndOrder( - sites, - AStrings::Vector{String}, - JAs::Vector{<:Number}, - BStrings::Vector{String}, - JBs::Vector{<:Number}, - CStrings::Vector{String}, - JCs::Vector{<:Number}, - DString::String, - JD::Number, - t::Number; -) - # Reference: Efficient higher-order matrix product operators for time evolution - # Maarten Van Damme, Jutho Haegeman, Ian McCulloch, Laurens Vanderstraeten - # SciPost Phys. 17, 135 (2024) · published 15 November 2024 - # https://scipost.org/SciPostPhys.17.5.135 - # Given an MPO W of a Hamiltonian H in form - # | Id | C | D | - # W = | 0 | A | B | - # | 0 | 0 | Id | - - # then we find that the time evolution opertator U(t) - # U(t) = exp(-ı t H) - # | Id + τD + τ^2/2 D^2 + τ^3/6 D^3 | C + τ^2/2 {CD} + τ^3/6 {CDD} | CC + τ/3 {CCD} - # ≈ | τB + τ^2/2 {B D} + τ^3/6 {BDD} | A + τ/2({BC} + {AD}) + τ^2/6( {CBD} + {ADD}) | {AC} + τ/3 ({ACD} + {CCB}) - # | τ^2/2 BB + τ^3/6 {BBD} | τ/2 {AB} + τ^2/6 ({ABD} + {BBC}) | AA + τ/3 ({ABC} +{AAD}) - - # NB: - # {AB} = AB + BA - # {ABB} = {BBA} = ABB + BAB + BBA - # {ABC} = ABC + ACB + BAC + BCA + CAB + CBA - # measure τ in imaginary units for actualy time evolution - τ = -1.0im * t - N = length(sites) - - if length(AStrings) != length(BStrings) || length(BStrings) != length(CStrings) - error( - "Input length of operator vectors is not the same!\n - Found length(AStrings)=$(length(AStrings)), length(BStrings)=$(length(BStrings)), length(CStrings)=$(length(CStrings))", - ) - end - if length(AStrings) != length(JAs) || - length(BStrings) != length(JBs) || - length(CStrings) != length(JCs) - error("Input length of couppling vectors is not the same as the operator vector!") - end - - linkindices = get_linkindices_timeEvo_MPO(sites, BStrings, CStrings) - - ## LEFT boundary - Ut_1 = Left_timeEvo_ITensor_2ndOrder(sites, linkindices, CStrings, JCs, DString, JD, τ;) - - ## BULK - # loop over BULK real space - bulk = map( - n -> bulk_timeEvo_ITensor_2ndOrder( - n, sites, linkindices, AStrings, JAs, BStrings, JBs, CStrings, JCs, DString, JD, τ; - ), - 2:(N - 1), - ) - - ## RIGHT boundary - Ut_N = Right_timeEvo_ITensor_2ndOrder(sites, linkindices, BStrings, JBs, DString, JD, τ;) - - return [Ut_1, bulk..., Ut_N] -end diff --git a/test/test-contraction-MPO-MPS.jl b/test/test-contraction-MPO-MPS.jl deleted file mode 100644 index edf9f65..0000000 --- a/test/test-contraction-MPO-MPS.jl +++ /dev/null @@ -1,185 +0,0 @@ -using ITensors, ITensorMPS -using ITensorInfiniteMPS - -includet("../examples/vumps/src/entropy.jl") -includet("../examples/vumps/src/vumps_subspace_expansion.jl") -includet("MPO-time-evo.jl") - -n_sites = 2 -initstate(n) = "↑" - -s = infsiteinds("S=1/2", n_sites; initstate, conserve_szparity=false) -# sqn = infsiteinds("S=1/2", n_sites; initstate, conserve_szparity=true) - -ψ = InfMPS(s, initstate) -# ψqn = InfMPS(sqn, initstate) - -Szs_infinite = [expect(ψ, "Z", n) for n in 1:length(s)] - -H = InfiniteSum{MPO}(Model("ising"), s; J=1.0, h=1.0) -H[1] - -dt = 0.05 - -AStrings = ["Id"] -BStrings = ["X"] -CStrings = ["X"] -DString = "Z" - -JAs = [0.6] -JBs = [1.0] -JCs = [-1.0] -JD = 0.0 - -T_Z = InfiniteMPO(CelledVector([op(s[1], "Z"), op(s[2], "Z")], translator(s))) - -# T_Z * ψ - -linkindices = get_linkindices_timeEvo_MPO( - infsiteinds("S=1/2", n_sites; initstate, conserve_szparity=true), BStrings, CStrings -) - -U_t = timeEvo_ITensor_2ndOrder( - s, AStrings, JAs, BStrings, JBs, CStrings, JCs, DString, JD, dt -) -U_tqn = timeEvo_ITensor_2ndOrder( - sqn, AStrings, JAs, BStrings, JBs, CStrings, JCs, DString, JD, dt -) - -dims(U_t[1]) -dims(ψ.AL[1]) - -dims(U_t[2]) -dims(ψ.AR[2]) - -apply(U_t[2], ψ.AR[2]) -apply(U_t[1], ψ.AL[1]) - -dims(U_tqn[1]) -dims(ψqn.AL[1]) - -dims(U_tqn[2]) -dims(ψqn.AR[2]) - -apply(U_tqn[2], ψqn.AR[2]) -apply(U_tqn[1], ψqn.AL[1]) - -apply(U_t[1], ψ.AL[1]) * ψ.C[1] * apply(U_t[2], ψ.AR[2]) - -apply(U_tqn[1], ψqn.AL[1]) * ψqn.C[1] * apply(U_tqn[2], ψqn.AR[2]) - -ψ.AL[1] * ψ.C[1] * ψ.AR[2] - -H_test_exp = InfiniteExpHTFI(s, λ, J, hz) - -################# -# Time evo -################# - -dt = 0.05 -tf = 5.0 - -timeline = dt:dt:tf - -ψ_tdvp = tdvp_subspace_expansion( - H_test_exp, - InfMPS(ComplexF64, s, initstate); - time_step=-im * dt, - outer_iters=length(timeline), - subspace_expansion_kwargs, - vumps_kwargs=( - tol=1e-10, maxiter=1, solver_tol=solver_tol, multisite_update_alg=multisite_update_alg - ), -) - -includet("MPO-time-evo.jl") - -As = ["Id"] -JAs = [λ] -Bs = ["X"] -JBs = [1.0] -Cs = ["X"] -JCs = [-J] -D = "Z" - -U_dt = timeEvo_ITensor_2ndOrder(s, As, JAs, Bs, JBs, Cs, JCs, D, hz, dt) - -ψ_trotter = InfMPS(ComplexF64, s, initstate) - -# for t_now in enumerate(timeline) -ACs = [ψ_trotter.C[1] * ψ_trotter.AL[1], ψ_trotter.AR[2] * dag(ψ_trotter.C[2])] -linkAC_l = [only(inds(ACs[1], "Link,c=0")), only(inds(ACs[2], "Link,c=1,l=1"))] -linkAC_r = [only(inds(ACs[1], "Link,c=1")), only(inds(ACs[2], "Link,c=1,l=2"))] - -linkU_l = [only(inds(U_dt[1], "Link,c=0")), only(inds(U_dt[2], "Link,c=1,l=1"))] -linkU_r = [only(inds(U_dt[1], "Link,c=1")), only(inds(U_dt[2], "Link,c=1,l=2"))] - -combiners_l = [ - combiner(linkAC_l[1], linkU_l[1]; tags="Link,c=0,l=2"),# dir=only(unique(dir.([linkAC_l[1], linkU_l[1]])))), - combiner(linkAC_l[2], linkU_l[2]; tags="Link,c=1,l=1"),# dir=only(unique(dir.([linkAC_l[2], linkU_l[2]])))), -] - -combiners_r = [ - # combiner(linkAC_r[1], linkU_r[1], tags="Link,c=1,l=1"), - dag(combiners_l[2]), - combiner(linkAC_r[2], linkU_r[2]; tags="Link,c=1,l=2"),# dir=only(unique(dir.([linkAC_r[2], linkU_r[2]])))), -] - -U_dt_ψC = map( - x -> combiners_l[x[1]] * apply(x[2][1], x[2][2]) * combiners_r[x[1]], - enumerate(zip(U_dt, ACs)), -) - -links_new_l = [only(inds(U_dt_ψC[1], "Link,c=0")), only(inds(U_dt_ψC[2], "Link,c=1,l=1"))] -links_new_r = [only(inds(U_dt_ψC[1], "Link,c=1")), only(inds(U_dt_ψC[2], "Link,c=1,l=2"))] - -UR1, SR1, VR1 = svd( - U_dt_ψC[1], - links_new_l[1]; - cutoff, - maxdim, - righttags=tags(links_new_l[1]), - lefttags=tags(links_new_r[1]), -) -UR2, SR2, VR2 = svd( - U_dt_ψC[2], - links_new_l[2]; - cutoff, - maxdim, - righttags=tags(links_new_l[2]), - lefttags=tags(links_new_r[2]), -) -UL1, SL1, VL1 = svd( - U_dt_ψC[1], - (links_new_l[1], only(inds(U_dt_ψC[1], "Site"))); - cutoff, - maxdim, - righttags=tags(links_new_l[1]), - lefttags=tags(links_new_r[1]), -) -UL2, SL2, VL2 = svd( - U_dt_ψC[2], - (links_new_l[2], only(inds(U_dt_ψC[2], "Site"))); - cutoff, - maxdim, - righttags=tags(links_new_l[2]), - lefttags=tags(links_new_r[2]), -) - -inds(SR1) -inds(SR2) -inds(SL1) -inds(SL2) - -VR1 -VR2 - -UL1 -UL2 - -Cs = [] -new_state = CelledVector([UL1, UL2], translator(s)) - -# end - -# ψ_tdvp[:] = tdvp(H_test_exp, ψ_tdvp; time_step=dt, vumps_kwargs...) From 844b9ca23f9b17a0588dc58b1c8216bcb8e93105 Mon Sep 17 00:00:00 2001 From: Jan Schneider Date: Tue, 30 Sep 2025 14:59:10 +0200 Subject: [PATCH 11/15] add trivial linkdim extension to product state --- src/infinitecanonicalmps.jl | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/src/infinitecanonicalmps.jl b/src/infinitecanonicalmps.jl index af974f8..6ee2dd2 100644 --- a/src/infinitecanonicalmps.jl +++ b/src/infinitecanonicalmps.jl @@ -191,7 +191,7 @@ function zero_qn(i::Index) return zero(qn(first(space(i)))) end -function insert_linkinds!(A; left_dir=ITensors.Out) +function insert_linkinds!(A; left_dir=ITensors.Out, extended_linkdim::Int=1) # TODO: use `celllength` here N = nsites(A) l = CelledVector{indtype(A)}(undef, N, translator(A)) @@ -200,19 +200,19 @@ function insert_linkinds!(A; left_dir=ITensors.Out) dim = if hasqns(s) kwargs = (; dir=left_dir) qn_ln = zero_qn(s) - [qn_ln => 1] #Default to 0 on the right + [qn_ln => extended_linkdim] #Default to 0 on the right else kwargs = () - 1 + extended_linkdim end l[N] = Index(dim, default_link_tags("l", n, 1); kwargs...) for n in 1:(N - 1) # TODO: is this correct? dim = if hasqns(s) qn_ln = flux(A[n]) * left_dir + qn_ln#Fixed a bug on flux conservation - [qn_ln => 1] + [qn_ln => extended_linkdim] else - 1 + extended_linkdim end l[n] = Index(dim, default_link_tags("l", n, 1); kwargs...) end @@ -226,7 +226,7 @@ function insert_linkinds!(A; left_dir=ITensors.Out) end function UniformMPS( - eltype::Type{<:Number}, s::CelledVector, f::Function; left_dir=ITensors.Out + eltype::Type{<:Number}, s::CelledVector, f::Function; left_dir=ITensors.Out, extended_linkdim::Int=1 ) sᶜ¹ = s[Cell(1)] A = InfiniteMPS([ITensor(eltype, sⁿ) for sⁿ in sᶜ¹], translator(s)) @@ -237,17 +237,17 @@ function UniformMPS( Aⁿ[indval(s[n] => f(n))] = 1.0 A[n] = Aⁿ end - insert_linkinds!(A; left_dir=left_dir) + insert_linkinds!(A; left_dir=left_dir, extended_linkdim=extended_linkdim) return A end -InfMPS(s::CelledVector, f::Function) = InfMPS(Float64, s::CelledVector, f::Function) +InfMPS(s::CelledVector, f::Function; kwargs) = InfMPS(Float64, s::CelledVector, f::Function; kwargs...) -function InfMPS(eltype::Type{<:Number}, s::CelledVector, f::Function) +function InfMPS(eltype::Type{<:Number}, s::CelledVector, f::Function; extended_linkdim::Int=1) # TODO: rename cell_length N = length(s) - ψL = UniformMPS(eltype, s, f; left_dir=ITensors.Out) - ψR = UniformMPS(eltype, s, f; left_dir=ITensors.In) + ψL = UniformMPS(eltype, s, f; left_dir=ITensors.Out, extended_linkdim) + ψR = UniformMPS(eltype, s, f; left_dir=ITensors.In, extended_linkdim) ψC = InfiniteMPS(N, translator(s)) l = linkinds(ψL) r = linkinds(ψR) From abfdd992661f9599d842395d1886327200f1062c Mon Sep 17 00:00:00 2001 From: Jan Schneider Date: Tue, 30 Sep 2025 15:05:39 +0200 Subject: [PATCH 12/15] forgot kwargs... --- src/infinitecanonicalmps.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/infinitecanonicalmps.jl b/src/infinitecanonicalmps.jl index 6ee2dd2..b3f5196 100644 --- a/src/infinitecanonicalmps.jl +++ b/src/infinitecanonicalmps.jl @@ -172,13 +172,13 @@ function ITensorMPS.linkinds(ψ::InfiniteMPS) end function InfMPS( - eltype::Type{<:Number}, s::Vector, f::Function, translator::Function=translatecelltags + eltype::Type{<:Number}, s::Vector, f::Function, translator::Function=translatecelltags; kwargs... ) - return InfMPS(eltype, infsiteinds(s, translator), f) + return InfMPS(eltype, infsiteinds(s, translator), f; kwargs...) end -function InfMPS(s::Vector, f::Function, translator::Function=translatecelltags) - return InfMPS(Float64, infsiteinds(s, translator), f) +function InfMPS(s::Vector, f::Function, translator::Function=translatecelltags; kwargs...) + return InfMPS(Float64, infsiteinds(s, translator), f; kwargs...) end function indval(iv::Pair) @@ -241,7 +241,7 @@ function UniformMPS( return A end -InfMPS(s::CelledVector, f::Function; kwargs) = InfMPS(Float64, s::CelledVector, f::Function; kwargs...) +InfMPS(s::CelledVector, f::Function; kwargs...) = InfMPS(Float64, s::CelledVector, f::Function; kwargs...) function InfMPS(eltype::Type{<:Number}, s::CelledVector, f::Function; extended_linkdim::Int=1) # TODO: rename cell_length From 597d540b294a83e6019b74140dc0a59c735dd750 Mon Sep 17 00:00:00 2001 From: Jan Schneider Date: Mon, 3 Nov 2025 15:04:35 +0100 Subject: [PATCH 13/15] up krylovkit --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 2357c4b..25ef896 100644 --- a/Project.toml +++ b/Project.toml @@ -23,7 +23,7 @@ ITensorMPS = "0.3" ITensors = "0.7, 0.8, 0.9" Infinities = "0.1" IterTools = "1" -KrylovKit = "0.5, 0.6, 0.7, 0.8, 0.9" +KrylovKit = "0.5, 0.6, 0.7, 0.8, 0.9, 0.10" OffsetArrays = "1" QuadGK = "2" SplitApplyCombine = "1.2.2" From a14f73bc72942193cef0ccb2e9d024c2f789bbf1 Mon Sep 17 00:00:00 2001 From: Jan Schneider Date: Thu, 27 Nov 2025 16:08:08 +0100 Subject: [PATCH 14/15] add kwarg to control imag part error --- src/orthogonalize.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/orthogonalize.jl b/src/orthogonalize.jl index e88cfef..5bde1ee 100644 --- a/src/orthogonalize.jl +++ b/src/orthogonalize.jl @@ -6,6 +6,7 @@ function right_orthogonalize( left_tags=ts"Left", right_tags=ts"Right", tol::Real=1e-12, + tol_imag::Real = 1e-15, eager=true, ishermitian_kwargs=(; rtol=tol * 100), ) @@ -32,7 +33,7 @@ function right_orthogonalize( @warn("Non-unique largest eigenvector of transfer matrix found") end - if imag(λ₁ᴿᴺ) / norm(λ₁ᴿᴺ) > 1e-15 + if imag(λ₁ᴿᴺ) / norm(λ₁ᴿᴺ) > tol_imag @show λ₁ᴿᴺ error( "Imaginary part of eigenvalue is large: imag(λ₁ᴿᴺ) / norm(λ₁ᴿᴺ) = $(imag(λ₁ᴿᴺ) / norm(λ₁ᴿᴺ))", @@ -45,9 +46,9 @@ function right_orthogonalize( @show norm(v₁ᴿᴺ - swapinds(dag(v₁ᴿᴺ), reverse(Pair(inds(v₁ᴿᴺ)...)))) @warn("v₁ᴿᴺ is not hermitian, passed kwargs: $ishermitian_kwargs") end - if norm(imag(v₁ᴿᴺ)) / norm(v₁ᴿᴺ) > 1e-13 + if norm(imag(v₁ᴿᴺ)) / norm(v₁ᴿᴺ) > tol_imag println( - "Norm of the imaginary part $(norm(imag(v₁ᴿᴺ))) is larger than the tolerance value 1e-15. Keeping as complex.", + "Norm of the imaginary part $(norm(imag(v₁ᴿᴺ))) is larger than the tolerance value $(tol_imag). Keeping as complex.", ) @show norm(v₁ᴿᴺ - swapinds(dag(v₁ᴿᴺ), reverse(Pair(inds(v₁ᴿᴺ)...)))) else From 19428d5010227141be8360eda6d6f5fc2490f9d6 Mon Sep 17 00:00:00 2001 From: Jan Schneider Date: Thu, 27 Nov 2025 16:09:32 +0100 Subject: [PATCH 15/15] add new kwargs to other functions --- src/orthogonalize.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/orthogonalize.jl b/src/orthogonalize.jl index 5bde1ee..8750623 100644 --- a/src/orthogonalize.jl +++ b/src/orthogonalize.jl @@ -102,10 +102,10 @@ function right_orthogonalize_polar( end function left_orthogonalize( - ψ::InfiniteMPS; left_tags=ts"Left", right_tags=ts"Right", tol::Real=1e-12 + ψ::InfiniteMPS; left_tags=ts"Left", right_tags=ts"Right", tol::Real=1e-12, tol_imag::Real = 1e-15, ) Cᴸ, ψᴸ, λᴸ = right_orthogonalize( - reverse(ψ); left_tags=right_tags, right_tags=left_tags, tol=tol + reverse(ψ); left_tags=right_tags, right_tags=left_tags, tol=tol, tol_imag=tol_imag, ) # Cᴸ has the unit cell shifted from what is expected Cᴸ = reverse(Cᴸ) @@ -119,10 +119,10 @@ end # TODO: rename to `orthogonalize(ψ)`? With no limit specified, it is like orthogonalizing to over point. # Alternatively, it could be called as `orthogonalize(ψ, :)` function mixed_canonical( - ψ::InfiniteMPS; left_tags=ts"Left", right_tags=ts"Right", tol::Real=1e-12 + ψ::InfiniteMPS; left_tags=ts"Left", right_tags=ts"Right", tol::Real=1e-12, tol_imag::Real = 1e-15, ) - _, ψᴿ, _ = right_orthogonalize(ψ; left_tags=ts"", right_tags) - ψᴸ, C, λ = left_orthogonalize(ψᴿ; left_tags, right_tags) + _, ψᴿ, _ = right_orthogonalize(ψ; left_tags=ts"", right_tags, tol, tol_imag) + ψᴸ, C, λ = left_orthogonalize(ψᴿ; left_tags, right_tags, tol, tol_imag) if λ ≉ one(λ) error("λ should be approximately 1 after orthogonalization, instead it is $λ") end