diff --git a/src/algorithms/contractions/local_patch/network_expr.jl b/src/algorithms/contractions/local_patch/network_expr.jl index 6120fa3b3..5f77d08bc 100644 --- a/src/algorithms/contractions/local_patch/network_expr.jl +++ b/src/algorithms/contractions/local_patch/network_expr.jl @@ -260,3 +260,40 @@ function operator_contraction_expr(::Type{<:AbstractTensorMap}, nsites) ), ] end + +# The MPO bond is small - the operator's Schmidt rank, e.g. 3 for Heisenberg XYZ - so it is +# labelled as a physical rather than a virtual dimension, which gives `@autoopt` the right +# order of magnitude when it searches for a contraction order. +mpolabel(args...) = physicallabel(:mpo, args...) + +function operator_contraction_expr(::Type{<:MPOTerm}, nsites) + # one factor per site, linked by a chain of bonds running left to right, rank-3 at the + # ends of the chain and rank-4 in the bulk: + # W₁ : bra₁ ← ket₁ ⊗ b₁ + # Wᵢ : bᵢ₋₁ ⊗ braᵢ ← ketᵢ ⊗ bᵢ + # W_N : b_{N-1} ⊗ bra_N ← ket_N + return map(1:nsites) do i + bra, ket = physicallabel(:O, 2, i), physicallabel(:O, 1, i) + out, ins = if nsites == 1 + (bra,), (ket,) + elseif i == 1 + (bra,), (ket, mpolabel(1)) + elseif i == nsites + (mpolabel(nsites - 1), bra), (ket,) + else + (mpolabel(i - 1), bra), (ket, mpolabel(i)) + end + return tensorexpr(:(operator[$i]), out, ins) + end +end + +# A tensor product term is the special case in which every bond has dimension 1, so its +# factors carry no bonds at all: each simply sits between the ket and bra index of its site. +# This method is more specific than the `MPOTerm` one above, so it takes precedence. +function operator_contraction_expr(::Type{<:TensorProductTerm}, nsites) + return map(1:nsites) do i + return tensorexpr( + :(operator[$i]), (physicallabel(:O, 2, i),), (physicallabel(:O, 1, i),) + ) + end +end diff --git a/src/algorithms/expectation_value/expectation_value.jl b/src/algorithms/expectation_value/expectation_value.jl index d2ea498f6..7c83aac90 100644 --- a/src/algorithms/expectation_value/expectation_value.jl +++ b/src/algorithms/expectation_value/expectation_value.jl @@ -54,6 +54,26 @@ function local_expectation_value(inds, state, operator::AbstractTensorMap, env) return trmul(operator, ρ) end +""" +$(SIGNATURES) + +Compute the contribution of an [`MPOTerm`](@ref), including a [`TensorProductTerm`](@ref), given as one tensor per site in `inds`. +For a PEPS or purified PEPO, compute ``⟨bra|O|ket⟩ / ⟨bra|ket⟩``; for a single-layer density matrix PEPO, compute ``tr(O * state) / tr(state)``. +Insert the factors into the patch contraction directly and divide by the local norm of the same patch. +""" +function local_expectation_value(inds, bra, operator::MPOTerm, ket, env) + return contract_local_operator(inds, operator, ket, bra, env) / + contract_local_norm(inds, ket, bra, env) +end +function local_expectation_value(inds, state, operator::MPOTerm, env) + return contract_local_operator(inds, operator, state, env) / + contract_local_norm(inds, state, env) +end + +# TODO: Implement boundary_contraction_expr for BPEnv to contract its messages directly. +local_expectation_value(inds, bra, operator::MPOTerm, ket, env::BPEnv) = + local_expectation_value(inds, bra, operator, ket, CTMRGEnv(env)) + # Expectation value of a local partition function tensor # ------------------------------------------------------ diff --git a/src/algorithms/expectation_value/patch_contractions.jl b/src/algorithms/expectation_value/patch_contractions.jl index 324d2ebb5..a1996743b 100644 --- a/src/algorithms/expectation_value/patch_contractions.jl +++ b/src/algorithms/expectation_value/patch_contractions.jl @@ -19,6 +19,24 @@ function contract_local_operator( static_inds = Tuple(Val.(inds)) return _contract_local_operator(static_inds, O, (ket, bra), env) end +function contract_local_operator( + inds::Vector{CartesianIndex{2}}, O::MPOTerm, + ket::InfinitePEPS, bra::InfinitePEPS, env + ) + length(inds) == length(O) || + throw(ArgumentError("Got $(length(inds)) sites but $(length(O)) MPO factors.")) + static_inds = Tuple(Val.(inds)) + return _contract_local_operator(static_inds, O, (ket, bra), env) +end +function contract_local_operator( + inds::Vector{CartesianIndex{2}}, O::TensorProductTerm, + ket::InfinitePEPS, bra::InfinitePEPS, env + ) + length(inds) == length(O) || + throw(ArgumentError("Got $(length(inds)) sites but $(length(O)) product factors.")) + static_inds = Tuple(Val.(inds)) + return _contract_local_operator(static_inds, O, (ket, bra), env) +end function contract_local_operator( inds::Vector{CartesianIndex{2}}, O, state::InfinitePEPO, env ) diff --git a/src/operators/localoperator.jl b/src/operators/localoperator.jl index c35e5ddb9..e4ac38180 100644 --- a/src/operators/localoperator.jl +++ b/src/operators/localoperator.jl @@ -50,10 +50,8 @@ end # Default to Any for eltype: needs to be abstract anyways so not that much to gain LocalOperator(lattice, terms) = LocalOperator{Any}(lattice, terms) LocalOperator(lattice, terms::Pair...) = LocalOperator(lattice, terms) -# TODO: add terms beyond AbstractTensorMap -# e.g. tensor product of 1-site operators, MPOs -add_term!(operator::LocalOperator, inds::Tuple, term::AbstractTensorMap) = add_term!(operator, collect(inds), term) -add_term!(operator::LocalOperator, inds::Vector, term::AbstractTensorMap) = add_term!(operator, map(CartesianIndex{2}, inds), term) +add_term!(operator::LocalOperator, inds::Tuple, term::AbstractTensorMap; kwargs...) = add_term!(operator, collect(inds), term; kwargs...) +add_term!(operator::LocalOperator, inds::AbstractVector, term::AbstractTensorMap; kwargs...) = add_term!(operator, CartesianIndex{2}[CartesianIndex{2}(ind) for ind in inds], term; kwargs...) function add_term!( operator::LocalOperator, inds::Vector{CartesianIndex{2}}, term::AbstractTensorMap; atol = zero(real(scalartype(term))), @@ -76,9 +74,11 @@ function add_term!( end # translate coordinates - _shift_into_unitcell!(inds, size(operator)) + inds = _shift_into_unitcell!(copy(inds), size(operator)) if haskey(operator.terms, inds) + operator.terms[inds] isa MPOTerm && + throw(ArgumentError("Accumulating terms with the same sites is not implemented for MPO or tensor product terms.")) operator.terms[inds] = VI.add!!(operator.terms[inds], term) else operator.terms[inds] = term @@ -87,6 +87,171 @@ function add_term!( return operator end +# Tensor product terms +# -------------------- +# A term given as one rank-2 operator per site, i.e. an explicit tensor product which is +# never actually formed. + +""" + TensorProductTerm + +A single term of a [`LocalOperator`](@ref) given as a tensor product of one rank-2 operator +per site acted on, rather than as the rank-`2N` tensor product itself. + +Only operators whose terms are individually tensor products can be represented this way. In +particular several such terms acting on the same set of sites cannot be combined, since a +sum of tensor products is not itself a tensor product. + +The factors are rank-2, carrying no bond indices, which is what distinguishes this from the +more general [`MPOTerm`](@ref): both are vectors of tensors, but a tensor product term is the +narrower type and therefore wins on dispatch wherever it applies. +""" +const TensorProductTerm{T} = AbstractVector{T} where {T <: AbstractTensorMap{<:Number, <:IndexSpace, 1, 1}} + +add_term!(operator::LocalOperator, inds::Tuple, term::TensorProductTerm; kwargs...) = + add_term!(operator, collect(inds), term; kwargs...) +add_term!(operator::LocalOperator, inds::AbstractVector, term::TensorProductTerm; kwargs...) = + add_term!(operator, CartesianIndex{2}[CartesianIndex{2}(ind) for ind in inds], term; kwargs...) +function add_term!( + operator::LocalOperator, inds::Vector{CartesianIndex{2}}, term::TensorProductTerm; + atol = 0, + ) + # input checks + isempty(term) && throw(ArgumentError("A tensor product term should contain at least one tensor.")) + length(inds) == length(term) || + throw(ArgumentError("Incompatible number of indices and tensor product factors")) + allunique(inds) || throw(ArgumentError("`inds` should not contain repeated coordinates.")) + for (i, ind) in enumerate(inds) + numin(term[i]) == numout(term[i]) == 1 || + throw(ArgumentError("Tensor product factors should be single-site operators")) + ind_translated = CartesianIndex(mod1.(Tuple(ind), size(operator))) + physicalspace(operator, ind_translated) == domain(term[i])[1] == codomain(term[i])[1] || + throw(SpaceMismatch("Incompatible physical spaces")) + end + prod(norm, term) <= atol && return operator # skip adding negligible terms + + # permute input, which for a product just reorders the factors along with their sites + if !issorted(inds) + I = sortperm(inds) + inds = inds[I] + term = term[I] + end + + # translate coordinates + inds = _shift_into_unitcell!(copy(inds), size(operator)) + + # a sum of tensor products is not a tensor product, so terms cannot be accumulated here + haskey(operator.terms, inds) && throw( + ArgumentError( + "A term acting on $inds is already present. Tensor product terms cannot be \ + summed, since a sum of tensor products is not itself a tensor product." + ) + ) + operator.terms[inds] = collect(term) + + return operator +end + +# MPO terms +# --------- +# A term given as a chain of MPO tensors, one per site, linked by virtual bonds which are +# contracted directly between neighbours rather than being fused into the PEPS bonds. This +# generalizes a tensor product term, which is the special case of bond dimension 1. + +""" + MPOTerm{T} + +A single term of a [`LocalOperator`](@ref) given as a matrix product operator with one +tensor per site acted on, rather than as the dense rank-`2N` tensor. + +The factors are ordered along the chain and follow the usual MPO convention, rank-3 at the ends and rank-4 in the bulk: + + W₁ : P₁ ← P₁ ⊗ B₁ + Wᵢ : Bᵢ₋₁ ⊗ Pᵢ ← Pᵢ ⊗ Bᵢ + W_N : B_{N-1} ⊗ P_N ← P_N + +so the physical index in the codomain is the bra index and the one in the domain is the ket +index, matching the convention of a dense term, and the bonds run left to right. This is what +[`gate_to_mpo`](@ref) produces, which is the way to obtain an `MPOTerm` from a dense operator. +A single-site term is just a rank-2 operator. + +A [`TensorProductTerm`](@ref) is the special case in which every bond is trivial, so its +factors carry no bond indices and are rank-2 throughout. That is exactly what separates the +two on dispatch: a vector of rank-2 operators is a tensor product term, anything else is an +MPO. +""" +const MPOTerm{T} = AbstractVector{T} where {T <: AbstractTensorMap} + + +""" +Validate an ordered MPO's sites, tensor partitions, physical spaces, and adjacent bond spaces. +""" +function _validate_mpo_term(inds, term::MPOTerm, lattice::AbstractMatrix{<:ElementarySpace}) + isempty(term) && throw(ArgumentError("An MPO term should contain at least one tensor.")) + length(inds) == length(term) || + throw(ArgumentError("Incompatible number of indices and MPO factors")) + allunique(inds) || throw(ArgumentError("`inds` should not contain repeated coordinates.")) + n = length(inds) + for (i, ind) in enumerate(inds) + # a factor carries its physical pair plus a bond towards each neighbour it has, so it + # is rank-2 for a lone site, rank-3 at the ends of a chain and rank-4 in the bulk + nout = (n == 1 || i == 1) ? 1 : 2 + nin = (n == 1 || i == n) ? 1 : 2 + (numout(term[i]) == nout && numin(term[i]) == nin) || throw( + ArgumentError( + "MPO factor $i of $n should have $nout index(es) out and $nin in, got \ + $(numout(term[i])) and $(numin(term[i]))." + ) + ) + # the bra index is the physical one in the codomain: last for a bulk factor, which + # carries the incoming bond first, and only for an end factor + bra = codomain(term[i])[nout] + ind_translated = CartesianIndex(mod1.(Tuple(ind), size(lattice))) + lattice[ind_translated] == bra == domain(term[i])[1] || + throw(SpaceMismatch("Incompatible physical spaces")) + end + for i in 1:(n - 1) + domain(term[i])[numin(term[i])] == codomain(term[i + 1])[1] || + throw(SpaceMismatch("Incompatible MPO bond spaces between tensors $i and $(i + 1).")) + end + return nothing +end + +add_term!(operator::LocalOperator, inds::Tuple, term::MPOTerm) = + add_term!(operator, collect(inds), term) +add_term!(operator::LocalOperator, inds::AbstractVector, term::MPOTerm) = + add_term!(operator, CartesianIndex{2}[CartesianIndex{2}(ind) for ind in inds], term) +function add_term!( + operator::LocalOperator, inds::Vector{CartesianIndex{2}}, term::MPOTerm + ) + _validate_mpo_term(inds, term, physicalspace(operator)) + _local_term_iszero(term) && return operator + + # NOTE: `inds` is deliberately *not* sorted here, unlike for dense and tensor product + # terms, since permuting MPOs is not straightforward. + + # translate coordinates + inds = _shift_into_unitcell!(copy(inds), size(operator)) + + # as for tensor products, a sum of MPOs of fixed bond dimension is not one of the same + # bond dimension, so terms are not accumulated here + haskey(operator.terms, inds) && throw( + ArgumentError( + "A term acting on $inds is already present. MPO terms cannot be summed in \ + place; combine the operators before splitting them." + ) + ) + operator.terms[inds] = collect(term) + + return operator +end + +""" +Detect an identically zero dense tensor or an MPO containing a zero factor. +""" +_local_term_iszero(term::AbstractTensorMap) = iszero(norm(term)) +_local_term_iszero(term::MPOTerm) = any(_local_term_iszero, term) + """ checklattice(Bool, args...) @@ -154,16 +319,37 @@ Base.eltype(::Type{LocalOperator{O, S}}) where {O, S} = O # Real and imaginary part # ----------------------- function Base.real(O::LocalOperator) + any(term -> term isa MPOTerm, values(O.terms)) && + throw(ArgumentError("Taking the real part is not implemented for MPO or tensor product terms.")) return LocalOperator(O.lattice, (sites => real(op) for (sites, op) in O.terms)...) end function Base.imag(O::LocalOperator) + any(term -> term isa MPOTerm, values(O.terms)) && + throw(ArgumentError("Taking the imaginary part is not implemented for MPO or tensor product terms.")) return LocalOperator(O.lattice, (sites => imag(op) for (sites, op) in O.terms)...) end # Linear Algebra # -------------- -Base.:*(α::Number, O::LocalOperator) = - LocalOperator(physicalspace(O), inds => α * operator for (inds, operator) in O.terms) +""" +Scale a single term of a [`LocalOperator`](@ref) by `α`. + +An [`MPOTerm`](@ref), including a [`TensorProductTerm`](@ref), is scaled by scaling its first factor only. +The factor container may widen its element type when the scalar type changes, while retaining tensor-product dispatch. +""" +_scale_local_term(term, α::Number) = α * term +function _scale_local_term(term::MPOTerm, α::Number) + return AbstractTensorMap[i == 1 ? α * tensor : tensor for (i, tensor) in enumerate(term)] +end +function _scale_local_term(term::TensorProductTerm, α::Number) + return AbstractTensorMap{<:Number, <:IndexSpace, 1, 1}[ + i == 1 ? α * tensor : tensor for (i, tensor) in enumerate(term) + ] +end + +Base.:*(α::Number, O::LocalOperator) = LocalOperator( + physicalspace(O), inds => _scale_local_term(operator, α) for (inds, operator) in O.terms +) Base.:*(O::LocalOperator, α::Number) = α * O Base.:/(O::LocalOperator, α::Number) = O * inv(α) @@ -171,7 +357,16 @@ Base.:\(α::Number, O::LocalOperator) = inv(α) * O function Base.:+(O1::LocalOperator, O2::LocalOperator) checklattice(O1, O2) - return LocalOperator(physicalspace(O1), mergewith(VI.add, O1.terms, O2.terms)) + return LocalOperator(physicalspace(O1), mergewith(_add_local_terms, O1.terms, O2.terms)) +end + +""" +Accumulate dense terms while rejecting addition involving an MPO or tensor product at the same sites. +""" +function _add_local_terms(term1, term2) + (term1 isa MPOTerm || term2 isa MPOTerm) && + throw(ArgumentError("Accumulating terms with the same sites is not implemented for MPO or tensor product terms.")) + return VI.add(term1, term2) end Base.:-(O::LocalOperator) = -1 * O @@ -182,9 +377,14 @@ Base.:-(O1::LocalOperator, O2::LocalOperator) = O1 + (-O2) # Since we allow abstract types in T, value and type domain might not match function VI.scalartype(operator::LocalOperator) - return promote_type((scalartype(term[2]) for term in operator.terms)...) + return promote_type((_local_term_scalartype(term) for term in values(operator.terms))...) end +""" +Return the promoted scalar type of a dense tensor or the factors of an MPO. +""" +_local_term_scalartype(term::AbstractTensorMap) = scalartype(term) +_local_term_scalartype(term::MPOTerm) = promote_type((scalartype(tensor) for tensor in term)...) # Equivalence # ----------- diff --git a/test/testsuite/TestSuite.jl b/test/testsuite/TestSuite.jl index 018a68a82..c121be3af 100644 --- a/test/testsuite/TestSuite.jl +++ b/test/testsuite/TestSuite.jl @@ -239,6 +239,19 @@ module ToolboxDensityMatrices end using .ToolboxDensityMatrices +module ToolboxTensorProductTerms + include("toolbox/tensorproduct_terms.jl") + export toolbox_tensorproduct_ising, toolbox_tensorproduct_bookkeeping +end +using .ToolboxTensorProductTerms + +module ToolboxMPOTerms + include("toolbox/mpo_terms.jl") + export toolbox_mpo_terms_dense, toolbox_mpo_bookkeeping, toolbox_mpo_validation + export toolbox_mpo_pepo +end +using .ToolboxMPOTerms + # Utility # ------- module UtilityCorrelator diff --git a/test/testsuite/toolbox/mpo_terms.jl b/test/testsuite/toolbox/mpo_terms.jl new file mode 100644 index 000000000..bb66488cb --- /dev/null +++ b/test/testsuite/toolbox/mpo_terms.jl @@ -0,0 +1,149 @@ +using Test +using Random +using TensorKit +using PEPSKit, Adapt +using PEPSKit: gate_to_mpo, add_term! +import TensorKitTensors.SpinOperators as SO + +const Dbond, χenv = 2, 8 + +const X = SO.σˣ(Float64, Trivial) +const Z = SO.σᶻ(Float64, Trivial) +const lattice11 = fill(ComplexSpace(2), 1, 1) + +const dense_terms = ( + ("NN horizontal", [(1, 1), (1, 2)], -(X ⊗ X) + Z ⊗ Z), + ("NN vertical", [(1, 1), (2, 1)], -(X ⊗ X) + Z ⊗ Z), + ("NNN diagonal", [(1, 1), (2, 2)], X ⊗ Z), + ("3 site line", [(1, 1), (1, 2), (1, 3)], X ⊗ X ⊗ X + Z ⊗ Z ⊗ Z), + # Asymmetric operators expose accidental reordering of sites or factors. + ("NNN anti-diagonal", [(1, 2), (2, 1)], X ⊗ Z), + ("reversed sites", [(1, 2), (1, 1)], X ⊗ Z), +) + +function toolbox_mpo_terms_dense(AT) + Random.seed!(2985721) + return @testset "MPO contractions ($AT)" begin + peps = adapt(AT, InfinitePEPS(ComplexSpace(2), ComplexSpace(Dbond))) + env = CTMRGEnv(randn, ComplexF64, peps, ComplexSpace(χenv)) + @testset "$name" for (name, inds, O) in dense_terms + O = adapt(AT, O) + H_dense = LocalOperator(lattice11, inds => O) + H_mpo = LocalOperator(lattice11, inds => gate_to_mpo(O)) + @test expectation_value(peps, H_mpo, env) ≈ expectation_value(peps, H_dense, env) rtol = 1.0e-9 + end + + # Complex scaling must promote real factors without scaling the coefficient twice. + _, inds, O = first(dense_terms) + O = adapt(AT, O) + H_mpo = LocalOperator(lattice11, inds => gate_to_mpo(O)) + α = 2 + 3im + @test expectation_value(peps, α * H_mpo, env) ≈ + α * expectation_value(peps, LocalOperator(lattice11, inds => O), env) rtol = 1.0e-9 + + # The BP fallback must agree with the existing dense BP contractions. + bp_env = BPEnv(peps) + x, z = adapt.(Ref(AT), (X, Z)) + @testset "BP ($name)" for (name, sites, dense, factors) in ( + ("one-site product", [(1, 1)], x, [x]), + ("two-site product", [(1, 1), (1, 2)], x ⊗ z, [x, z]), + ("MPO", inds, O, gate_to_mpo(O)), + ) + H_dense = LocalOperator(lattice11, sites => dense) + H_factors = LocalOperator(lattice11, sites => factors) + @test expectation_value(peps, H_factors, bp_env) ≈ + expectation_value(peps, H_dense, bp_env) rtol = 1.0e-9 + end + end +end + +function toolbox_mpo_pepo(AT) + return @testset "Factored PEPO observables ($p, purified=$purified) ($AT)" for p in (ℂ^2, Vect[FermionParity](0 => 1, 1 => 1)), purified in (false, true) + Random.seed!(425) + rho = adapt(AT, InfinitePEPO(p, p; unitcell = (1, 1, 1))) + env = CTMRGEnv(purified ? InfinitePEPS(rho) : InfinitePartitionFunction(rho), p) + # The purified form supplies rho as both bra and ket. + args = purified ? (rho, env) : (env,) + a, b = ntuple(_ -> adapt(AT, randn(ComplexF64, p ← p)), 2) + gate = adapt(AT, randn(ComplexF64, p ⊗ p ← p ⊗ p)) + @testset "$name" for (name, sites, dense, factors) in ( + ("one-site product", [(1, 1)], a, [a]), + ("two-site product", [(1, 1), (1, 2)], a ⊗ b, [a, b]), + ("MPO", [(1, 1), (1, 2)], gate, gate_to_mpo(gate; trunc = notrunc())), + ) + H_dense = LocalOperator(physicalspace(rho), sites => dense) + H_factors = LocalOperator(physicalspace(rho), sites => factors) + @test expectation_value(rho, H_factors, args...) ≈ + expectation_value(rho, H_dense, args...) rtol = 1.0e-9 + end + end +end + +function toolbox_mpo_bookkeeping(AT) + Random.seed!(2985721) + return @testset "MPO bookkeeping ($AT)" begin + d = ℂ^2 + lattice = fill(d, 2, 2) + dense = adapt(AT, randn(Float64, d ⊗ d ← d ⊗ d)) + mpo = gate_to_mpo(dense; trunc = notrunc()) + sites = CartesianIndex.([(3, 4), (3, 3)]) + original_sites, original_mpo = copy(sites), copy(mpo) + operator = LocalOperator(lattice, sites => mpo) + @test sites == original_sites + reverse!(sites) + reverse!(mpo) + @test only(operator.terms) == (CartesianIndex.([(1, 2), (1, 1)]) => original_mpo) + + snapshot = deepcopy(operator) + scaled = (2 + 3im) * operator + @test operator == snapshot + onsite = adapt(AT, randn(ComplexF64, d ← d)) + @test scalartype(operator + LocalOperator(lattice, ((2, 1),) => onsite)) == ComplexF64 + @test scalartype(scaled) == ComplexF64 + @test_throws ArgumentError real(scaled) + @test_throws ArgumentError imag(scaled) + + # Both insertion and addition must reject accumulation involving an MPO. + inds = CartesianIndex.([(1, 1), (1, 2)]) + dense_operator = LocalOperator(lattice, inds => dense) + mpo_operator = LocalOperator(lattice, inds => original_mpo) + for (left, right) in ((mpo_operator, mpo_operator), (mpo_operator, dense_operator), (dense_operator, mpo_operator)) + @test_throws ArgumentError add_term!(deepcopy(left), inds, only(values(right.terms))) + @test_throws ArgumentError left + right + end + @test isempty(LocalOperator(lattice, inds => [zero(first(original_mpo)), last(original_mpo)]).terms) + + # Unequal physical spaces expose factor reordering during coordinate transformations. + lattice = [ℂ^2 ℂ^3 ℂ^4; ℂ^5 ℂ^6 ℂ^7] + dense = adapt(AT, randn(Float64, ℂ^6 ⊗ ℂ^2 ← ℂ^6 ⊗ ℂ^2)) + operator = LocalOperator(lattice, ((2, 2), (1, 1)) => gate_to_mpo(dense; trunc = notrunc())) + @test rotr90(rotl90(operator)) == operator + sites, term = only(operator.terms) + @test repeat(operator, 2, 1).terms == Dict(sites => term, (sites .+ CartesianIndex(2, 0)) => term) + end +end + +function toolbox_mpo_validation(AT) + Random.seed!(2985721) + return @testset "MPO validation ($AT)" begin + d, b1, b2 = ℂ^2, ℂ^3, ℂ^4 + lattice = fill(d, 2, 2) + first_tensor = adapt(AT, randn(Float64, d ← d ⊗ b1)) + last_tensor = adapt(AT, randn(Float64, b1 ⊗ d ← d)) + mpo = [first_tensor, last_tensor] + sites = ((1, 1), (1, 2)) + + @test_throws ArgumentError LocalOperator(lattice, CartesianIndex{2}[] => AbstractTensorMap[]) + @test_throws ArgumentError LocalOperator(lattice, ((1, 1),) => mpo) + @test_throws ArgumentError LocalOperator(lattice, ((1, 1), (1, 1)) => mpo) + @test_throws ArgumentError LocalOperator(lattice, sites => reverse(mpo)) + + # Bond matching and the two physical legs are independent constraints. + wrong_bond = adapt(AT, randn(Float64, b2 ⊗ d ← d)) + wrong_input = adapt(AT, randn(Float64, b1 ⊗ d ← ℂ^3)) + wrong_output = adapt(AT, randn(Float64, b1 ⊗ ℂ^3 ← d)) + @test_throws SpaceMismatch LocalOperator(lattice, sites => [first_tensor, wrong_bond]) + @test_throws SpaceMismatch LocalOperator(lattice, sites => [first_tensor, wrong_input]) + @test_throws SpaceMismatch LocalOperator(lattice, sites => [first_tensor, wrong_output]) + end +end diff --git a/test/testsuite/toolbox/tensorproduct_terms.jl b/test/testsuite/toolbox/tensorproduct_terms.jl new file mode 100644 index 000000000..58d7e3616 --- /dev/null +++ b/test/testsuite/toolbox/tensorproduct_terms.jl @@ -0,0 +1,71 @@ +using Test +using Random +using TensorKit +using PEPSKit, Adapt +using PEPSKit: add_term! +import TensorKitTensors.SpinOperators as SO + +const Dbond, χenv, g = 2, 8, 3.1 +const unitcells = ((1, 1), (2, 2)) + +""" +Build the transverse-field Ising Hamiltonian from real single-site factors. +""" +function transverse_field_ising_tensorproduct(AT, lattice::InfiniteSquare; g = 1.0) + Z, X = adapt.(Ref(AT), (SO.σᶻ(Float64, Trivial), SO.σˣ(Float64, Trivial))) + spaces = fill(domain(X)[1], (lattice.Nrows, lattice.Ncols)) + return LocalOperator( + spaces, + (neighbor => [-Z, Z] for neighbor in nearest_neighbours(lattice))..., + ([idx] => [-g * X] for idx in vertices(lattice))..., + ) +end + +function toolbox_tensorproduct_ising(AT) + return @testset "Tensor product Ising Hamiltonian ($uc) ($AT)" for uc in unitcells + Random.seed!(2985721) + H = transverse_field_ising(InfiniteSquare(uc...); g) + H_prod = transverse_field_ising_tensorproduct(AT, InfiniteSquare(uc...); g) + peps = adapt(AT, InfinitePEPS(ComplexSpace(2), ComplexSpace(Dbond); unitcell = uc)) + env, = leading_boundary( + CTMRGEnv(peps, ComplexSpace(χenv)), peps; tol = 1.0e-8, verbosity = 0 + ) + E = expectation_value(peps, H, env) + @test expectation_value(peps, H_prod, env) ≈ E rtol = 1.0e-9 + # Real factors must retain tensor-product dispatch after complex scaling. + α = 2 + 3im + @test expectation_value(peps, α * H_prod, env) ≈ α * E rtol = 1.0e-9 + end +end + +function toolbox_tensorproduct_bookkeeping(AT) + return @testset "Tensor product bookkeeping ($AT)" begin + lattice = fill(ComplexSpace(2), 1, 1) + X, Z = adapt.(Ref(AT), (SO.σˣ(Float64, Trivial), SO.σᶻ(Float64, Trivial))) + sites = CartesianIndex.([(3, 4), (3, 5)]) + factors = [X, Z] + product = LocalOperator(lattice, sites => factors) + @test sites == CartesianIndex.([(3, 4), (3, 5)]) + reverse!(sites) + reverse!(factors) + @test only(product.terms) == (CartesianIndex.([(1, 1), (1, 2)]) => [X, Z]) + @test product == LocalOperator(lattice, [(1, 2), (1, 1)] => [Z, X]) + + snapshot = deepcopy(product) + scaled = (2 + 3im) * product + @test product == snapshot + @test scalartype(scaled) == ComplexF64 + @test_throws ArgumentError add_term!(product, [(1, 1), (1, 2)], [Z, X]) + @test_throws ArgumentError product + product + + # Factorwise real/imaginary parts are incorrect even for tensor products. + imaginary = LocalOperator(lattice, [(1, 1), (1, 2)] => [im * Z, im * Z]) + @test_throws ArgumentError real(imaginary) + @test_throws ArgumentError imag(imaginary) + @test_throws ArgumentError LocalOperator(lattice, CartesianIndex{2}[] => typeof(X)[]) + @test_throws ArgumentError LocalOperator(lattice, [(1, 1), (1, 2)] => [X]) + @test_throws SpaceMismatch LocalOperator( + lattice, [(1, 1)] => [SO.S_x(Float64, Trivial; spin = 1)] + ) + end +end diff --git a/test/toolbox/mpo_terms.jl b/test/toolbox/mpo_terms.jl new file mode 100644 index 000000000..1a8117c62 --- /dev/null +++ b/test/toolbox/mpo_terms.jl @@ -0,0 +1,14 @@ +using Test +using PEPSKit + +@isdefined(TestSuite) || include("../testsuite/TestSuite.jl") +using .TestSuite + +is_buildkite = get(ENV, "BUILDKITE", "false") == "true" + +if !is_buildkite + TestSuite.toolbox_mpo_terms_dense(Vector) + TestSuite.toolbox_mpo_pepo(Vector) + TestSuite.toolbox_mpo_bookkeeping(Vector) + TestSuite.toolbox_mpo_validation(Vector) +end diff --git a/test/toolbox/tensorproduct_terms.jl b/test/toolbox/tensorproduct_terms.jl new file mode 100644 index 000000000..a9b5a8d9f --- /dev/null +++ b/test/toolbox/tensorproduct_terms.jl @@ -0,0 +1,12 @@ +using Test +using PEPSKit + +@isdefined(TestSuite) || include("../testsuite/TestSuite.jl") +using .TestSuite + +is_buildkite = get(ENV, "BUILDKITE", "false") == "true" + +if !is_buildkite + TestSuite.toolbox_tensorproduct_ising(Vector) + TestSuite.toolbox_tensorproduct_bookkeeping(Vector) +end