Skip to content

Commit

Permalink
add KondoNConserved
Browse files Browse the repository at this point in the history
  • Loading branch information
@tmisawa
tmisawa committed Apr 14, 2024
1 parent cbec6cb commit 262091b
Show file tree
Hide file tree
Showing 2 changed files with 94 additions and 114 deletions.
28 changes: 28 additions & 0 deletions src/check.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -157,6 +157,33 @@ int check(struct BindStruct *X){
comb_sum += comb_1*comb_2*comb_3;
}
break;
case KondoNConserved:
//idim_max
// calculation of dimension
// Nup = u_loc+u_cond
// Ndown = d_loc+d_cond
// NLocSpn = u_loc+d_loc
// Ncond = Nsite-NLocSpn
// idim_max = \sum_{u_loc=0}^{u_loc=Nup}
// Binomial(NLocSpn,u_loc)
// *Binomial(NCond,Nup-u_loc)
// *Binomial(NCond,Ndown+u_loc-NLocSpn)
//comb_1 = Binomial(NLocSpn,u_loc)
//comb_2 = Binomial(NCond,Nup-u_loc)
//comb_3 = Binomial(NCond,Ndown+u_loc-NLocSpn)
NLocSpn = X->Def.NLocSpn;
NCond = X->Def.Nsite-NLocSpn;
comb_1 = pow(2,NLocSpn);//Tpow is not defined
comb_sum = 0;
for(int tmp_Nup=0;tmp_Nup<=NCond;tmp_Nup++){
comb_2 = Binomial(NCond,tmp_Nup,comb,NCond);
comb_3 = Binomial(NCond,NCond-tmp_Nup,comb,NCond);
comb_sum += comb_1*comb_2*comb_3;
printf("tmp %ld %ld %ld\n",comb_sum,comb_2,comb_3);
}
printf("Ne=%d NLocSpn=%d NCond=%d comb_1=%ld comb_2=%ld comb_3=%ld comb_sum=%ld\n",X->Def.Ne,NLocSpn,NCond,comb_1,comb_2,comb_3,comb_sum);
break;

case KondoGC:
comb_sum = 1;
NCond = X->Def.Nsite-X->Def.NLocSpn;
Expand Down Expand Up @@ -373,6 +400,7 @@ int check(struct BindStruct *X){
case HubbardNConserved:
case Hubbard:
case Kondo:
case KondoNConserved:
for(i=1;i<=2*X->Def.Nsite-1;i++){
u_tmp=u_tmp*2;
X->Def.Tpow[i]=u_tmp;
Expand Down
180 changes: 66 additions & 114 deletions src/sz.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -219,7 +219,7 @@ int sz(
icnt = X->Def.Tpow[X->Def.Nsite-1]*X->Def.SiteToBit[X->Def.Nsite-1];
}
break;
case Hubbard:
case Hubbard:
hacker = X->Def.read_hacker;
if(hacker==0){
calculate_jb_Hubbard(X,list_jb,ihfbit,N2);
Expand Down Expand Up @@ -297,15 +297,17 @@ int sz(
calculate_jb_KondoNConserved(X, list_jb,ihfbit);
TimeKeeper(X, cFileNameSzTimeKeep, cOMPSzMid, "a");
TimeKeeper(X, cFileNameTimeKeep, cOMPSzMid, "a");


icnt = 0;
#pragma omp parallel for default(none) \
reduction(+:icnt) private(ib) firstprivate(ihfbit, N2, X) \
shared(list_1_, list_2_1_, list_2_2_, list_jb)
for(ib=0;ib<X->Check.sdim;ib++){
icnt+=omp_sz_KondoNConserved(ib,ihfbit, X, list_1_, list_2_1_, list_2_2_, list_jb);
//printf("ib = %ld icnt =%ld \n",ib,icnt);
}
//printf("AAA icnt=%ld\n",icnt);
break;
case Kondo:
calculate_jb_Kondo(X, list_jb,ihfbit);
TimeKeeper(X, cFileNameSzTimeKeep, cOMPSzMid, "a");
Expand Down Expand Up @@ -383,6 +385,7 @@ int sz(
exitMPI(-1);
}
i_max=icnt;
//printf("BBB i_max %ld icnt=%ld\n",i_max,icnt);
//fprintf(stdoutMPI, "Debug: Xicnt=%ld \n",icnt);
TimeKeeper(X, cFileNameSzTimeKeep, cOMPSzFinish, "a");
TimeKeeper(X, cFileNameTimeKeep, cOMPSzFinish, "a");
Expand All @@ -399,6 +402,7 @@ int sz(

//Error message
if(i_max!=X->Check.idim_max){
//printf("DDD %lu %lu \n",i_max, X->Check.idim_max);
fprintf(stderr, "%s", cErrSz);
fprintf(stderr, cErrSz_ShowDim, i_max, X->Check.idim_max);
strcpy(sdt_err,cFileNameErrorSz);
Expand Down Expand Up @@ -712,81 +716,55 @@ int omp_sz_KondoNConserved(
long unsigned int num_up,num_down;
long unsigned int tmp_num_up,tmp_num_down;
int icheck_loc;

jb = list_jb_[ib];
i = ib*ihfbit;

num_up = 0;
num_down = 0;
icheck_loc = 1;
for(j=X->Def.Nsite/2; j< X->Def.Nsite ;j++){
div_up = i & X->Def.Tpow[2*j];
div_up = div_up/X->Def.Tpow[2*j];
div_down = i & X->Def.Tpow[2*j+1];
div_down = div_down/X->Def.Tpow[2*j+1];

if(X->Def.LocSpn[j] == ITINERANT){
num_up += div_up;
num_down += div_down;
}else{
num_up += div_up;
num_down += div_down;
if(X->Def.Nsite%2==1 && j==(X->Def.Nsite/2)){
icheck_loc= icheck_loc;
}else{
icheck_loc = icheck_loc*(div_up^div_down);// exclude doubllly ocupited site
}
}
}

ja=1;
tmp_num_up = num_up;
tmp_num_down = num_down;
if(icheck_loc ==1){
for(ia=0;ia<X->Check.sdim;ia++){
i = ia;
num_up = tmp_num_up;
num_down = tmp_num_down;
icheck_loc=1;
for(j=0;j<(X->Def.Nsite+1)/2;j++){
div_up = i & X->Def.Tpow[2*j];
div_up = div_up/X->Def.Tpow[2*j];
div_down = i & X->Def.Tpow[2*j+1];
div_down = div_down/X->Def.Tpow[2*j+1];

if(X->Def.LocSpn[j] == ITINERANT){
num_up += div_up;
num_down += div_down;
}else{
num_up += div_up;
num_down += div_down;
if(X->Def.Nsite%2==1 && j==(X->Def.Nsite/2)){
icheck_loc= icheck_loc;
}else{
icheck_loc = icheck_loc*(div_up^div_down);// exclude doubllly ocupited site
}
}
}
if(icheck_loc == 1 && X->Def.LocSpn[X->Def.Nsite/2] != ITINERANT && X->Def.Nsite%2==1){
div_up = ia & X->Def.Tpow[X->Def.Nsite-1];
div_up = div_up/X->Def.Tpow[X->Def.Nsite-1];
div_down = (ib*ihfbit) & X->Def.Tpow[X->Def.Nsite];
div_down = div_down/X->Def.Tpow[X->Def.Nsite];
icheck_loc= icheck_loc*(div_up^div_down);
}

//if(num_up == X->Def.Nup && num_down == X->Def.Ndown && icheck_loc==1){
if(icheck_loc==1){
list_1_[ja+jb] = ia+ib*ihfbit;
list_2_1_[ia] = ja+1;
list_2_2_[ib] = jb+1;
//printf("DEBUG: rank=%d, list_1[%d]=%d, list_2_1_[%d]=%d, list_2_2_[%d]=%d\n", myrank, ja+jb, list_1_[ja+jb], ia, list_2_1[ia], ib, list_2_2[ib]);
ja += 1;
}
}
jb = list_jb_[ib];
i = ib*ihfbit;
num_up = 0;
num_down = 0;
icheck_loc = 1;
for(j=X->Def.Nsite/2; j< X->Def.Nsite ;j++){
div_up = i & X->Def.Tpow[2*j];
div_up = div_up/X->Def.Tpow[2*j];
div_down = i & X->Def.Tpow[2*j+1];
div_down = div_down/X->Def.Tpow[2*j+1];

if(X->Def.LocSpn[j] == ITINERANT){
num_up += div_up;
num_down += div_down;
}
ja = ja-1;
return ja;
}

ja=1;
//printf("DEBUG: num_up %d num_down %d Ne %d NLocSpn %d \n",num_up,num_down,X->Def.Ne,X->Def.NLocSpn);
if(num_up+num_down == X->Def.Ne-X->Def.NLocSpn){
for(ia=0;ia<X->Check.sdim;ia++){
i=ia;
num_up = 0;
num_down = 0;
icheck_loc=1;
//for(j=0;j<(X->Def.Nsite+1)/2;j++){ /* is this correct?*/
for(j=0;j<(X->Def.Nsite)/2;j++){ /*we assume local spins exist i<X->Def.Nsite/2, X->Def.Nsite should be a even number.*/
div_up = i & X->Def.Tpow[2*j];
div_up = div_up/X->Def.Tpow[2*j];
div_down = i & X->Def.Tpow[2*j+1];
div_down = div_down/X->Def.Tpow[2*j+1];

if (div_up^div_down == 0){
icheck_loc = 0; /* if double occupied site exists, break!*/
}
}
if (icheck_loc ==1){
list_1_[ja+jb] = ia+ib*ihfbit;
list_2_1_[ia] = ja+1;
list_2_2_[ib] = jb+1;
//printf("DEBUG: ja+jb %lu ja %lu jb %lu ia %lu ib %lu\n",ja+jb,ja,jb,ia,ib);
ja+=1;
}
}
}
//printf("XXX ib %lu ja %lu\n",ib,ja);
ja=ja-1;
return ja;
}


Expand Down Expand Up @@ -1541,7 +1519,7 @@ void calculate_jb_Kondo(struct BindStruct *X, long unsigned int *list_jb, long u
comb = li_2d_allocate(X->Def.Nsite+1,X->Def.Nsite+1);
N_all_up = X->Def.Nup;
N_all_down = X->Def.Ndown;
printf("DEBUG N_all_up N_all_down Ne %d %d %d\n",N_all_up,N_all_down,X->Def.Ne);
//printf("DEBUG N_all_up N_all_down Ne %d %d %d\n",N_all_up,N_all_down,X->Def.Ne);
fprintf(stdoutMPI, cStateNupNdown, N_all_up,N_all_down);

jb = 0;
Expand All @@ -1551,14 +1529,15 @@ void calculate_jb_Kondo(struct BindStruct *X, long unsigned int *list_jb, long u
num_loc += 1;
}
}
printf("DEBUG num_loc %d \n",num_loc);
//printf("DEBUG num_loc %d \n",num_loc);
for(ib=0;ib<X->Check.sdim;ib++){ //sdim = 2^(N/2)
list_jb[ib] = jb;
//printf("DEBUG ib %lu jb %lu %lu\n",ib,jb,X->Large.SizeOflistjb);
i = ib*ihfbit; // ihfbit=pow(2,((Nsite+1)/2))
num_up = 0;
num_down = 0;
icheck_loc = 1;
printf("DEBUG num_loc %lu %d i=%lu %lu\n",ib,X->Def.Nsite,i,ihfbit);
//printf("DEBUG num_loc %lu %d i=%lu %lu\n",ib,X->Def.Nsite,i,ihfbit);

for(j=X->Def.Nsite/2; j< X->Def.Nsite ;j++){
//printf("DEBUG j=%lu %d %d\n",j,X->Def.LocSpn[j],ITINERANT);
Expand Down Expand Up @@ -1594,7 +1573,7 @@ void calculate_jb_Kondo(struct BindStruct *X, long unsigned int *list_jb, long u
all_up = (X->Def.Nsite)/2-all_loc;
all_down = (X->Def.Nsite)/2-all_loc;
}
printf("all_loc %d all_up %d all_down %d \n",all_loc,all_up,all_down);
//printf("all_loc %d all_up %d all_down %d \n",all_loc,all_up,all_down);

for(num_loc_up=0; num_loc_up <= all_loc; num_loc_up++){
tmp_1 = Binomial(all_loc, num_loc_up, comb, all_loc);
Expand Down Expand Up @@ -1624,12 +1603,10 @@ void calculate_jb_KondoNConserved(struct BindStruct *X, long unsigned int *list_
int num_up, num_down, icheck_loc;
int all_up, all_down, all_loc,tmp_res,num_loc,num_loc_up,num_loc_down;
long unsigned int tmp_1,tmp_2,tmp_3;
long int **comb;
int all_cond,all_cond_up,all_cond_down,num_cond;

comb = li_2d_allocate(X->Def.Nsite+1,X->Def.Nsite+1);
Ne = X->Def.Ne;
printf("DEBUG Ne %d\n",X->Def.Ne);
//printf("DEBUG Ne %d\n",X->Def.Ne);

jb = 0;
num_loc = 0;
Expand All @@ -1638,15 +1615,16 @@ void calculate_jb_KondoNConserved(struct BindStruct *X, long unsigned int *list_
num_loc += 1;
}
}
printf("DEBUG num_loc %d \n",num_loc);
//printf("DEBUG num_loc %d \n",num_loc);
all_loc = X->Def.NLocSpn;
for(ib=0;ib<X->Check.sdim;ib++){ //sdim = 2^(N/2)
list_jb[ib] = jb;
//printf("DEBUG ib %lu jb %lu %lu Ne %d\n",ib,jb,X->Large.SizeOflistjb,X->Def.Ne);
i = ib*ihfbit; // ihfbit=pow(2,((Nsite+1)/2))
//printf("DEBUG num_loc %lu %d i=%lu %lu\n",ib,X->Def.Nsite,i,ihfbit);

num_up = 0;
num_down = 0;
icheck_loc = 1;
printf("DEBUG num_loc %lu %d i=%lu %lu\n",ib,X->Def.Nsite,i,ihfbit);

for(j=X->Def.Nsite/2; j< X->Def.Nsite ;j++){
//printf("DEBUG j=%lu %d %d\n",j,X->Def.LocSpn[j],ITINERANT);
div_up = i & X->Def.Tpow[2*j];
Expand All @@ -1656,38 +1634,12 @@ void calculate_jb_KondoNConserved(struct BindStruct *X, long unsigned int *list_
if(X->Def.LocSpn[j] == ITINERANT){
num_up += div_up;
num_down += div_down;
}else{
/*this part may not be used*/
num_up += div_up;
num_down += div_down;
if(X->Def.Nsite%2==1 && j==(X->Def.Nsite/2)){ // odd site
icheck_loc = icheck_loc;
ihfSpinDown = div_down;
if(div_down ==0){
num_up += 1;
}
}else{
icheck_loc = icheck_loc*(div_up^div_down);// exclude empty or doubly occupied site
}
}
}

if(icheck_loc == 1){ // itinerant of local spins without holon or doublon
tmp_res = X->Def.Nsite%2; // even Ns-> 0, odd Ns -> 1
all_loc = X->Def.NLocSpn-num_loc; // # of local spins
all_cond = X->Def.Nsite-all_loc;
all_cond_up = X->Def.Nsite/2-num_up;
all_cond_down = X->Def.Nsite/2-num_down;

tmp_1 = X->Def.Tpow[all_loc];
for(num_cond=0; num_cond <= all_cond; num_cond++){
tmp_2 = Binomial(all_cond_up,all_cond, comb, all_cond);
tmp_3 = Binomial(all_cond_down,all_cond, comb, all_cond);
jb += tmp_1*tmp_2*tmp_3;
}
if (num_up+num_down==X->Def.Ne-all_loc){
jb += X->Def.Tpow[all_loc];
}
}
free_li_2d_allocate(comb);
}

void calculate_jb_KondoGC(struct BindStruct *X, int num_loc, long unsigned int *list_jb, long unsigned int ihfbit){
Expand Down

0 comments on commit 262091b

Please sign in to comment.