From c567950df62948aa20e57e368d07260237b20691 Mon Sep 17 00:00:00 2001 From: afeuerpfeil Date: Wed, 5 Nov 2025 19:23:32 -0500 Subject: [PATCH 1/7] refactor IDMRG by following the IterativeSolver interface and by sharing code between IDMRG and IDMRG2 --- src/algorithms/groundstate/idmrg.jl | 346 ++++++++++++++-------------- 1 file changed, 172 insertions(+), 174 deletions(-) diff --git a/src/algorithms/groundstate/idmrg.jl b/src/algorithms/groundstate/idmrg.jl index 693fbb790..1163fc86c 100644 --- a/src/algorithms/groundstate/idmrg.jl +++ b/src/algorithms/groundstate/idmrg.jl @@ -24,81 +24,6 @@ $(TYPEDFIELDS) alg_eigsolve::A = Defaults.alg_eigsolve() end -function find_groundstate(ost::InfiniteMPS, H, alg::IDMRG, envs = environments(ost, H)) - ϵ::Float64 = calc_galerkin(ost, H, ost, envs) - ψ = copy(ost) - log = IterLog("IDMRG") - local iter, E_current - - LoggingExtras.withlevel(; alg.verbosity) do - @infov 2 begin - E_current = expectation_value(ψ, H, envs) - loginit!(log, ϵ, E_current) - end - for outer iter in 1:(alg.maxiter) - alg_eigsolve = updatetol(alg.alg_eigsolve, iter, ϵ) - C_current = ψ.C[0] - - # left to right sweep - for pos in 1:length(ψ) - h = AC_hamiltonian(pos, ψ, H, ψ, envs) - _, ψ.AC[pos] = fixedpoint(h, ψ.AC[pos], :SR, alg_eigsolve) - if pos == length(ψ) - # AC needed in next sweep - ψ.AL[pos], ψ.C[pos] = left_orth(ψ.AC[pos]) - else - ψ.AL[pos], ψ.C[pos] = left_orth!(ψ.AC[pos]) - end - transfer_leftenv!(envs, ψ, H, ψ, pos + 1) - end - - # right to left sweep - for pos in length(ψ):-1:1 - h = AC_hamiltonian(pos, ψ, H, ψ, envs) - _, ψ.AC[pos] = fixedpoint(h, ψ.AC[pos], :SR, alg_eigsolve) - - ψ.C[pos - 1], temp = right_orth!(_transpose_tail(ψ.AC[pos]; copy = (pos == 1))) - ψ.AR[pos] = _transpose_front(temp) - - transfer_rightenv!(envs, ψ, H, ψ, pos - 1) - end - - ϵ = norm(C_current - ψ.C[0]) - - if ϵ < alg.tol - @infov 2 begin - E_next = expectation_value(ψ, H, envs) - ΔE = E_next - E_current - E_current = E_next - logfinish!(log, iter, ϵ, ΔE) - end - break - end - if iter == alg.maxiter - @warnv 1 begin - E_next = expectation_value(ψ, H, envs) - ΔE = E_next - E_current - E_current = E_next - logcancel!(log, iter, ϵ, ΔE) - end - else - @infov 3 begin - E_next = expectation_value(ψ, H, envs) - ΔE = E_next - E_current - E_current = E_next - logiter!(log, iter, ϵ, ΔE) - end - end - end - end - - alg_gauge = updatetol(alg.alg_gauge, iter, ϵ) - ψ′ = InfiniteMPS(ψ.AR[1:end]; alg_gauge.tol, alg_gauge.maxiter) - recalculate!(envs, ψ′, H, ψ′) - - return ψ′, envs, ϵ -end - """ $(TYPEDEF) @@ -131,121 +56,194 @@ $(TYPEDFIELDS) trscheme::TruncationStrategy end -function find_groundstate(ost::InfiniteMPS, H, alg::IDMRG2, envs = environments(ost, H)) - length(ost) < 2 && throw(ArgumentError("unit cell should be >= 2")) - ϵ::Float64 = calc_galerkin(ost, H, ost, envs) - ψ = copy(ost) - log = IterLog("IDMRG2") - local iter +# Internal state of the IDMRG algorithm +struct IDMRGState{S, O, E, T <: Number} + mps::S + operator::O + envs::E + iter::Int + ϵ::Float64 + E_current::T +end + +function find_groundstate(mps::AbstractMPS, operator, alg::alg_type, envs = environments(mps, operator)) where {alg_type <: Union{<:IDMRG, <:IDMRG2}} + # isfinite(mps) && throw(ArgumentError("mps should be an 'InfiniteMPS'")) + (length(mps) ≤ 1 && alg isa IDMRG2) && throw(ArgumentError("unit cell should be >= w")) + log = alg isa IDMRG ? IterLog("IDMRG") : IterLog("IDMRG2") + mps = copy(mps) + iter = 0 + ϵ = calc_galerkin(mps, operator, mps, envs) + E_current = 0 LoggingExtras.withlevel(; alg.verbosity) do - @infov 2 loginit!(log, ϵ) - for outer iter in 1:(alg.maxiter) - alg_eigsolve = updatetol(alg.alg_eigsolve, iter, ϵ) - C_current = ψ.C[0] - - # sweep from left to right - for pos in 1:(length(ψ) - 1) - ac2 = AC2(ψ, pos; kind = :ACAR) - h_ac2 = AC2_hamiltonian(pos, ψ, H, ψ, envs) - _, ac2′ = fixedpoint(h_ac2, ac2, :SR, alg_eigsolve) - - al, c, ar = svd_trunc!(ac2′; trunc = alg.trscheme, alg = alg.alg_svd) - normalize!(c) - - ψ.AL[pos] = al - ψ.C[pos] = complex(c) - ψ.AR[pos + 1] = _transpose_front(ar) - ψ.AC[pos + 1] = _transpose_front(c * ar) - - transfer_leftenv!(envs, ψ, H, ψ, pos + 1) - transfer_rightenv!(envs, ψ, H, ψ, pos) - end + @infov 2 begin + E_current = expectation_value(mps, operator, envs) + loginit!(log, ϵ, E_current) + end + end - # update the edge - ψ.AL[end] = ψ.AC[end] / ψ.C[end] - ψ.AC[1] = _mul_tail(ψ.AL[1], ψ.C[1]) - ac2 = AC2(ψ, 0; kind = :ALAC) - h_ac2 = AC2_hamiltonian(0, ψ, H, ψ, envs) - _, ac2′ = fixedpoint(h_ac2, ac2, :SR, alg_eigsolve) + state = IDMRGState(mps, operator, envs, iter, ϵ, E_current) + it = IterativeSolver(alg, state) - al, c, ar = svd_trunc!(ac2′; trunc = alg.trscheme, alg = alg.alg_svd) - normalize!(c) + return LoggingExtras.withlevel(; alg.verbosity) do + for (mps, envs, ϵ, ΔE) in it + if ϵ ≤ alg.tol + @infov 2 logfinish!(log, it.iter, ϵ, ΔE) + break + end + if it.iter ≥ alg.maxiter + @warnv 1 logcancel!(log, it.iter, ϵ, ΔE) + break + end + @infov 3 logiter!(log, it.iter, ϵ, ΔE) + end - ψ.AL[end] = al - ψ.C[end] = complex(c) - ψ.AR[1] = _transpose_front(ar) + alg_gauge = updatetol(alg.alg_gauge, it.state.iter, it.state.ϵ) + ψ′ = InfiniteMPS(it.state.mps.AR[1:end]; alg_gauge.tol, alg_gauge.maxiter) + envs = recalculate!(it.state.envs, ψ′, it.state.operator, ψ′) + return ψ′, envs, it.state.ϵ + end +end - ψ.AC[end] = _mul_tail(al, c) - ψ.AC[1] = _transpose_front(c * ar) - ψ.AL[1] = ψ.AC[1] / ψ.C[1] +function Base.iterate(it::IterativeSolver{alg_type}, state = it.state) where {alg_type <: Union{<:IDMRG, <:IDMRG2}} + mps, envs, C_old = localupdate_step!(it, state) - C_current = complex(c) + # error criterion + C = mps.C[0] + smallest = infimum(_firstspace(C_old), _firstspace(C)) + e1 = isometry(_firstspace(C_old), smallest) + e2 = isometry(_firstspace(C), smallest) + ϵ = norm(e2' * C * e2 - e1' * C_old * e1) - # update environments - transfer_leftenv!(envs, ψ, H, ψ, 1) - transfer_rightenv!(envs, ψ, H, ψ, 0) + # New energy + E_new = expectation_value(mps, state.operator, envs) + ΔE = E_new - state.E_current - # sweep from right to left - for pos in (length(ψ) - 1):-1:1 - ac2 = AC2(ψ, pos; kind = :ALAC) - h_ac2 = AC2_hamiltonian(pos, ψ, H, ψ, envs) - _, ac2′ = fixedpoint(h_ac2, ac2, :SR, alg_eigsolve) + # update state + it.state = IDMRGState(mps, state.operator, envs, state.iter + 1, ϵ, E_new) - al, c, ar = svd_trunc!(ac2′; trunc = alg.trscheme, alg = alg.alg_svd) - normalize!(c) + return (mps, envs, ϵ, ΔE), it.state +end - ψ.AL[pos] = al - ψ.AC[pos] = _mul_tail(al, c) - ψ.C[pos] = complex(c) - ψ.AR[pos + 1] = _transpose_front(ar) - ψ.AC[pos + 1] = _transpose_front(c * ar) +function MPSKit.localupdate_step!( + it::IterativeSolver{<:Union{IDMRG, IDMRG2}}, state + ) + alg_eigsolve = updatetol(it.alg_eigsolve, state.iter, state.ϵ) + mps, envs, C_old = _localupdate_sweep_idmrg!(state.mps, state.operator, state.envs, alg_eigsolve, it.alg) - transfer_leftenv!(envs, ψ, H, ψ, pos + 1) - transfer_rightenv!(envs, ψ, H, ψ, pos) - end + return mps, envs, C_old +end - # update the edge - ψ.AC[end] = _mul_front(ψ.C[end - 1], ψ.AR[end]) - ψ.AR[1] = _transpose_front(ψ.C[end] \ _transpose_tail(ψ.AC[1])) - ac2 = AC2(ψ, 0; kind = :ACAR) - h_ac2 = AC2_hamiltonian(0, ψ, H, ψ, envs) - _, ac2′ = fixedpoint(h_ac2, ac2, :SR, alg_eigsolve) - al, c, ar = svd_trunc!(ac2′; trunc = alg.trscheme, alg = alg.alg_svd) - normalize!(c) - - ψ.AL[end] = al - ψ.C[end] = complex(c) - ψ.AR[1] = _transpose_front(ar) - - ψ.AR[end] = _transpose_front(ψ.C[end - 1] \ _transpose_tail(al * c)) - ψ.AC[1] = _transpose_front(c * ar) - - transfer_leftenv!(envs, ψ, H, ψ, 1) - transfer_rightenv!(envs, ψ, H, ψ, 0) - - # update error - smallest = infimum(_firstspace(C_current), _firstspace(c)) - e1 = isometry(_firstspace(C_current), smallest) - e2 = isometry(_firstspace(c), smallest) - ϵ = norm(e2' * c * e2 - e1' * C_current * e1) - - if ϵ < alg.tol - @infov 2 logfinish!(log, iter, ϵ) - break - end - if iter == alg.maxiter - @warnv 1 logcancel!(log, iter, ϵ) - else - @infov 3 logiter!(log, iter, ϵ) - end +function _localupdate_sweep_idmrg!(ψ::AbstractMPS, H, envs, alg_eigsolve, ::IDMRG) + C_old = ψ.C[0] + # left to right sweep + for pos in 1:length(ψ) + h = AC_hamiltonian(pos, ψ, H, ψ, envs) + _, ψ.AC[pos] = fixedpoint(h, ψ.AC[pos], :SR, alg_eigsolve) + if pos == length(ψ) + # AC needed in next sweep + ψ.AL[pos], ψ.C[pos] = left_orth(ψ.AC[pos]) + else + ψ.AL[pos], ψ.C[pos] = left_orth!(ψ.AC[pos]) end + transfer_leftenv!(envs, ψ, H, ψ, pos + 1) end - alg_gauge = updatetol(alg.alg_gauge, iter, ϵ) - ψ′ = InfiniteMPS(ψ.AR[1:end]; alg_gauge.tol, alg_gauge.maxiter) - recalculate!(envs, ψ′, H, ψ′) + # right to left sweep + for pos in length(ψ):-1:1 + h = AC_hamiltonian(pos, ψ, H, ψ, envs) + _, ψ.AC[pos] = fixedpoint(h, ψ.AC[pos], :SR, alg_eigsolve) - return ψ′, envs, ϵ + ψ.C[pos - 1], temp = right_orth!(_transpose_tail(ψ.AC[pos]; copy = (pos == 1))) + ψ.AR[pos] = _transpose_front(temp) + + transfer_rightenv!(envs, ψ, H, ψ, pos - 1) + end + return ψ, envs, C_old end + + +function _localupdate_sweep_idmrg!(ψ::AbstractMPS, H, envs, alg_eigsolve, alg::IDMRG2) + # sweep from left to right + for pos in 1:(length(ψ) - 1) + ac2 = AC2(ψ, pos; kind = :ACAR) + h_ac2 = AC2_hamiltonian(pos, ψ, H, ψ, envs) + _, ac2′ = fixedpoint(h_ac2, ac2, :SR, alg_eigsolve) + + al, c, ar = svd_trunc!(ac2′; trunc = alg.trscheme, alg = alg.alg_svd) + normalize!(c) + + ψ.AL[pos] = al + ψ.C[pos] = complex(c) + ψ.AR[pos + 1] = _transpose_front(ar) + ψ.AC[pos + 1] = _transpose_front(c * ar) + + transfer_leftenv!(envs, ψ, H, ψ, pos + 1) + transfer_rightenv!(envs, ψ, H, ψ, pos) + end + + # update the edge + ψ.AL[end] = ψ.AC[end] / ψ.C[end] + ψ.AC[1] = _mul_tail(ψ.AL[1], ψ.C[1]) + ac2 = AC2(ψ, 0; kind = :ALAC) + h_ac2 = AC2_hamiltonian(0, ψ, H, ψ, envs) + _, ac2′ = fixedpoint(h_ac2, ac2, :SR, alg_eigsolve) + + al, c, ar = svd_trunc!(ac2′; trunc = alg.trscheme, alg = alg.alg_svd) + normalize!(c) + + ψ.AL[end] = al + ψ.C[end] = complex(c) + ψ.AR[1] = _transpose_front(ar) + + ψ.AC[end] = _mul_tail(al, c) + ψ.AC[1] = _transpose_front(c * ar) + ψ.AL[1] = ψ.AC[1] / ψ.C[1] + + C_old = complex(c) + + # update environments + transfer_leftenv!(envs, ψ, H, ψ, 1) + transfer_rightenv!(envs, ψ, H, ψ, 0) + + # sweep from right to left + for pos in (length(ψ) - 1):-1:1 + ac2 = AC2(ψ, pos; kind = :ALAC) + h_ac2 = AC2_hamiltonian(pos, ψ, H, ψ, envs) + _, ac2′ = fixedpoint(h_ac2, ac2, :SR, alg_eigsolve) + + al, c, ar = svd_trunc!(ac2′; trunc = alg.trscheme, alg = alg.alg_svd) + normalize!(c) + + ψ.AL[pos] = al + ψ.AC[pos] = _mul_tail(al, c) + ψ.C[pos] = complex(c) + ψ.AR[pos + 1] = _transpose_front(ar) + ψ.AC[pos + 1] = _transpose_front(c * ar) + + transfer_leftenv!(envs, ψ, H, ψ, pos + 1) + transfer_rightenv!(envs, ψ, H, ψ, pos) + end + + # update the edge + ψ.AC[end] = _mul_front(ψ.C[end - 1], ψ.AR[end]) + ψ.AR[1] = _transpose_front(ψ.C[end] \ _transpose_tail(ψ.AC[1])) + ac2 = AC2(ψ, 0; kind = :ACAR) + h_ac2 = AC2_hamiltonian(0, ψ, H, ψ, envs) + _, ac2′ = fixedpoint(h_ac2, ac2, :SR, alg_eigsolve) + al, c, ar = svd_trunc!(ac2′; trunc = alg.trscheme, alg = alg.alg_svd) + normalize!(c) + + ψ.AL[end] = al + ψ.C[end] = complex(c) + ψ.AR[1] = _transpose_front(ar) + + ψ.AR[end] = _transpose_front(ψ.C[end - 1] \ _transpose_tail(al * c)) + ψ.AC[1] = _transpose_front(c * ar) + + transfer_leftenv!(envs, ψ, H, ψ, 1) + transfer_rightenv!(envs, ψ, H, ψ, 0) + + return ψ, envs, C_old +end \ No newline at end of file From 3095c0f34071e7ee5cbb376e5ce41b102cee677e Mon Sep 17 00:00:00 2001 From: afeuerpfeil Date: Wed, 5 Nov 2025 19:30:28 -0500 Subject: [PATCH 2/7] fix format --- src/algorithms/groundstate/idmrg.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/algorithms/groundstate/idmrg.jl b/src/algorithms/groundstate/idmrg.jl index 1163fc86c..80df1e1e7 100644 --- a/src/algorithms/groundstate/idmrg.jl +++ b/src/algorithms/groundstate/idmrg.jl @@ -246,4 +246,4 @@ function _localupdate_sweep_idmrg!(ψ::AbstractMPS, H, envs, alg_eigsolve, alg:: transfer_rightenv!(envs, ψ, H, ψ, 0) return ψ, envs, C_old -end \ No newline at end of file +end From 06f02279f9fd3a5beef9bd4c9974d5f52f0fe9e2 Mon Sep 17 00:00:00 2001 From: afeuerpfeil Date: Thu, 6 Nov 2025 11:49:49 -0500 Subject: [PATCH 3/7] Fix energy definition for IDMRG2 --- src/algorithms/groundstate/idmrg.jl | 30 +++++++++++++++-------------- 1 file changed, 16 insertions(+), 14 deletions(-) diff --git a/src/algorithms/groundstate/idmrg.jl b/src/algorithms/groundstate/idmrg.jl index 80df1e1e7..5047bab66 100644 --- a/src/algorithms/groundstate/idmrg.jl +++ b/src/algorithms/groundstate/idmrg.jl @@ -82,6 +82,7 @@ function find_groundstate(mps::AbstractMPS, operator, alg::alg_type, envs = envi loginit!(log, ϵ, E_current) end end + @show E_current state = IDMRGState(mps, operator, envs, iter, ϵ, E_current) it = IterativeSolver(alg, state) @@ -106,8 +107,8 @@ function find_groundstate(mps::AbstractMPS, operator, alg::alg_type, envs = envi end end -function Base.iterate(it::IterativeSolver{alg_type}, state = it.state) where {alg_type <: Union{<:IDMRG, <:IDMRG2}} - mps, envs, C_old = localupdate_step!(it, state) +function Base.iterate(it::IterativeSolver{alg_type}, state::IDMRGState{S, O, E, T} = it.state) where {alg_type <: Union{<:IDMRG, <:IDMRG2}, S, O, E, T <: Number} + mps, envs, C_old, E_new = localupdate_step!(it, state) # error criterion C = mps.C[0] @@ -117,11 +118,11 @@ function Base.iterate(it::IterativeSolver{alg_type}, state = it.state) where {al ϵ = norm(e2' * C * e2 - e1' * C_old * e1) # New energy - E_new = expectation_value(mps, state.operator, envs) - ΔE = E_new - state.E_current - + # E_new = expectation_value(mps, state.operator, envs) + ΔE = (E_new - state.E_current)/(length(mps)) + alg_type <: IDMRG2 && (ΔE /= 2) # update state - it.state = IDMRGState(mps, state.operator, envs, state.iter + 1, ϵ, E_new) + it.state = IDMRGState(mps, state.operator, envs, state.iter + 1, ϵ, T(E_new)) return (mps, envs, ϵ, ΔE), it.state end @@ -130,12 +131,13 @@ function MPSKit.localupdate_step!( it::IterativeSolver{<:Union{IDMRG, IDMRG2}}, state ) alg_eigsolve = updatetol(it.alg_eigsolve, state.iter, state.ϵ) - mps, envs, C_old = _localupdate_sweep_idmrg!(state.mps, state.operator, state.envs, alg_eigsolve, it.alg) + mps, envs, C_old, E = _localupdate_sweep_idmrg!(state.mps, state.operator, state.envs, alg_eigsolve, it.alg) - return mps, envs, C_old + return mps, envs, C_old, E end function _localupdate_sweep_idmrg!(ψ::AbstractMPS, H, envs, alg_eigsolve, ::IDMRG) + E=0 C_old = ψ.C[0] # left to right sweep for pos in 1:length(ψ) @@ -153,18 +155,19 @@ function _localupdate_sweep_idmrg!(ψ::AbstractMPS, H, envs, alg_eigsolve, ::IDM # right to left sweep for pos in length(ψ):-1:1 h = AC_hamiltonian(pos, ψ, H, ψ, envs) - _, ψ.AC[pos] = fixedpoint(h, ψ.AC[pos], :SR, alg_eigsolve) + E, ψ.AC[pos] = fixedpoint(h, ψ.AC[pos], :SR, alg_eigsolve) ψ.C[pos - 1], temp = right_orth!(_transpose_tail(ψ.AC[pos]; copy = (pos == 1))) ψ.AR[pos] = _transpose_front(temp) transfer_rightenv!(envs, ψ, H, ψ, pos - 1) end - return ψ, envs, C_old + return ψ, envs, C_old, E end function _localupdate_sweep_idmrg!(ψ::AbstractMPS, H, envs, alg_eigsolve, alg::IDMRG2) + E=0 # sweep from left to right for pos in 1:(length(ψ) - 1) ac2 = AC2(ψ, pos; kind = :ACAR) @@ -231,7 +234,7 @@ function _localupdate_sweep_idmrg!(ψ::AbstractMPS, H, envs, alg_eigsolve, alg:: ψ.AR[1] = _transpose_front(ψ.C[end] \ _transpose_tail(ψ.AC[1])) ac2 = AC2(ψ, 0; kind = :ACAR) h_ac2 = AC2_hamiltonian(0, ψ, H, ψ, envs) - _, ac2′ = fixedpoint(h_ac2, ac2, :SR, alg_eigsolve) + E, ac2′ = fixedpoint(h_ac2, ac2, :SR, alg_eigsolve) al, c, ar = svd_trunc!(ac2′; trunc = alg.trscheme, alg = alg.alg_svd) normalize!(c) @@ -244,6 +247,5 @@ function _localupdate_sweep_idmrg!(ψ::AbstractMPS, H, envs, alg_eigsolve, alg:: transfer_leftenv!(envs, ψ, H, ψ, 1) transfer_rightenv!(envs, ψ, H, ψ, 0) - - return ψ, envs, C_old -end + return ψ, envs, C_old, E +end \ No newline at end of file From 8161a2cc28fb179e69758a4195bc6fd1b014baae Mon Sep 17 00:00:00 2001 From: afeuerpfeil Date: Thu, 6 Nov 2025 11:50:09 -0500 Subject: [PATCH 4/7] fix energy of IDMRG2 --- src/algorithms/groundstate/idmrg.jl | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/src/algorithms/groundstate/idmrg.jl b/src/algorithms/groundstate/idmrg.jl index 5047bab66..cec186f7c 100644 --- a/src/algorithms/groundstate/idmrg.jl +++ b/src/algorithms/groundstate/idmrg.jl @@ -58,13 +58,13 @@ end # Internal state of the IDMRG algorithm -struct IDMRGState{S, O, E, T <: Number} +struct IDMRGState{S, O, E} mps::S operator::O envs::E iter::Int ϵ::Float64 - E_current::T + E_current::Number end function find_groundstate(mps::AbstractMPS, operator, alg::alg_type, envs = environments(mps, operator)) where {alg_type <: Union{<:IDMRG, <:IDMRG2}} @@ -82,7 +82,6 @@ function find_groundstate(mps::AbstractMPS, operator, alg::alg_type, envs = envi loginit!(log, ϵ, E_current) end end - @show E_current state = IDMRGState(mps, operator, envs, iter, ϵ, E_current) it = IterativeSolver(alg, state) @@ -107,7 +106,7 @@ function find_groundstate(mps::AbstractMPS, operator, alg::alg_type, envs = envi end end -function Base.iterate(it::IterativeSolver{alg_type}, state::IDMRGState{S, O, E, T} = it.state) where {alg_type <: Union{<:IDMRG, <:IDMRG2}, S, O, E, T <: Number} +function Base.iterate(it::IterativeSolver{alg_type}, state= it.state) where {alg_type <: Union{<:IDMRG, <:IDMRG2}} mps, envs, C_old, E_new = localupdate_step!(it, state) # error criterion @@ -118,11 +117,11 @@ function Base.iterate(it::IterativeSolver{alg_type}, state::IDMRGState{S, O, E, ϵ = norm(e2' * C * e2 - e1' * C_old * e1) # New energy - # E_new = expectation_value(mps, state.operator, envs) - ΔE = (E_new - state.E_current)/(length(mps)) - alg_type <: IDMRG2 && (ΔE /= 2) + ΔE = (E_new - state.E_current)/2 + (alg_type <: IDMRG2 && length(mps) == 2) && (ΔE /= 2) + # update state - it.state = IDMRGState(mps, state.operator, envs, state.iter + 1, ϵ, T(E_new)) + it.state = IDMRGState(mps, state.operator, envs, state.iter + 1, ϵ, E_new) return (mps, envs, ϵ, ΔE), it.state end From 883d244f3ec095537c1f320d57916f4e91be50c1 Mon Sep 17 00:00:00 2001 From: AFeuerpfeil Date: Mon, 1 Dec 2025 20:59:23 -0500 Subject: [PATCH 5/7] make energy in IDMRGState type stable and implement suggestion by Lukas --- src/algorithms/groundstate/idmrg.jl | 53 +++++++++++++++++------------ src/utility/iterativesolvers.jl | 2 +- 2 files changed, 32 insertions(+), 23 deletions(-) diff --git a/src/algorithms/groundstate/idmrg.jl b/src/algorithms/groundstate/idmrg.jl index cec186f7c..506f3132c 100644 --- a/src/algorithms/groundstate/idmrg.jl +++ b/src/algorithms/groundstate/idmrg.jl @@ -58,32 +58,34 @@ end # Internal state of the IDMRG algorithm -struct IDMRGState{S, O, E} +struct IDMRGState{S, O, E, T} mps::S operator::O envs::E iter::Int ϵ::Float64 - E_current::Number + energy::T +end +function IDMRGState{T}(mps::S, operator::O, envs::E, iter::Int, ϵ::Float64, energy) where {S, O, E, T} + return IDMRGState{S, O, E, T}(mps, operator, envs, iter, ϵ, T(energy)) end -function find_groundstate(mps::AbstractMPS, operator, alg::alg_type, envs = environments(mps, operator)) where {alg_type <: Union{<:IDMRG, <:IDMRG2}} - # isfinite(mps) && throw(ArgumentError("mps should be an 'InfiniteMPS'")) - (length(mps) ≤ 1 && alg isa IDMRG2) && throw(ArgumentError("unit cell should be >= w")) +function find_groundstate(mps, operator, alg::alg_type, envs = environments(mps, operator)) where {alg_type <: Union{<:IDMRG, <:IDMRG2}} + (length(mps) ≤ 1 && alg isa IDMRG2) && throw(ArgumentError("unit cell should be >= 2")) log = alg isa IDMRG ? IterLog("IDMRG") : IterLog("IDMRG2") mps = copy(mps) iter = 0 ϵ = calc_galerkin(mps, operator, mps, envs) - E_current = 0 + E = zero(Base.promote_type(scalartype(mps), scalartype(operator))) LoggingExtras.withlevel(; alg.verbosity) do @infov 2 begin - E_current = expectation_value(mps, operator, envs) - loginit!(log, ϵ, E_current) + E = expectation_value(mps, operator, envs) + loginit!(log, ϵ, E) end end - state = IDMRGState(mps, operator, envs, iter, ϵ, E_current) + state = IDMRGState(mps, operator, envs, iter, ϵ, E) it = IterativeSolver(alg, state) return LoggingExtras.withlevel(; alg.verbosity) do @@ -106,22 +108,30 @@ function find_groundstate(mps::AbstractMPS, operator, alg::alg_type, envs = envi end end -function Base.iterate(it::IterativeSolver{alg_type}, state= it.state) where {alg_type <: Union{<:IDMRG, <:IDMRG2}} +function Base.iterate( + it::IterativeSolver{alg_type}, state::IDMRGState{<:Any, <:Any, <:Any, T}= it.state + ) where {alg_type <: Union{<:IDMRG, <:IDMRG2}, T} mps, envs, C_old, E_new = localupdate_step!(it, state) # error criterion C = mps.C[0] - smallest = infimum(_firstspace(C_old), _firstspace(C)) - e1 = isometry(_firstspace(C_old), smallest) - e2 = isometry(_firstspace(C), smallest) - ϵ = norm(e2' * C * e2 - e1' * C_old * e1) + space_C_old = _firstspace(C_old) + space_C = _firstspace(C) + if space_C != space_C_old + smallest = infimum(space_C_old, space_C) + e1 = isometry(space_C_old, smallest) + e2 = isometry(space_C, smallest) + ϵ = norm(e2' * C * e2 - e1' * C_old * e1) + else + ϵ = norm(C - C_old) + end # New energy - ΔE = (E_new - state.E_current)/2 - (alg_type <: IDMRG2 && length(mps) == 2) && (ΔE /= 2) + ΔE = (E_new - state.energy)/2 + (alg_type <: IDMRG2 && length(mps) == 2) && (ΔE /= 2) # This extra factor gives the correct energy per unit cell. I have no clue why right now. # update state - it.state = IDMRGState(mps, state.operator, envs, state.iter + 1, ϵ, E_new) + it.state = IDMRGState{T}(mps, state.operator, envs, state.iter + 1, ϵ, E_new) return (mps, envs, ϵ, ΔE), it.state end @@ -135,8 +145,8 @@ function MPSKit.localupdate_step!( return mps, envs, C_old, E end -function _localupdate_sweep_idmrg!(ψ::AbstractMPS, H, envs, alg_eigsolve, ::IDMRG) - E=0 +function _localupdate_sweep_idmrg!(ψ, H, envs, alg_eigsolve, ::IDMRG) + local E C_old = ψ.C[0] # left to right sweep for pos in 1:length(ψ) @@ -165,8 +175,7 @@ function _localupdate_sweep_idmrg!(ψ::AbstractMPS, H, envs, alg_eigsolve, ::IDM end -function _localupdate_sweep_idmrg!(ψ::AbstractMPS, H, envs, alg_eigsolve, alg::IDMRG2) - E=0 +function _localupdate_sweep_idmrg!(ψ, H, envs, alg_eigsolve, alg::IDMRG2) # sweep from left to right for pos in 1:(length(ψ) - 1) ac2 = AC2(ψ, pos; kind = :ACAR) @@ -247,4 +256,4 @@ function _localupdate_sweep_idmrg!(ψ::AbstractMPS, H, envs, alg_eigsolve, alg:: transfer_leftenv!(envs, ψ, H, ψ, 1) transfer_rightenv!(envs, ψ, H, ψ, 0) return ψ, envs, C_old, E -end \ No newline at end of file +end diff --git a/src/utility/iterativesolvers.jl b/src/utility/iterativesolvers.jl index e74218544..8eb863d12 100644 --- a/src/utility/iterativesolvers.jl +++ b/src/utility/iterativesolvers.jl @@ -7,7 +7,7 @@ mutable struct IterativeSolver{A, B} end function Base.getproperty(it::IterativeSolver{A, B}, name::Symbol) where {A, B} - name === :alg || name === :state && return getfield(it, name) + (name === :alg || name === :state) && return getfield(it, name) alg = getfield(it, :alg) name in propertynames(alg) && return getproperty(alg, name) From 0aa0ecb45e9a320bc96f4ea0628506826b894830 Mon Sep 17 00:00:00 2001 From: AFeuerpfeil Date: Mon, 1 Dec 2025 21:00:26 -0500 Subject: [PATCH 6/7] format --- src/algorithms/groundstate/idmrg.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/algorithms/groundstate/idmrg.jl b/src/algorithms/groundstate/idmrg.jl index 506f3132c..762aaf9f5 100644 --- a/src/algorithms/groundstate/idmrg.jl +++ b/src/algorithms/groundstate/idmrg.jl @@ -66,7 +66,7 @@ struct IDMRGState{S, O, E, T} ϵ::Float64 energy::T end -function IDMRGState{T}(mps::S, operator::O, envs::E, iter::Int, ϵ::Float64, energy) where {S, O, E, T} +function IDMRGState{T}(mps::S, operator::O, envs::E, iter::Int, ϵ::Float64, energy) where {S, O, E, T} return IDMRGState{S, O, E, T}(mps, operator, envs, iter, ϵ, T(energy)) end @@ -109,7 +109,7 @@ function find_groundstate(mps, operator, alg::alg_type, envs = environments(mps, end function Base.iterate( - it::IterativeSolver{alg_type}, state::IDMRGState{<:Any, <:Any, <:Any, T}= it.state + it::IterativeSolver{alg_type}, state::IDMRGState{<:Any, <:Any, <:Any, T} = it.state ) where {alg_type <: Union{<:IDMRG, <:IDMRG2}, T} mps, envs, C_old, E_new = localupdate_step!(it, state) @@ -127,7 +127,7 @@ function Base.iterate( end # New energy - ΔE = (E_new - state.energy)/2 + ΔE = (E_new - state.energy) / 2 (alg_type <: IDMRG2 && length(mps) == 2) && (ΔE /= 2) # This extra factor gives the correct energy per unit cell. I have no clue why right now. # update state From 510be3e7e8cba4f2de7a12cb6cd9039cc34553a2 Mon Sep 17 00:00:00 2001 From: AFeuerpfeil Date: Tue, 2 Dec 2025 08:47:42 -0500 Subject: [PATCH 7/7] apply Lukas comments --- src/algorithms/groundstate/idmrg.jl | 29 +++++++++++++++++------------ 1 file changed, 17 insertions(+), 12 deletions(-) diff --git a/src/algorithms/groundstate/idmrg.jl b/src/algorithms/groundstate/idmrg.jl index 762aaf9f5..6cb71afd3 100644 --- a/src/algorithms/groundstate/idmrg.jl +++ b/src/algorithms/groundstate/idmrg.jl @@ -63,7 +63,7 @@ struct IDMRGState{S, O, E, T} operator::O envs::E iter::Int - ϵ::Float64 + ϵ::Float64 # TODO: Could be any <:Real energy::T end function IDMRGState{T}(mps::S, operator::O, envs::E, iter::Int, ϵ::Float64, energy) where {S, O, E, T} @@ -76,7 +76,7 @@ function find_groundstate(mps, operator, alg::alg_type, envs = environments(mps, mps = copy(mps) iter = 0 ϵ = calc_galerkin(mps, operator, mps, envs) - E = zero(Base.promote_type(scalartype(mps), scalartype(operator))) + E = zero(TensorOperations.promote_contract(scalartype(mps), scalartype(operator))) LoggingExtras.withlevel(; alg.verbosity) do @infov 2 begin @@ -136,16 +136,21 @@ function Base.iterate( return (mps, envs, ϵ, ΔE), it.state end -function MPSKit.localupdate_step!( - it::IterativeSolver{<:Union{IDMRG, IDMRG2}}, state +function localupdate_step!( + it::IterativeSolver{<:IDMRG}, state ) alg_eigsolve = updatetol(it.alg_eigsolve, state.iter, state.ϵ) - mps, envs, C_old, E = _localupdate_sweep_idmrg!(state.mps, state.operator, state.envs, alg_eigsolve, it.alg) + return _localupdate_sweep_idmrg!(state.mps, state.operator, state.envs, alg_eigsolve) +end - return mps, envs, C_old, E +function localupdate_step!( + it::IterativeSolver{<:IDMRG2}, state + ) + alg_eigsolve = updatetol(it.alg_eigsolve, state.iter, state.ϵ) + return _localupdate_sweep_idmrg2!(state.mps, state.operator, state.envs, alg_eigsolve, it.trscheme, it.alg_svd) end -function _localupdate_sweep_idmrg!(ψ, H, envs, alg_eigsolve, ::IDMRG) +function _localupdate_sweep_idmrg!(ψ, H, envs, alg_eigsolve) local E C_old = ψ.C[0] # left to right sweep @@ -175,14 +180,14 @@ function _localupdate_sweep_idmrg!(ψ, H, envs, alg_eigsolve, ::IDMRG) end -function _localupdate_sweep_idmrg!(ψ, H, envs, alg_eigsolve, alg::IDMRG2) +function _localupdate_sweep_idmrg2!(ψ, H, envs, alg_eigsolve, alg_trscheme, alg_svd) # sweep from left to right for pos in 1:(length(ψ) - 1) ac2 = AC2(ψ, pos; kind = :ACAR) h_ac2 = AC2_hamiltonian(pos, ψ, H, ψ, envs) _, ac2′ = fixedpoint(h_ac2, ac2, :SR, alg_eigsolve) - al, c, ar = svd_trunc!(ac2′; trunc = alg.trscheme, alg = alg.alg_svd) + al, c, ar = svd_trunc!(ac2′; trunc = alg_trscheme, alg = alg_svd) normalize!(c) ψ.AL[pos] = al @@ -201,7 +206,7 @@ function _localupdate_sweep_idmrg!(ψ, H, envs, alg_eigsolve, alg::IDMRG2) h_ac2 = AC2_hamiltonian(0, ψ, H, ψ, envs) _, ac2′ = fixedpoint(h_ac2, ac2, :SR, alg_eigsolve) - al, c, ar = svd_trunc!(ac2′; trunc = alg.trscheme, alg = alg.alg_svd) + al, c, ar = svd_trunc!(ac2′; trunc = alg_trscheme, alg = alg_svd) normalize!(c) ψ.AL[end] = al @@ -224,7 +229,7 @@ function _localupdate_sweep_idmrg!(ψ, H, envs, alg_eigsolve, alg::IDMRG2) h_ac2 = AC2_hamiltonian(pos, ψ, H, ψ, envs) _, ac2′ = fixedpoint(h_ac2, ac2, :SR, alg_eigsolve) - al, c, ar = svd_trunc!(ac2′; trunc = alg.trscheme, alg = alg.alg_svd) + al, c, ar = svd_trunc!(ac2′; trunc = alg_trscheme, alg = alg_svd) normalize!(c) ψ.AL[pos] = al @@ -243,7 +248,7 @@ function _localupdate_sweep_idmrg!(ψ, H, envs, alg_eigsolve, alg::IDMRG2) ac2 = AC2(ψ, 0; kind = :ACAR) h_ac2 = AC2_hamiltonian(0, ψ, H, ψ, envs) E, ac2′ = fixedpoint(h_ac2, ac2, :SR, alg_eigsolve) - al, c, ar = svd_trunc!(ac2′; trunc = alg.trscheme, alg = alg.alg_svd) + al, c, ar = svd_trunc!(ac2′; trunc = alg_trscheme, alg = alg_svd) normalize!(c) ψ.AL[end] = al