I'd like to follow up on your explanation in github:
I tried evolving with ExpMPOEvolution in the first step, and then applying TDVP works, as you said. But does this mean that, at some later time TDVP will saturate the new bond dimension, and I'll have to do another ExpMPOEvolution step, because the MPO I use is not nearest-neighbor w.r.t. the MPO index?This issue is not a bug of TeNPy, it's a limitation of TDVP trying to keep the MPS bond dimension fixed and starting from a product state.
While SingleSite TDVP strictly keeps the MPS bond dimension fixed, two-site TDVP can increase it only if you get a non-trivial evolution (increasing the bond dimension on the given bond you evolve). This is not the case here, because you don't have nearest-neighbor interactions:
> print(H.all_coupling_terms().to_TermList())
1.00000 * Sx_2 Sx_4 +
1.00000 * Sx_3 Sx_5
e.g. when trying to evolve sites 2-3, the "Sx_4" and "Sx_5" are projected out, so the total effective H on sites 2-3 is zero.
The solution is probably to use another time evolution algorithm to initially increase the MPS bond dimension; after that TDVP should work as expected as well. Indeed, after one time step with the ExpMPOEvolution, I find that TDVP also changes total Sz.