Manually constructing MPO from set of ITensors

I have a hamiltonian

H = -\sum_i H_i

where

H_i\ket{ijk}

is a 3-site operator. For various reasons, I am finding it difficult to express this hamiltonian as a sum of one-site operators, and rather prefer to calculate the matrix elements directly. I wish to construct an MPO for this hamiltonian.
What I have done so far is manually construct the 3-site ITensor using the appropriate matrix elements., and performed successive SVDs to break it down to 3 ITensors - say T1, T2, T3.

My question is, how do I form an MPO of N-sites using these three ITensors?

something like

s = siteinds("S=1/2",10)

s3 = [3,4,5]
Tops = [T1,T2,T3]
H = MPO([op(si,"I") for si in s])
for (si,Ti) in zip(s3,Tops)
  H[si] = Ti
end

could work

Ryan’s suggestion above would work for turning a single (3-operator) term in your sum into an MPO. Converting an entire sum of such (non-product) terms into an MPO is more technical, and beyond the ability of our OpSum system currently.

But here is something you can try which might work decently well, at least for smaller systems: you could make an array of MPOs using Ryan’s approach then pass this array into DMRG (assuming that’s what you’re doing) because our DMRG code does accept an array of MPO’s which it treats implicitly as a sum. This will likely have more memory overhead compared to using a single MPO but might not be too much worse. Please give it a try, or let us know if it doesn’t work or you need more details.