Hi, thanks for explaining some more. To me, though, those sound like slightly different or incompatible goals. We can definitely help you to make a Hamiltonian MPO that represents the sum over 4 indices. The other goal about an MPO that has (perhaps?) an extra non-standard index for the kind of Pauli matrix might be a different goal and I’m not sure that the result would technically be an MPO? But perhaps I don’t understand your goal in that case.
To make an MPO which makes Heisenberg interactions connecting sites i and j with interaction coefficients J_ij you can use the OpSum and AutoMPO system to do something like this:
os = OpSum()
for i=1:N, j=1:N
os += J[i,j], “Sz”, i, “Sz”, j
os += J[i,j]/2, “S+”, i , “S-“, j
os += J[i,j]/2, “S-“, i, “S+”, j
and then turn the OpSum into an MPO as in the ITensor documentation and code examples.
Here in the code above J[i,j] is just a regular Julia array with two dimensions that you would make earlier in the code.
You can generalize the above pattern to 4 indices by just summing over the 4 indices and including four operators in each term if that is what you want. Please let me know if you have a question about that.
If by alpha and beta, though, you just mean the indices which run over x,y,z, for the spin operators to make a Heisenberg S.S interaction, then the above code already includes these in the form of the three terms I put which are (schematically) these terms: Sz*SZ, S+*S-, S-*S+