Correlation Function for three Band Hubbard (Cuprate)

Polat · June 21, 2023, 10:27pm

Hello together,

I try to write the correlation function of the three Band Hubbard model in c++ and have some problems to find the right way.

As example:
One part of the correlation function looks like this <c_{i,2,\downarrow} c_{i,1,\uparrow} c^{\dag}_{j,1,\uparrow} c^{\dag}_{j,2,\downarrow} >

Assume I calc. bevor psi with DMRG.
My first try out with the correlation function in c++ looks like this:

        auto si = 1;
        auto sj = 2;
        
        psi.position(d_list1[si]);
        auto psidag = dag(psi);
        psidag.prime("Link");
        
        
        auto li_1 = leftLinkIndex(psi,d_list2[si]);
        auto lj = rightLinkIndex(psi,d_list2[sj]);

        auto op_up = op(sites,"Cup",d_list1[si]);
        auto op_dup = op(sites,"Cdagup",d_list1[sj]);
        auto op_dn = op(sites,"Cdn",d_list2[si]);
        auto op_ddn = op(sites,"Cdagdn",d_list2[sj]);
        
        auto C = prime(psi(d_list2[si]),li_1)*op_up;
        
        C *= psi(d_list1[sj])*op_dn;
        C *= prime(psidag(d_list1[sj]), "Site");
        
        C *= psi(d_list1[sj])*op_ddn;
        C *= prime(psidag(d_list1[sj]), "Site");
        
        C *= prime(psi(sj),lj)*op_dup;
        C *= prime(psidag(sj),"Site");
        
        auto result = elt(C); //or eltC(C) if expecting complex
        cout<<"Correlation Value: "<<result<<endl;

After running the Code I got “zsh: killed”.

Has someone of you a similar correlation function code in c++ that I can use as orientation for my problem?

Thank you an BR

Gökmen Polat

miles · June 22, 2023, 5:03pm

Hi Gökmen,
Regarding the error you are getting, probably the contractions are not happening in the code in quite the way you are planning on paper (i.e. there is a bug in your code). I would suggest printing out the indices of each of the intermediate tensors and examining them to see if the contractions are succeeding. I like to use the following pattern in C++ to do this:

#include "itensor/util/print_macro.h" //add this line to the top of your code
//...
Print(C);
PAUSE;

where the Print macro is a nice shortcut that is like doing println("C = ",C); but saves you from having to write the name of the variable twice. And the PAUSE; macro pauses your code at that line until you hit the Enter key.

Regarding similar code, you might find this tutorial page helpful:
http://itensor.org/docs.cgi?vers=cppv3&page=tutorials/correlations

Some of the details may be different from yours and the code slightly out of date, but it still could help you.

Polat · July 19, 2023, 3:13pm

Hi Miles, thank you for answering.

I try like the tutorial with this Code:

        auto si = 0;
        auto sj = 2;
        
        psi.position(d_list2[si]);
        auto psidag = dag(psi);
        psidag.prime("Link");
        
        
        auto li = leftLinkIndex(psi,d_list1[si]);
        auto lj = rightLinkIndex(psi,d_list2[sj]);


        auto op_up = op(sites,"Cup",d_list1[si]);
        auto op_dup = op(sites,"Cdagup",d_list1[sj]);
        auto op_dn = op(sites,"Cdn",d_list2[si]);
        auto op_ddn = op(sites,"Cdagdn",d_list2[sj]);
        
        auto C = prime(psi(d_list1[si]),li);
        C *= op_up;
        C *= prime(psidag(d_list1[si]),"Site");
        Print(C);
        cout<<"#######_1_####### k="<<d_list1[si]<<endl;
        
        C *= psi(d_list1[si]+1);
        C *= psidag(d_list1[si]+1);
        Print(C);
        cout<<"#######_2_####### k="<<d_list1[si]+1<<endl;
        
        C *= psi(d_list2[si]);
        C *= op_dn;
        C *= psidag(d_list2[si]);
        Print(C);
        cout<<"#######_3_####### k="<<d_list2[si]<<endl;

        
        for(int k=d_list2[si]+1; k<d_list1[sj]; k++)
        {
            C *= psi(k);
            C *= psidag(k);
            Print(C);
            cout<<"#######_4_####### k="<<k<<endl;
        }
        
        C *= psi(d_list1[sj]);
        C *= op_dup;
        C *= psidag(d_list1[sj]);
        
        
        Print(C);
        cout<<"#######_5_####### k="<<d_list1[sj]<<endl;
        
        C *= psi(d_list1[sj]+1);
        C *= psidag(d_list1[sj]+1);
        Print(C);
        cout<<"#######_6_####### k="<<d_list1[sj]+1<<endl;
        
        C *= prime(psi(d_list2[sj]),lj);
        Print(C);
        cout<<"#######_7.1_####### k="<<d_list2[sj]<<endl;
        C *= op_ddn;
        Print(C);
        cout<<"#######_7.2_####### k="<<d_list2[sj]<<endl;
        C *= prime(psidag(d_list2[sj]),"Site");
        
        Print(C);
        cout<<"#######_7.3_####### k="<<d_list2[sj]<<endl;
        
        
        auto result = elt(C); //or eltC(C) if expecting complex
        cout<<"Correlation Value: "<<result<<endl;

And I got this Output: “Wrong number of IndexVals passed to elt/eltC (expected 0, got 4)”
I tried a lot of different ways, do you see some failures in my code?

Tank you and BR

Gökmen

The last Outputs…
#######7.1####### k=13
C =
ITensor ord=7:
(dim=4|id=408|“n=3,Elec,Site”)’
1: 1 QN({“Nf”,0,-1},{“Sz”,0})
2: 1 QN({“Nf”,1,-1},{“Sz”,1})
3: 1 QN({“Nf”,1,-1},{“Sz”,-1})
4: 1 QN({“Nf”,2,-1},{“Sz”,0})
(dim=4|id=408|“n=3,Elec,Site”)
1: 1 QN({“Nf”,0,-1},{“Sz”,0})
2: 1 QN({“Nf”,1,-1},{“Sz”,1})
3: 1 QN({“Nf”,1,-1},{“Sz”,-1})
4: 1 QN({“Nf”,2,-1},{“Sz”,0})
(dim=4|id=438|“n=11,Elec,Site”)’
1: 1 QN({“Nf”,0,-1},{“Sz”,0})
2: 1 QN({“Nf”,1,-1},{“Sz”,1})
3: 1 QN({“Nf”,1,-1},{“Sz”,-1})
4: 1 QN({“Nf”,2,-1},{“Sz”,0})
(dim=4|id=438|“n=11,Elec,Site”)
1: 1 QN({“Nf”,0,-1},{“Sz”,0})
2: 1 QN({“Nf”,1,-1},{“Sz”,1})
3: 1 QN({“Nf”,1,-1},{“Sz”,-1})
4: 1 QN({“Nf”,2,-1},{“Sz”,0})
(dim=200|id=667|“l=12,Link”)’
1: 1 QN({“Nf”,0,-1},{“Sz”,0})
2: 3 QN({“Nf”,1,-1},{“Sz”,1})
3: 3 QN({“Nf”,1,-1},{“Sz”,-1})
4: 14 QN({“Nf”,2,-1},{“Sz”,0})
5: 6 QN({“Nf”,2,-1},{“Sz”,2})
6: 22 QN({“Nf”,3,-1},{“Sz”,1})
7: 6 QN({“Nf”,2,-1},{“Sz”,-2})
8: 22 QN({“Nf”,3,-1},{“Sz”,-1})
9: 30 QN({“Nf”,4,-1},{“Sz”,0})
10: 5 QN({“Nf”,3,-1},{“Sz”,3})
11: 16 QN({“Nf”,4,-1},{“Sz”,2})
12: 15 QN({“Nf”,5,-1},{“Sz”,1})
13: 5 QN({“Nf”,3,-1},{“Sz”,-3})
14: 15 QN({“Nf”,4,-1},{“Sz”,-2})
15: 15 QN({“Nf”,5,-1},{“Sz”,-1})
16: 7 QN({“Nf”,6,-1},{“Sz”,0})
17: 1 QN({“Nf”,4,-1},{“Sz”,4})
18: 2 QN({“Nf”,5,-1},{“Sz”,3})
19: 3 QN({“Nf”,6,-1},{“Sz”,2})
20: 1 QN({“Nf”,4,-1},{“Sz”,-4})
21: 4 QN({“Nf”,5,-1},{“Sz”,-3})
22: 2 QN({“Nf”,6,-1},{“Sz”,-2})
23: 1 QN({“Nf”,7,-1},{“Sz”,1})
24: 1 QN({“Nf”,7,-1},{“Sz”,-1})
(dim=200|id=328|“l=13,Link”)’
1: 1 QN({“Nf”,0,-1},{“Sz”,0})
2: 4 QN({“Nf”,1,-1},{“Sz”,1})
3: 4 QN({“Nf”,1,-1},{“Sz”,-1})
4: 16 QN({“Nf”,2,-1},{“Sz”,0})
5: 8 QN({“Nf”,2,-1},{“Sz”,2})
6: 27 QN({“Nf”,3,-1},{“Sz”,1})
7: 8 QN({“Nf”,2,-1},{“Sz”,-2})
8: 27 QN({“Nf”,3,-1},{“Sz”,-1})
9: 29 QN({“Nf”,4,-1},{“Sz”,0})
10: 5 QN({“Nf”,3,-1},{“Sz”,3})
11: 13 QN({“Nf”,4,-1},{“Sz”,2})
12: 13 QN({“Nf”,5,-1},{“Sz”,1})
13: 5 QN({“Nf”,3,-1},{“Sz”,-3})
14: 14 QN({“Nf”,4,-1},{“Sz”,-2})
15: 13 QN({“Nf”,5,-1},{“Sz”,-1})
16: 5 QN({“Nf”,6,-1},{“Sz”,0})
17: 1 QN({“Nf”,4,-1},{“Sz”,4})
18: 2 QN({“Nf”,5,-1},{“Sz”,3})
19: 1 QN({“Nf”,6,-1},{“Sz”,2})
20: 1 QN({“Nf”,4,-1},{“Sz”,-4})
21: 2 QN({“Nf”,5,-1},{“Sz”,-3})
22: 1 QN({“Nf”,6,-1},{“Sz”,-2})
(dim=4|id=1|“n=13,Elec,Site”)’
1: 1 QN({“Nf”,0,-1},{“Sz”,0})
2: 1 QN({“Nf”,1,-1},{“Sz”,1})
3: 1 QN({“Nf”,1,-1},{“Sz”,-1})
4: 1 QN({“Nf”,2,-1},{“Sz”,0})
{norm=0.01 (QDense Real)}

#######7.2####### k=13
C =
ITensor ord=4:
(dim=4|id=408|“n=3,Elec,Site”)’
1: 1 QN({“Nf”,0,-1},{“Sz”,0})
2: 1 QN({“Nf”,1,-1},{“Sz”,1})
3: 1 QN({“Nf”,1,-1},{“Sz”,-1})
4: 1 QN({“Nf”,2,-1},{“Sz”,0})
(dim=4|id=408|“n=3,Elec,Site”)
1: 1 QN({“Nf”,0,-1},{“Sz”,0})
2: 1 QN({“Nf”,1,-1},{“Sz”,1})
3: 1 QN({“Nf”,1,-1},{“Sz”,-1})
4: 1 QN({“Nf”,2,-1},{“Sz”,0})
(dim=4|id=438|“n=11,Elec,Site”)’
1: 1 QN({“Nf”,0,-1},{“Sz”,0})
2: 1 QN({“Nf”,1,-1},{“Sz”,1})
3: 1 QN({“Nf”,1,-1},{“Sz”,-1})
4: 1 QN({“Nf”,2,-1},{“Sz”,0})
(dim=4|id=438|“n=11,Elec,Site”)
1: 1 QN({“Nf”,0,-1},{“Sz”,0})
2: 1 QN({“Nf”,1,-1},{“Sz”,1})
3: 1 QN({“Nf”,1,-1},{“Sz”,-1})
4: 1 QN({“Nf”,2,-1},{“Sz”,0})
{norm=0.00 (QDense Real)}

#######7.3####### k=13
inds() =
(dim=4|id=408|“n=3,Elec,Site”)’
1: 1 QN({“Nf”,0,-1},{“Sz”,0})
2: 1 QN({“Nf”,1,-1},{“Sz”,1})
3: 1 QN({“Nf”,1,-1},{“Sz”,-1})
4: 1 QN({“Nf”,2,-1},{“Sz”,0})
(dim=4|id=408|“n=3,Elec,Site”)
1: 1 QN({“Nf”,0,-1},{“Sz”,0})
2: 1 QN({“Nf”,1,-1},{“Sz”,1})
3: 1 QN({“Nf”,1,-1},{“Sz”,-1})
4: 1 QN({“Nf”,2,-1},{“Sz”,0})
(dim=4|id=438|“n=11,Elec,Site”)’
1: 1 QN({“Nf”,0,-1},{“Sz”,0})
2: 1 QN({“Nf”,1,-1},{“Sz”,1})
3: 1 QN({“Nf”,1,-1},{“Sz”,-1})
4: 1 QN({“Nf”,2,-1},{“Sz”,0})
(dim=4|id=438|“n=11,Elec,Site”)
1: 1 QN({“Nf”,0,-1},{“Sz”,0})
2: 1 QN({“Nf”,1,-1},{“Sz”,1})
3: 1 QN({“Nf”,1,-1},{“Sz”,-1})
4: 1 QN({“Nf”,2,-1},{“Sz”,0})

From line 89, file itensor.cc

Wrong number of IndexVals passed to elt/eltC (expected 0, got 4)

miles · July 20, 2023, 2:59am

Hi, so if you are getting that kind of error, then it usually means there was a mistake in your logic of how you managed “prime levels” of ITensor indices, or similar issues like skipping a step in a contraction (leaving out a tensor you meant to contract), or thinking that an index would stay uncontracted whereas it got contracted automatically on a step you did not expect.

For these kinds of issues, a good approach can be to print out the indices of the tensor you are making (here the C tensor) and the other tensors you are contracting with it to confirm that each one has exactly the indices and prime levels you are expecting, and draw diagrams on paper to visualize each step and match it with the code.

Polat · July 27, 2023, 2:15pm

Hi Miles, I got it, and my program calculates know without error, but I think there is still a mistake, because I get implausible values.

I discovered the correlationmatrix function:

My question: is it possible to use the correlationmatrix function for i,i+1,j,j+1 → in total 4 operators on different sites ?

In my case, ladder look like this (d&p orbitals/sites):

I have to calculate this correlation function P(r)=0.5 (<\Delta^+_i \Delta_j>+<\Delta_i \Delta^+_j>) with \Delta^+_i=2^{-1/2} (d^+_{i,1,\uparrow}d^+_{i,2,\downarrow}-d^+_{i,1,\downarrow}d^+_{i,2,\uparrow}) and \Delta_i=2^{-1/2} (d_{i,2,\downarrow}d_{i,1,\uparrow}-d_{i,2,\uparrow}d_{i,1,\downarrow}).

That’s why I’m asking about the possibility to use the correlationmatrix function for my problem, it would save me time, but I don’t see a way to use it like this

BR

Gökmen

Topic		Replies	Views
Unstable behavior of the self-defined MPO's DMRG and Numerical Methods	1	185	February 3, 2023
correlation matrix of custom operators DMRG and Numerical Methods julia , dmrg , fermions	3	329	December 13, 2023
Computing out-of-time correlation function ITensor C++ Questions	0	207	March 20, 2023
question about asymmetric hubbard model ITensor Julia Questions dmrg	4	59	October 14, 2024
Problems with on site interactions Bose-Hubbard Julia DMRG and Numerical Methods julia , dmrg	2	700	March 10, 2023

Correlation Function for three Band Hubbard (Cuprate)

Related topics