# Correlation Function for three Band Hubbard (Cuprate)

**URL:** <https://itensor.discourse.group/t/correlation-function-for-three-band-hubbard-cuprate/957>\
**Category:** DMRG and Numerical Methods\
**Created:** [June 21, 2023, 10:27pm UTC](https://itensor.discourse.group/t/correlation-function-for-three-band-hubbard-cuprate/957 "2023-06-21T22:27:32Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![Polat](https://avatars.discourse-cdn.com/v4/letter/p/ee7513/32.png) [@Polat](https://itensor.discourse.group/u/Polat)\
**Post date:** [June 21, 2023, 10:27pm UTC](https://itensor.discourse.group/t/correlation-function-for-three-band-hubbard-cuprate/957/1 "2023-06-21T22:27:32Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![miles](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/miles/32/6_2.png) [@miles](https://itensor.discourse.group/u/miles)\
**Post date:** [June 22, 2023, 5:03pm UTC](https://itensor.discourse.group/t/correlation-function-for-three-band-hubbard-cuprate/957/2 "2023-06-22T17:03:30Z")

</div>

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:

```auto
#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](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.

---

<div class="post-metadata">

**Author:** ![Polat](https://avatars.discourse-cdn.com/v4/letter/p/ee7513/32.png) [@Polat](https://itensor.discourse.group/u/Polat)\
**Post date:** [July 19, 2023, 3:13pm UTC](https://itensor.discourse.group/t/correlation-function-for-three-band-hubbard-cuprate/957/3 "2023-07-19T15:13:18Z")

</div>

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)

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

---

<div class="post-metadata">

**Author:** ![miles](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/miles/32/6_2.png) [@miles](https://itensor.discourse.group/u/miles)\
**Post date:** [July 20, 2023, 2:59am UTC](https://itensor.discourse.group/t/correlation-function-for-three-band-hubbard-cuprate/957/4 "2023-07-20T02:59:53Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Polat](https://avatars.discourse-cdn.com/v4/letter/p/ee7513/32.png) [@Polat](https://itensor.discourse.group/u/Polat)\
**Post date:** [July 27, 2023, 2:15pm UTC](https://itensor.discourse.group/t/correlation-function-for-three-band-hubbard-cuprate/957/5 "2023-07-27T14:15:38Z")

</div>

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:

> [@Measuring correlations](https://itensor.discourse.group/t/measuring-correlations/51):
>
> There’s a couple of links on measuring correlations in ITensor. [http://itensor.org/docs.cgi?vers=cppv3&page=formulas/correlator\_mps](http://itensor.org/docs.cgi?vers=cppv3&page=formulas/correlator_mps)[http://itensor.org/docs.cgi?vers=cppv3&page=formulas/spinless\_correlator\_mps](http://itensor.org/docs.cgi?vers=cppv3&page=formulas/spinless_correlator_mps) Those recipes work when i \< j. What should I do to measure Cdagup \* Cup or Cdagup\*Cdagdn\*Cup\*Cdn for Spinful fermions and i=j (on the same site)?

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):

d-p-d-p-d -p -d \\ | \quad \qquad \ |\qquad \quad \ |\qquad \quad \ |\\ p \quad \qquad p \qquad \quad p \qquad \quad p \\ | \quad \qquad \ |\qquad \quad \ |\qquad \quad \ |\\ d-p-d-p-d -p -d

Site numbers:  
3-4-8-9-\ \ 13 -14 -18 \\ | \quad \qquad \ |\qquad \quad \ \ \ \ |\qquad \quad \ \ \ \ \ |\\ 2 \quad \qquad 7 \qquad \quad \ \ 12 \qquad \quad \ \ 17 \\ | \quad \qquad \ |\qquad \quad \ \ \ \ |\qquad \quad \ \ \ \ \ |\\ 1-5-6-10-11 -15 -16

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