Fix variance SpaceMismatch for mixed-scalartype MPOHamiltonians (#524)

src/algorithms/toolbox.jl

@@ -115,16 +115,13 @@ function variance(
             state.AC[i], envs.GLs[i], H[i][:, :, :, end], envs.GRs[i][end]
         )
     end
-    lattice = physicalspace(H)
-    H_renormalized = InfiniteMPOHamiltonian(
-        lattice, i => e * id(storagetype(eltype(H)), lattice[i]) for (i, e) in enumerate(e_local)
-    )
-    return real(expectation_value(state, (H - H_renormalized)^2))
+    return real(expectation_value(state, (H - e_local)^2))
 end
 
 function variance(state::FiniteMPS, H::FiniteMPOHamiltonian, envs = environments(state, H, state))
-    H2 = H * H
-    return real(expectation_value(state, H2) - expectation_value(state, H, envs)^2)
+    E = expectation_value(state, H, envs)
+    λs = fill(E / length(H), length(H))
+    return real(expectation_value(state, (H - λs)^2))
 end
 
 function variance(state::FiniteQP, H::FiniteMPOHamiltonian, args...)

src/operators/jordanmpotensor.jl

@@ -162,6 +162,14 @@ end
 function jordanmpotensortype(::Type{O}) where {O <: AbstractTensorMap}
     return jordanmpotensortype(spacetype(O), storagetype(O))
 end
+function Base.promote_rule(::Type{O1}, ::Type{O2}) where {O1 <: JordanMPOTensor, O2 <: JordanMPOTensor}
+    O1 === O2 && return O1
+    spacetype(O1) === spacetype(O2) ||
+        throw(ArgumentError("cannot promote JordanMPOTensor types with different spacetypes"))
+    T = promote_type(scalartype(O1), scalartype(O2))
+    A = TensorKit.similarstoragetype(storagetype(O1), T)
+    return jordanmpotensortype(spacetype(O1), A)
+end
 function Base.similar(W::JordanMPOTensor, ::Type{T}) where {T <: Number}
     TE = TensorKit.similarstoragetype(TensorKit.storagetype(W), T)
     return jordanmpotensortype(spacetype(W), TE)(undef, space(W))

src/operators/mpohamiltonian.jl

@@ -903,6 +903,17 @@ function Base.convert(
     return InfiniteMPOHamiltonian(convert.(O1, parent(H)))
 end
 
+function Base.promote_rule(
+        ::Type{FiniteMPOHamiltonian{O1}}, ::Type{FiniteMPOHamiltonian{O2}}
+    ) where {O1 <: JordanMPOTensor, O2 <: JordanMPOTensor}
+    return FiniteMPOHamiltonian{promote_type(O1, O2)}
+end
+function Base.promote_rule(
+        ::Type{InfiniteMPOHamiltonian{O1}}, ::Type{InfiniteMPOHamiltonian{O2}}
+    ) where {O1 <: JordanMPOTensor, O2 <: JordanMPOTensor}
+    return InfiniteMPOHamiltonian{promote_type(O1, O2)}
+end
+
 function add_physical_charge(H::MPOHamiltonian, charges::AbstractVector{<:Sector})
     W = map(add_physical_charge, parent(H), charges)
     if isfinite(H)
@@ -930,8 +941,9 @@ Base.circshift(H::InfiniteMPOHamiltonian, shift::Integer) = InfiniteMPOHamiltoni
 # Linear Algebra
 # --------------
 function Base.:+(
-        H₁::FiniteMPOHamiltonian{O}, H₂::FiniteMPOHamiltonian{O}
-    ) where {O <: JordanMPOTensor}
+        H₁::FiniteMPOHamiltonian{O1}, H₂::FiniteMPOHamiltonian{O2}
+    ) where {O1 <: JordanMPOTensor, O2 <: JordanMPOTensor}
+    O1 === O2 || return +(promote(H₁, H₂)...)
     N = check_length(H₁, H₂)
     H = similar(parent(H₁))
     # same as rightunitspace (asserted within construction FiniteMPOHamiltonian)
@@ -954,9 +966,10 @@ function Base.:+(
     return FiniteMPOHamiltonian(H)
 end
 function Base.:+(
-        H₁::InfiniteMPOHamiltonian{O},
-        H₂::InfiniteMPOHamiltonian{O}
-    ) where {O <: JordanMPOTensor}
+        H₁::InfiniteMPOHamiltonian{O1},
+        H₂::InfiniteMPOHamiltonian{O2}
+    ) where {O1 <: JordanMPOTensor, O2 <: JordanMPOTensor}
+    O1 === O2 || return +(promote(H₁, H₂)...)
     N = check_length(H₁, H₂)
     H = similar(parent(H₁))
     # same as rightunitspace (asserted within construction of InfiniteMPOHamiltonian)
@@ -979,7 +992,7 @@ end
 function Base.:+(H::FiniteMPOHamiltonian, λs::AbstractVector{<:Number})
     check_length(H, λs)
     lattice = [physicalspace(H, i) for i in 1:length(H)]
-    M = storagetype(H)
+    M = TensorKit.similarstoragetype(storagetype(H), promote_type(scalartype(H), eltype(λs)))
     Hλ = FiniteMPOHamiltonian(
         lattice,
         i => scale!(id(M, lattice[i]), λs[i]) for i in 1:length(H)
@@ -989,7 +1002,7 @@ end
 function Base.:+(H::InfiniteMPOHamiltonian, λs::AbstractVector{<:Number})
     check_length(H, λs)
     lattice = [physicalspace(H, i) for i in 1:length(H)]
-    M = storagetype(H)
+    M = TensorKit.similarstoragetype(storagetype(H), promote_type(scalartype(H), eltype(λs)))
     Hλ = InfiniteMPOHamiltonian(
         lattice,
         i => scale!(id(M, lattice[i]), λs[i]) for i in 1:length(H)

test/groundstate/groundstate.jl

@@ -203,11 +203,9 @@ end
     ψ = InfiniteMPS(ℙ^2, ℙ^D)
     v₀ = variance(ψ, H_ref)
 
-    # VUMPS spawns over the unit cell, so it is run under both schedulers: which allocator serves
-    # the local updates follows from the scheduler, and must not change any number
-    # NOTE: the starting state is built from scratch rather than from this block's `ψ`. The testsets
-    # around this one rebind `ψ` to their own (already repeated) result, and a `@testset for` body is
-    # a single scope, so reading it here would feed a length-3 state back into `repeat`.
+    H_ref_realT = force_planar(transverse_field_ising(Float64; g))
+    @test variance(ψ, H_ref_realT) ≈ v₀ atol = 1.0e-10
+
     @testset "VUMPS (unit cell $unit_cell_size, $schedname)" for unit_cell_size in [1, 3],
             (schedname, scheduler) in SCHEDULERS
 

test/hamiltonian/infinite.jl

@@ -51,4 +51,20 @@ end
 
     H4 = H1 + H3
     @test real(expectation_value(ψ2, H4)) >= 0
+
+    O1_real = project_hermitian!(randn(Float64, pspace, pspace))
+    O2_real = project_hermitian!(randn(Float64, pspace ⊗ pspace, pspace ⊗ pspace))
+    Hreal_long_range = InfiniteMPOHamiltonian(O1_real) + InfiniteMPOHamiltonian(O2_real)
+    Hcomplex_local = InfiniteMPOHamiltonian(operators[1])
+
+    @test scalartype(Hreal_long_range) == Float64
+    @test scalartype(Hcomplex_local) == ComplexF64
+    for H in (
+            Hreal_long_range + Hcomplex_local, Hcomplex_local + Hreal_long_range,
+            Hreal_long_range - Hcomplex_local, Hcomplex_local - Hreal_long_range,
+        )
+        @test scalartype(H) == ComplexF64
+    end
+    @test expectation_value(ψ1, Hreal_long_range + Hcomplex_local) ≈
+        expectation_value(ψ1, Hreal_long_range) + expectation_value(ψ1, Hcomplex_local) atol = 1.0e-10
 end