[ITensors] [ENHANCEMENT] Add support for defining MPOs from operators…

… represented as matrices (#904)

src/mps/mpo.jl

@@ -66,7 +66,7 @@ MPO(N::Int) = MPO(Vector{ITensor}(undef, N))
 Make an MPO with pairs of sites `s[i]` and `s[i]'`
 and operators `ops` on each site.
 """
-function MPO(::Type{ElT}, sites::Vector{<:Index}, ops::Vector{String}) where {ElT<:Number}
+function MPO(::Type{ElT}, sites::Vector{<:Index}, ops::Vector) where {ElT<:Number}
   N = length(sites)
   ampo = OpSum() + [ops[n] => n for n in 1:N]
   M = MPO(ampo, sites)
@@ -101,6 +101,15 @@ end
 
 MPO(sites::Vector{<:Index}, op::String) = MPO(Float64, sites, op)
 
+function MPO(::Type{ElT}, sites::Vector{<:Index}, op::Matrix{<:Number}) where {ElT<:Number}
+  # return MPO(ElT, sites, fill(op, length(sites)))
+  return error(
+    "Not defined on purpose because of potential ambiguity with `MPO(A::Array, sites::Vector)`. Pass the on-site matrices as functions like `MPO(sites, n -> [1 0; 0 1])` instead.",
+  )
+end
+
+MPO(sites::Vector{<:Index}, op::Matrix{ElT}) where {ElT<:Number} = MPO(ElT, sites, op)
+
 function randomMPO(sites::Vector{<:Index}, m::Int=1)
   M = MPO(sites, "Id")
   for i in eachindex(sites)

src/physics/autompo.jl

@@ -137,7 +137,7 @@ function MPOTerm(op1::Union{String,AbstractArray}, ops...)
   return MPOTerm(one(Float64), op1, ops...)
 end
 
-function MPOTerm(ops::Vector{Pair{String,Int}})
+function MPOTerm(ops::Vector{<:Pair})
   return MPOTerm(Iterators.flatten(ops)...)
 end
 
@@ -301,7 +301,7 @@ function (ampo::OpSum + term::MPOTerm)
 end
 
 (ampo::OpSum + term::Tuple) = ampo + MPOTerm(term...)
-(ampo::OpSum + term::Vector{Pair{String,Int64}}) = ampo + MPOTerm(term)
+(ampo::OpSum + term::Vector{<:Pair}) = ampo + MPOTerm(term)
 
 function (ampo::OpSum - term::Tuple)
   ampo_plus_term = copy(ampo)

test/mpo.jl

@@ -399,7 +399,7 @@ end
     l = [Index(3, "left_$n") for n in 1:2]
     r = [Index(3, "right_$n") for n in 1:2]
 
-    sis = IndexSet.(prime.(s), s)
+    sis = [[sₙ', sₙ] for sₙ in s]
 
     A = randomITensor(s..., prime.(s)...)
     ψ = MPO(A, sis; orthocenter=4)
@@ -429,6 +429,10 @@ end
     @test ITensors.orthocenter(ψ) == 3
     @test maxlinkdim(ψ) == 1
 
+    # Use matrix
+    @test_throws ErrorException MPO(s, [1/2 0; 0 1/2])
+    @test MPO(s, _ -> [1/2 0; 0 1/2]) ≈ MPO(s, "Id") ./ 2
+
     ψ0 = MPO(s, "Id")
     A = prod(ψ0)
     ψ = MPO(A, s; cutoff=1e-15, orthocenter=3)