From f1c877d409f8de763f101593c8d1d0e93df5558a Mon Sep 17 00:00:00 2001 From: tmisawa Date: Mon, 15 Apr 2024 15:02:20 +0900 Subject: [PATCH] add tJ --- src/check.c | 34 ++++- src/include/DefCommon.h | 25 +-- src/include/sz.h | 119 ++++++++------- src/sz.c | 331 +++++++++++++++++++++++++++++++++++++++- 4 files changed, 441 insertions(+), 68 deletions(-) diff --git a/src/check.c b/src/check.c index a6e07c2b..7339e87b 100644 --- a/src/check.c +++ b/src/check.c @@ -51,7 +51,7 @@ int check(struct BindStruct *X){ FILE *fp; - long unsigned int i,tmp_sdim; + long unsigned int i,j,tmp_sdim; int NLocSpn,NCond,Nup,Ndown; long unsigned int u_tmp; long unsigned int tmp; @@ -130,6 +130,37 @@ int check(struct BindStruct *X){ comb_sum +=comb_up*comb_down; } break; + + case tJ: + comb_up = Binomial(Ns, X->Def.Nup, comb, Ns); + comb_down = Binomial(Ns-X->Def.Nup, X->Def.Ndown, comb, Ns); + comb_sum = comb_up*comb_down; + break; + + case tJNConserved: + comb_sum=0; + if(X->Def.Ne > X->Def.Nsite){ + iMinup = X->Def.Ne-X->Def.Nsite; + iAllup = X->Def.Nsite; + } + + for(i=iMinup; i<= iAllup; i++){ + comb_up = Binomial(Ns, i, comb, Ns); + comb_down = Binomial(Ns-i, X->Def.Ne-i, comb, Ns-i); + comb_sum += comb_up*comb_down; + } + break; + + case tJGC: + comb_sum=0; + for(i=0; i<= X->Def.Ne; i++){ + comb_up = Binomial(Ns, i, comb, Ns); + for(j=0; j<= X->Def.Ne; j++){ + comb_down = Binomial(Ns-i, j, comb, Ns-i); + comb_sum += comb_up*comb_down; + } + } + break; case Kondo: //idim_max @@ -195,7 +226,6 @@ int check(struct BindStruct *X){ } break; case Spin: - if(X->Def.iFlgGeneralSpin ==FALSE){ if(X->Def.Nup+X->Def.Ndown != X->Def.Nsite){ fprintf(stderr, " 2Sz is incorrect.\n"); diff --git a/src/include/DefCommon.h b/src/include/DefCommon.h index ead33b56..af4c30e2 100644 --- a/src/include/DefCommon.h +++ b/src/include/DefCommon.h @@ -27,17 +27,20 @@ #define cTPQ 5 /*!< CalcType is canonical TPQ*/ /*!< CalcModel */ -#define NUM_CALCMODEL 6 /*!< Number of model types defined by CalcModel in calcmodfile. Note: HubbardNConserved is not explicitly defined in calcmod file and thus not counted. SpinlessFermion and SpinlessFermionGC are not yet supported*/ -#define Hubbard 0 /*!< CalcModel is Hubbard model.*/ -#define Spin 1 /*!< CalcModel is Spin system.*/ -#define Kondo 2 /*!< CalcModel is Kondo model.*/ -#define HubbardGC 3 /*!< CalcModel is GrandCanonical Hubbard model.*/ -#define SpinGC 4 /*!< CalcModel is GrandCanonical Spin system.*/ -#define KondoGC 5 /*!< CalcModel is GrandCanonical Kondo model.*/ -#define HubbardNConserved 6 /*!< CalcModel is Hubbard model under particle number conserved. This symmetry is automatically introduced by not defining 2Sz in a modpara file.*/ -#define SpinlessFermion 7 /*!< CalcModel is GrandCanonical Spinless fermion model.*/ -#define SpinlessFermionGC 8 /*!< CalcModel is GrandCanonical Spinless fermionGC model.*/ -#define KondoNConserved 20 /*!< CalcModel is Kondo model under particle number conserved.*/ +#define NUM_CALCMODEL 12 /*!< Number of model types defined by CalcModel in calcmodfile. Note: HubbardNConserved is not explicitly defined in calcmod file and thus not counted. SpinlessFermion and SpinlessFermionGC are not yet supported*/ +#define Hubbard 0 /*!< CalcModel is Hubbard model.*/ +#define Spin 1 /*!< CalcModel is Spin system.*/ +#define Kondo 2 /*!< CalcModel is Kondo model.*/ +#define HubbardGC 3 /*!< CalcModel is GrandCanonical Hubbard model.*/ +#define SpinGC 4 /*!< CalcModel is GrandCanonical Spin system.*/ +#define KondoGC 5 /*!< CalcModel is GrandCanonical Kondo model.*/ +#define HubbardNConserved 6 /*!< CalcModel is Hubbard model under particle number conserved. This symmetry is automatically introduced by not defining 2Sz in a modpara file.*/ +#define SpinlessFermion 7 /*!< CalcModel is GrandCanonical Spinless fermion model.*/ +#define SpinlessFermionGC 8 /*!< CalcModel is GrandCanonical Spinless fermionGC model.*/ +#define tJ 9 /*!< CalcModel is Canonical tJ model.*/ +#define tJGC 10 /*!< CalcModel is GrandCanonical tJ model.*/ +#define tJNConserved 11 /*!< CalcModel is Canonical w/o Sz conservation tJ model*/ +#define KondoNConserved 12 /*!< CalcModel is Kondo model under particle number conserved.*/ /*!< OutputMode */ #define NUM_OUTPUTMODE 2 /*!< Number of output mode.*/ diff --git a/src/include/sz.h b/src/include/sz.h index 95616b67..d3512ec4 100644 --- a/src/include/sz.h +++ b/src/include/sz.h @@ -26,15 +26,25 @@ int omp_sz( long unsigned int *list_jb_ ); +int omp_tJ( + long unsigned int ib, + long unsigned int ihfbit, + struct BindStruct *X, + long unsigned int *list_1_, + long unsigned int *list_2_1_, + long unsigned int *list_2_2_, + long unsigned int *list_jb_ +); + int omp_sz_hacker( - long unsigned int ib, - long unsigned int ihfbit, - struct BindStruct *X, - long unsigned int *list_1_, - long unsigned int *list_2_1_, - long unsigned int *list_2_2_, - long unsigned int *list_jb_ - ); + long unsigned int ib, + long unsigned int ihfbit, + struct BindStruct *X, + long unsigned int *list_1_, + long unsigned int *list_2_1_, + long unsigned int *list_2_2_, + long unsigned int *list_jb_ +); int omp_sz_Kondo( long unsigned int ib, @@ -57,35 +67,35 @@ int omp_sz_KondoNConserved( ); int omp_sz_KondoGC( - long unsigned int ib, - long unsigned int ihfbit, - struct BindStruct *X, - long unsigned int *list_1_, - long unsigned int *list_2_1_, - long unsigned int *list_2_2_, - long unsigned int *list_jb_ - ); + long unsigned int ib, + long unsigned int ihfbit, + struct BindStruct *X, + long unsigned int *list_1_, + long unsigned int *list_2_1_, + long unsigned int *list_2_2_, + long unsigned int *list_jb_ +); int omp_sz_spin( - long unsigned int ib, - long unsigned int ihfbit, - unsigned int N, - struct BindStruct *X, - long unsigned int *list_1_, - long unsigned int *list_2_1_, - long unsigned int *list_2_2_, - long unsigned int *list_jb_ - ); + long unsigned int ib, + long unsigned int ihfbit, + unsigned int N, + struct BindStruct *X, + long unsigned int *list_1_, + long unsigned int *list_2_1_, + long unsigned int *list_2_2_, + long unsigned int *list_jb_ +); int omp_sz_spin_hacker( - long unsigned int ib, - long unsigned int ihfbit, - unsigned int N, - struct BindStruct *X, - long unsigned int *list_1_, - long unsigned int *list_2_1_, - long unsigned int *list_2_2_, - long unsigned int *list_jb_ + long unsigned int ib, + long unsigned int ihfbit, + unsigned int N, + struct BindStruct *X, + long unsigned int *list_1_, + long unsigned int *list_2_1_, + long unsigned int *list_2_2_, + long unsigned int *list_jb_ ); int omp_sz_Kondo_hacker( @@ -100,30 +110,30 @@ int omp_sz_Kondo_hacker( int omp_sz_GeneralSpin( - long unsigned int ib, - long unsigned int ihfbit, - struct BindStruct *X, - long unsigned int *list_1_, - long unsigned int *list_2_1_, - long unsigned int *list_2_2_, - long int *list_2_1_Sz_, - long int *list_2_2_Sz_, - long unsigned int *list_jb_ - ); + long unsigned int ib, + long unsigned int ihfbit, + struct BindStruct *X, + long unsigned int *list_1_, + long unsigned int *list_2_1_, + long unsigned int *list_2_2_, + long int *list_2_1_Sz_, + long int *list_2_2_Sz_, + long unsigned int *list_jb_ +); long int Binomial( - int n, - int k, - long int **comb, - int Nsite - ); + int n, + int k, + long int **comb, + int Nsite +); int sz( - struct BindStruct *X, - long unsigned int *list_1_, - long unsigned int *list_2_1_, - long unsigned int *list_2_2_ - ); + struct BindStruct *X, + long unsigned int *list_1_, + long unsigned int *list_2_1_, + long unsigned int *list_2_2_ +); int Read_sz( struct BindStruct *X, @@ -147,4 +157,7 @@ void calculate_jb_Hubbard(struct BindStruct *X,long unsigned int *list_jb, long void calculate_jb_Hubbard_Hacker(struct BindStruct *X,long unsigned int *list_jb, long unsigned int ihfbit, unsigned int N2); void calculate_jb_HubbardNCoserved(struct BindStruct *X,long unsigned int *list_jb, long unsigned int ihfbit, unsigned int N2); void calculate_jb_HubbardNCoserved_Hacker(struct BindStruct *X,long unsigned int *list_jb, long unsigned int ihfbit, unsigned int N2); +void calculate_jb_tJ(struct BindStruct *X,long unsigned int *list_jb, long unsigned int ihfbit, unsigned int N2); +void calculate_jb_tJNConserved(struct BindStruct *X,long unsigned int *list_jb, long unsigned int ihfbit, unsigned int N2); +void calculate_jb_tJGC(struct BindStruct *X,long unsigned int *list_jb, long unsigned int ihfbit, unsigned int N2); /*[e] func. for refactoring */ diff --git a/src/sz.c b/src/sz.c index 551abd07..75b48aa9 100644 --- a/src/sz.c +++ b/src/sz.c @@ -129,6 +129,9 @@ int sz( case HubbardGC: case HubbardNConserved: case Hubbard: + case tJGC: + case tJNConserved: + case tJ: N2=2*X->Def.Nsite; idim = pow(2.0,N2); break; @@ -164,6 +167,9 @@ int sz( switch(X->Def.iCalcModel){ case HubbardNConserved: case Hubbard: + case tJGC: + case tJNConserved: + case tJ: case KondoGC: case Kondo: case KondoNConserved: @@ -332,6 +338,33 @@ int sz( } } break; + case tJ: + calculate_jb_tJ(X,list_jb,ihfbit,N2); + TimeKeeper(X, cFileNameSzTimeKeep, cOMPSzMid, "a"); + TimeKeeper(X, cFileNameTimeKeep, cOMPSzMid, "a"); + icnt = 0; + for(ib=0;ibCheck.sdim;ib++){ + icnt += omp_sz(ib,ihfbit, X, list_1_, list_2_1_, list_2_2_, list_jb); + } + break; + case tJNConserved: + calculate_jb_tJNConserved(X,list_jb,ihfbit,N2); + TimeKeeper(X, cFileNameSzTimeKeep, cOMPSzMid, "a"); + TimeKeeper(X, cFileNameTimeKeep, cOMPSzMid, "a"); + icnt = 0; + for(ib=0;ibCheck.sdim;ib++){ + icnt += omp_sz(ib,ihfbit, X, list_1_, list_2_1_, list_2_2_, list_jb); + } + break; + case tJGC: + calculate_jb_tJGC(X,list_jb,ihfbit,N2); + TimeKeeper(X, cFileNameSzTimeKeep, cOMPSzMid, "a"); + TimeKeeper(X, cFileNameTimeKeep, cOMPSzMid, "a"); + icnt = 0; + for(ib=0;ibCheck.sdim;ib++){ + icnt += omp_sz(ib,ihfbit, X, list_1_, list_2_1_, list_2_2_, list_jb); + } + break; case Spin: if(X->Def.iFlgGeneralSpin==FALSE){ hacker = X->Def.read_hacker; @@ -379,7 +412,6 @@ int sz( icnt+=omp_sz_GeneralSpin(ib,ihfbit,X, list_1_, list_2_1_, list_2_2_, list_2_1_Sz, list_2_2_Sz,list_jb); } } - break; default: exitMPI(-1); @@ -472,6 +504,145 @@ long int Binomial(int n,int k,long int **comb,int Nsite){ return comb[n][k]; } +/** + * @brief calculating restricted Hilbert space for the tJ systems + * + * @param[in] ib upper half bit of i + * @param[in] ihfbit 2^(Ns/2) + * @param[in] X + * @param[out] list_1_ list_1_[icnt] = i : i is divided into ia and ib (i=ib*ihfbit+ia) + * @param[out] list_2_1_ list_2_1_[ib] = jb + * @param[out] list_2_2_ list_2_2_[ia] = ja : icnt=jb+ja + * @param[in] list_jb_ list_jb_[ib] = jb + * + * @return + * @author Takahiro Misawa (The University of Tokyo) + */ +int omp_sz_tJ( + long unsigned int ib, //!<[in] + long unsigned int ihfbit, //!<[in] + struct BindStruct *X, //!<[in] + long unsigned int *list_1_, //!<[out] + long unsigned int *list_2_1_,//!<[out] + long unsigned int *list_2_2_,//!<[out] + long unsigned int *list_jb_ //!<[in] + ) +{ + long unsigned int i,j; + long unsigned int ia,ja,jb; + long unsigned int div_down, div_up; + long unsigned int num_up,num_down; + long unsigned int tmp_num_up,tmp_num_down; + int check_doublon_1,check_doublon_2; + + jb = list_jb_[ib]; + i = ib*ihfbit; + + num_up = 0; + num_down = 0; + check_doublon_1 = 1; + ja = 1; + for(j=0;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]; + check_doublon_1 = div_up^div_down; + if (check_doublon_1==0){ + break; + } + num_up += div_up; + num_down += div_down; + } + + if(check_doublon_1==1){ + tmp_num_up = num_up; + tmp_num_down = num_down; + + if(X->Def.iCalcModel==tJ){ + for(ia=0;iaCheck.sdim;ia++){ + i = ia; + num_up = tmp_num_up; + num_down = tmp_num_down; + check_doublon_2 = 1; + for(j=0;jDef.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]; + check_doublon_2 = div_up^div_down; + if (check_doublon_2==0){ + break; + } + num_up += div_up; + num_down += div_down; + if(check_doublon_2==1){ + if(num_up == X->Def.Nup && num_down == X->Def.Ndown ){ + list_1_[ja+jb]=ia+ib*ihfbit; + list_2_1_[ia]=ja+1; + list_2_2_[ib]=jb+1; + ja+=1; + } + } + } + } + }else if(X->Def.iCalcModel==tJNConserved){ + for(ia=0;iaCheck.sdim;ia++){ + i = ia; + num_up = tmp_num_up; + num_down = tmp_num_down; + check_doublon_2 = 1; + for(j=0;jDef.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 (check_doublon_2==0){ + break; + } + num_up += div_up; + num_down += div_down; + } + if(check_doublon_2==1){ + if( (num_up+num_down) == X->Def.Ne){ + list_1_[ja+jb]=ia+ib*ihfbit; + list_2_1_[ia]=ja+1; + list_2_2_[ib]=jb+1; + ja+=1; + } + } + } + }else if(X->Def.iCalcModel==tJGC){ + for(ia=0;iaCheck.sdim;ia++){ + i = ia; + num_up = tmp_num_up; + num_down = tmp_num_down; + check_doublon_2 = 1; + for(j=0;jDef.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 (check_doublon_2==0){ + break; + } + num_up += div_up; + num_down += div_down; + } + if(check_doublon_2==1){ + list_1_[ja+jb]=ia+ib*ihfbit; + list_2_1_[ia]=ja+1; + list_2_2_[ib]=jb+1; + ja+=1; + } + } + } + } + ja=ja-1; + return ja; +} + + /** * @brief calculating restricted Hilbert space for Hubbard systems * @@ -1903,4 +2074,160 @@ void calculate_jb_HubbardNCoserved_Hacker(struct BindStruct *X,long unsigned int } }//omp parallel free_lui_1d_allocate(jbthread); -} \ No newline at end of file +} + +void calculate_jb_tJ(struct BindStruct *X,long unsigned int *list_jb, long unsigned int ihfbit, unsigned int N2){ + /*[s] this part can not be parallelized*/ + long unsigned int jb = 0,div_up,div_down,i,j,tmp_1,tmp_2; + long int **comb; + int num_up,num_down,check_doublon; + int tmp_res,all_up,all_down; + comb = li_2d_allocate(X->Def.Nsite+1,X->Def.Nsite+1); + + for(long unsigned ib=0;ibCheck.sdim;ib++){ // sdim = 2^(N/2) + list_jb[ib] = jb; + i = ib*ihfbit; + check_doublon = 1; + //[s] counting # of up and down electrons + num_up = 0; + num_down = 0; + for(j=0;j<=N2-2;j+=2){ + div_up = i & X->Def.Tpow[2*j];// even -> up spin + div_up = div_up/X->Def.Tpow[2*j]; + div_down = i & X->Def.Tpow[2*j+1];//odd -> down spin + div_down = div_down/X->Def.Tpow[2*j+1]; + check_doublon = div_up^div_down; + if (check_doublon==0){ + break; + } + num_up += div_up; + num_down += div_down; + } + //[e] counting # of up and down electrons + + /* eg of even sites: 4site-> |DU|DU || |DU|DU|*/ + /* eg of odd sites: 3site-> |DU|D || |U|DU|*/ + /* all_up -> # of up sites in the lower half of bits*/ + /* all_down -> # of down sites in the lower half of bits*/ + if (check_doublon==1){ + tmp_res = X->Def.Nsite%2; // even Ns-> 0, odd Ns -> 1 + all_up = (X->Def.Nsite+tmp_res)/2; + all_down = (X->Def.Nsite-tmp_res)/2; + tmp_1 = Binomial(all_up,X->Def.Nup-num_up,comb,all_up); + tmp_2 = Binomial(all_down-(X->Def.Nup-num_up),X->Def.Ndown-num_down,comb,all_down-(X->Def.Nup-num_up));/* tJ all_down-(X->Def.Nup-num_up)*/ + jb += tmp_1*tmp_2; + } + } + free_li_2d_allocate(comb); + /*[e] this part can not be parallelized*/ +} + +void calculate_jb_tJNConserved(struct BindStruct *X,long unsigned int *list_jb, long unsigned int ihfbit, unsigned int N2){ + /*[s] this part can not be parallelized*/ + long unsigned int jb = 0,div_up,div_down,i,j,tmp_1,tmp_2; + long int **comb; + int num_up,num_down,check_doublon; + int tmp_res,all_up,all_down; + int iSpnup, iMinup,iAllup; + + comb = li_2d_allocate(X->Def.Nsite+1,X->Def.Nsite+1); + iMinup = 0; + iAllup = X->Def.Ne; + if(X->Def.Ne > X->Def.Nsite){ + iMinup = X->Def.Ne-X->Def.Nsite; + iAllup = X->Def.Nsite; + } + + for(long unsigned ib=0;ibCheck.sdim;ib++){ // sdim = 2^(N/2) + list_jb[ib] = jb; + i = ib*ihfbit; + check_doublon = 1; + //[s] counting # of up and down electrons + num_up = 0; + num_down = 0; + for(j=0;j<=N2-2;j+=2){ + div_up = i & X->Def.Tpow[2*j];// even -> up spin + div_up = div_up/X->Def.Tpow[2*j]; + div_down = i & X->Def.Tpow[2*j+1];//odd -> down spin + div_down = div_down/X->Def.Tpow[2*j+1]; + check_doublon = div_up^div_down; + if (check_doublon==0){ + break; + } + num_up += div_up; + num_down += div_down; + } + //[e] counting # of up and down electrons + + /* eg of even sites: 4site-> |DU|DU || |DU|DU|*/ + /* eg of odd sites: 3site-> |DU|D || |U|DU|*/ + /* all_up -> # of up sites in the lower half of bits*/ + /* all_down -> # of down sites in the lower half of bits*/ + if (check_doublon==1){ + tmp_res = X->Def.Nsite%2; // even Ns-> 0, odd Ns -> 1 + all_up = (X->Def.Nsite+tmp_res)/2; + all_down = (X->Def.Nsite-tmp_res)/2; + + for(iSpnup=iMinup; iSpnup<= iAllup; iSpnup++){ + tmp_1 = Binomial(all_up,iSpnup-num_up,comb,all_up); + tmp_2 = Binomial(all_down-(iSpnup-num_up),X->Def.Ne-(iSpnup+num_down),comb,all_down-(iSpnup-num_up));/* tJ all_down-(iSpnup-num_up)*/ + jb += tmp_1*tmp_2; + } + } + } + free_li_2d_allocate(comb); + /*[e] this part can not be parallelized*/ +} + +void calculate_jb_tJGC(struct BindStruct *X,long unsigned int *list_jb, long unsigned int ihfbit, unsigned int N2){ + /*[s] this part can not be parallelized*/ + long unsigned int jb = 0,div_up,div_down,i,j,tmp_1,tmp_2; + long int **comb; + int num_up,num_down,check_doublon; + int tmp_res,all_up,all_down; + int iSpnup,iSpndown,iMinup,iAllup; + + comb = li_2d_allocate(X->Def.Nsite+1,X->Def.Nsite+1); + + for(long unsigned ib=0;ibCheck.sdim;ib++){ // sdim = 2^(N/2) + list_jb[ib] = jb; + i = ib*ihfbit; + check_doublon = 1; + //[s] counting # of up and down electrons + num_up = 0; + num_down = 0; + for(j=0;j<=N2-2;j+=2){ + div_up = i & X->Def.Tpow[2*j];// even -> up spin + div_up = div_up/X->Def.Tpow[2*j]; + div_down = i & X->Def.Tpow[2*j+1];//odd -> down spin + div_down = div_down/X->Def.Tpow[2*j+1]; + check_doublon = div_up^div_down; + if (check_doublon==0){ + break; + } + num_up += div_up; + num_down += div_down; + } + //[e] counting # of up and down electrons + + /* eg of even sites: 4site-> |DU|DU || |DU|DU|*/ + /* eg of odd sites: 3site-> |DU|D || |U|DU|*/ + /* all_up -> # of up sites in the lower half of bits*/ + /* all_down -> # of down sites in the lower half of bits*/ + if (check_doublon==1){ + tmp_res = X->Def.Nsite%2; // even Ns-> 0, odd Ns -> 1 + all_up = (X->Def.Nsite+tmp_res)/2; + all_down = (X->Def.Nsite-tmp_res)/2; + for(iSpnup=0; iSpnup<= all_up; iSpnup++){ + tmp_1 = Binomial(all_up,iSpnup,comb,all_up); + for(iSpndown=0; iSpndown<= all_down; iSpndown++){ + tmp_2 = Binomial(all_down-iSpnup,iSpndown,comb,all_down-iSpnup); + jb += tmp_1*tmp_2; + } + } + } + } + free_li_2d_allocate(comb); + /*[e] this part can not be parallelized*/ +} +