I am interested in developing on ITensors.jl the lattice Schwinger model with alternating (anti)fermion sites and links (representing electric flux). Each spatial site has 4 ITensors sites--first a fermion site, then a link, then an antifermion site, and then another link. It incorporates a local constraint given by Gauss's Law, and I've gotten it to mostly work (I think) by using a few tricks.

However, I am still running into one problem. Basically, the smallest "unit" that you can change without changing anything else is a unit consisting of two (anti)fermion sites and the link between them. For example, I can change sites 1-3 without changing anything else; similar for 3-5, 5-7, etc. This is why I run 3-site DMRG (this isn't natively in ITensors; I developed this code and have put a pull request for it here: https://github.com/ITensor/ITensors.jl/pull/434). However, this spells trouble when you reach the end. One of the "units" that you should be able to change is that consisting of sites 4N-1, 4N, and 1, where N is the number of spatial sites. However, DMRG doesn't do such a thing, because it only works on a bundle of adjacent sites. This means that, due to the Gauss's Law constraint enforced by sites 4N, 1, and 2, site 4N is stuck unless one also changes either site 1 or site 2, which only happens if one uses k-site DMRG for k=4N-1 or 4N (i.e. adjusting the whole state, minus at most one site, in each step). For example, for 2 spatial sites, I need 7-site DMRG to be able to change site 8 at all, and for 3 spatial sites, I need 11-site DMRG. This is simply infeasible for any lattice larger than 3 spatial sites.

I am wondering whether anyone has any ideas about how to handle this or has dealt with something similar before. One idea I had was "rotating" the state between or during sweeps in the DMRG process, so that every site has at least some time not being the one at the end, which is stuck in place. As an example, one could "rotate" the state so that the last two sites are now the first two, and then sites 4N, 1, and 2 (which are now sites 2-4), are together and can be optimized. And then one can rotate the state back and complete the rest of the sweep. But I'm not sure how one would do such a thing, or whether it would really produce the desired results. Does anyone know of a feasible fix for this, whether by rotating the state or by any other means?

It is true that this error should in principle become relatively small for a sufficiently large number of spatial sites, but in my research, I will mostly be working with fairly small lattices, for which such limits won't apply. Also, I need my code to be accurate in finding the ground state so that I can accurately evaluate excited states, as well as assess my code against existing results on the Schwinger model produced using other computational methods (e.g. exact diagonalization). Any help on this would be greatly appreciated.