If a spin Hamiltonian contains a term which is the product of ALL local site operators, say \prod_i \sigma^z_i, how to use OpSum to add this term into the Hamiltonian MPO? Since it becomes impractical to write down this term in the usual way as other local terms…

Is it possible to construct such term into an MPO manually in ITensor, and then add it into the Hamiltonian MPO which is constructed separately. Can someone elaborate on this in details.

s = siteinds("S=1/2",N)
pSz = MPO(s,"Sz")
os = OpSum()
for j=1:N-1
# add terms to OpSum
end
H = MPO(os,s)
#...
energy,psi = dmrg([H,pSz],psi0; nsweeps, maxdim, cutoff)

where the key points are:

you can make a product of all “Sz” operators by calling MPO(s,"Sz")

you can pass two or more MPOs to DMRG in a vector like [H,pSz] and DMRG will treat them as if they are summed for the purpose of the ground state calculation. It will not actually sum them, but loop over them internally in an efficient way when doing each core step of the DMRG algorithm.