diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl index 331d4ab36..44a9fbb79 100644 --- a/src/algorithms/toolbox.jl +++ b/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...) diff --git a/src/operators/jordanmpotensor.jl b/src/operators/jordanmpotensor.jl index d7133f8b6..d1c34e57e 100644 --- a/src/operators/jordanmpotensor.jl +++ b/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)) diff --git a/src/operators/mpohamiltonian.jl b/src/operators/mpohamiltonian.jl index 4d26b165e..6a8c6e055 100644 --- a/src/operators/mpohamiltonian.jl +++ b/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) diff --git a/test/groundstate/groundstate.jl b/test/groundstate/groundstate.jl index dd42b951b..816c04809 100644 --- a/test/groundstate/groundstate.jl +++ b/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 diff --git a/test/hamiltonian/infinite.jl b/test/hamiltonian/infinite.jl index 864cef5d3..f44ed42b4 100644 --- a/test/hamiltonian/infinite.jl +++ b/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