I would like to convert an OpSum to a list of gates. I’m doing so by converting the OpSum terms to MPO then contracting with prod. The problem is that this yields pauli strings with all the identities. How can I drop the identities to only keep non trivial indices ?
For example, in the following code, A is the gate that I want, but I want to construct it like B, however B also includes identity on sites 3 and 4. How can I drop these ?
N = 4
sites = siteinds("S=1/2", N)
A = op("Sx", sites[1]) * op("Sx", sites[2])
println(A)
B = OpSum()
B += 1, "Sx", 1, "Sx", 2
B = prod(MPO(B[1], sites))
println(B)
Hi Nichola,
The direct answer to your question is that we don’t have a totally polished or “production ready” system right now to convert a whole OpSum to a product of gates. Also for that case, we’d anyway not call it an OpSum but more like an OpProduct or similar, and have a multiplication syntax for adding gates to it.
We plan to expand support greatly for that sort of thing in the future, but we are more in the planning stages of that right now.
Meanwhile, I’d have two suggestions for you (or questions):
what was the limitation you found with your “A” pattern? Even though the syntax there isn’t as nice as with OpSum, you can make fairly arbitrary gates that way and push them into a Vector of ITensors. Is there a reason you really needed to use an OpSum or the OpSum syntax?
assuming you do want to convert terms of an OpSum to gates, you are on the right track with your prod(MPO... workaround but a better version of it Matt Fishman suggested to me is to do prod(MPO(B[1], [sites[1], sites[2]])) and the generalization thereof to other sites or more sites. I think you know already, but the MPO constructor there is passed an array of sites on which it will be “supported” so you don’t have to pass all the sites but instead can pass the ones you want the MPO to be supported on, and they don’t even have to be contiguous sites especially if afterward you’re just going to product all the tensors together to make a single ITensor.
I’m doing TEBD for a bunch of different Hamiltonians, and I also need the same Hamiltonians as MPO’s for other applications. So now, for each H, I have a function that generate list of gates for TEBD and another function that construct the MPO. It would be much clearer to have a single function that construct each H as an OpSum, then convert them to the appropriate representation for each application. Also, “A” is to my opinion harder to read than the way the OpSum are built.
Yes passing the support to the MPO constructor seems to be the solution, however I run into this issue :
N = 4
sites = siteinds("S=1/2", N)
O = OpSum()
O += 1, "Sx", 2, "Sx", 3
B = prod(MPO(O[1], [sites[2], sites[3]]))
ERROR: LoadError: The OpSum contains a term 1.0 Sx(2,) Sx(3,) that extends beyond the number of sites 2.