Hi Miles,

Thanks for your offer to help and look through my code. I have tried to clarify it as far as possible and take out all the pieces that are not (yet) relevant.

Attached is the preliminary DMRG code to calculate ground state energy and two-particle correlator Sz*i Sz*(i+1) of spin-1/2 zig-zag chain. The system is defined as follows:

(1) The hamiltonian for the system consists of NN and NNN hopping and interaction terms as well as the magnetic field term. The total number of particles in the system is always half the total number of sites. The NN hopping terms, however, are set to zero for the presumably easiest case.

(2) The interactions are dipole interactions which are functions of the geometry. As a result, there are two different NN interactions: V1odd is the NN interaction for odd site index and V1even is the NN interaction for even site index - typically they are not just different but can have different signs. (V2 is the NNN interaction, which is set to zero in the simplest case.)

The DMRG code is as follows:

## include "dmrg.h"

## include "sites/spinhalf.h"

## include "autompo.h"

## include iomanip> //To use setw()

using namespace itensor;

using namespace std;

int main()

{

for(int N = 4; N <= 20; N = N + 2)

{

auto sites = SpinHalf(N);

```
auto ampo = AutoMPO(sites);
#define pi 3.1415926535
double gamma = 120.0;
double theta = 22.0;
double V1odd = 1 - 3 * pow(cos(pi - (gamma * pi/180)/2 - theta * pi/180 ), 2);
double V1even = 1 - 3 * pow(cos((gamma * pi/180)/2 - theta * pi/180), 2);
double V2 = 0.0;
double h = 0.0;
double J1 = 0.001;
double J2 = 0.0;
for(int j = 1; j < N; ++j)
{
ampo += (-J1),"S+",j,"S-",j+1;
ampo += (-J1),"S-",j,"S+",j+1;
if(j%2 != 0)
{
ampo += V1odd,"Sz",j,"Sz",j+1;
}
else
{
ampo += V1even,"Sz",j,"Sz",j+1;
}
}
for(int j = 1; j < N-1; j = j+2)
{
ampo += J2,"S+",j,"S-",j+2;
ampo += J2,"S-",j,"S+",j+2;
ampo += V2,"Sz",j,"Sz",j+2;
}
for(int j = 1; j <= N; ++j)
{
ampo += h,"Sz",j;
}
auto H = IQMPO(ampo);
// Set the initial wavefunction matrix product state
// to be a classical Neel state
auto initState = InitState(sites);
for(int i = 1; i <= N; ++i)
{
if(i%2 == 1)
initState.set(i,"Up");
else
initState.set(i,"Dn");
}
/*
//Set the initial wavefunction matrix product state to be a ferromagnetic state
auto initState = InitState(sites);
for(int i = 1; i <= N; ++i)
{
if(i <= N/2)
initState.set(i,"Up");
else
initState.set(i,"Dn");
}
*/
auto psi = IQMPS(initState);
auto sweeps = Sweeps(5);
sweeps.maxm() = 10,20,100,100,200;
sweeps.cutoff() = 1E-10;
sweeps.noise() = 1E-7,1E-8,0;
println(sweeps);
auto energy = dmrg(psi,H,sweeps,"Quiet");
printfln("\nGround State Energy = %.10f",energy);
printfln("\nUsing psiHphi = %.10f", psiHphi(psi,H,psi) );
println("\nTotal QN of Ground State = ",totalQN(psi));
//
//To calculate the two particle correlator Sz_i Sz_j where j = i+1
//using the instructions on http://itensor.org/docs.cgi?page=formulas/correlator_mps
Print(N);
Print(V1odd);
Print(V1even);
double szi_szip1_array[N] = {0};
for(int i = 1; i < N; ++i)
{
int j = i + 1;
auto Sz_i = sites.op("Sz",i);
auto Sz_j = sites.op("Sz",j);
psi.position(i);
auto ir = commonIndex(psi.A(i),psi.A(i+1),Link);
auto C = psi.A(i)*Sz_i*dag(prime(prime(psi.A(i),Site),ir));
for(int k = i+1; k < j; ++k)
{
C *= psi.A(k);
C *= dag(prime(psi.A(k),Link));
}
C *= psi.A(j);
C *= Sz_j;
auto jl = commonIndex(psi.A(j),psi.A(j-1),Link);
//Original line found on http://itensor.org/docs.cgi?page=formulas/correlator_mps DIDN'T WORK
//C *= dag(prime(psi.A(j),jl,Site));
C *= dag(prime(prime(psi.A(j),Site),jl));
auto result = toReal(C);
szi_szip1_array[i] = result;
}
cout << left << setw(15) << "Site index i" << left << setw(15) << "<Sz_i Sz_(i+1)>" << endl;
for(int i = 1; i < N; ++i)
{
cout << left << setw(15) << i << left << setw(15) << szi_szip1_array[i] << endl;
}
}
return 0;
}
```

A small list of things that I have found so far and that might be of relevance:

(1) There is a problem in the description of interaction terms, V1odd and V1even. Each time I run the code, the value of ground state energy for a given value of N, theta, gamma and h changes, even if none of the variables (N, theta, gamma, h) changes.

(2) For gamma = 120 degrees, theta = 22 degrees and J1 = J2 = V2 =0, the two particle correlator Sz*i Sz*(i+1) should report values -0.25 for odd i and +0.25 for even i (because in this case V1odd and V1even have different signs) but that's not what my code reports. What it reports is as follows (Number on left is site index and that on right is Sz*i Sz*(i+1)):

N = 4

1 -0.25

2 0.25

3 -0.25

N = 6

1 -0.25

2 -0.25

3 -0.25

4 -0.25

5 -0.25

N = 8

1 -0.25

2 0.25

3 -0.25

4 0.25

5 -0.25

6 -0.25

7 -0.25

N = 10

1 -0.25

2 0.25

3 -0.25

4 0.25

5 -0.25

6 0.25

7 -0.25

8 0.25

9 -0.25

N = 12

1 -0.25

2 0.25

3 -0.25

4 0.25

5 -0.25

6 0.25

7 -0.25

8 -0.238734

9 -0.25

10 0.238734

11 -0.25

N = 14

1 -0.25

2 -0.25

3 -0.25

4 -0.25

5 -0.25

6 -0.25

7 -0.25

8 -0.25

9 -0.25

10 -0.25

11 -0.25

12 -0.25

13 -0.25

(3) If I turn on very small hopping (J1 = 0.001), keeping everything else exactly the same as before (i.e., gamma = 120 degrees, theta = 22 degrees, J2 = V2 = 0, etc.), the code reports correct values for Sz*i Sz*(i+1) (-0.25 for odd i and +0.25 for even i).

(4) Everything before was done setting initial wavefunction MPS to be a classical Neel state. If I set it to a ferromagnetic state and run the code with gamma = 120 degrees, theta = 22 degrees and J1 = J2 = V2 = 0, the values of Sz*i Sz*(i+1) I obtain imply that the state remains unchanged, i.e., it stays ferromagnetic (in this case, all up-spins on one side, all down-spins on the other). Even if I set J1 = 0.001, it produces no significant difference.