I am trying the define a three local Hamiltonian H= \sum_i h_i with h_j= X_j \; CZ_{j-1,j+1} but keep getting an error due to the two local CZ gate at the end. Is there an straightforward way to get around the issue?

using ITensors
N=16
acc=1E-13
bond_dim=50
sites = siteinds("Qubit",N)
psi0 = randomMPS(sites,bond_dim)
ampo = AutoMPO()
for j=2:N-1
# Want to define h_j= X_j CZ_{j-1,j+1}
ampo += -2,"CZ",j-1,j+1,"σx",j;
end
H = MPO(ampo,sites)
sweeps = Sweeps(15)
maxdim!(sweeps,10,20,50,100,200,400,500)
cutoff!(sweeps,acc)
energy,T = dmrg(H,psi0,sweeps,outputlevel=0)
end

Output: ArgumentError: Tuple contains 2 elements, must contain exactly 1 element

It’s an interesting thing to do, and would be a good feature for us to have, but unfortunately right now the AutoMPO system (since renamed the OpSum system) doesn’t support two-site operators like “CZ” currently. Internally what it would have to do is to compute the SVD of this operator to write it as a sum of products of single-site operators, and then put these into the sum of terms (like a rewriting rule). So it’s something we could add but would take a bit of development.

As a workaround, you could work out a factorization of CZ like this for yourself, either numerically or probably by hand if you worked on it a bit. Basically you could define it as:

P_0 \otimes I + P_1 \otimes Z

where P_0 projects the control site onto the 0 state and P_1 projects the control site onto the 1 state. That sum of terms would then be something our system can read in.

Currently for the “Qubit” site type, we have the operators P_0 and P_1 defined as “ProjUp” and “ProjDn” respectively, but it would be good to add the aliases “Proj0” and “Proj1” so I’ll do that and it will be in the next version. For now you could use “ProjUp” and “ProjDn” or define your own custom operators.