Here is an idea that I had (written using the Julia code). I still highly doubt it will work, for some reasons I will explain, but I have some hope that something like it may work. In particular, I think understanding why it doesn't work would still be helpful. Anyway, here is the code:

```
ferm_full = QN("Sz",-1) => 1
ferm_vac = QN("Sz",0) => 1
antiferm_full = QN("Sz",+1) => 1
antiferm_vac = QN("Sz",0) => 1
even_link_p1 = QN("Sz",+2) => 1 # p1 = "positive 1"
even_link_0 = QN("Sz",0) => 1
even_link_n1 = QN("Sz",-2) => 1 # n1 = "negative 1"
odd_link_p1 = QN("Sz",-2) => 1
odd_link_0 = QN("Sz",0) => 1
odd_link_n1 = QN("Sz",+2) => 1
function Kogut_Susskind_sites(num_spatial_sites)
sites = Index[]
for i = 1:num_spatial_sites
push!(sites, Index(ferm_full,ferm_vac; tags="Site,S=1/2,n=$(4*i-3)"))
push!(sites, Index(odd_link_p1,odd_link_0,odd_link_n1; tags="Site,S=1,n=$(4*i-2)"))
push!(sites, Index(antiferm_full,antiferm_vac; tags="Site,S=1/2,n=$(4*i-1)"))
push!(sites, Index(even_link_p1,even_link_0,even_link_n1; tags="Site,S=1,n=$(4*i)"))
end
return sites
end
function Schwinger_Hamiltonian(sites, num_spatial_sites, x, μ)
N = num_spatial_sites * 2
ampo_H = AutoMPO()
for j=0:N-1
add!(ampo_H, 1, "Sz2", 2*j+2)
add!(ampo_H, (μ/2)*(-1)^j * 2, "Sz", 2*j+1)
if j != N-1
add!(ampo_H, x * 2*2/sqrt(2), "S+", 2*j+1, "S-", 2*j+2, "S-", 2*j+3)
add!(ampo_H, x * 2*2/sqrt(2), "S+", 2*j+3, "S+", 2*j+2, "S-", 2*j+1)
end
end
add!(ampo_H, x * 2*2/sqrt(2), "S+", 2*N-1, "S-", 2*N, "S-", 1)
add!(ampo_H, x * 2*2/sqrt(2), "S+", 1, "S+", 2*N, "S-", 2*N-1)
H = MPO(ampo_H, sites)
return H
end
function Schwinger_DMRG(num_spatial_sites)
sites = Kogut_Susskind_sites(num_spatial_sites)
H = Schwinger_Hamiltonian(sites, num_spatial_sites, x, μ)
sweeps = Sweeps(10)
maxdim!(sweeps, 10,20,100,100,200)
cutoff!(sweeps, 1E-10)
# construct initial state as completely empty
init_state = ["Up" for n = 1:4*num_spatial_sites] # just to initialize
for n = 1:4*num_spatial_sites
if n % 4 == 1
init_state[n] = "Dn" # empty fermion sites
elseif n % 4 == 3
init_state[n] = "Up" # empty antifermion states
else # n % 4 == 2 || n % 4 == 0
init_state[n] = "Up" # flux = 0 at all links
end
end
psi0 = randomMPS(sites, init_state, 100)
energy, psi = dmrg(H, psi0, sweeps)
println("Ground state energy = $energy")
return energy
end
```

In case you're wondering why I haven't already tested it, basically it's because I ran into some other more fundamental issue with the randomMPS function in Julia. This issue is mentioned on the Support Q&A (http://itensor.org/support/2186/questions-random-initial-conservation-randommps-linkdim) and on the GitHub page (https://github.com/ITensor/ITensors.jl/issues/424), so I think it will be fixed soon.

The idea is that, by assigning these particular quantum numbers, we have now made it so that each term in the Hamiltonian indeed has zero QN flux. So I think this means that, when the DMRG process varies the state, it will only end up doing things that respect the gauge symmetry? But maybe I'm misunderstanding.

The first problem I see with this is that I have still only defined one QN, so I'm pretty sure this QN conservation is not actually capturing the symmetry at each site. I believe there would be ways for DMRG to change the state in a way that respects this QN but doesn't respect the local constraints, in some way that I don't really understand. If this is the case, then I wonder whether setting up this type of quantum number for every site could be the way to go.

The second problem I see with this is that it is unclear to me how DMRG will work properly. Thanks to the gauge symmetries, the smallest "legal" operation is creating a fermion-antifermion pair on adjacent sites and changing the flux between the two particles by 1 (I'm assuming the fermion is negative, so +1 if the antifermion comes first and -1 if the fermion comes first), or the opposite (destroying a fermion-antifermion pair and making the corresponding flux change for the site between them). This means that one needs to change three adjacent sites (two particle sites and the link between them) to make a legal change to the state, and DMRG only isolates two sites at a time when performing the local optimization. This seems to me like something that would lead to trouble.

Anyway, I'm sorry to disappoint if you were expecting a definitive answer, but I thought this was extensive enough to put as an answer to at least get some thoughts going. I would love to hear people's thoughts on this.