Why does DMRG give different results for spin-1/2 and spin-1 cases?

Question

Dear forum,
I have a 1D chain where half of the sites are singly occupied while the other half are vacant and double occupancy on any site is not allowed. I modeled this system with spin-1/2 and spin-1 site sets to compare my results. For both the cases, the values of the parameters in the Hamiltonian are such that the the ground state is expected to be antiferromagnetic, which is the superposition of the two Neel states |101010...> and |010101...>.

Spin-1/2 model: The spin operators are defined such that S^z |1> = +1 |1> and S^z |0> = -1 |0> where |1> = |spin up> refers to a singly occupied site while |0> = |spin down> refers to a vacant site.

Spin-1 model: The spin operators are defined such that S^z |2> = +1 |2>, S^z |1> = 0 |1> and S^z|0> = -1 |0> where |2> = |spin up>, |1> = |spin 0> and |0> = |spin down> refer to doubly occupied, singly occupied and vacant sites respectively. There is an onsite potential U with a value large enough to avoid double occupancy on any site.

Here are the results I obtained:

For spin-1/2, ITensor gives < S^z _ j > = +-1 alternating along the 1D chain while < S^z _ j S^z_(j+1) > = -1 for all j, where j is the site index. This implies that DMRG returns only one of the Neel states as the ground state.

For spin-1, ITensor gives < S^z _ j > = -0.5 and < S^z _ j S^z_(j+1) > = -0.25 for all j. This implies that the ground state returned by DMRG is a superposition of the two Neel states |101010...> and |010101...>.

Both the spin-1/2 and spin-1 cases seem to make sense, but while the former returns only one of the two states in the superposition, the latter seems to give an average. Could you please explain why DMRG works differently for the two cases?

Thanks in advance,
Niraj

Comments

Hi Niraj,
Is this with IQMPS/IQMPO quantum number conserving code or with MPS/MPO non-QN code? That may clarify what's going on. Thanks -

---

Hi Miles,
          Thanks for the question and sorry for the late reply. I used IQMPS/IQMPO quantum number conserving code for both the spin-1/2 and spin-1 models.

MOREOVER, HERE'S A SMALL EDIT TO MY QUESTION:

For the AFM ground state in the spin-1/2 model, the values of < S^z _ j S^z _ (j+1)> are nearly equal to  -1 and for the spin-1 model, they are nearly equal to -0.25.

---

Hi Niraj,
One other question: what is the initial state you use for the two cases? That might be behind what's happening here.

Ultimately, to answer your question I may need you to send me a copy of your code so I can see more details and run it myself.

But meanwhile, one thing you can try is this: try putting a small potential/field term on the first or the last site (or both) that favors one of the two ground states. Doing so for the S=1 case may resolve the issue.

Also, one other comment is that if this system can be thought of as being a "hard core" boson system, then we have a site set for that called "Spinless" which defines operators such as "Adag" and "A" and "N" so that might be a more natural Hilbert space to use.

---

Hi Miles,
       Thanks again for the reply. Here are my answers:

(1) For the initial state, instead of using the Neel state |10101010...>, I used the one where singly occupied site (denoted by 1) and vacant site (denoted by 0) are randomly placed in the lattice.

(2) As you suggested, I added a small field term on the first site for the spin-1 case. The only difference I found was in the values of the the correlation between the first few sites.

(3) Likewise, I modeled this system with "Spinless" siteset, I obtained results nearly identical to the spin-1/2 case.

(4) Here's the code for the spin-1 case:

double pi = M_PI;
       
int N = 100;

double c = 0.1;
        
double h = 1;
   
double J_factor = -0.5;

int gamma = 180;
   
int theta = 0;

double B_o = 1;
   
double B_e = B_o;
       
double B_2 = B_o * exp(-c*(2 * sin ((gamma*pi/180)/2) - 1) );
   
double U = 30;
   
double V_e = 1 - 3 * pow(cos(pi - (gamma*pi/180)/2 - theta*pi/180), 2);
   
double V_o = 1 - 3 * pow(cos((gamma*pi/180)/2 - theta*pi/180), 2);
   
double V_2 = (1/pow((2*(1 - cos(gamma*pi/180))), 1.5)) * (1 - 3 * pow(cos(pi/2 - theta*pi/180), 2));
       
double a_o = V_o/B_o;
   
double a_e = V_e/B_e;
   
double a_2 = V_2/B_2;
         
auto sites = SpinOne(N);
   
auto ampo = AutoMPO(sites);

for(int j = 1; j <= N; ++j)
    {
    ampo += U, "Sz2",j;
    ampo += U,"Sz",j;
    }               
       
for(int j = 1; j < N; ++j)
    {
       if(j%2 != 0)
           {
           ampo += B_o * J_factor,"S+",j,"S-",j+1;
        ampo += B_o * J_factor,"S-",j,"S+",j+1;
           ampo += B_o * a_o,"Sz",j,"Sz",j+1;
           }
       else
           {
           ampo += B_e * J_factor,"S+",j,"S-",j+1;
        ampo += B_e * J_factor,"S-",j,"S+",j+1;
           ampo += B_e * a_e,"Sz",j,"Sz",j+1;
           }
    }

for(int j = 1; j < N-1; ++j)
    {
    ampo += B_2 * J_factor,"S+",j,"S-",j+2;
    ampo += B_2 * J_factor,"S-",j,"S+",j+2;
        ampo += B_2 * a_2,"Sz",j,"Sz",j+2;
    }        

for(int j = 1; j <= N; ++j)
    {
    ampo += h,"Sz",j;
    }
                
auto H = IQMPO(ampo);
   
auto initState = InitState(sites);
       
cout << endl << "#Vector v(N):" << endl;
       
int v[N] = {1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0};
       
cout << endl << endl << "#Initial wavefunction MPS mapped to vector v(N):" << endl;
   
for(int i = 1; i <= N; ++i)
       {
       if(v[i-1] == 1)
           {
           initState.set(i,"Z0");
        cout << "0 ";
        }
    else if(v[i-1] == 0)
           {
           initState.set(i,"Dn");
        cout << "-1 ";
        }   
       }
      
auto psi = IQMPS(initState);
           
printfln("Initial energy = %.5f", overlap(psi,H,psi) );

auto sweeps = Sweeps(7);
sweeps.maxm() = 10,20,80,200,300,400,400;
sweeps.minm() = 1,1,1,1,1,1,1;
sweeps.cutoff() = 1E-5,1E-6,1E-7,1E-8,1E-8,1E-8;
sweeps.niter() = 4,3,3,2,2,2,2;
sweeps.noise() = 1E-5,1E-5,1E-8,1E-9,1E-10,1E-10,1E-10;
println(sweeps);

auto energy = dmrg(psi,H,sweeps,{"Quiet=",true});