Hello,

I am running some calculations on a 4-site square chain with a tight-binding Hamiltonian, and have found an issue when applying the DMRG function to my system. Depending on the state I initialize the system, the DMRG can either run normally and return the expected value for the energy, or almost instantly return a (wrong) value of 0 for the ground state energy. Here is a sample code to demonstrate the problem:

```
#include "itensor/all.h"
using namespace itensor;
int main()
{
int N = 4;
double t = 1.0;
int nsweeps = 5;
int maxm = 200;
double cutoff = 1E-15;
int niter = 2;
double noise = 0;
auto sites = Hubbard(N);
auto ampo = AutoMPO(sites);
ampo += -t,"Cdagup",1,"Cup",2;
ampo += -t,"Cdagup",2,"Cup",1;
ampo += -t,"Cdagdn",1,"Cdn",2;
ampo += -t,"Cdagdn",2,"Cdn",1;
ampo += -t,"Cdagup",1,"Cup",3;
ampo += -t,"Cdagup",3,"Cup",1;
ampo += -t,"Cdagdn",1,"Cdn",3;
ampo += -t,"Cdagdn",3,"Cdn",1;
ampo += -t,"Cdagup",2,"Cup",4;
ampo += -t,"Cdagup",4,"Cup",2;
ampo += -t,"Cdagdn",2,"Cdn",4;
ampo += -t,"Cdagdn",4,"Cdn",2;
ampo += -t,"Cdagup",3,"Cup",4;
ampo += -t,"Cdagup",4,"Cup",3;
ampo += -t,"Cdagdn",3,"Cdn",4;
ampo += -t,"Cdagdn",4,"Cdn",3;
//Use IQMPO or MPO to generate the hamiltonian
auto H = IQMPO(ampo);
Print(ampo);
//Define initial wavefunction MPS for IQMPS
auto state1 = InitState(sites);
state1.set(1,"Up");
state1.set(2,"Up");
state1.set(3,"0");
state1.set(4,"0");
//Create IQMPS
auto psi = IQMPS(state1);
//Overlap test
auto state2 = InitState(sites);
state2.set(1,"0");
state2.set(2,"Up");
state2.set(3,"Up");
state2.set(4,"0");
auto psi2 = IQMPS(state2);
println(overlap(psi2,H,psi));
//Define DMRG sweeps
auto sweeps = Sweeps(nsweeps);
sweeps.maxm() = maxm/20,maxm/10,maxm/2,maxm/2,maxm;
sweeps.cutoff() = cutoff;
sweeps.niter() = niter;
sweeps.noise() = noise;
auto energy = dmrg(psi,H,sweeps,{"Quiet",true});
return 0;
} //End of main
```

When I run this code exactly as it is, the DMRG returns a value of 0 for the ground state energy. The same happens if the two electrons are initialized in sites 3 and 4. However, if I change the initial conditions to initialize the electrons on other combinations of sites (for instance sites 1 and 3), I get the correct value for the energy (-2).

I've made some tests and found out this problem does not happen if I change the ordering of the sites to a cyclic chain, that is, with interactions 1-2, 2-3, 3-4, 4-1. However, since my goal is to work with larger chains, I figured I should number the sites the way I am doing in the sample code, as to avoid creating super-long-range interactions when working with more sites. Therefore, I think having a fix for this problem while still using my ordering would be ideal.

I believe the cause of this issue is the way the ITensor library handles the fermionic operators "Cdagup" and "Cup". Depending on the order they are applied when initializing the system and how many anticommutations are made when the Hamiltonian is applied, the final sign of the tight-binding interaction might become positive when it should be negative, thus "fooling" the DMRG algorithm into believing the system is already at the lowest energy level. To exemplify, if one calculates the overlap <0 Up Up 0| H |Up Up 0 0> in the sample code, instead of the expected result of -1 (given the sign of the hopping term), the code returns +1.

This overlap test made me think a simple change of signs on the Hamiltonian could solve everything. However, I discovered that changing the sign of the hopping term on the specific part of the Hamiltonian that annihilates one electron on the first site and creates it on the third site doesn't seem to completely fix the problem; although the overlap described above now returns a value of -1, the DMRG surprisingly still returns a value of 0 for the energy.

So, I wanted to ask if there is anyway I can have consistent (a.k.a. not dependent on the initial conditions) DMRG results for my system while keeping the ordering of the sites intact. I would appreciate any input into this. Thanks!