Issues with my DRMG code

Question

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 Szi 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 Szi 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 Szi 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 Szi 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 Szi 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.

Answer 1

Hi Niraj,
Sorry about the slow reply! I looked over your code and found a few things that could be causing the problems you report. However, when I ran it I always got the same energy value so I was not able to reproduce that particular issue. What exact parameters were you using when you saw the energy come out different for each run of the code?

Regarding your other issues / questions:

1. Did you intend for the J2 and V2 interactions to skip every other bond in addition to being next-nearest-neighbor? Maybe you don't actually want to increase j by 2 (j += 2) for every step of the loop. Maybe you do; just wanted to point this out.

2. I would not recommend combining different sizes in a single loop in your code. This may be an ok thing to do later once you have thoroughly debugged your code and studied the best DMRG accuracy parameters to use to ensure convergence. But since you are having issues, it is best to fix a single system size and figure out the best way to use DMRG for this system size. For longer systems you may need to do more DMRG sweeps or keep more states, etc.

3. The reason the code C *= dag(prime(psi.A(j),jl,Site)) didn't compile is that you are using ITensor version 1.x and that particular style of code only works for version 2.x. The website examples mostly refer to the new version now, since it would be difficult to maintain the site for both versions simultanously. Please upgrade to version 2.x if you can since it is faster and has some nicer features and is better maintained.

Ok those are my main recommendations for now. If you are still seeing different energies for each run, please simplify your code to use only one size; check that your Hamiltonian is really correct and what you want; and then tell me exactly which Hamiltonian and sweeping/accuracy parameters you are using and we can discuss more.

Best,
Miles

Comments

Hi Miles,
       Thank you so much for the reply. I apologize for a slow reply - I was on a vacation. As you suggested I upgraded to version 2.x of ITensor, fixed the system size to a single value (N = 10) and ran the code. Here are the details:

(1) The Hamiltonian is:
$$ H =  –J_1 \sum_{j= 1}^{N-1} (S^+_{j} S^-_{j+1} + S^-_{j} S^+_{j+1}) \ \
             – J_2 \sum_{j= 1}^{N-2} (S^+_{j} S^-_{j+2} + S^-_{j} S^+_{j+2}) \\
             + \sum_{j=1}^{N-1}((V_{1even} + V_{1odd})/2 + (-1)^j (V_{1even} –  V_{1odd} )/2) S^z_{j} S^z_{j+1} \\
            + V_2 \sum_{j=1}^{N-2}S^z_{j} S^z_{j+2} – h \sum_{j=1}^{N} S^z_j  $$
 
(2) The sweeping parameters are:
     auto sweeps = Sweeps(5);
    sweeps.maxm() = 10,20,100,100,200;
    sweeps.cutoff() = 1E-10;
    sweeps.noise() = 1E-7,1E-8,0;

(3) When I set J1 to 0.001, I don't see any change in the ground state energy for N = 10,  J2 = 0.0, V2 = 0.0 and h = 0.0 (or h = 10.0 or h = 100.0). The value I obtain each time I run the code is -2.0402566357. But when I set J1 to zero keeping all the other parameters the same, the energy changes. The values I got on 3 different occasions were -0.7459242582, -1.6088071017 and -2.0402485234 respectively. Why does the code mess around when I turn off the nearest-neighbor hopping?

(4) The V1odd and V1even are  nearest-neighbor dipole interactions for odd and even site indices respectively while V2 is the next-nearest-neighbor dipole interaction. When all of them are turned on, every dipole in the chain must interact with every neighbor dipole and next-neighbor dipole. Is there any problem if I want to increase j by 2 for every step of the loop? If so, is there a better way to do this?

Thanks again,
Niraj

---

Hi Niraj,
My guess is the problem with changing energies is connected to the incorrect definition of the J2 terms in your code. Please carefully step through the for-loop for the J2 terms that you wrote and you will see the issue. It only couples the odd-numbered sites the way you wrote it in the code. In the formula you wrote both even and odd numbered sites are connected by J2 terms. When you turn J1 to zero, then your Hamiltonian becomes somewhat pathological, especially with no J2 terms connecting to the even sites. So probably the different energies you're seeing are sensitivity to initial conditions that can occur in DMRG for Hamiltonians which are "classical" or otherwise pathological, non-typical quantum Hamiltonians.

See for example this other recent question: http://itensor.org/support/105/sweeps-system-size-clarification-choice-intial-state-dmrg

---

Hi Miles,
        Thank you so much for pointing out the mistake in my code regarding the J2 term, I fixed it. Even after this,  the ground state energy reported by DMRG changes each time I run the code as long as J1 is set equal to zero. However, if J1 is set equal to a very small value (for instance 1E-10, which can be considered to be zero for practical purposes), the ground state energy remains unchanged. Like you suggested, this might be due to the nature of the Hamiltonian. I will continue working on it and come back if I have any more questions. I truly appreciate all your help.

Regards,
Niraj