where \phi is just a phase.
Without the exponential factor, it is clear how to construct MPO via autoMPO. Is there any clever way to include the exponential factor via autoMPO routines or it is better to construct it manually?

I thought of a simple method you can use. You should be able to show that

e^{i \phi S^z} = \cos( \phi/2) I + i \sin( \phi/2) S^z

but please check me on this. To prove this you can just use the fact that S^z = \frac{1}{2} \sigma^z and that (\sigma^z)^2 = 1, then just use the Taylor series for e^x, \cos(x), and \sin(x).

I had originally written a longer answer below. I’ll leave it here if it’s helpful, but you should be able to just check my formula above then plug it into your Hamiltonian and then get something that AutoMPO can accept as input.

Other answer:

It’s a good question. The best way to do this with the AutoMPO system is to define a new local operator. And the easiest way to do that with the C++ version of ITensor is to make your own custom SiteType that is based on the SpinHalf site type but which you customize to include other local operators beyond the standard set.

If you copy the file itensor/mps/sites/spinhalf.h to another file custom_spinhalf.h and edit the class name to be something like CustomSpinHalf then you could add a new operator inside the op method with a name like "ZPhase" say. You’d probably have to pass the value of \phi through this site set class as well, since currently there’s not a good way to pass parameters that go inside operators through the AutoMPO interface.

(Note that all of the above is much, much easier to do in the Julia version of ITensor if that’s a possbility for you to use for your project.)

Thank you very much for your detailed answer. Regarding you additional point, the true Hamiltonian is for higher spin, I just wrote the question to get the general concept. Although, I got nice expressions in this case also. Let me try!