diff --git a/src/algorithms/derivatives/hamiltonian_derivatives.jl b/src/algorithms/derivatives/hamiltonian_derivatives.jl index 58eecca43..3fa5ab263 100644 --- a/src/algorithms/derivatives/hamiltonian_derivatives.jl +++ b/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 diff --git a/src/algorithms/timestep/tdvp.jl b/src/algorithms/timestep/tdvp.jl index 84a6b6202..125c7a286 100644 --- a/src/algorithms/timestep/tdvp.jl +++ b/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 diff --git a/src/operators/jordanmpotensor.jl b/src/operators/jordanmpotensor.jl index 916809038..4166fc1fc 100644 --- a/src/operators/jordanmpotensor.jl +++ b/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 diff --git a/src/states/finitemps.jl b/src/states/finitemps.jl index d6c309638..3d306d3f9 100644 --- a/src/states/finitemps.jl +++ b/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) diff --git a/test/algorithms.jl b/test/algorithms.jl index a7606d1eb..caa6bdda0 100644 --- a/test/algorithms.jl +++ b/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]) diff --git a/test/other.jl b/test/other.jl index 0306f1888..1b9b84bb4 100644 --- a/test/other.jl +++ b/test/other.jl @@ -70,13 +70,19 @@ 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 @@ -84,13 +90,13 @@ end 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 diff --git a/test/setup.jl b/test/setup.jl index e5b007684..4689a2aff 100644 --- a/test/setup.jl +++ b/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) diff --git a/test/states.jl b/test/states.jl index 5c694da01..24d1690cb 100644 --- a/test/states.jl +++ b/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