Thanks for the good question. The short answer is yes it can be done with ITensor, though with the support for different algorithms currently differing a bit from the C++ versus Julia version.
The current status is:
- there is a good quality TDVP package for the C++ version
- the function "apply" which performs TEBD in Julia can handle three and even more site operators very easily
Ideally we'd have both algorithms in both language - that's happening this year for TDVP in Julia but isn't done yet.
But since you said you're working in C++, TDVP would be a really great choice for this kind of interaction, which is quite straightforward to make into an MPO using the AutoMPO system. (Truly long range is trickier in both versions of ITensor, though works well on modestly sized systems or with more custom techniques to make MPOs "by hand" as it were.)
The main tricky thing to know about TDVP is that it can need a little help in the early part of a time evolution, where the state may change very rapidly. Just lowering the time step is one thing that can help, but can be insufficient. Other techniques are to use a more reliable, if slower, technique like TEBD for the few handful of time steps or to use "basis expansion" methods like those of Yang and White to help with convergence. (M. Yang is the person who actually wrote the C++ TDVP ITensor code!)