# how to measure correlator of two point with 4 operators in two component bose-hubbard model(c++)

+1 vote

Dear ITensor users

I have tried to obtain two correlation functions in a two component bose-hubbard model. But I got nothing about correlation function in my output files and the output file reported errors as the following:

Mismatched QN Index from set 1 (dim=6|id=345|"n=32,Site,Boson")' <Out>
1: 1 QN({"Nb",0})
2: 1 QN({"Nb",1})
3: 1 QN({"Nb",2})
4: 1 QN({"Nb",3})
5: 1 QN({"Nb",4})
6: 1 QN({"Nb",5})
Mismatched QN Index from set 2 (dim=6|id=345|"n=32,Site,Boson")' <Out>
1: 1 QN({"Nb",0})
2: 1 QN({"Nb",1})
3: 1 QN({"Nb",2})
4: 1 QN({"Nb",3})
5: 1 QN({"Nb",4})
6: 1 QN({"Nb",5})

Mismatched QN Index arrows
From line 912, file itensor.cc


You can obtain the Hamiltonian of two component bose-hubbard model from equation (1) in this paper: https://doi.org/10.1103/PhysRevA.80.023619
Or you can see it in the question posted before in the community：http://itensor.org/support/2973/two-component-bose-hubbard-model?show=2973#q2973

I want to measure two different correlation functions of two point with 4 operators as the followings.
1. $$\braket{\hat{a}^{\dagger }i\hat{b}^{\dagger}i\hat{a}j\hat{b}j }$$
2. $$\braket{\hat{a}^{\dagger }i\hat{b}i\hat{a}j\hat{b}^{\dagger }j }$$
They were also presented in the paper I mentioned above.(equation (13) (14))

In my codes, I made a copy of boson.h and called it aboson.h, and I renamed all operators in aboson.h as "Aa", "Aadag", "Na". Also, I rename its class as "ABosonSite".

Here is my codes for groundstate

float t=atof(argv);
float U_12=atof(argv);
int L = atof(argv);
int max = atof(argv);
using twocom=MixedSiteSet<ABosonSite,BosonSite>;
autosites=twocom(L,{"ConserveQNs",true,"ConserveQNs",true,
"ConserveNb",true,"ConserveNb",true,"MaxOcc=",max});
auto sweeps = Sweeps(8);
sweeps.maxdim() = 100,200,400,400,800,800,1000,1000;
sweeps.cutoff() = 1E-12;
float U = 1;
float t_1=t;
float t_2=t;
auto ampo = AutoMPO(sites);
for (int i = 1;i <= L-3; i += 2)
{
ampo += -U/2, "Na", i;
ampo += U/2, "Na", i,"Na",i;
ampo += -t_1, "Aa", i, "Aadag", i + 2;
ampo += -t_1, "Aadag", i, "Aa", i + 2;
}
ampo += -U/2, "Na", L-1;
ampo += U/2, "Na", L-1,"Na",L-1;
for (int i = 2;i <= L-2; i += 2)
{
ampo += -U/2, "N", i;
ampo += U/2, "N", i,"N",i;
ampo += -t_2, "A", i, "Adag", i + 2;
ampo += -t_2, "Adag", i, "A", i + 2;
}
ampo += -U/2, "N", L;
ampo += U/2, "N", L,"N",L;
for (int i = 1;i < L; i += 2)
{
ampo +=U_12, "Na", i,'N',i+1;
ampo += V, "N", i, "N", i + 1;
}
auto H = toMPO(ampo);
auto state = InitState(sites);
for (int i = 1; i <= L; i += 4)
{
state.set(i,"1");
}
for (int i = 4; i <= L; i += 4)
{
state.set(i,"1");
}
PrintData(sweeps);
auto psi0 = randomMPS(state);
auto [energy, psi] = dmrg(H, psi0, sweeps, {"Quiet=",true})


And here is my codes for one of two correlation functions.(sorry for my long codes)
To minimize the effect of boundary condition, I place i and j symmetrically around the center.

int j_p;
j_p=L/2-1;
int i_p;
i_p=L/2-1;
for (auto r : range1(1,L))
{
auto corf_op3=op(sites,"Aa",j_p);
auto corf_op4=op(sites,"A",j_p+1);
if (r%2==1)
{
j_p=j_p+2;
psi.position(j_p+1);
ITensor C=psi(j_p+1);
C *= corf_op4;
double corf=0;
if (r==1)
{
C *=dag(prime(prime(psi(j_p+1),"Site"),il));
C *=psi(j_p);
C *=corf_op3;
C *=dag(prime(prime(psi(j_p),"Site"),ill));
C *=psi(i_p+1);
C *=corf_op2;
C *=dag(prime(prime(psi(i_p+1),"Site"),illl));
C *=psi(i_p);
C *=corf_op1;
C *=dag(prime(prime(psi(i_p),"Site"),ir));
corf=elt(C);
}
if (r>=2)
{
C *=dag(prime(prime(psi(j_p+1),"Site"),il));
C *=psi(j_p);
C *=corf_op3;
C *=dag(prime(prime(psi(j_p),"Site"),ill));
for(int k=1;k<(j_p-i_p-1);k++)
{
C *= psi(j_p-k);
}
C *=psi(i_p+1);
C *=corf_op2;
C *=dag(prime(prime(psi(i_p+1),"Site"),illl));
C *=psi(i_p);
C *=corf_op1;
C *=dag(prime(prime(psi(i_p),"Site"),ir));
corf=elt(C);
}
fprintf(fp2,"%i,%i,%i,%f\n",r,i_p,j_p,corf);
}
else
{
i_p=i_p-2;
psi.position(j_p+1);
ITensor C=psi(j_p+1);
C *= corf_op4;
double corf=0;
C *=dag(prime(prime(psi(j_p+1),"Site"),il));
C *=psi(j_p);
C *=corf_op3;
C *=dag(prime(prime(psi(j_p),"Site"),ill));
for(int k=1;k<(j_p-i_p-1);k++)
{
C *= psi(j_p-k);
}
C *=psi(i_p+1);
C *=corf_op2;
C *=dag(prime(prime(psi(i_p+1),"Site"),illl));
C *=psi(i_p);
C *=corf_op1;
C *=dag(prime(prime(psi(i_p),"Site"),ir));
corf=elt(C);
fprintf(fp2,"%i,%i,%i,%f\n",r,i_p,j_p,corf);
}
}


Another correlation function is the same thing but with different operators.
After calculation, I got energies and density distributions but without correlation functions.
Are there some mistakes I make in my codes?

Looking forward to your response and suggestions.
Thank you very much!

Hi Wang,
Thanks for the question. Yes, this kind of error means there is a mistake in your code, where you are trying to contract two ITensors that have a matching Index, but where the Index has the same arrow direction. As you may know, in order for two indices to contract, they have to have opposite arrows. Usually this kind of error results from not putting a "dag" somewhere in your code where it is needed. So please look for any places that you should need to take the Hermitian conjugate (or "dag") of an ITensor that may be missing.

Also it would be very helpful in this case to either use a debugger or a printing statement to identify which line of your code is causing the error to happen.

If you spend some more time looking for the error and can't find it, or have a question about how to work with dag and Index arrows, please just ask by posting a comment below.

Best regards,
Miles

commented by (180 points)
OK, thank you, Miles.
One of my questions is how to measure correlators in a mixed site or are there some examples and guides to  do this.
Additionally, as I mentioned above, I made a copy of boson.h and rename its operators, so except renaming these, are there other steps I need to do to measure correlators?
commented by (70.1k points)
Hi Wang,
We have recently added a function called correlationMatrix to ITensor that can help you with measuring two-point correlation functions. It should work for mixed site sets as well as regular ones.

Here is an example of how to use it:

the returned value corr will be a vector<vector<Real>>.