Fix derivatives at the edges of finite operators (#296)

* Fix edge logic for AC_hamiltonian

* Fix edge logic for AC2_hamiltonian

* Fix conversion logic for `SparseBlockTensorMap`

* Better type detection in `complex`

* Handle real/complex environments

* update storagetype determination

* Additional tests

* fix tests?

* Fix tests!

* Fix braille checks

* Update tests for sparsity improvements in BlockTensorKit

Author: lkdvos

src/algorithms/derivatives/hamiltonian_derivatives.jl

@@ -13,10 +13,11 @@ struct JordanMPO_AC_Hamiltonian{O1,O2,O3} <: DerivativeOperator
     A::O3                # continuing
     function JordanMPO_AC_Hamiltonian(onsite, not_started, finished, starting, ending,
                                       continuing)
-        tensor = coalesce(onsite, not_started, finished, starting, ending)
-        ismissing(tensor) && throw(ArgumentError("unable to determine type"))
-        S = spacetype(tensor)
-        M = storagetype(tensor)
+        # 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,
@@ -46,10 +47,11 @@ struct JordanMPO_AC2_Hamiltonian{O1,O2,O3,O4} <: DerivativeOperator
     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)
-        tensor = coalesce(II, IC, ID, CB, CA, AB, AA, BE, DE, EE)
-        ismissing(tensor) && throw(ArgumentError("unable to determine type"))
-        S = spacetype(tensor)
-        M = storagetype(tensor)
+        # 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)
@@ -107,7 +109,9 @@ function AC_hamiltonian(site::Int, below::_HAM_MPS_TYPES,
     end
 
     # not_started
-    if (!isfinite(operator) || site < length(operator)) && !ismissing(starting)
+    if isfinite(operator) && site == length(operator)
+        not_started = missing
+    elseif !ismissing(starting)
         I = id(storagetype(GR[1]), physicalspace(W))
         @plansor starting[-1 -2; -3 -4] += I[-1; -3] * removeunit(GR[1], 2)[-4; -2]
         not_started = missing
@@ -116,7 +120,9 @@ function AC_hamiltonian(site::Int, below::_HAM_MPS_TYPES,
     end
 
     # finished
-    if (!isfinite(operator) || site > 1) && !ismissing(ending)
+    if isfinite(operator) && site == 1
+        finished = missing
+    elseif !ismissing(ending)
         I = id(storagetype(GL[end]), physicalspace(W))
         @plansor ending[-1 -2; -3 -4] += removeunit(GL[end], 2)[-1; -3] * I[-2; -4]
         finished = missing
@@ -241,7 +247,9 @@ function AC2_hamiltonian(site::Int, below::_HAM_MPS_TYPES,
     end
 
     # finished
-    if !ismissing(IC)
+    if isfinite(operator) && site + 1 == length(operator)
+        II = missing
+    elseif !ismissing(IC)
         I = id(storagetype(GR[1]), physicalspace(W2))
         @plansor IC[-1 -2; -3 -4] += I[-1; -3] * removeunit(GR[1], 2)[-4; -2]
         II = missing
@@ -254,7 +262,9 @@ function AC2_hamiltonian(site::Int, below::_HAM_MPS_TYPES,
     end
 
     # unstarted
-    if !ismissing(BE)
+    if isfinite(operator) && site == 1
+        EE = missing
+    elseif !ismissing(BE)
         I = id(storagetype(GL[end]), physicalspace(W1))
         @plansor BE[-1 -2; -3 -4] += removeunit(GL[end], 2)[-1; -3] * I[-2; -4]
         EE = missing

src/algorithms/timestep/tdvp.jl

@@ -178,7 +178,18 @@ end
 
 #copying version
 function timestep(ψ::AbstractFiniteMPS, H, time::Number, timestep::Number,
-                  alg::Union{TDVP,TDVP2}, envs::AbstractMPSEnvironments=environments(ψ, H);
+                  alg::Union{TDVP,TDVP2}, envs::AbstractMPSEnvironments...;
                   kwargs...)
-    return timestep!(copy(complex(ψ)), H, time, timestep, alg, envs; kwargs...)
+    isreal = scalartype(ψ) <: Real
+    ψ′ = isreal ? complex(ψ) : copy(ψ)
+    if length(envs) != 0 && isreal
+        @warn "Currently cannot reuse real environments for complex evolution"
+        envs′ = environments(ψ′, H)
+    elseif length(envs) == 1
+        envs′ = only(envs)
+    else
+        @assert length(envs) == 0 "Invalid signature"
+        envs′ = environments(ψ′, H)
+    end
+    return timestep!(ψ′, H, time, timestep, alg, envs′; kwargs...)
 end

src/operators/jordanmpotensor.jl

@@ -155,10 +155,10 @@ function SparseBlockTensorMap(W::JordanMPOTensor)
     τ = BraidingTensor{scalartype(W)}(eachspace(W)[1])
     W′ = SparseBlockTensorMap{AbstractTensorMap{scalartype(W),spacetype(W),2,2}}(undef_blocks,
                                                                                  space(W))
-    if size(W, 1) > 1
+    if size(W, 4) > 1
         W′[1, 1, 1, 1] = τ
     end
-    if size(W, 4) > 1
+    if size(W, 1) > 1
         W′[end, 1, 1, end] = τ
     end
 
@@ -174,7 +174,9 @@ function SparseBlockTensorMap(W::JordanMPOTensor)
     for (I, v) in nonzero_pairs(W.C)
         W′[1, I + Ic] = insertleftunit(v, 1)
     end
-    W′[1, 1, 1, end] = insertrightunit(insertleftunit(only(W.D), 1), 3)
+    if nonzero_length(W.D) > 0
+        W′[1, 1, 1, end] = insertrightunit(insertleftunit(only(W.D), 1), 3)
+    end
 
     return W′
 end

src/states/finitemps.jl

@@ -325,10 +325,12 @@ function Base.complex(mps::FiniteMPS)
     ARs = _complex_if_not_missing.(mps.ARs)
     Cs = _complex_if_not_missing.(mps.Cs)
     ACs = _complex_if_not_missing.(mps.ACs)
-    return FiniteMPS(collect(Union{Missing,eltype(ALs)}, ALs),
-                     collect(Union{Missing,eltype(ARs)}, ARs),
-                     collect(Union{Missing,eltype(ACs)}, ACs),
-                     collect(Union{Missing,eltype(Cs)}, Cs))
+    TA = Base.promote_op(complex, site_type(mps))
+    TB = Base.promote_op(complex, bond_type(mps))
+    return FiniteMPS(collect(Union{Missing,TA}, ALs),
+                     collect(Union{Missing,TA}, ARs),
+                     collect(Union{Missing,TA}, ACs),
+                     collect(Union{Missing,TB}, Cs))
 end
 
 @inline function Base.getindex(ψ::FiniteMPS, I::AbstractUnitRange)

test/algorithms.jl

@@ -326,14 +326,17 @@ end
     algs = [TDVP(), TDVP2(; trscheme=truncdim(10))]
     L = 10
 
-    H = force_planar(heisenberg_XXX(; spin=1 // 2, L))
-    ψ₀ = FiniteMPS(L, ℙ^2, ℙ^1)
+    H = force_planar(heisenberg_XXX(Trivial, Float64; spin=1 // 2, L))
+    ψ = FiniteMPS(rand, Float64, L, ℙ^2, ℙ^4)
+    E = expectation_value(ψ, H)
+    ψ₀, = find_groundstate(ψ, H)
     E₀ = expectation_value(ψ₀, H)
 
     @testset "Finite $(alg isa TDVP ? "TDVP" : "TDVP2")" for alg in algs
-        ψ, envs = timestep(ψ₀, H, 0.0, dt, alg)
-        E = expectation_value(ψ, H, envs)
-        @test E₀ ≈ E atol = 1e-2
+        ψ1, envs = timestep(ψ₀, H, 0.0, dt, alg)
+        E1 = expectation_value(ψ1, H, envs)
+        @test E₀ ≈ E1 atol = 1e-2
+        @test dot(ψ1, ψ₀) ≈ exp(im * dt * E₀) atol = 1e-4
     end
 
     Hlazy = LazySum([3 * H, 1.55 * H, -0.1 * H])

test/other.jl

@@ -70,27 +70,33 @@ end
     buffer = IOBuffer()
     braille(buffer, H)
     output = String(take!(buffer))
-    check = "... 🭻⎡⠉⢀⎤🭻 ...\n     ⎣⠀⢀⎦ \n"
+    check = """
+            ... 🭻⎡⠉⢈⎤🭻 ...
+                 ⎣⠀⢀⎦ 
+            """
     @test output == check
 
     O = make_time_mpo(H, 1.0, TaylorCluster(3, false, false))
     braille(buffer, O)
     output = String(take!(buffer))
-    check = "... 🭻⎡⡏⠉⠒⠔⎤🭻 ...\n     ⎣⡇⠀⠀⡂⎦ \n"
+    check = """
+            ... 🭻⎡⡏⠉⠛⠟⎤🭻 ...
+                 ⎣⡇⠀⠀⡂⎦ 
+            """
     @test output == check
 
     # Finite Hamiltonians and MPOs
     # ----------------------------
     H = transverse_field_ising(; L=4)
     braille(buffer, H)
     output = String(take!(buffer))
-    check = " ⎡⠉⠀⎤🭻🭻⎡⠉⢀⎤🭻🭻⎡⠉⢀⎤🭻🭻⎡⡀⠀⎤ \n ⎣⠀⠀⎦  ⎣⠀⢀⎦  ⎣⠀⢀⎦  ⎣⡀⠀⎦ \n"
+    check = " ⎡⠉⠈⎤🭻🭻⎡⠉⢈⎤🭻🭻⎡⠉⢈⎤🭻🭻⎡⡁⠀⎤ \n ⎣⠀⠀⎦  ⎣⠀⢀⎦  ⎣⠀⢀⎦  ⎣⡀⠀⎦ \n"
     @test output == check
 
     O = make_time_mpo(H, 1.0, TaylorCluster(3, false, false))
     braille(buffer, O)
     output = String(take!(buffer))
-    check = " ⎡⠉⠉⎤🭻🭻⎡⠍⠉⠤⠠⎤🭻🭻⎡⡏⠉⠒⠔⎤🭻🭻⎡⡇⠀⎤ \n ⎣⠀⠀⎦  ⎣⡂⠀⠀⠂⎦  ⎣⡇⠀⠀⡂⎦  ⎣⡇⠀⎦ \n"
+    check = " ⎡⠉⠉⠉⠉⎤🭻🭻⎡⡏⠉⠛⠟⎤🭻🭻⎡⡏⠉⠛⠟⎤🭻🭻⎡⡇⠀⎤ \n ⎣⠀⠀⠀⠀⎦  ⎣⡇⠀⠀⡂⎦  ⎣⡇⠀⠀⡂⎦  ⎣⡇⠀⎦ \n"
     @test output == check
 end
 end

test/setup.jl

@@ -123,33 +123,20 @@ function S_zz(::Type{Trivial}=Trivial, ::Type{T}=ComplexF64; spin=1 // 2) where
     return S_z(Trivial, T; spin) ⊗ S_z(Trivial, T; spin)
 end
 
-function transverse_field_ising(; g=1.0, L=Inf)
-    X = S_x(; spin=1 // 2)
-    ZZ = S_zz(; spin=1 // 2)
-    E = TensorMap(ComplexF64[1 0; 0 1], ℂ^2 ← ℂ^2)
+function transverse_field_ising(::Type{T}=ComplexF64; g=1.0, L=Inf) where {T<:Number}
+    X = S_x(Trivial, T; spin=1 // 2)
+    ZZ = S_zz(Trivial, T; spin=1 // 2)
 
-    # lattice = L == Inf ? PeriodicVector([ℂ^2]) : fill(ℂ^2, L)
     if L == Inf
         lattice = PeriodicArray([ℂ^2])
-        return InfiniteMPOHamiltonian(lattice,
-                                      (i, i + 1) => -(ZZ + (g / 2) * (X ⊗ E + E ⊗ X))
-                                      for i in 1:1)
-        # return MPOHamiltonian(-ZZ - (g / 2) * (X ⊗ E + E ⊗ X))
+        H₁ = InfiniteMPOHamiltonian(lattice, i => -g * X for i in 1:1)
+        H₂ = InfiniteMPOHamiltonian(lattice, (i, i + 1) => -ZZ for i in 1:1)
     else
         lattice = fill(ℂ^2, L)
-        return FiniteMPOHamiltonian(lattice,
-                                    (i, i + 1) => -(ZZ + (g / 2) * (X ⊗ E + E ⊗ X))
-                                    for i in 1:(L - 1)) #+
-        # FiniteMPOHamiltonian(lattice, (i,) => -g * X for i in 1:L)
+        H₁ = FiniteMPOHamiltonian(lattice, i => -g * X for i in 1:L)
+        H₂ = FiniteMPOHamiltonian(lattice, (i, i + 1) => -ZZ for i in 1:(L - 1))
     end
-
-    H = S_zz(; spin=1 // 2) + (g / 2) * (X ⊗ E + E ⊗ X)
-    return if L == Inf
-        MPOHamiltonian(H)
-    else
-        FiniteMPOHamiltonian(fill(ℂ^2, L), (i, i + 1) => H for i in 1:(L - 1))
-    end
-    return MPOHamiltonian(-H)
+    return H₁ + H₂
 end
 
 function heisenberg_XXX(::Type{SU2Irrep}; spin=1, L=Inf)
@@ -169,8 +156,9 @@ function heisenberg_XXX(::Type{SU2Irrep}; spin=1, L=Inf)
     end
 end
 
-function heisenberg_XXX(; spin=1, L=Inf)
-    h = ones(ComplexF64, SU2Space(spin => 1)^2 ← SU2Space(spin => 1)^2)
+function heisenberg_XXX(::Type{Trivial}=Trivial, ::Type{T}=ComplexF64; spin=1,
+                        L=Inf) where {T<:Number}
+    h = ones(T, SU2Space(spin => 1)^2 ← SU2Space(spin => 1)^2)
     for (c, b) in blocks(h)
         S = (dim(c) - 1) / 2
         b .= S * (S + 1) / 2 - spin * (spin + 1)

test/states.jl

@@ -168,18 +168,18 @@ end
     @test 9 * 9 ≈ norm(window)^2
     normalize!(window)
 
-    e1 = expectation_value(window, (1, 2) => O)
+    e1 = expectation_value(window, (2, 3) => O)
 
     window, envs, _ = find_groundstate(window, ham, DMRG(; verbosity=0))
 
-    e2 = expectation_value(window, (1, 2) => O)
+    e2 = expectation_value(window, (2, 3) => O)
 
     @test real(e2) ≤ real(e1)
 
     window, envs = timestep(window, ham, 0.1, 0.0, TDVP2(; trscheme=truncdim(20)), envs)
     window, envs = timestep(window, ham, 0.1, 0.0, TDVP(), envs)
 
-    e3 = expectation_value(window, (1, 2) => O)
+    e3 = expectation_value(window, (2, 3) => O)
 
     @test e2 ≈ e3 atol = 1e-4
 end