Thanks for the extra details.
So actually you’re in luck and this is the simplest kind of MPO to construct. You could still do it using OpSum but it would be overkill in this case. All you have to do in this case is the following. Make an MPO for your N qubits and a collection of site indices you are going to use
N = 10
s = siteinds("Qubit",N)
M = MPO(s)
then just set the MPO tensors one by one to the values you want, following the ITensor convention that the j’th MPO tensor has the indices s[j], s[j]'
.
So like this:
for j=1:N
T = ITensor(s[j],s[j]')
T[1,1] = ...
T[2,1] = ...
T[1,2] = ...
T[2,2] = ...
M[j] = T
end
You could either work out the matrix elements of T ahead of time and then set them one by one like I am suggesting in the code snippet above. Or you could use some built in ITensor facilities for making each MPO tensor. What I mean is, say you want to make the operator
3|x+\rangle\langle x\!+\!| - I
You could do this with the code
T = 3*state("X+",s[j]')*state("X+",s[j]) - op("I",s[j])
Here I haven’t included the U’s explicity, but based on what you said they are Pauli matrices so you could work out what U^\dagger_j |b_j\rangle is in advance without numerically acting U^\dagger_j onto a state.
For a list of predefined states and operators for “Qubit” sites please see this page:
https://itensor.github.io/ITensors.jl/dev/IncludedSiteTypes.html#“Qubit”-SiteType
You can also define your own custom overloads of the state
and op
functions if you wish:
https://itensor.github.io/ITensors.jl/dev/examples/Physics.html