diff --git a/src/algorithms/derivatives/hamiltonian_derivatives.jl b/src/algorithms/derivatives/hamiltonian_derivatives.jl index b922205f8..b35bf85d2 100644 --- a/src/algorithms/derivatives/hamiltonian_derivatives.jl +++ b/src/algorithms/derivatives/hamiltonian_derivatives.jl @@ -10,25 +10,13 @@ struct JordanMPO_AC_Hamiltonian{O1, O2, O3} <: DerivativeOperator E::Union{O1, Missing} # finished C::Union{O2, Missing} # starting B::Union{O2, Missing} # ending - A::O3 # continuing - function JordanMPO_AC_Hamiltonian( - onsite, not_started, finished, starting, ending, continuing - ) - # obtaining storagetype of environments since these should have already mixed - # the types of the operator and state - gl = continuing[1] - S = spacetype(gl) - M = storagetype(gl) - O1 = tensormaptype(S, 1, 1, M) - O2 = tensormaptype(S, 2, 2, M) - return new{O1, O2, typeof(continuing)}( - onsite, not_started, finished, starting, ending, continuing - ) - end + A::Union{O3, Missing} # continuing + + # need inner constructor to prohibit no-type-param constructor with unbound vars function JordanMPO_AC_Hamiltonian{O1, O2, O3}( - onsite, not_started, finished, starting, ending, continuing + D, I, E, C, B, A, ) where {O1, O2, O3} - return new{O1, O2, O3}(onsite, not_started, finished, starting, ending, continuing) + return new{O1, O2, O3}(D, I, E, C, B, A) end end @@ -45,21 +33,12 @@ struct JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4} <: DerivativeOperator CB::Union{O2, Missing} # starting left - ending right CA::Union{O3, Missing} # starting left - continuing right AB::Union{O3, Missing} # continuing left - ending right - AA::O4 # continuing left - continuing right + AA::Union{O4, Missing} # continuing left - continuing right BE::Union{O2, Missing} # ending left DE::Union{O1, Missing} # onsite left EE::Union{O1, Missing} # finished - function JordanMPO_AC2_Hamiltonian(II, IC, ID, CB, CA, AB, AA, BE, DE, EE) - # obtaining storagetype of environments since these should have already mixed - # the types of the operator and state - gl = AA[1] - S = spacetype(gl) - M = storagetype(gl) - O1 = tensormaptype(S, 1, 1, M) - O2 = tensormaptype(S, 2, 2, M) - O3 = tensormaptype(S, 3, 3, M) - return new{O1, O2, O3, typeof(AA)}(II, IC, ID, CB, CA, AB, AA, BE, DE, EE) - end + + # need inner constructor to prohibit no-type-param constructor with unbound vars function JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4}( II, IC, ID, CB, CA, AB, AA, BE, DE, EE ) where {O1, O2, O3, O4} @@ -69,12 +48,6 @@ end # Constructors # ------------ -const _HAM_MPS_TYPES = Union{ - FiniteMPS{<:MPSTensor}, - WindowMPS{<:MPSTensor}, - InfiniteMPS{<:MPSTensor}, -} - function AC_hamiltonian( site::Int, below::_HAM_MPS_TYPES, operator::MPOHamiltonian{<:JordanMPOTensor}, above::_HAM_MPS_TYPES, envs @@ -83,9 +56,11 @@ function AC_hamiltonian( GR = rightenv(envs, site, below) W = operator[site] + GR_2 = GR[2:(end - 1)] + GL_2 = GL[2:(end - 1)] + # starting if nonzero_length(W.C) > 0 - GR_2 = GR[2:(end - 1)] @plansor starting_[-1 -2; -3 -4] ≔ W.C[-1; -3 1] * GR_2[-4 1; -2] starting = only(starting_) else @@ -94,7 +69,6 @@ function AC_hamiltonian( # ending if nonzero_length(W.B) > 0 - GL_2 = GL[2:(end - 1)] @plansor ending_[-1 -2; -3 -4] ≔ GL_2[-1 1; -3] * W.B[1 -2; -4] ending = only(ending_) else @@ -138,11 +112,19 @@ function AC_hamiltonian( finished = removeunit(GL[end], 2) end - # continuing - A = W.A - continuing = (GL[2:(end - 1)], A, GR[2:(end - 1)]) + if nonzero_length(W.A) > 0 + continuing = AC_hamiltonian(GL_2, W.A, GR_2) + else + continuing = missing + end + + S = spacetype(GL) + M = storagetype(GL) + O1 = tensormaptype(S, 1, 1, M) + O2 = tensormaptype(S, 2, 2, M) + O3 = Core.Compiler.return_type(AC_hamiltonian, typeof((GL_2, W.A, GR_2))) - return JordanMPO_AC_Hamiltonian( + return JordanMPO_AC_Hamiltonian{O1, O2, O3}( onsite, not_started, finished, starting, ending, continuing ) end @@ -156,10 +138,13 @@ function AC2_hamiltonian( W1 = operator[site] W2 = operator[site + 1] + GR_2 = GR[2:(end - 1)] + GL_2 = GL[2:(end - 1)] + # starting left - continuing right if nonzero_length(W1.C) > 0 && nonzero_length(W2.A) > 0 @plansor CA_[-1 -2 -3; -4 -5 -6] ≔ W1.C[-1; -4 2] * W2.A[2 -2; -5 1] * - GR[2:(end - 1)][-6 1; -3] + GR_2[-6 1; -3] CA = only(CA_) else CA = missing @@ -167,7 +152,7 @@ function AC2_hamiltonian( # continuing left - ending right if nonzero_length(W1.A) > 0 && nonzero_length(W2.B) > 0 - @plansor AB_[-1 -2 -3; -4 -5 -6] ≔ GL[2:(end - 1)][-1 2; -4] * W1.A[2 -2; -5 1] * + @plansor AB_[-1 -2 -3; -4 -5 -6] ≔ GL_2[-1 2; -4] * W1.A[2 -2; -5 1] * W2.B[1 -3; -6] AB = only(AB_) else @@ -197,10 +182,11 @@ function AC2_hamiltonian( if !ismissing(CA) I = id(storagetype(GR[1]), physicalspace(W1)) @plansor CA[-1 -2 -3; -4 -5 -6] += (I[-1; -4] * W2.C[-2; -5 1]) * - GR[2:(end - 1)][-6 1; -3] + GR_2[-6 1; -3] IC = missing else - @plansor IC[-1 -2; -3 -4] ≔ W2.C[-1; -3 1] * GR[2:(end - 1)][-4 1; -2] + @plansor IC[-1 -2; -3 -4] ≔ W2.C[-1; -3 1] * GR_2[-4 1; -2] + IC = only(IC) end else IC = missing @@ -210,11 +196,12 @@ function AC2_hamiltonian( if nonzero_length(W1.B) > 0 if !ismissing(AB) I = id(storagetype(GL[end]), physicalspace(W2)) - @plansor AB[-1 -2 -3; -4 -5 -6] += GL[2:(end - 1)][-1 1; -4] * + @plansor AB[-1 -2 -3; -4 -5 -6] += GL_2[-1 1; -4] * (W1.B[1 -2; -5] * I[-3; -6]) BE = missing else - @plansor BE[-1 -2; -3 -4] ≔ GL[2:(end - 1)][-1 2; -3] * W1.B[2 -2; -4] + @plansor BE[-1 -2; -3 -4] ≔ GL_2[-1 2; -3] * W1.B[2 -2; -4] + BE = only(BE) end else BE = missing @@ -286,33 +273,38 @@ function AC2_hamiltonian( end # continuing - continuing - # TODO: MPODerivativeOperator code reuse + optimization - AA = (GL[2:(end - 1)], W1.A, W2.A, GR[2:(end - 1)]) + ## TODO: Think about how one could and whether one should store these objects and use them for (a) advancing environments in iDMRG, (b) reuse ind backwards-sweep in IDMRG, (c) subspace expansion + if nonzero_length(W1.A) > 0 && nonzero_length(W2.A) > 0 + AA = AC2_hamiltonian(GL_2, W1.A, W2.A, GR_2) + else + AA = missing + end - return JordanMPO_AC2_Hamiltonian(II, IC, ID, CB, CA, AB, AA, BE, DE, EE) + S = spacetype(GL) + M = storagetype(GL) + O1 = tensormaptype(S, 1, 1, M) + O2 = tensormaptype(S, 2, 2, M) + O3 = tensormaptype(S, 3, 3, M) + O4 = Core.Compiler.return_type(AC2_hamiltonian, typeof((GL_2, W1.A, W2.A, GR_2))) + + return JordanMPO_AC2_Hamiltonian{O1, O2, O3, O4}(II, IC, ID, CB, CA, AB, AA, BE, DE, EE) end # Actions # ------- function (H::JordanMPO_AC_Hamiltonian)(x::MPSTensor) - y = zerovector(x) - + y = ismissing(H.A) ? zerovector(x) : H.A(x) ismissing(H.D) || @plansor y[-1 -2; -3] += x[-1 1; -3] * H.D[-2; 1] ismissing(H.E) || @plansor y[-1 -2; -3] += H.E[-1; 1] * x[1 -2; -3] ismissing(H.I) || @plansor y[-1 -2; -3] += x[-1 -2; 1] * H.I[1; -3] ismissing(H.C) || @plansor y[-1 -2; -3] += x[-1 2; 1] * H.C[-2 -3; 2 1] ismissing(H.B) || @plansor y[-1 -2; -3] += H.B[-1 -2; 1 2] * x[1 2; -3] - GL, A, GR = H.A - if nonzero_length(A) > 0 - @plansor y[-1 -2; -3] += GL[-1 5; 4] * x[4 2; 1] * A[5 -2; 2 3] * GR[1 3; -3] - end - return y end function (H::JordanMPO_AC2_Hamiltonian)(x::MPOTensor) - y = zerovector(x) + y = ismissing(H.AA) ? zerovector(x) : H.AA(x) ismissing(H.II) || @plansor y[-1 -2; -3 -4] += x[-1 -2; 1 -4] * H.II[-3; 1] ismissing(H.IC) || @plansor y[-1 -2; -3 -4] += x[-1 -2; 1 2] * H.IC[-4 -3; 2 1] ismissing(H.ID) || @plansor y[-1 -2; -3 -4] += x[-1 -2; -3 1] * H.ID[-4; 1] @@ -323,12 +315,5 @@ function (H::JordanMPO_AC2_Hamiltonian)(x::MPOTensor) ismissing(H.DE) || @plansor y[-1 -2; -3 -4] += x[-1 1; -3 -4] * H.DE[-2; 1] ismissing(H.EE) || @plansor y[-1 -2; -3 -4] += x[1 -2; -3 -4] * H.EE[-1; 1] - GL, A1, A2, GR = H.AA - if nonzero_length(A1) > 0 && nonzero_length(A2) > 0 - # TODO: there are too many entries here, this could be further optimized - @plansor y[-1 -2; -3 -4] += GL[-1 7; 6] * x[6 5; 1 3] * A1[7 -2; 5 4] * - A2[4 -4; 3 2] * GR[1 2; -3] - end - return y end diff --git a/src/algorithms/derivatives/mpo_derivatives.jl b/src/algorithms/derivatives/mpo_derivatives.jl index 12b030b1e..02ec00b34 100644 --- a/src/algorithms/derivatives/mpo_derivatives.jl +++ b/src/algorithms/derivatives/mpo_derivatives.jl @@ -9,7 +9,13 @@ struct MPODerivativeOperator{L, O <: Tuple, R} <: DerivativeOperator rightenv::R end +struct MPOContractedDerivativeOperator{L, R, N} <: DerivativeOperator + leftenv::L + rightenv::R +end + Base.length(H::MPODerivativeOperator) = length(H.operators) +Base.length(::MPOContractedDerivativeOperator{L, R, N}) where {L, R, N} = N const MPO_C_Hamiltonian{L, R} = MPODerivativeOperator{L, Tuple{}, R} MPO_C_Hamiltonian(GL, GR) = MPODerivativeOperator(GL, (), GR) @@ -20,8 +26,20 @@ MPO_AC_Hamiltonian(GL, O, GR) = MPODerivativeOperator(GL, (O,), GR) const MPO_AC2_Hamiltonian{L, O₁, O₂, R} = MPODerivativeOperator{L, Tuple{O₁, O₂}, R} MPO_AC2_Hamiltonian(GL, O1, O2, GR) = MPODerivativeOperator(GL, (O1, O2), GR) +const MPO_Contracted_AC_Hamiltonian{L, R} = MPOContractedDerivativeOperator{L, R, 1} +MPO_Contracted_AC_Hamiltonian(GL::L, GR::R) where {L, R} = MPOContractedDerivativeOperator{L,R,1}(GL, GR) + +const MPO_Contracted_AC2_Hamiltonian{L, R} = MPOContractedDerivativeOperator{L, R, 2} +MPO_Contracted_AC2_Hamiltonian(GL::L, GR::R) where {L, R} = MPOContractedDerivativeOperator{L,R,2}(GL, GR) + # Constructors # ------------ +const _HAM_MPS_TYPES = Union{ + FiniteMPS{<:MPSTensor}, + WindowMPS{<:MPSTensor}, + InfiniteMPS{<:MPSTensor}, +} + function C_hamiltonian(site::Int, below, operator, above, envs) return MPO_C_Hamiltonian(leftenv(envs, site + 1, below), rightenv(envs, site, below)) end @@ -35,6 +53,27 @@ function AC2_hamiltonian(site::Int, below, operator, above, envs) leftenv(envs, site, below), O1, O2, rightenv(envs, site + 1, below) ) end +function AC_hamiltonian(site::Int, below::_HAM_MPS_TYPES, operator::MPO{<:MPOTensor}, above::_HAM_MPS_TYPES, envs) + O = operator[site] + GL = leftenv(envs, site, below) + return AC_hamiltonian(GL, O, rightenv(envs, site, below)) +end +function AC_hamiltonian(GL::MPSTensor, O::MPOTensor, GR::MPSTensor) + @plansor GLW[-1 -2 -3; -4 -5] ≔ GL[-1 1; -4] * O[1 -2; -5 -3] + return MPO_Contracted_AC_Hamiltonian(GLW, GR) +end +function AC2_hamiltonian(site::Int, below::_HAM_MPS_TYPES, operator::MPO{<:MPOTensor}, above::_HAM_MPS_TYPES, envs) + O1 = operator[site] + O2 = operator[site + 1] + GL = leftenv(envs, site, below) + GR = rightenv(envs, site + 1, below) + return AC2_hamiltonian(GL, O1, O2, GR) +end +function AC2_hamiltonian(GL::MPSTensor, O1::MPOTensor, O2::MPOTensor, GR::MPSTensor) + @plansor GLW[-1 -2 -3; -4 -5] ≔ GL[-1 1; -4] * O1[1 -2; -5 -3] + @plansor GWR[-1 -2 -3; -4 -5] ≔ O2[-3 -5; -2 1] * GR[-1 1; -4] + return MPO_Contracted_AC2_Hamiltonian(GLW, GWR) +end # Properties # ---------- @@ -50,6 +89,20 @@ function TensorKit.codomain(H::MPODerivativeOperator) V_o = prod(physicalspace, H.O; init = one(V_l)) return V_l ⊗ V_o ⊗ V_r end +function TensorKit.domain(H::MPOContractedDerivativeOperator) + V_l = right_virtualspace(H.leftenv) + V_r = left_virtualspace(H.rightenv) + ## TODO: How to deal with the H.O here? + V_o = prod(physicalspace, H.O; init = one(V_l)) + return V_l ⊗ V_o ⊗ V_r +end +function TensorKit.codomain(H::MPOContractedDerivativeOperator) + V_l = left_virtualspace(H.leftenv) + V_r = right_virtualspace(H.rightenv) + ## TODO: How to deal with the H.O here? + V_o = prod(physicalspace, H.O; init = one(V_l)) + return V_l ⊗ V_o ⊗ V_r +end # Actions # ------- @@ -105,3 +158,15 @@ function (h::MPO_AC2_Hamiltonian{<:MPSTensor, <:MPOTensor, <:MPOTensor, <:MPSTen h.operators[2][7 -6; 4 5] * τ[5 -5; 2 3] return y isa AbstractBlockTensorMap ? only(y) : y end +function (h::MPO_Contracted_AC_Hamiltonian)( + x::AbstractTensorMap{<:Any, <:Any, 2, 1} + ) + @plansor y[-1 -2; -3] ≔ h.leftenv[-1 -2 3; 1 2] * x[1 2; 4] * h.rightenv[4 3; -3] + return y isa AbstractBlockTensorMap ? only(y) : y +end +function (h::MPO_Contracted_AC2_Hamiltonian)( + x::MPOTensor + ) + @plansor y[-1 -2; -3 -4] ≔ h.leftenv[-1 -2 5; 1 2] * x[1 2; 3 4] * h.rightenv[3 4 5; -3 -4] + return y isa AbstractBlockTensorMap ? only(y) : y +end \ No newline at end of file diff --git a/src/algorithms/statmech/idmrg.jl b/src/algorithms/statmech/idmrg.jl index 5aa7b0b74..4f4204733 100644 --- a/src/algorithms/statmech/idmrg.jl +++ b/src/algorithms/statmech/idmrg.jl @@ -19,7 +19,7 @@ function leading_boundary( for row in 1:size(ψ, 1) ac = ψ.AC[row, col] (col == size(ψ, 2)) && (ac = copy(ac)) # needed in next sweep - ψ.AL[row, col], ψ.C[row, col] = leftorth!(ac) + ψ.AL[row, col], ψ.C[row, col] = qr_compact!(ac) end transfer_leftenv!(envs, ψ, operator, ψ, col + 1) @@ -31,7 +31,7 @@ function leading_boundary( _, ψ.AC[:, col] = fixedpoint(Hac, ψ.AC[:, col], :LM, alg_eigsolve) for row in 1:size(ψ, 1) - ψ.C[row, col - 1], temp = rightorth!(_transpose_tail(ψ.AC[row, col])) + ψ.C[row, col - 1], temp = lq_compact!(_transpose_tail(ψ.AC[row, col])) ψ.AR[row, col] = _transpose_front(temp) end diff --git a/test/algorithms.jl b/test/algorithms.jl index c1d0ad9e9..c91be0656 100644 --- a/test/algorithms.jl +++ b/test/algorithms.jl @@ -190,7 +190,7 @@ module TestAlgorithms @testset "LazySum FiniteMPS groundstate" verbose = true begin tol = 1.0e-8 - D = 15 + D = 32 atol = 1.0e-2 L = 10