Hi everyone,

I am trying to observe a phase transition in the non interacting SSH model, described by the following Hamiltonian( for reference see https://arxiv.org/pdf/1810.13286.pdf):

The SSH model is formulated on a one-dimensional lattice, divided into two sub-lattices: A = {1, 3, . . . , N − 1}, involving odd lattice sites, and B = {2, 4, . . . , N}, with an even number of sites N and staggered hopping of particles between A and B ( here I considered one hopping strength equal to 1 and the other equal to J). Depending on the value of J <1 or >1 , I should observe a phase transition, from a trivial state with to a topological state, characterized by the presence of edge state (see the following image for context).

I calculated the density operator for each site ( since I mapped the model to a spin half chain I used the following density operator)

and printed the result and I always get 0,5 for every site, no matter the value of J and I don’t see why the code it is not working as it should. I attach the full code:

```
#include "itensor/all.h"
#include "myclass.h"
#include <complex>
#include <cmath>
#include<iostream>
#include<stdlib.h>
#include<time.h>
#include <vector>
using namespace itensor;
using std::vector;
int main()
{
int N = 100;
auto sites = SpinHalf(N,{"ConserveQNs=",false});
auto state = InitState(sites);
auto psi0 = randomMPS(state);
auto ampo = AutoMPO(sites);
std ::vector<double> J{0.1, 0.2, 0.3, 0.5, 0.7, 0.8, 1, 2, 3, 4, 5, 10};
for(int l=0; l < J.size(); l++)
{
for(int j = 1; j <= N-1; ++j)
{
if (j%2 == 1)
{
ampo += -1,"S+",j,"S-",j+1;
ampo += -1,"S+",j+1,"S-",j;
}
else (j%2 ==0);
{
ampo += -J[l],"S+",j+1,"S-",j;
ampo += -J[l],"S+",j,"S-",j+1;
}
}
auto H = toMPO(ampo);
auto sweeps = Sweeps(5);
sweeps.maxdim() = 10,40,50,60,100;
sweeps.cutoff() = 1E-10;
auto [energy,psi] = dmrg(H,psi0,sweeps,{"Quiet",true});
auto psidag = dag(prime(psi,"Site"));
for(int i = 1; i <= N; ++i)
{
auto k = sites(i);
auto n = ITensor(k,prime(k));
n.set(k(2),prime(k)(2), 1);
psi.position(i);
auto phi = psi(i);
auto bra_phi = dag(prime(phi,"Site"));
auto exp_n = elt(bra_phi * n * phi);
print(" ",exp_n);
}
}
return 0;
}
```

I always get the following:

Thanks a lot for your time and patience to go through all of this !

Francesco