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[1]);
float U_12=atof(argv[2]);
int L = atof(argv[3]);
int max = atof(argv[4]);
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_op1=op(sites,"Aadag",i_p);
auto corf_op2=op(sites,"Adag",i_p+1);
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)
{
auto il = commonIndex(psi(j_p+1),psi(j_p),"Link");
C *=dag(prime(prime(psi(j_p+1),"Site"),il));
C *=psi(j_p);
C *=corf_op3;
auto ill = commonIndex(psi(j_p),psi(i_p+1),"Link");
C *=dag(prime(prime(psi(j_p),"Site"),ill));
C *=psi(i_p+1);
C *=corf_op2;
auto illl = commonIndex(psi(i_p+1),psi(i_p),"Link");
C *=dag(prime(prime(psi(i_p+1),"Site"),illl));
C *=psi(i_p);
C *=corf_op1;
auto ir =commonIndex(psi(i_p),psi(i_p+1),"Link");
C *=dag(prime(prime(psi(i_p),"Site"),ir));
corf=elt(C);
}
if (r>=2)
{
auto il = commonIndex(psi(j_p+1),psi(j_p),"Link");
C *=dag(prime(prime(psi(j_p+1),"Site"),il));
C *=psi(j_p);
C *=corf_op3;
auto ill = commonIndex(psi(j_p),psi(j_p-1),"Link");
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 *= dag(prime(psi(j_p-k),"Link"));
}
C *=psi(i_p+1);
C *=corf_op2;
auto illl = commonIndex(psi(i_p+1),psi(i_p),"Link");
C *=dag(prime(prime(psi(i_p+1),"Site"),illl));
C *=psi(i_p);
C *=corf_op1;
auto ir =commonIndex(psi(i_p),psi(i_p+1),"Link");
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;
auto il = commonIndex(psi(j_p+1),psi(j_p),"Link");
C *=dag(prime(prime(psi(j_p+1),"Site"),il));
C *=psi(j_p);
C *=corf_op3;
auto ill = commonIndex(psi(j_p),psi(j_p-1),"Link");
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 *= dag(prime(psi(j_p-k),"Link"));
}
C *=psi(i_p+1);
C *=corf_op2;
auto illl = commonIndex(psi(i_p+1),psi(i_p),"Link");
C *=dag(prime(prime(psi(i_p+1),"Site"),illl));
C *=psi(i_p);
C *=corf_op1;
auto ir =commonIndex(psi(i_p),psi(i_p+1),"Link");
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!