diff --git a/src/CalcByLOBPCG.c b/src/CalcByLOBPCG.c index b2616869..3f43fccb 100644 --- a/src/CalcByLOBPCG.c +++ b/src/CalcByLOBPCG.c @@ -650,12 +650,14 @@ int CalcByLOBPCG( switch (X->Bind.Def.iCalcModel) { case HubbardGC: case SpinGC: - case KondoGC: case SpinlessFermionGC: initial_mode = 1; // 1 -> random initial vector break; case Hubbard: case Kondo: + case KondoGC: + case tJ: + case tJGC: case Spin: case SpinlessFermion: diff --git a/src/CalcByLanczos.c b/src/CalcByLanczos.c index 63a99088..ab601f89 100644 --- a/src/CalcByLanczos.c +++ b/src/CalcByLanczos.c @@ -69,12 +69,14 @@ int CalcByLanczos( switch(X->Bind.Def.iCalcModel){ case HubbardGC: case SpinGC: - case KondoGC: case SpinlessFermionGC: initial_mode = 1; // 1 -> random initial vector break; case Hubbard: case Kondo: + case KondoGC: + case tJ: + case tJGC: case Spin: case SpinlessFermion: diff --git a/src/CalcSpectrum.c b/src/CalcSpectrum.c index e6d59226..55a251c2 100644 --- a/src/CalcSpectrum.c +++ b/src/CalcSpectrum.c @@ -455,9 +455,12 @@ int MakeExcitedList( case HubbardGC: break; case HubbardNConserved: - case KondoGC: case Hubbard: case Kondo: + case KondoGC: + case tJ: + case tJNConserved: + case tJGC: *iFlgListModifed = TRUE; break; case Spin: @@ -466,13 +469,16 @@ int MakeExcitedList( } } else if (X->Def.NPairExcitationOperator > 0) { switch (X->Def.iCalcModel) { - case HubbardGC: case SpinGC: case HubbardNConserved: - break; + case HubbardGC: case KondoGC: + case tJNConserved: + case tJGC: + break; case Hubbard: case Kondo: + case tJ: case Spin: if (X->Def.PairExcitationOperator[0][1] != X->Def.PairExcitationOperator[0][3]) { *iFlgListModifed = TRUE; @@ -519,6 +525,8 @@ int MakeExcitedList( case HubbardGC: break; case HubbardNConserved: + case KondoNConserved:/*To be confirmed*/ + case tJNConserved:/*To be confirmed*/ if (X->Def.SingleExcitationOperator[0][2] == 1) { //cis X->Def.Ne = X->Def.NeMPI + 1; } @@ -526,9 +534,11 @@ int MakeExcitedList( X->Def.Ne = X->Def.NeMPI - 1; } break; - case KondoGC: case Hubbard: case Kondo: + case KondoGC: + case tJ: + case tJGC: if (X->Def.SingleExcitationOperator[0][2] == 1) { //cis X->Def.Ne = X->Def.NeMPI + 1; if (X->Def.SingleExcitationOperator[0][1] == 0) {//up @@ -557,13 +567,17 @@ int MakeExcitedList( } else if (X->Def.NPairExcitationOperator > 0) { X->Def.Ne=X->Def.NeMPI; switch (X->Def.iCalcModel) { - case HubbardGC: case SpinGC: case HubbardNConserved: + case HubbardGC: + case KondoNConserved:/*To be confirmed*/ + case tJNConserved:/*To be confirmed*/ break; - case KondoGC: case Hubbard: case Kondo: + case KondoGC: + case tJ: + case tJGC: if (X->Def.PairExcitationOperator[0][1] != X->Def.PairExcitationOperator[0][3]) { if (X->Def.PairExcitationOperator[0][1] == 0) {//up X->Def.Nup = X->Def.NupOrg + 1; diff --git a/src/CheckMPI.c b/src/CheckMPI.c index e95bb9f8..6aab2bb6 100644 --- a/src/CheckMPI.c +++ b/src/CheckMPI.c @@ -39,7 +39,11 @@ int CheckMPI(struct BindStruct *X/**< [inout] */) case Hubbard: case HubbardNConserved: case Kondo: + case KondoNConserved: case KondoGC: + case tJ: + case tJNConserved: + case tJGC: /**@brief
  • For Hubbard & Kondo @@ -104,6 +108,32 @@ int CheckMPI(struct BindStruct *X/**< [inout] */) break;/*case Hubbard:*/ + case tJ: + /**@brief +
  • For canonical Hubbard + DefineList::Nup, DefineList::Ndown, and DefineList::Ne should be + differerent in each PE.
    • + */ + SmallDim = myrank; + for (isite = X->Def.Nsite; isite < X->Def.NsiteMPI; isite++) { + SpinNum = SmallDim % 4; + SmallDim /= 4; + if (SpinNum == 1 /*01*/) { + X->Def.Nup -= 1; + X->Def.Ne -= 1; + } + else if (SpinNum == 2 /*10*/) { + X->Def.Ndown -= 1; + X->Def.Ne -= 1; + } + /*else if (SpinNum == 3 //11){ + X->Def.Nup -= 1; + X->Def.Ndown -= 1; + X->Def.Ne -= 2; + }*/ + } /*for (isite = X->Def.Nsite; isite < X->Def.NsiteMPI; isite++)*/ + break;/*case tJ:*/ + case HubbardNConserved: /**@brief
  • For N-conserved canonical Hubbard @@ -119,8 +149,24 @@ int CheckMPI(struct BindStruct *X/**< [inout] */) break; /*case HubbardNConserved:*/ - case KondoGC: + case tJNConserved: + case tJGC:/*is this correct?*/ + /**@brief +
  • For N-conserved canonical Hubbard + DefineList::Ne should be differerent in each PE.
    • + */ + SmallDim = myrank; + for (isite = X->Def.Nsite; isite < X->Def.NsiteMPI; isite++) { + SpinNum = SmallDim % 4; + SmallDim /= 4; + if (SpinNum == 1 /*01*/ || SpinNum == 2 /*10*/) X->Def.Ne -= 1; + //else if (SpinNum == 3 /*11*/) X->Def.Ne -= 2; + } /*for (isite = X->Def.Nsite; isite < X->Def.NsiteMPI; isite++)*/ + break; /*case HubbardNConserved:*/ + case Kondo: + case KondoNConserved: + case KondoGC: /**@brief
  • For canonical Kondo system DefineList::Nup, DefineList::Ndown, and DefineList::Ne should be @@ -155,7 +201,32 @@ int CheckMPI(struct BindStruct *X/**< [inout] */) } }/*for (isite = X->Def.Nsite; isite < X->Def.NsiteMPI; isite++)*/ } /*if (X->Def.iCalcModel == Kondo)*/ - + else if(X->Def.iCalcModel == KondoNConserved){ + SmallDim = myrank; + for (isite = X->Def.Nsite; isite < X->Def.NsiteMPI; isite++) { + SpinNum = SmallDim % 4; + SmallDim /= 4; + if (X->Def.LocSpn[isite] == ITINERANT) { + if (SpinNum == 1 /*01*/) { + //X->Def.Nup -= 1; + X->Def.Ne -= 1; + } + else if (SpinNum == 2 /*10*/) { + //X->Def.Ndown -= 1; + X->Def.Ne -= 1; + } + else if (SpinNum == 3 /*11*/) { + //X->Def.Nup -= 1; + //X->Def.Ndown -= 1; + X->Def.Ne -= 2; + } + } + else { + fprintf(stdoutMPI, "\n Stop because local spin in the inter process region\n"); + return FALSE; + } + }/*for (isite = X->Def.Nsite; isite < X->Def.NsiteMPI; isite++)*/ + } break; /*case KondoGC, Kondo*/ } /*switch (X->Def.iCalcModel) 2(inner)*/ @@ -318,11 +389,15 @@ void CheckMPI_Summary(struct BindStruct *X/**< [inout] */) { fprintf(stdoutMPI, " Site Bit\n"); for (isite = 0; isite < X->Def.Nsite; isite++) { switch (X->Def.iCalcModel) { - case HubbardGC: case Hubbard: case HubbardNConserved: + case HubbardGC: case Kondo: + case KondoNConserved: case KondoGC: + case tJ: + case tJNConserved: + case tJGC: fprintf(stdoutMPI, " %4d %4d\n", isite, 4); break; @@ -346,11 +421,15 @@ void CheckMPI_Summary(struct BindStruct *X/**< [inout] */) { fprintf(stdoutMPI, " Site Bit\n"); for (isite = X->Def.Nsite; isite < X->Def.NsiteMPI; isite++) { switch (X->Def.iCalcModel) { - case HubbardGC: case Hubbard: case HubbardNConserved: + case HubbardGC: case Kondo: + case KondoNConserved: case KondoGC: + case tJ: + case tJNConserved: + case tJGC: fprintf(stdoutMPI, " %4d %4d\n", isite, 4); break; @@ -404,11 +483,15 @@ void CheckMPI_Summary(struct BindStruct *X/**< [inout] */) { as a binary (excepting general spin) format. */ switch (X->Def.iCalcModel) { - case HubbardGC: /****************************************************/ case Hubbard: case HubbardNConserved: + case HubbardGC: /****************************************************/ case Kondo: + case KondoNConserved: case KondoGC: + case tJ: + case tJNConserved: + case tJGC: SmallDim = iproc; for (isite = X->Def.Nsite; isite < X->Def.NsiteMPI; isite++) { @@ -466,11 +549,15 @@ void CheckMPI_Summary(struct BindStruct *X/**< [inout] */) { affected by the number of processes. */ switch (X->Def.iCalcModel) { - case HubbardGC: /****************************************************/ case Hubbard: case HubbardNConserved: + case HubbardGC: /****************************************************/ case Kondo: + case KondoNConserved: case KondoGC: + case tJ: + case tJNConserved: + case tJGC: X->Def.Tpow[2 * X->Def.Nsite] = 1; for (isite = 2 * X->Def.Nsite + 1; isite < 2 * X->Def.NsiteMPI; isite++) diff --git a/src/PairEx.c b/src/PairEx.c index 5d64e29b..6b872e4a 100644 --- a/src/PairEx.c +++ b/src/PairEx.c @@ -76,9 +76,11 @@ int GetPairExcitedState iret=GetPairExcitedStateHubbardGC(X, tmp_v0, tmp_v1); break; - case KondoGC: case Hubbard: case Kondo: + case KondoGC: + case tJ: + case tJGC: iret=GetPairExcitedStateHubbard(X, tmp_v0, tmp_v1); break; diff --git a/src/SingleEx.c b/src/SingleEx.c index 4f79d340..ea56d0bb 100644 --- a/src/SingleEx.c +++ b/src/SingleEx.c @@ -41,9 +41,11 @@ int GetSingleExcitedState( iret = GetSingleExcitedStateHubbardGC(X, tmp_v0, tmp_v1); break; - case KondoGC: case Hubbard: case Kondo: + case KondoGC: + case tJ: + case tJGC: iret = GetSingleExcitedStateHubbard(X, tmp_v0, tmp_v1); break; diff --git a/src/bitcalc.c b/src/bitcalc.c index 3183aea0..89527c74 100644 --- a/src/bitcalc.c +++ b/src/bitcalc.c @@ -85,11 +85,15 @@ int GetSplitBitByModel( { int tmpNsite=Nsite; switch(iCalcModel){ + case Hubbard: case HubbardGC: - case KondoGC: case HubbardNConserved: - case Hubbard: case Kondo: + case KondoNConserved: + case KondoGC: + case tJ: + case tJNConserved: + case tJGC: tmpNsite *= 2; break; case Spin: diff --git a/src/check.c b/src/check.c index 040a50dd..2a79a034 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; @@ -129,6 +129,39 @@ int check(struct BindStruct *X){ comb_down= Binomial(Ns, X->Def.Ne-i, comb, Ns); comb_sum +=comb_up*comb_down; } + //printf("Ns %ld iMinup %d iAllup %d Ne %d comb_sum %ld\n",Ns,iMinup,iAllup,X->Def.Ne,comb_sum); + 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); + comb_sum += comb_up*comb_down; + } + //printf("Ns %ld iMinup %d iAllup %d Ne %d comb_sum %ld\n",Ns,iMinup,iAllup,X->Def.Ne,comb_sum); + 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); + comb_sum += comb_up*comb_down; + } + } break; case Kondo: @@ -157,6 +190,34 @@ 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.Ne-NLocSpn; + int NsCond = 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(NsCond,tmp_Nup,comb,NsCond); + comb_3 = Binomial(NsCond,NCond-tmp_Nup,comb,NsCond); + 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 NsCond %d NCond=%d comb_1=%ld comb_2=%ld comb_3=%ld comb_sum=%ld\n",X->Def.Ne,NLocSpn,NsCond,NCond,comb_1,comb_2,comb_3,comb_sum); + break; + case KondoGC: comb_sum = 1; NCond = X->Def.Nsite-X->Def.NLocSpn; @@ -167,7 +228,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"); @@ -213,7 +273,11 @@ int check(struct BindStruct *X){ case Hubbard: case HubbardNConserved: case Kondo: + case KondoNConserved: case KondoGC: + case tJ: + case tJNConserved: + case tJGC: case Spin: X->Check.max_mem = 5.5 * X->Check.idim_max * 8.0 / (pow(10, 9)); break; @@ -228,7 +292,11 @@ int check(struct BindStruct *X){ case Hubbard: case HubbardNConserved: case Kondo: + case KondoNConserved: case KondoGC: + case tJ: + case tJNConserved: + case tJGC: case Spin: X->Check.max_mem = (6 * X->Def.k_exct + 2) * X->Check.idim_max * 16.0 / (pow(10, 9)); break; @@ -244,7 +312,11 @@ int check(struct BindStruct *X){ case Hubbard: case HubbardNConserved: case Kondo: + case KondoNConserved: case KondoGC: + case tJ: + case tJNConserved: + case tJGC: case Spin: if (X->Def.iFlgCalcSpec != CALCSPEC_NOT) { X->Check.max_mem = (2) * X->Check.idim_max * 16.0 / (pow(10, 9)); @@ -299,10 +371,14 @@ int check(struct BindStruct *X){ switch(X->Def.iCalcModel){ case HubbardGC: - case KondoGC: case HubbardNConserved: case Hubbard: case Kondo: + case KondoNConserved: + case KondoGC: + case tJ: + case tJNConserved: + case tJGC: while(tmp <= X->Def.Nsite){ tmp_sdim=tmp_sdim*2; tmp+=1; @@ -338,6 +414,10 @@ int check(struct BindStruct *X){ case HubbardNConserved: case Hubbard: case Kondo: + case KondoNConserved: + case tJ: + case tJNConserved: + case tJGC: //fprintf(stdoutMPI, "sdim=%ld =2^%d\n",X->Check.sdim,X->Def.Nsite); fprintf(fp,"sdim=%ld =2^%d\n",X->Check.sdim,X->Def.Nsite); break; @@ -359,15 +439,19 @@ int check(struct BindStruct *X){ switch(X->Def.iCalcModel){ case HubbardGC: case KondoGC: + case tJGC: for(i=1;i<=2*X->Def.Nsite;i++){ u_tmp=u_tmp*2; X->Def.Tpow[i]=u_tmp; fprintf(fp,"%ld %ld \n",i,u_tmp); } break; - case HubbardNConserved: case Hubbard: + case HubbardNConserved: case Kondo: + case KondoNConserved: + case tJ: + case tJNConserved: for(i=1;i<=2*X->Def.Nsite-1;i++){ u_tmp=u_tmp*2; X->Def.Tpow[i]=u_tmp; diff --git a/src/diagonalcalc.c b/src/diagonalcalc.c index 6518cc19..31612fa6 100644 --- a/src/diagonalcalc.c +++ b/src/diagonalcalc.c @@ -278,9 +278,11 @@ int SetDiagonalCoulombIntra switch (X->Def.iCalcModel) { case HubbardGC: - case KondoGC: case Hubbard: case Kondo: + case KondoGC: + case tJ: + case tJGC: is1_up = X->Def.Tpow[2 * isite1 - 2]; is1_down = X->Def.Tpow[2 * isite1 - 1]; @@ -292,7 +294,7 @@ int SetDiagonalCoulombIntra for (j = 1; j <= i_max; j++) list_Diagonal[j] += dtmp_V; } - break; /*case HubbardGC, KondoGC, Hubbard, Kondo:*/ + break; /*case HubbardGC, KondoGC, Hubbard, Kondo*/ case Spin: case SpinGC: @@ -326,9 +328,11 @@ int SetDiagonalCoulombIntra } break; - case KondoGC: case Hubbard: case Kondo: + case KondoGC: + case tJ: + case tJGC: is1_up = X->Def.Tpow[2*isite1-2]; is1_down = X->Def.Tpow[2*isite1-1]; is=is1_up+is1_down; @@ -393,10 +397,12 @@ int SetDiagonalChemi switch (X->Def.iCalcModel) { - case HubbardGC: - case KondoGC: case Hubbard: + case HubbardGC: case Kondo: + case KondoGC: + case tJ: + case tJGC: if (spin == 0) { is1 = X->Def.Tpow[2 * isite1 - 2]; @@ -410,7 +416,7 @@ int SetDiagonalChemi firstprivate(i_max, dtmp_V, num1) private(j) for (j = 1; j <= i_max; j++) list_Diagonal[j] += num1*dtmp_V; - break;/*case HubbardGC, case KondoGC, Hubbard, Kondo:*/ + break;/*case HubbardGC, case KondoGC, Hubbard, Kondo*/ case SpinGC: case Spin: @@ -460,9 +466,11 @@ firstprivate(i_max, dtmp_V) private(j) list_Diagonal[j]+=num1*dtmp_V; } break; - case KondoGC: case Hubbard: case Kondo: + case KondoGC: + case tJ: + case tJGC: if(spin==0){ is1 = X->Def.Tpow[2*isite1-2]; }else{ @@ -578,10 +586,12 @@ int SetDiagonalCoulombInter switch (X->Def.iCalcModel) { - case HubbardGC: - case KondoGC: case Hubbard: + case HubbardGC: case Kondo: + case KondoGC: + case tJ: + case tJGC: is1_up = X->Def.Tpow[2 * isite1 - 2]; is1_down = X->Def.Tpow[2 * isite1 - 1]; @@ -605,7 +615,7 @@ int SetDiagonalCoulombInter firstprivate(i_max, dtmp_V, num1, num2) private(j) for (j = 1; j <= i_max; j++) list_Diagonal[j] += num1*num2*dtmp_V; - break;/*case HubbardGC, KondoGC, Hubbard, Kondo:*/ + break;/*case HubbardGC, KondoGC, Hubbard, Kondo*/ case Spin: case SpinGC: @@ -628,9 +638,11 @@ int SetDiagonalCoulombInter switch(X->Def.iCalcModel){ case HubbardGC: - case KondoGC: case Hubbard: case Kondo: + case KondoGC: + case tJ: + case tJGC: is1_up = X->Def.Tpow[2 * isite1 - 2]; is1_down = X->Def.Tpow[2 * isite1 - 1]; is2_up = X->Def.Tpow[2 * isite2 - 2]; @@ -670,9 +682,11 @@ private(num1, ibit1_up, ibit1_down, j) break;/*case HubbardGC*/ - case KondoGC: case Hubbard: case Kondo: + case KondoGC: + case tJ: + case tJGC: #pragma omp parallel for default(none) shared(list_1, list_Diagonal) \ firstprivate(i_max, dtmp_V, is1_up, is1_down, num2) \ @@ -686,7 +700,7 @@ private(num1, ibit1_up, ibit1_down, j) list_Diagonal[j] += num1*num2*dtmp_V; } - break;/*case KondoGC, Hubbard, Kondo:*/ + break;/*case KondoGC, Hubbard, Kondo*/ case Spin: case SpinGC: @@ -729,9 +743,11 @@ private(num1, ibit1_up, ibit1_down, j) list_Diagonal[j]+=num1*num2*dtmp_V; } break; - case KondoGC: case Hubbard: case Kondo: + case KondoGC: + case tJ: + case tJGC: is1_up = X->Def.Tpow[2*isite1-2]; is1_down = X->Def.Tpow[2*isite1-1]; is2_up = X->Def.Tpow[2*isite2-2]; @@ -819,10 +835,12 @@ int SetDiagonalHund switch (X->Def.iCalcModel) { - case HubbardGC: - case KondoGC: case Hubbard: + case HubbardGC: case Kondo: + case KondoGC: + case tJ: + case tJGC: is1_up = X->Def.Tpow[2 * isite1 - 2]; is1_down = X->Def.Tpow[2 * isite1 - 1]; @@ -849,7 +867,7 @@ int SetDiagonalHund for (j = 1; j <= i_max; j++) list_Diagonal[j] += dtmp_V*(num1_up*num2_up + num1_down*num2_down); - break;/*case HubbardGC, KondoGC, Hubbard, Kondo:*/ + break;/*case HubbardGC, KondoGC, Hubbard, Kondo*/ case SpinGC: case Spin: @@ -908,9 +926,11 @@ private(num1_up, num1_down, ibit1_up, ibit1_down, j) } break;/*case HubbardGC:*/ - case KondoGC: case Hubbard: case Kondo: + case KondoGC: + case tJ: + case tJGC: is1_up = X->Def.Tpow[2 * isite1 - 2]; is1_down = X->Def.Tpow[2 * isite1 - 1]; @@ -939,7 +959,7 @@ private(num1_up, num1_down, ibit1_up, ibit1_down, j) list_Diagonal[j] += dtmp_V*(num1_up*num2_up + num1_down*num2_down); } - break;/*case KondoGC, Hubbard, Kondo:*/ + break;/*case KondoGC, Hubbard, Kondo*/ case SpinGC: is1_up = X->Def.Tpow[isite1 - 1]; @@ -1032,9 +1052,11 @@ firstprivate(i_max, dtmp_V, is1_up) private(j, ibit1_up) list_Diagonal[j]+=dtmp_V*(num1_up*num2_up+num1_down*num2_down); } break; - case KondoGC: case Hubbard: case Kondo: + case KondoGC: + case tJ: + case tJGC: is1_up = X->Def.Tpow[2*isite1-2]; is1_down = X->Def.Tpow[2*isite1-1]; is2_up = X->Def.Tpow[2*isite2-2]; @@ -1152,10 +1174,12 @@ int SetDiagonalInterAll switch (X->Def.iCalcModel) { - case HubbardGC: - case KondoGC: case Hubbard: + case HubbardGC: case Kondo: + case KondoGC: + case tJ: + case tJGC: is1_spin = X->Def.Tpow[2 * isite1 - 2 + isigma1]; is2_spin = X->Def.Tpow[2 * isite2 - 2 + isigma2]; @@ -1172,7 +1196,7 @@ int SetDiagonalInterAll firstprivate(i_max, dtmp_V, num2, num1) private(ibit1_spin, j) for (j = 1; j <= i_max; j++) list_Diagonal[j] += num1*num2*dtmp_V; - break;/*case HubbardGC, KondoGC, Hubbard, Kondo:*/ + break;/*case HubbardGC, KondoGC, Hubbard, Kondo*/ case SpinGC: case Spin: @@ -1235,9 +1259,11 @@ firstprivate(i_max, dtmp_V, is1_spin, num2) private(num1, ibit1_spin, j) } break;/*case HubbardGC:*/ - case KondoGC: case Hubbard: case Kondo: + case KondoGC: + case tJ: + case tJGC: is1_spin = X->Def.Tpow[2 * isite1 - 2 + isigma1]; is2_spin = X->Def.Tpow[2 * isite2 - 2 + isigma2]; @@ -1254,7 +1280,7 @@ firstprivate(i_max, dtmp_V, is1_spin, num2) private(num1, ibit1_spin, j) num1 += ibit1_spin / is1_spin; list_Diagonal[j] += num1*num2*dtmp_V; } - break;/*case KondoGC, Hubbard, Kondo:*/ + break;/*case KondoGC, Hubbard, Kondo*/ case SpinGC: @@ -1339,9 +1365,11 @@ firstprivate(i_max, dtmp_V, isite1, isigma1, X) private(j, num1) list_Diagonal[j]+=num1*num2*dtmp_V; } break; - case KondoGC: case Hubbard: case Kondo: + case KondoGC: + case tJ: + case tJGC: is1_spin = X->Def.Tpow[2*isite1-2+isigma1]; is2_spin = X->Def.Tpow[2*isite2-2+isigma2]; @@ -1476,10 +1504,12 @@ int SetDiagonalTEInterAll( switch (X->Def.iCalcModel) { + case Hubbard: case HubbardGC: + case Kondo: case KondoGC: - case Hubbard: - case Kondo: + case tJ: + case tJGC: is1_spin = X->Def.Tpow[2 * isite1 - 2 + isigma1]; is2_spin = X->Def.Tpow[2 * isite2 - 2 + isigma2]; num1 = 0; @@ -1488,7 +1518,7 @@ int SetDiagonalTEInterAll( num2 = 0; ibit2_spin = (unsigned long int)myrank&is2_spin; num2 += ibit2_spin / is2_spin; - break;/*case HubbardGC, KondoGC, Hubbard, Kondo:*/ + break;/*case HubbardGC, KondoGC, Hubbard, Kondo*/ case SpinGC: case Spin: @@ -1549,9 +1579,11 @@ firstprivate(i_max, dtmp_V) private(j) } break;/*case HubbardGC:*/ - case KondoGC: case Hubbard: case Kondo: + case KondoGC: + case tJ: + case tJGC: is1_spin = X->Def.Tpow[2 * isite1 - 2 + isigma1]; is2_spin = X->Def.Tpow[2 * isite2 - 2 + isigma2]; @@ -1570,7 +1602,7 @@ firstprivate(i_max, dtmp_V) private(j) dam_pr += dtmp_V * num1*conj(tmp_v1[j]) * tmp_v1[j]; } } - break;/*case KondoGC, Hubbard, Kondo:*/ + break;/*case KondoGC, Hubbard, Kondo*/ case SpinGC: @@ -1664,9 +1696,11 @@ firstprivate(i_max, dtmp_V, isite1, isigma1, X) private(j, num1) dam_pr += dtmp_V * num1*num2*conj(tmp_v1[j]) * tmp_v1[j]; } break; - case KondoGC: case Hubbard: case Kondo: + case KondoGC: + case tJ: + case tJGC: is1_spin = X->Def.Tpow[2*isite1-2+isigma1]; is2_spin = X->Def.Tpow[2*isite2-2+isigma2]; @@ -1787,10 +1821,12 @@ int SetDiagonalTEChemi( switch (X->Def.iCalcModel) { - case HubbardGC: - case KondoGC: case Hubbard: + case HubbardGC: case Kondo: + case KondoGC: + case tJ: + case tJGC: if (spin == 0) { is1 = X->Def.Tpow[2 * isite1 - 2]; @@ -1800,7 +1836,7 @@ int SetDiagonalTEChemi( } ibit1 = (unsigned long int)myrank & is1; num1 = ibit1 / is1; - break;/*case HubbardGC, case KondoGC, Hubbard, Kondo:*/ + break;/*case HubbardGC, case KondoGC, Hubbard, Kondo*/ case SpinGC: case Spin: @@ -1850,9 +1886,11 @@ firstprivate(i_max, dtmp_V) private(j) dam_pr += dtmp_V * num1*conj(tmp_v1[j]) * tmp_v1[j]; } break; - case KondoGC: case Hubbard: case Kondo: + case KondoGC: + case tJ: + case tJGC: if(spin==0){ is1 = X->Def.Tpow[2*isite1-2]; }else{ @@ -1966,10 +2004,12 @@ int SetDiagonalTETransfer switch (X->Def.iCalcModel) { - case HubbardGC: - case KondoGC: case Hubbard: + case HubbardGC: case Kondo: + case KondoGC: + case tJ: + case tJGC: if (spin == 0) { is1 = X->Def.Tpow[2 * isite1 - 2]; } @@ -1978,7 +2018,7 @@ int SetDiagonalTETransfer } ibit1 = (unsigned long int)myrank & is1; num1 = ibit1 / is1; - break;/*case HubbardGC, case KondoGC, Hubbard, Kondo:*/ + break;/*case HubbardGC, case KondoGC, Hubbard, Kondo*/ case SpinGC: case Spin: @@ -2025,9 +2065,11 @@ int SetDiagonalTETransfer } break; - case KondoGC: case Hubbard: case Kondo: + case KondoGC: + case tJ: + case tJGC: if (spin == 0) { is1 = X->Def.Tpow[2 * isite1 - 2]; } else { diff --git a/src/expec_cisajs.c b/src/expec_cisajs.c index c2b2f789..22b089bd 100644 --- a/src/expec_cisajs.c +++ b/src/expec_cisajs.c @@ -133,9 +133,11 @@ int expec_cisajs(struct BindStruct *X,double complex *vec){ } break; - case KondoGC: case Hubbard: case Kondo: + case KondoGC: + case tJ: + case tJGC: if(expec_cisajs_Hubbard(X, vec, &fp)!=0){ return -1; } diff --git a/src/expec_cisajscktaltdc.c b/src/expec_cisajscktaltdc.c index 2b760349..2aec9a7b 100644 --- a/src/expec_cisajscktaltdc.c +++ b/src/expec_cisajscktaltdc.c @@ -180,8 +180,6 @@ int expec_cisajscktaltdc fp_4 = fp; } - - switch(X->Def.iCalcModel){ case HubbardGC: if(expec_cisajscktalt_HubbardGC(X, vec, &fp)!=0){ @@ -189,9 +187,11 @@ int expec_cisajscktaltdc } break; - case KondoGC: case Hubbard: case Kondo: + case KondoGC: + case tJ: + case tJGC: if(expec_cisajscktalt_Hubbard(X, vec, &fp)!=0){ return -1; } diff --git a/src/expec_energy_flct.c b/src/expec_energy_flct.c index d022d250..eec53003 100644 --- a/src/expec_energy_flct.c +++ b/src/expec_energy_flct.c @@ -86,9 +86,11 @@ int expec_energy_flct(struct BindStruct *X){ case HubbardGC: expec_energy_flct_HubbardGC(X); break; - case KondoGC: case Hubbard: case Kondo: + case KondoGC: + case tJ: + case tJGC: expec_energy_flct_Hubbard(X); break; diff --git a/src/expec_totalspin.c b/src/expec_totalspin.c index d7dd5c3f..aca66e20 100644 --- a/src/expec_totalspin.c +++ b/src/expec_totalspin.c @@ -61,12 +61,14 @@ int expec_totalspin case SpinGC: totalspin_SpinGC(X,vec); break; - case Hubbard: - case Kondo: + case Hubbard: + case Kondo: + case KondoGC: + case tJ: + case tJGC: totalspin_Hubbard(X,vec); break; case HubbardGC: - case KondoGC: totalspin_HubbardGC(X,vec); break; default: diff --git a/src/include/DefCommon.h b/src/include/DefCommon.h index 41f3f6f0..af4c30e2 100644 --- a/src/include/DefCommon.h +++ b/src/include/DefCommon.h @@ -27,16 +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 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/readdef.h b/src/include/readdef.h index 5c01491d..bdc628f7 100644 --- a/src/include/readdef.h +++ b/src/include/readdef.h @@ -91,6 +91,17 @@ int CheckInterAllHermite const int iCalcModel ); +int CheckInterAllHermite_simple + ( + int **InterAll, + double complex* ParaInterAll, + int **InterAllOffDiagonal, + double complex*ParaInterAllOffDiagonal, + const int NInterAllOffDiagonal, + const int iCalcModel + ); + + /* int GetDiagonalInterAll ( @@ -114,6 +125,30 @@ int GetDiagonalInterAll const int iCalcModel ); +int GetDiagonalInterAll_simple + ( + int **InterAll, + complex double *ParaInterAll, + const int NInterAll, + int **InterAllDiagonal, + double *ParaInterAllDiagonal, + int **InterAllOffDiagonal, + complex double *ParaInterAllOffDiagonal, + int *Chemi, + int *SpinChemi, + double *ParaChemi, + unsigned int *NChemi, + const int iCalcModel + ); + +int ArrangeInterAllOffDiagonal +( + const int NInterAllOffDiagonal, + int **InterAllOffDiagonal, + complex double *ParaInterAllOffDiagonal, + const int iCalcModel +); + int JudgeDefType ( const int argc, diff --git a/src/include/sz.h b/src/include/sz.h index 22382f0f..51c21d33 100644 --- a/src/include/sz.h +++ b/src/include/sz.h @@ -17,100 +17,147 @@ #include "Common.h" int omp_sz( - 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_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, - 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_KondoNConserved( + 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_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( + 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_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_ - ); - -int Read_sz -( - struct BindStruct *X, - const long unsigned int irght, - const long unsigned int ilft, - const long unsigned int ihfbit, - long unsigned int *i_max - ); + 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, + const long unsigned int irght, + const long unsigned int ilft, + const long unsigned int ihfbit, + long unsigned int *i_max +); + +/*[s] func. for refactoring */ +int count_localized_spins(struct BindStruct *X); +void calculate_jb_GeneralSpin(struct BindStruct *X, long unsigned int *list_jb, long int *list_2_1_Sz,long int *list_2_2_Sz, long unsigned int ihfbit,long unsigned int ilftdim,unsigned int N); +void calculate_jb_Spin_m1(struct BindStruct *X, long unsigned int *list_jb, long unsigned int *list_1_, long unsigned int *list_2_1,long unsigned int *list_2_2,\ +long unsigned int ihfbit,long unsigned int irght,long unsigned int ilft,long unsigned int ibpatn, unsigned int N); +void calculate_jb_Spin_Old(struct BindStruct *X, long unsigned int *list_jb, long unsigned int ihfbit,unsigned int N); +void calculate_jb_Spin_Hacker(struct BindStruct *X, long unsigned int *list_jb, long unsigned int ihfbit,unsigned int N); +void calculate_jb_Kondo(struct BindStruct *X, long unsigned int *list_jb, long unsigned int ihfbit); +void calculate_jb_KondoGC(struct BindStruct *X, int num_loc, long unsigned int *list_jb, long unsigned int ihfbit); +void calculate_jb_KondoNConserved(struct BindStruct *X, long unsigned int *list_jb, long unsigned int ihfbit); +void calculate_jb_Hubbard(struct BindStruct *X,long unsigned int *list_jb, long unsigned int ihfbit, unsigned int N2); +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/makeHam.c b/src/makeHam.c index d492bce1..24fad5c9 100644 --- a/src/makeHam.c +++ b/src/makeHam.c @@ -210,9 +210,11 @@ int makeHam(struct BindStruct *X) { } } break; - case KondoGC: case Hubbard: case Kondo: + case KondoGC: + case tJ: + case tJGC: //Transfer for (i = 0; i < X->Def.EDNTransfer / 2; i++) { for (ihermite = 0; ihermite < 2; ihermite++) { diff --git a/src/mltply.c b/src/mltply.c index bd72ef1d..2f851a80 100644 --- a/src/mltply.c +++ b/src/mltply.c @@ -108,9 +108,11 @@ int mltply(struct BindStruct *X, double complex *tmp_v0,double complex *tmp_v1) mltplyHubbardGC(X, tmp_v0, tmp_v1); break; - case KondoGC: case Hubbard: case Kondo: + case KondoGC: + case tJ: + case tJGC: mltplyHubbard(X, tmp_v0, tmp_v1); break; diff --git a/src/output.c b/src/output.c index 64b027fe..4b2e4923 100644 --- a/src/output.c +++ b/src/output.c @@ -36,11 +36,13 @@ int output(struct BindStruct *X) { case Spin: case Hubbard: case Kondo: + case tJ: sprintf(sdt, cFileNamePhys_FullDiag, X->Def.CDataFileHead, X->Def.Nup, X->Def.Ndown); break; case SpinGC: case HubbardGC: case KondoGC: + case tJGC: sprintf(sdt, cFileNamePhys_FullDiag_GC, X->Def.CDataFileHead); break; default: diff --git a/src/output_list.c b/src/output_list.c index da3f8550..504bae85 100644 --- a/src/output_list.c +++ b/src/output_list.c @@ -48,6 +48,8 @@ int output_list(struct BindStruct *X){ switch(X->Def.iCalcModel){ case HubbardGC: case Hubbard: + case tJGC: + case tJ: case Spin: case SpinGC: sprintf(sdt, cFileNameListModel, X->Def.Nsite,X->Def.Nup,X->Def.Ndown); diff --git a/src/readdef.c b/src/readdef.c index 80919d19..4b47c71b 100644 --- a/src/readdef.c +++ b/src/readdef.c @@ -848,6 +848,7 @@ int ReadDefFileNInt( switch(X->iCalcModel){ case Spin: case Hubbard: + case tJ: case Kondo: case SpinlessFermion: @@ -869,17 +870,29 @@ int ReadDefFileNInt( X->Ndown=X->NLocSpn+X->NCond-X->Total2Sz; X->Nup/=2; X->Ndown/=2; - } - else{ - if(X->iCalcModel == Hubbard){ + }else{ + if(X->iCalcModel == Hubbard ){ X->Ne=X->NCond; if(X->Ne <1){ fprintf(stdoutMPI, "Ncond is incorrect.\n"); return(-1); } X->iCalcModel=HubbardNConserved; - } - else if(X->iCalcModel ==SpinlessFermion){ + }else if(X->iCalcModel == tJ ){ + X->Ne=X->NCond; + if(X->Ne <1){ + fprintf(stdoutMPI, "Ncond is incorrect.\n"); + return(-1); + } + X->iCalcModel=tJNConserved; + }else if(X->iCalcModel == Kondo){ + X->Ne=X->NCond + X->NLocSpn; + if(X->Ne <1){ + fprintf(stdoutMPI, "Ncond is incorrect.\n"); + return(-1); + } + X->iCalcModel=KondoNConserved; + }else if(X->iCalcModel ==SpinlessFermion){ X->Ne=X->NCond; X->Nup=X->NCond; X->Ndown=0; @@ -931,6 +944,7 @@ int ReadDefFileNInt( case SpinGC: case KondoGC: case HubbardGC: + case tJGC: case SpinlessFermionGC: if(iReadNCond == TRUE || X->iFlgSzConserved ==TRUE){ fprintf(stdoutMPI, "\n Warning: For GC, both Ncond and 2Sz should not be defined.\n"); @@ -1416,26 +1430,60 @@ int ReadDefFileIdxPara( X->InterAll_OffDiagonal, X->ParaInterAll_OffDiagonal, X->InterAll_Diagonal, X->ParaInterAll_Diagonal, NInterAllSet); */ - if(GetDiagonalInterAll( - X->InterAll, X->ParaInterAll, X->NInterAll, - X->InterAll_Diagonal, X->ParaInterAll_Diagonal, - X->InterAll_OffDiagonal, X->ParaInterAll_OffDiagonal, - X->EDChemi, X->EDSpinChemi, X->EDParaChemi, &X->EDNChemi, - X->iCalcModel - )!=0){ - fclose(fp); - return(-1); - } - if(CheckInterAllHermite( - X->InterAll, X->ParaInterAll, - X->InterAll_OffDiagonal, X->ParaInterAll_OffDiagonal, - X->NInterAll_OffDiagonal, X->iCalcModel - )!=0) { - fprintf(stdoutMPI, "%s", cErrNonHermiteInterAllForAll); - fclose(fp); - return (-1); - } +// if(GetDiagonalInterAll( +// X->InterAll, X->ParaInterAll, X->NInterAll, +// X->InterAll_Diagonal, X->ParaInterAll_Diagonal, +// X->InterAll_OffDiagonal, X->ParaInterAll_OffDiagonal, +// X->EDChemi, X->EDSpinChemi, X->EDParaChemi, &X->EDNChemi, +// X->iCalcModel +// )!=0){ +// fclose(fp); +// return(-1); +// } +// +// if(CheckInterAllHermite( +// X->InterAll, X->ParaInterAll, +// X->InterAll_OffDiagonal, X->ParaInterAll_OffDiagonal, +// X->NInterAll_OffDiagonal, X->iCalcModel +// )!=0) { +// fprintf(stdoutMPI, "%s", cErrNonHermiteInterAllForAll); +// fclose(fp); +// return (-1); +// } + + if(GetDiagonalInterAll_simple( + X->InterAll, X->ParaInterAll, X->NInterAll, + X->InterAll_Diagonal, X->ParaInterAll_Diagonal, + X->InterAll_OffDiagonal, X->ParaInterAll_OffDiagonal, + X->EDChemi, X->EDSpinChemi, X->EDParaChemi, &X->EDNChemi, + X->iCalcModel + )!=0){ + fclose(fp); + return(-1); + } + + if(CheckInterAllHermite_simple( + X->InterAll, X->ParaInterAll, + X->InterAll_OffDiagonal, X->ParaInterAll_OffDiagonal, + X->NInterAll_OffDiagonal, X->iCalcModel + )!=0) { + fprintf(stdoutMPI, "%s", cErrNonHermiteInterAllForAll); + fclose(fp); + return (-1); + } + + if(ArrangeInterAllOffDiagonal( + X->NInterAll_OffDiagonal, + X->InterAll_OffDiagonal, X->ParaInterAll_OffDiagonal, + X->iCalcModel + )!=0){ + fclose(fp); + return(-1); + } + + + break; case KWOneBodyG: @@ -1827,24 +1875,53 @@ int ReadDefFileIdxPara( X->NTEInterAll[idx] = icnt_interall; X->NTEInterAllDiagonal[idx] = icnt_diagonal; X->NTEInterAllOffDiagonal[idx] = icnt_interall - icnt_diagonal; + //Diagonal -> OffDiagonal -> search pair -> hermite - if (GetDiagonalInterAll(X->TEInterAll[idx], X->ParaTEInterAll[idx], X->NTEInterAll[idx], X->TEInterAllDiagonal[idx], X->ParaTEInterAllDiagonal[idx], - X->TEInterAllOffDiagonal[idx], X->ParaTEInterAllOffDiagonal[idx], X->TEChemi[idx], X->SpinTEChemi[idx], X->ParaTEChemi[idx], &X->NTEChemi[idx], X->iCalcModel) != 0) - { - fclose(fp); - return (-1); - } +// if (GetDiagonalInterAll(X->TEInterAll[idx], X->ParaTEInterAll[idx], X->NTEInterAll[idx], X->TEInterAllDiagonal[idx], X->ParaTEInterAllDiagonal[idx], +// X->TEInterAllOffDiagonal[idx], X->ParaTEInterAllOffDiagonal[idx], X->TEChemi[idx], X->SpinTEChemi[idx], X->ParaTEChemi[idx], &X->NTEChemi[idx], X->iCalcModel) != 0) +// { +// fclose(fp); +// return (-1); +// } +// +// if(CheckInterAllHermite( +// X->TEInterAll[idx], X->ParaTEInterAll[idx], +// X->TEInterAllOffDiagonal[idx], X->ParaTEInterAllOffDiagonal[idx], +// X->NTEInterAllOffDiagonal[idx], X->iCalcModel +// )!=0) { +// fprintf(stdoutMPI, "%s", cErrNonHermiteInterAllForAll); +// fclose(fp); +// return (-1); +// } +// idx++; +// } + + if (GetDiagonalInterAll_simple(X->TEInterAll[idx], X->ParaTEInterAll[idx], X->NTEInterAll[idx], X->TEInterAllDiagonal[idx], X->ParaTEInterAllDiagonal[idx], + X->TEInterAllOffDiagonal[idx], X->ParaTEInterAllOffDiagonal[idx], X->TEChemi[idx], X->SpinTEChemi[idx], X->ParaTEChemi[idx], &X->NTEChemi[idx], X->iCalcModel) != 0) + { + fclose(fp); + return (-1); + } - if(CheckInterAllHermite( - X->TEInterAll[idx], X->ParaTEInterAll[idx], - X->TEInterAllOffDiagonal[idx], X->ParaTEInterAllOffDiagonal[idx], - X->NTEInterAllOffDiagonal[idx], X->iCalcModel - )!=0) { - fprintf(stdoutMPI, "%s", cErrNonHermiteInterAllForAll); - fclose(fp); - return (-1); - } - idx++; + if(CheckInterAllHermite_simple( + X->TEInterAll[idx], X->ParaTEInterAll[idx], + X->TEInterAllOffDiagonal[idx], X->ParaTEInterAllOffDiagonal[idx], + X->NTEInterAllOffDiagonal[idx], X->iCalcModel + )!=0) { + fprintf(stdoutMPI, "%s", cErrNonHermiteInterAllForAll); + fclose(fp); + return (-1); + } + + if(ArrangeInterAllOffDiagonal( + X->NTEInterAllOffDiagonal[idx], + X->TEInterAllOffDiagonal[idx], X->ParaTEInterAllOffDiagonal[idx], + X->iCalcModel + )!=0){ + fclose(fp); + return(-1); + } + idx++; } if(idx!=X->NTETimeSteps){ @@ -2160,8 +2237,6 @@ int CheckTransferHermite isite2=X->GeneralTransfer[i][2]; isigma2=X->GeneralTransfer[i][3]; icheckHermiteCount=FALSE; - // fprintf(stdoutMPI, "Debug: isite1=%d, isigma1=%d, isite2=%d, isigma2=%d, reTrans=%lf, imTrans = %lf\n", - // isite1, isigma1, isite2, isigma2, creal(X->ParaGeneralTransfer[i]), cimag((X->ParaGeneralTransfer[i]))); for(j=0; jNTransfer; j++){ itmpsite1=X->GeneralTransfer[j][0]; itmpsigma1=X->GeneralTransfer[j][1]; @@ -2273,6 +2348,8 @@ int CheckInterAllHermite double complex ddiff_intall; icntincorrect = 0; icntHermite = 0; + + for (i = 0; i < NInterAllOffDiagonal; i++) { itmpret = 0; isite1 = InterAllOffDiagonal[i][0]; @@ -2284,7 +2361,6 @@ int CheckInterAllHermite isite4 = InterAllOffDiagonal[i][6]; isigma4 = InterAllOffDiagonal[i][7]; icheckHermiteCount = FALSE; - for (j = 0; j < NInterAllOffDiagonal; j++) { itmpsite1 = InterAllOffDiagonal[j][0]; itmpsigma1 = InterAllOffDiagonal[j][1]; @@ -2395,6 +2471,96 @@ int CheckInterAllHermite return 0; } +int CheckInterAllHermite_simple + ( + int **InterAll, + double complex* ParaInterAll, + int **InterAllOffDiagonal, + double complex*ParaInterAllOffDiagonal, + const int NInterAllOffDiagonal, + const int iCalcModel + ) { + unsigned int i, j, icntincorrect, itmpret; + int isite1, isite2, isite3, isite4; + int isigma1, isigma2, isigma3, isigma4; + int itmpsite1, itmpsite2, itmpsite3, itmpsite4; + int itmpsigma1, itmpsigma2, itmpsigma3, itmpsigma4; + unsigned int itmpIdx, icntHermite; + int icheckHermiteCount = FALSE; + double complex ddiff_intall; + icntincorrect = 0; + icntHermite = 0; + + for (i = 0; i < NInterAllOffDiagonal; i++) { + itmpret = 0; + isite1 = InterAllOffDiagonal[i][0]; + isigma1 = InterAllOffDiagonal[i][1]; + isite2 = InterAllOffDiagonal[i][2]; + isigma2 = InterAllOffDiagonal[i][3]; + isite3 = InterAllOffDiagonal[i][4]; + isigma3 = InterAllOffDiagonal[i][5]; + isite4 = InterAllOffDiagonal[i][6]; + isigma4 = InterAllOffDiagonal[i][7]; + icheckHermiteCount = FALSE; + for (j = 0; j < NInterAllOffDiagonal; j++) { + itmpsite1 = InterAllOffDiagonal[j][0]; + itmpsigma1 = InterAllOffDiagonal[j][1]; + itmpsite2 = InterAllOffDiagonal[j][2]; + itmpsigma2 = InterAllOffDiagonal[j][3]; + itmpsite3 = InterAllOffDiagonal[j][4]; + itmpsigma3 = InterAllOffDiagonal[j][5]; + itmpsite4 = InterAllOffDiagonal[j][6]; + itmpsigma4 = InterAllOffDiagonal[j][7]; + if (isite1 == itmpsite4 && isite2 == itmpsite3 && isite3 == itmpsite2 && isite4 == itmpsite1) { + if (isigma1 == itmpsigma4 && isigma2 == itmpsigma3 && isigma3 == itmpsigma2 && isigma4 == itmpsigma1) { + ddiff_intall = cabs(ParaInterAllOffDiagonal[i] - conj(ParaInterAllOffDiagonal[j])); + if (cabs(ddiff_intall) < eps_CheckImag0) { + itmpret = 1; + if (icheckHermiteCount == FALSE) { + if (i <= j) { + icheckHermiteCount = TRUE; //for not double counting + if (2 * icntHermite >= NInterAllOffDiagonal) { + fprintf(stdoutMPI, "Elements of InterAll are incorrect.\n"); + return (-1); + } + for (itmpIdx = 0; itmpIdx < 8; itmpIdx++) { + InterAll[2 * icntHermite][itmpIdx] = InterAllOffDiagonal[i][itmpIdx]; + InterAll[2 * icntHermite + 1][itmpIdx] = InterAllOffDiagonal[j][itmpIdx]; + } + ParaInterAll[2 * icntHermite] = ParaInterAllOffDiagonal[i]; + ParaInterAll[2 * icntHermite + 1] = ParaInterAllOffDiagonal[j]; + icntHermite++; + break; + } + } + } + } + } + } + //if counterpart for satisfying hermite conjugate does not exist. + if (itmpret != 1) { + fprintf(stdoutMPI, cErrNonHermiteInterAll, isite1, isigma1, isite2, isigma2, isite3, isigma3, isite4, isigma4, + creal(ParaInterAllOffDiagonal[i]), cimag(ParaInterAllOffDiagonal[i])); + icntincorrect++; + } + } + + //if (icntincorrect != 0) { + if (icntincorrect != 0 || NInterAllOffDiagonal != 2*icntHermite) { + fprintf(stdoutMPI,"DEBUG: icntincorrect=%d, NInterAllOffDiagonal=%d, icntHermite=%d\n",icntincorrect,NInterAllOffDiagonal,icntHermite); + return (-1); + } + + for (i = 0; i < NInterAllOffDiagonal; i++) { + for (itmpIdx = 0; itmpIdx < 8; itmpIdx++) { + InterAllOffDiagonal[i][itmpIdx] = InterAll[i][itmpIdx]; + } + ParaInterAllOffDiagonal[i] = ParaInterAll[i]; + } + + return 0; +} + /// \brief function of getting diagonal components /// \param InterAll arrays of information of interall interactions /// \param ParaInterAll arrays of values of interall interactions @@ -2478,9 +2644,13 @@ int GetDiagonalInterAll switch(iCalcModel){ case Hubbard: case HubbardNConserved: + case HubbardGC: case Kondo: + case KondoNConserved: case KondoGC: - case HubbardGC: + case tJ: + case tJNConserved: + case tJGC: if(isigma1 == isigma2 && isigma3 == isigma4){ for(tmp_i=0; tmp_i<8; tmp_i++){ InterAllOffDiagonal[icnt_offdiagonal][tmp_i]=InterAll[i][tmp_i]; @@ -2548,6 +2718,152 @@ int GetDiagonalInterAll return 0; } +int GetDiagonalInterAll_simple + ( + int **InterAll, + complex double *ParaInterAll, + const int NInterAll, + int **InterAllDiagonal, + double *ParaInterAllDiagonal, + int **InterAllOffDiagonal, + complex double *ParaInterAllOffDiagonal, + int *Chemi, + int *SpinChemi, + double *ParaChemi, + unsigned int *NChemi, + const int iCalcModel + ) +{ + unsigned int i,icnt_diagonal, icnt_offdiagonal, tmp_i; + int isite1, isite2, isite3, isite4; + int isigma1, isigma2, isigma3, isigma4; + icnt_diagonal=0; + icnt_offdiagonal=0; + + for(i=0; iNsite; i++){ @@ -2739,6 +3058,7 @@ int CheckLocSpin break; case Kondo: + case KondoNConserved: case KondoGC: for(i=0; iNsite; i++){ if(X->LocSpn[i]>LOCSPIN){ @@ -3020,7 +3340,8 @@ int CheckInterAllCondition( return -1; } } - else if(iCalcModel == Kondo){ + //else if(iCalcModel == Kondo){ + else if(iCalcModel == Kondo || iCalcModel == KondoNConserved){ if(CheckFormatForKondoInt(isite1, isite2, isite3, isite4, iLocInfo)!=0){ return -1; } diff --git a/src/sz.c b/src/sz.c index 028f5491..29e8899c 100644 --- a/src/sz.c +++ b/src/sz.c @@ -36,15 +36,6 @@ * */ -int omp_sz_Kondo_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_ -); /** * @@ -59,776 +50,430 @@ int omp_sz_Kondo_hacker( * @author Takahiro Misawa (The University of Tokyo) * @author Kazuyoshi Yoshimi (The University of Tokyo) */ -int sz -( - struct BindStruct *X, - long unsigned int *list_1_, - long unsigned int *list_2_1_, - long unsigned int *list_2_2_ - ) -{ - FILE *fp,*fp_err; - char sdt[D_FileNameMax],sdt_err[D_FileNameMax]; - long unsigned int *HilbertNumToSz; - long unsigned int i,icnt; - long unsigned int ib,jb,ib_start,ib_end, sdim_div, sdim_rest; +int sz( + struct BindStruct *X, + long unsigned int *list_1_, + long unsigned int *list_2_1_, + long unsigned int *list_2_2_ +){ + FILE *fp,*fp_err; + char sdt[D_FileNameMax],sdt_err[D_FileNameMax]; + long unsigned int *HilbertNumToSz; + long unsigned int i,icnt; + long unsigned int ib,jb,ib_start,ib_end, sdim_div, sdim_rest; - long unsigned int j; - long unsigned int div; - long unsigned int num_up,num_down; - long unsigned int irght,ilft,ihfbit; - long unsigned int *jbthread; - - //*[s] for omp parall - int mythread; - unsigned int all_up,all_down,tmp_res,num_threads; - long unsigned int tmp_1,tmp_2,tmp_3; - long int **comb, **comb2; - //*[e] for omp parall - - // [s] for Kondo - unsigned int N_all_up, N_all_down; - unsigned int all_loc; - long unsigned int num_loc, div_down; - unsigned int num_loc_up; - int icheck_loc; - int ihfSpinDown=0; - // [e] for Kondo - - long unsigned int i_max=0; - double idim=0.0; - long unsigned int div_up; - - // [s] for general spin - long int *list_2_1_Sz = NULL; - long int *list_2_2_Sz = NULL; - if(X->Def.iFlgGeneralSpin==TRUE){ - list_2_1_Sz = li_1d_allocate(X->Check.sdim+2); - list_2_2_Sz = li_1d_allocate((X->Def.Tpow[X->Def.Nsite-1]*X->Def.SiteToBit[X->Def.Nsite-1]/X->Check.sdim)+2); - for(j=0; jCheck.sdim+2;j++){ - list_2_1_Sz[j]=0; - } - for(j=0; j< (X->Def.Tpow[X->Def.Nsite-1]*X->Def.SiteToBit[X->Def.Nsite-1]/X->Check.sdim)+2; j++){ - list_2_2_Sz[j]=0; - } - } - // [e] for general spin - - long unsigned int *list_jb; - list_jb = lui_1d_allocate(X->Large.SizeOflistjb); - for(i=0; iLarge.SizeOflistjb; i++){ - list_jb[i]=0; - } - -//hacker - int hacker; - long unsigned int tmp_i,tmp_j,tmp_pow,max_tmp_i; - long unsigned int ia,ja; - long unsigned int ibpatn=0; -//hacker - - int iSpnup, iMinup,iAllup; - unsigned int N2=0; - unsigned int N=0; - fprintf(stdoutMPI, "%s", cProStartCalcSz); - TimeKeeper(X, cFileNameSzTimeKeep, cInitalSz, "w"); - TimeKeeper(X, cFileNameTimeKeep, cInitalSz, "a"); - - if(X->Check.idim_max!=0){ - switch(X->Def.iCalcModel){ - case HubbardGC: - case HubbardNConserved: - case Hubbard: - N2=2*X->Def.Nsite; - idim = pow(2.0,N2); - break; - case KondoGC: - case Kondo: - N2 = 2*X->Def.Nsite; - N = X->Def.Nsite; - idim = pow(2.0,N2); - for(j=0;jDef.LocSpn[j]); - } - break; - case SpinGC: - case Spin: - N=X->Def.Nsite; - if(X->Def.iFlgGeneralSpin==FALSE){ - idim = pow(2.0, N); - } - else{ - idim=1; - for(j=0; jDef.SiteToBit[j]; - } - } - break; - } - comb = li_2d_allocate(X->Def.Nsite+1,X->Def.Nsite+1); - i_max=X->Check.idim_max; - - switch(X->Def.iCalcModel){ - case HubbardNConserved: - case Hubbard: - case KondoGC: - case Kondo: - case Spin: - if(X->Def.iFlgGeneralSpin==FALSE){ - if(GetSplitBitByModel(X->Def.Nsite, X->Def.iCalcModel, &irght, &ilft, &ihfbit)!=0){ - exitMPI(-1); - } - X->Large.irght=irght; - X->Large.ilft=ilft; - X->Large.ihfbit=ihfbit; - //fprintf(stdoutMPI, "idim=%lf irght=%ld ilft=%ld ihfbit=%ld \n",idim,irght,ilft,ihfbit); - } - else{ - ihfbit=X->Check.sdim; - //fprintf(stdoutMPI, "idim=%lf ihfbit=%ld \n",idim, ihfbit); - } - break; - default: - break; - } - - icnt=1; - jb=0; - - if(X->Def.READ==1){ - if(Read_sz(X, irght, ilft, ihfbit, &i_max)!=0){ - exitMPI(-1); - } - } - else{ - sprintf(sdt, cFileNameSzTimeKeep, X->Def.CDataFileHead); -#ifdef _OPENMP - num_threads = omp_get_max_threads(); -#else - num_threads=1; -#endif - childfopenMPI(sdt,"a", &fp); - fprintf(fp, "num_threads==%d\n",num_threads); - fclose(fp); - - //*[s] omp parallel - - TimeKeeper(X, cFileNameSzTimeKeep, cOMPSzStart, "a"); - TimeKeeper(X, cFileNameTimeKeep, cOMPSzStart, "a"); - switch(X->Def.iCalcModel){ - case HubbardGC: - icnt = X->Def.Tpow[2*X->Def.Nsite-1]*2+0;/*Tpow[2*X->Def.Nsit]=1*/ - break; + long unsigned int j; + long unsigned int div; + long unsigned int num_up,num_down; + long unsigned int irght,ilft,ihfbit; + long unsigned int *jbthread; + + /*[s] for omp parall*/ + int mythread; + unsigned int all_up,all_down,tmp_res,num_threads; + long unsigned int tmp_1,tmp_2,tmp_3; + long int **comb, **comb2; + /*[e] for omp parall*/ + + /*[s] for Kondo*/ + unsigned int N_all_up, N_all_down; + unsigned int all_loc; + long unsigned int num_loc, div_down; + unsigned int num_loc_up; + int icheck_loc; + int ihfSpinDown=0; + /*[e] for Kondo*/ - case SpinGC: - if(X->Def.iFlgGeneralSpin==FALSE){ - icnt = X->Def.Tpow[X->Def.Nsite-1]*2+0;/*Tpow[X->Def.Nsit]=1*/ - } - else{ - icnt = X->Def.Tpow[X->Def.Nsite-1]*X->Def.SiteToBit[X->Def.Nsite-1]; - } - break; - - case KondoGC: - // this part can not be parallelized - jb = 0; - num_loc=0; - for(j=X->Def.Nsite/2; j< X->Def.Nsite ;j++){ // counting # of localized spins - if(X->Def.LocSpn[j] != ITINERANT){ // //ITINERANT ==0 -> itinerant - num_loc += 1; - } - } - - for(ib=0;ibCheck.sdim;ib++){ - list_jb[ib]=jb; - i=ib*ihfbit; - icheck_loc=1; - for(j=(X->Def.Nsite+1)/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){ - 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){ - if(X->Def.Nsite%2==1 && X->Def.LocSpn[X->Def.Nsite/2] != ITINERANT){ - jb +=X->Def.Tpow[X->Def.Nsite-1-(X->Def.NLocSpn-num_loc)]; - }else{ - jb +=X->Def.Tpow[X->Def.Nsite-(X->Def.NLocSpn-num_loc)]; - } - } - } - - 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;ibCheck.sdim;ib++){ - icnt+=omp_sz_KondoGC(ib, ihfbit, X, list_1_, list_2_1_, list_2_2_, list_jb); - } - break; + long unsigned int i_max=0; + double idim=0.0; + long unsigned int div_up; - case Hubbard: - hacker = X->Def.read_hacker; - if(hacker==0){ - // this part can not be parallelized - jb = 0; - for(ib=0;ibCheck.sdim;ib++){ // sdim = 2^(N/2) - list_jb[ib] = jb; - i = ib*ihfbit; - //[s] counting # of up and down electrons - num_up = 0; - for(j=0;j<=N2-2;j+=2){ // even -> up spin - div=i & X->Def.Tpow[j]; - div=div/X->Def.Tpow[j]; - num_up+=div; - } - num_down=0; - for(j=1;j<=N2-1;j+=2){ // odd -> down spin - div=i & X->Def.Tpow[j]; - div=div/X->Def.Tpow[j]; - num_down+=div; - } - //[e] counting # of up and down electrons - 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.Ndown-num_down,comb,all_down); - jb += tmp_1*tmp_2; + /*[s] for general spin*/ + long int *list_2_1_Sz = NULL; + long int *list_2_2_Sz = NULL; + if(X->Def.iFlgGeneralSpin==TRUE){ + list_2_1_Sz = li_1d_allocate(X->Check.sdim+2); + list_2_2_Sz = li_1d_allocate((X->Def.Tpow[X->Def.Nsite-1]*X->Def.SiteToBit[X->Def.Nsite-1]/X->Check.sdim)+2); + for(j=0; jCheck.sdim+2;j++){ + list_2_1_Sz[j]=0; } - - //#pragma omp barrier - 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); + for(j=0; j< (X->Def.Tpow[X->Def.Nsite-1]*X->Def.SiteToBit[X->Def.Nsite-1]/X->Check.sdim)+2; j++){ + list_2_2_Sz[j]=0; } - break; - }else if(hacker==1){ - jbthread = lui_1d_allocate(nthreads); - #pragma omp parallel default(none) \ - shared(X,list_jb,ihfbit,N2,nthreads,jbthread) \ - private(ib,i,j,num_up,num_down,div,tmp_res,tmp_1,tmp_2,jb,all_up,all_down, \ - comb2,mythread,sdim_div,sdim_rest,ib_start,ib_end) - { - jb = 0; -#ifdef _OPENMP - mythread = omp_get_thread_num(); -#else - mythread = 0; -#endif - comb2 = li_2d_allocate(X->Def.Nsite+1,X->Def.Nsite+1); - // - // explict loop decomposition is nessesary to fix the asignment to each thread - // - sdim_div = X->Check.sdim / nthreads; - sdim_rest = X->Check.sdim % nthreads; - if(mythread < sdim_rest){ - ib_start = sdim_div*mythread + mythread; - ib_end = ib_start + sdim_div + 1; - } - else{ - ib_start = sdim_div*mythread + sdim_rest; - ib_end = ib_start + sdim_div; - } - // - for(ib=ib_start;ibDef.Tpow[j]; - div=div/X->Def.Tpow[j]; - num_up+=div; - } - num_down=0; - for(j=1;j<=N2-1;j+=2){ - div=i & X->Def.Tpow[j]; - div=div/X->Def.Tpow[j]; - num_down+=div; - } - - 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,comb2,all_up); - tmp_2 = Binomial(all_down,X->Def.Ndown-num_down,comb2,all_down); - jb += tmp_1*tmp_2; - } - free_li_2d_allocate(comb2); - if(mythread != nthreads-1) jbthread[mythread+1] = jb; - #pragma omp barrier - #pragma omp single - { - jbthread[0] = 0; - for(j=1;jLarge.SizeOflistjb); + for(i=0; iLarge.SizeOflistjb; i++){ + list_jb[i]=0; + } - icnt = 0; -#pragma omp parallel for default(none) reduction(+:icnt) private(ib) firstprivate(ihfbit, X) shared(list_1_, list_2_1_, list_2_2_, list_jb) - for(ib=0;ibCheck.sdim;ib++){ - icnt+=omp_sz_hacker(ib,ihfbit,X,list_1_, list_2_1_, list_2_2_, list_jb); + /*hacker*/ + int hacker; + long unsigned int tmp_i,tmp_j,tmp_pow,max_tmp_i; + long unsigned int ia,ja; + long unsigned int ibpatn=0; + /*hacker*/ + + int iSpnup, iMinup,iAllup; + unsigned int N2 = 0; + unsigned int N = 0; + fprintf(stdoutMPI, "%s", cProStartCalcSz); + TimeKeeper(X, cFileNameSzTimeKeep, cInitalSz, "w"); + TimeKeeper(X, cFileNameTimeKeep, cInitalSz, "a"); + if(X->Check.idim_max!=0){ + /*[s] calculating the maximum size of Hilbert dimensions*/ + switch(X->Def.iCalcModel){ + case HubbardGC: + case HubbardNConserved: + case Hubbard: + case tJGC: + case tJNConserved: + case tJ: + N2=2*X->Def.Nsite; + idim = pow(2.0,N2); + break; + case KondoGC: + case Kondo: + case KondoNConserved: + N2 = 2*X->Def.Nsite; + N = X->Def.Nsite; + idim = pow(2.0,N2); + for(j=0;jDef.LocSpn[j]); + } + break; + case SpinGC: + case Spin: + N=X->Def.Nsite; + if(X->Def.iFlgGeneralSpin==FALSE){ + idim = pow(2.0, N); + }else{ + idim=1; + for(j=0; jDef.SiteToBit[j]; + } + } + break; } - break; - } - else{ - fprintf(stderr, "Error: CalcHS in ModPara file must be 0 or 1 for Hubbard model."); - return -1; - } - - case HubbardNConserved: - hacker = X->Def.read_hacker; - if(hacker==0){ - // this part can not be parallelized - jb = 0; - iSpnup=0; - 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(ib=0;ibCheck.sdim;ib++){ - list_jb[ib]=jb; - i=ib*ihfbit; - num_up=0; - for(j=0;j<=N2-2;j+=2){ - div=i & X->Def.Tpow[j]; - div=div/X->Def.Tpow[j]; - num_up+=div; - } - num_down=0; - for(j=1;j<=N2-1;j+=2){ - div=i & X->Def.Tpow[j]; - div=div/X->Def.Tpow[j]; - num_down+=div; - } - 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; + /*[e] calculating the maximum size of Hilbert dimensions*/ - for(iSpnup=iMinup; iSpnup<= iAllup; iSpnup++){ - tmp_1 = Binomial(all_up, iSpnup-num_up,comb,all_up); - tmp_2 = Binomial(all_down, X->Def.Ne-iSpnup-num_down,comb,all_down); - jb += tmp_1*tmp_2; - } - } - //#pragma omp barrier - 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;ibCheck.sdim;ib++){ - icnt+=omp_sz(ib,ihfbit, X,list_1_, list_2_1_, list_2_2_, list_jb); - } - break; - } - else if(hacker==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; + comb = li_2d_allocate(X->Def.Nsite+1,X->Def.Nsite+1); /*comb=Binomial coefficient */ + i_max = X->Check.idim_max; /*i_max, idim_max : Hilbert dimensions*/ + + /*[s] calculating irght, ilht, ihfbit*/ + switch(X->Def.iCalcModel){ + case HubbardNConserved: + case Hubbard: + case tJGC: + case tJNConserved: + case tJ: + case KondoGC: + case Kondo: + case KondoNConserved: + case Spin: + if(X->Def.iFlgGeneralSpin==FALSE){ + if(GetSplitBitByModel(X->Def.Nsite, X->Def.iCalcModel, &irght, &ilft, &ihfbit)!=0){ + exitMPI(-1); + } + X->Large.irght=irght; + X->Large.ilft=ilft; + X->Large.ihfbit=ihfbit; + //fprintf(stdoutMPI, "idim=%lf irght=%ld ilft=%ld ihfbit=%ld \n",idim,irght,ilft,ihfbit); + }else{ + ihfbit=X->Check.sdim; + //fprintf(stdoutMPI, "idim=%lf ihfbit=%ld \n",idim, ihfbit); + } + break; + default: + break; } - jbthread = lui_1d_allocate(nthreads); - #pragma omp parallel default(none) \ - shared(X,iMinup,iAllup,list_jb,ihfbit,N2,nthreads,jbthread) \ - private(ib,i,j,num_up,num_down,div,tmp_res,tmp_1,tmp_2,iSpnup,jb,all_up,all_down,comb2, \ - mythread,sdim_rest,sdim_div,ib_start,ib_end) - { - jb = 0; - iSpnup=0; -#ifdef _OPENMP - mythread = omp_get_thread_num(); -#else - mythread = 0; -#endif - comb2 = li_2d_allocate(X->Def.Nsite+1,X->Def.Nsite+1); - // - // explict loop decomposition is nessesary to fix the asignment to each thread - // - sdim_div = X->Check.sdim / nthreads; - sdim_rest = X->Check.sdim % nthreads; - if(mythread < sdim_rest){ - ib_start = sdim_div*mythread + mythread; - ib_end = ib_start + sdim_div + 1; - } - else{ - ib_start = sdim_div*mythread + sdim_rest; - ib_end = ib_start + sdim_div; - } - // - for(ib=ib_start;ibDef.Tpow[j]; - div=div/X->Def.Tpow[j]; - num_up+=div; - } - num_down=0; - for(j=1;j<=N2-1;j+=2){ - div=i & X->Def.Tpow[j]; - div=div/X->Def.Tpow[j]; - num_down+=div; - } - 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; + /*[e] calculating irght, ilht, ihfbit*/ + + icnt=1; + jb=0; - for(iSpnup=iMinup; iSpnup<= iAllup; iSpnup++){ - tmp_1 = Binomial(all_up, iSpnup-num_up,comb2,all_up); - tmp_2 = Binomial(all_down, X->Def.Ne-iSpnup-num_down,comb2,all_down); - jb += tmp_1*tmp_2; + if(X->Def.READ==1){ + if(Read_sz(X, irght, ilft, ihfbit, &i_max)!=0){ + exitMPI(-1); } + }else{ + sprintf(sdt, cFileNameSzTimeKeep, X->Def.CDataFileHead); + #ifdef _OPENMP + num_threads = omp_get_max_threads(); + #else + num_threads = 1; + #endif + childfopenMPI(sdt,"a", &fp); + fprintf(fp, "num_threads==%d\n",num_threads); + fclose(fp); + + //*[s] omp parallel + + TimeKeeper(X, cFileNameSzTimeKeep, cOMPSzStart, "a"); + TimeKeeper(X, cFileNameTimeKeep, cOMPSzStart, "a"); + switch(X->Def.iCalcModel){ + case HubbardGC: + icnt = X->Def.Tpow[2*X->Def.Nsite-1]*2+0;/*Tpow[2*X->Def.Nsit]=1*/ + break; + case SpinGC: + if(X->Def.iFlgGeneralSpin==FALSE){ + icnt = X->Def.Tpow[X->Def.Nsite-1]*2+0;/*Tpow[X->Def.Nsit]=1*/ + }else{ + icnt = X->Def.Tpow[X->Def.Nsite-1]*X->Def.SiteToBit[X->Def.Nsite-1]; + } + break; + case Hubbard: + hacker = X->Def.read_hacker; + if(hacker==0){ + calculate_jb_Hubbard(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; + }else if(hacker==1){ + calculate_jb_Hubbard_Hacker(X,list_jb,ihfbit,N2); + TimeKeeper(X, cFileNameSzTimeKeep, cOMPSzMid, "a"); + TimeKeeper(X, cFileNameTimeKeep, cOMPSzMid, "a"); + + icnt = 0; + #pragma omp parallel for default(none) \ + reduction(+:icnt) private(ib) firstprivate(ihfbit, X) \ + shared(list_1_, list_2_1_, list_2_2_, list_jb) + for(ib=0;ibCheck.sdim;ib++){ + icnt += omp_sz_hacker(ib,ihfbit,X,list_1_, list_2_1_, list_2_2_, list_jb); + } + break; + }else{ + fprintf(stderr, "Error: CalcHS in ModPara file must be 0 or 1 for Hubbard model."); + return -1; + } + case HubbardNConserved: + hacker = X->Def.read_hacker; + if(hacker==0){ + calculate_jb_HubbardNCoserved(X,list_jb,ihfbit,N2); + 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;ibCheck.sdim;ib++){ + icnt+=omp_sz(ib,ihfbit, X,list_1_, list_2_1_, list_2_2_, list_jb); + } + break; + }else if(hacker==1){ + calculate_jb_HubbardNCoserved_Hacker(X,list_jb,ihfbit,N2); + 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;ibCheck.sdim;ib++){ + icnt+=omp_sz_hacker(ib,ihfbit, X,list_1_, list_2_1_, list_2_2_, list_jb); + } + break; + }else{ + fprintf(stderr, "Error: CalcHS in ModPara file must be 0 or 1 for Hubbard model."); + return -1; + } + case KondoGC: + /*[s] this part can not be simply parallelized*/ + num_loc = count_localized_spins(X); + calculate_jb_KondoGC(X,num_loc,list_jb,ihfbit); + /*[e] this part can not be simply parallelized*/ + 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;ibCheck.sdim;ib++){ + icnt+=omp_sz_KondoGC(ib, ihfbit, X, list_1_, list_2_1_, list_2_2_, list_jb); + } + break; + case KondoNConserved: + 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;ibCheck.sdim;ib++){ + icnt+=omp_sz_KondoNConserved(ib,ihfbit, X, list_1_, list_2_1_, list_2_2_, list_jb); + } + BarrierMPI(); + //printf("AAA icnt=%ld\n",icnt); + break; + case Kondo: + calculate_jb_Kondo(X, list_jb,ihfbit); + TimeKeeper(X, cFileNameSzTimeKeep, cOMPSzMid, "a"); + TimeKeeper(X, cFileNameTimeKeep, cOMPSzMid, "a"); + + hacker = X->Def.read_hacker; + if(hacker==0){ + 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;ibCheck.sdim;ib++){ + icnt+=omp_sz_Kondo(ib,ihfbit, X, list_1_, list_2_1_, list_2_2_, list_jb); + } + }else if(hacker==1){ + 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;ibCheck.sdim;ib++){ + icnt+=omp_sz_Kondo_hacker(ib,ihfbit, X, list_1_, list_2_1_, list_2_2_, list_jb); + } + } + 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_tJ(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_tJ(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_tJ(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; + if(hacker == -1){ + calculate_jb_Spin_m1(X,list_jb,list_1_,list_2_1_,list_2_2_,ihfbit,irght,ilft,ibpatn, N); + }else if(hacker == 1){ + calculate_jb_Spin_Hacker(X,list_jb,ihfbit,N); + TimeKeeper(X, cFileNameSzTimeKeep, cOMPSzMid, "a"); + TimeKeeper(X, cFileNameTimeKeep, cOMPSzMid, "a"); + icnt = 0; + #pragma omp parallel for default(none)\ + reduction(+:icnt) private(ib) \ + firstprivate(ihfbit, N, X, list_1_, list_2_1_, list_2_2_, list_jb) + for(ib=0;ibCheck.sdim;ib++){ + icnt+=omp_sz_spin_hacker(ib,ihfbit,N,X, list_1_, list_2_1_, list_2_2_, list_jb); + } + }else if(hacker == 0){ + calculate_jb_Spin_Old(X,list_jb,ihfbit,N); + //#pragma omp barrier + TimeKeeper(X, cFileNameSzTimeKeep, cOMPSzMid, "a"); + TimeKeeper(X, cFileNameTimeKeep, cOMPSzMid, "a"); + + icnt = 0; + #pragma omp parallel for default(none) reduction(+:icnt)\ + private(ib) firstprivate(ihfbit, N, X)\ + shared(list_1_, list_2_1_, list_2_2_, list_jb) + for(ib=0;ibCheck.sdim;ib++){ + icnt+=omp_sz_spin(ib,ihfbit,N,X,list_1_, list_2_1_, list_2_2_, list_jb); + } + }else{ + fprintf(stderr, "Error: CalcHS in ModPara file must be -1 or 0 or 1 for Spin model."); + return -1; + } + }else{ + long unsigned int ilftdim = (X->Def.Tpow[X->Def.Nsite-1]*X->Def.SiteToBit[X->Def.Nsite-1])/ihfbit; + calculate_jb_GeneralSpin(X,list_jb,list_2_1_Sz,list_2_2_Sz,ihfbit,ilftdim,N); + TimeKeeper(X, cFileNameSzTimeKeep, cOMPSzMid, "a"); + TimeKeeper(X, cFileNameTimeKeep, cOMPSzMid, "a"); + + icnt = 0; + #pragma omp parallel for default(none)\ + reduction(+:icnt) private(ib) firstprivate(ilftdim, ihfbit, X)\ + shared(list_1_, list_2_1_, list_2_2_, list_2_1_Sz, list_2_2_Sz,list_jb) + for(ib=0;ibCheck.sdim;ib++){ - icnt+=omp_sz_hacker(ib,ihfbit, X,list_1_, list_2_1_, list_2_2_, list_jb); - } - - break; - } - else{ - fprintf(stderr, "Error: CalcHS in ModPara file must be 0 or 1 for Hubbard model."); - return -1; - } - - case Kondo: - // this part can not be parallelized - N_all_up = X->Def.Nup; - N_all_down = X->Def.Ndown; - fprintf(stdoutMPI, cStateNupNdown, N_all_up,N_all_down); - - jb = 0; - num_loc=0; - for(j=X->Def.Nsite/2; j< X->Def.Nsite ;j++){// counting localized # of spins - if(X->Def.LocSpn[j] != ITINERANT){ - num_loc += 1; - } - } - - for(ib=0;ibCheck.sdim;ib++){ //sdim = 2^(N/2) - list_jb[ib] = jb; - i = ib*ihfbit; // ihfbit=pow(2,((Nsite+1)/2)) - 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)){ // odd site - icheck_loc= icheck_loc; - ihfSpinDown=div_down; - if(div_down ==0){ - num_up += 1; + /* NConserved -> Normal */ + // this part is move to the ene of sz.c so that this change is made even + // when idim_max=0 + /* + if(X->Def.iFlgCalcSpec == CALCSPEC_NOT){ + if(X->Def.iCalcModel == HubbardNConserved){ + X->Def.iCalcModel = Hubbard; } - }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_up = (X->Def.Nsite+tmp_res)/2-all_loc; - all_down = (X->Def.Nsite-tmp_res)/2-all_loc; - if(X->Def.Nsite%2==1 && X->Def.LocSpn[X->Def.Nsite/2] != ITINERANT){ - all_up = (X->Def.Nsite)/2-all_loc; - all_down = (X->Def.Nsite)/2-all_loc; - } - - for(num_loc_up=0; num_loc_up <= all_loc; num_loc_up++){ - tmp_1 = Binomial(all_loc, num_loc_up, comb, all_loc); - if( X->Def.Nsite%2==1 && X->Def.LocSpn[X->Def.Nsite/2] != ITINERANT){ - if(ihfSpinDown !=0){ - tmp_2 = Binomial(all_up, X->Def.Nup-num_up-num_loc_up, comb, all_up); - tmp_3 = Binomial(all_down, X->Def.Ndown-num_down-(all_loc-num_loc_up), comb, all_down); + if(X->Def.iCalcModel == KondoNConserved){ + X->Def.iCalcModel = Kondo; } - else{ - tmp_2 = Binomial(all_up, X->Def.Nup-num_up-num_loc_up, comb, all_up); - tmp_3 = Binomial(all_down, X->Def.Ndown-num_down-(all_loc-num_loc_up), comb, all_down); + if(X->Def.iCalcModel == tJNConserved){ + X->Def.iCalcModel = tJ; } - } - else{ - tmp_2 = Binomial(all_up, X->Def.Nup-num_up-num_loc_up, comb, all_up); - tmp_3 = Binomial(all_down, X->Def.Ndown-num_down-(all_loc-num_loc_up), comb, all_down); - } - jb += tmp_1*tmp_2*tmp_3; - } - } - - } - //#pragma omp barrier - TimeKeeper(X, cFileNameSzTimeKeep, cOMPSzMid, "a"); - TimeKeeper(X, cFileNameTimeKeep, cOMPSzMid, "a"); - - hacker = X->Def.read_hacker; - if(hacker==0){ - 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;ibCheck.sdim;ib++){ - icnt+=omp_sz_Kondo(ib,ihfbit, X, list_1_, list_2_1_, list_2_2_, list_jb); - } - }else if(hacker==1){ - 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;ibCheck.sdim;ib++){ - icnt+=omp_sz_Kondo_hacker(ib,ihfbit, X, list_1_, list_2_1_, list_2_2_, list_jb); - } - } - break; - - case Spin: - // this part can not be parallelized - if(X->Def.iFlgGeneralSpin==FALSE){ - hacker = X->Def.read_hacker; - //printf(" rank=%d:Ne=%ld ihfbit=%ld sdim=%ld\n", myrank,X->Def.Ne,ihfbit,X->Check.sdim); - // using hacker's delight only + no open mp - if(hacker == -1){ - icnt = 1; - tmp_pow = 1; - tmp_i = 0; - jb = 0; - ja = 0; - while(tmp_pow < X->Def.Tpow[X->Def.Ne]){ - tmp_i += tmp_pow; - tmp_pow = tmp_pow*2; - } - //printf("DEBUG: %ld %ld %ld %ld\n",tmp_i,X->Check.sdim,X->Def.Tpow[X->Def.Ne],X->Def.Nsite); - if(X->Def.Nsite%2==0){ - max_tmp_i = X->Check.sdim*X->Check.sdim; - }else{ - max_tmp_i = X->Check.sdim*X->Check.sdim*2-1; - } - while(tmp_iCheck.sdim;ib++){ - list_jb[ib]=jb; - i=ib*ihfbit; - num_up=0; - for(j=0;jDef.Tpow[j]; - div_up = div_up/X->Def.Tpow[j]; - num_up +=div_up; - } - all_up = (X->Def.Nsite+1)/2; - tmp_1 = Binomial(all_up,X->Def.Ne-num_up,comb,all_up); - jb += tmp_1; - } - //#pragma omp barrier - TimeKeeper(X, cFileNameSzTimeKeep, cOMPSzMid, "a"); - TimeKeeper(X, cFileNameTimeKeep, cOMPSzMid, "a"); - - icnt = 0; -#pragma omp parallel for default(none) reduction(+:icnt) private(ib) firstprivate(ihfbit, N, X, list_1_, list_2_1_, list_2_2_, list_jb) - for(ib=0;ibCheck.sdim;ib++){ - icnt+=omp_sz_spin_hacker(ib,ihfbit,N,X, list_1_, list_2_1_, list_2_2_, list_jb); } - //printf(" rank=%d ib=%ld:Ne=%d icnt=%ld :idim_max=%ld N=%d\n", myrank,ib,X->Def.Ne,icnt,X->Check.idim_max,N); - // old version - }else if(hacker == 0){ - jb = 0; - for(ib=0;ibCheck.sdim;ib++){ - list_jb[ib]=jb; - i=ib*ihfbit; - num_up=0; - for(j=0;jDef.Tpow[j]; - div_up = div_up/X->Def.Tpow[j]; - num_up +=div_up; + */ + //Error message + if(i_max!=X->Check.idim_max){ + //printf("DDD %d %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); + if(childfopenMPI(sdt_err,"a",&fp_err)!=0){ + exitMPI(-1); } - all_up = (X->Def.Nsite+1)/2; - tmp_1 = Binomial(all_up,X->Def.Ne-num_up,comb,all_up); - jb += tmp_1; - } - //#pragma omp barrier - TimeKeeper(X, cFileNameSzTimeKeep, cOMPSzMid, "a"); - TimeKeeper(X, cFileNameTimeKeep, cOMPSzMid, "a"); - - icnt = 0; -#pragma omp parallel for default(none) reduction(+:icnt) private(ib) firstprivate(ihfbit, N, X) shared(list_1_, list_2_1_, list_2_2_, list_jb) - for(ib=0;ibCheck.sdim;ib++){ - icnt+=omp_sz_spin(ib,ihfbit,N,X,list_1_, list_2_1_, list_2_2_, list_jb); - } - } - else{ - fprintf(stderr, "Error: CalcHS in ModPara file must be -1 or 0 or 1 for Spin model."); - return -1; - } - }else{ - unsigned int Max2Sz=0; - unsigned int irghtsite=1; - long unsigned int itmpSize=1; - int i2Sz=0; - for(j=0; jDef.Nsite; j++){ - itmpSize *= X->Def.SiteToBit[j]; - if(itmpSize==ihfbit){ - break; + fprintf(fp_err, "%s",cErrSz_OutFile); + fclose(fp_err); + exitMPI(-1); } - irghtsite++; - } - for(j=0; jDef.Nsite; j++){ - Max2Sz += X->Def.LocSpn[j]; - } + free_li_2d_allocate(comb); + } + fprintf(stdoutMPI, "%s", cProEndCalcSz); + free(list_jb); - HilbertNumToSz = lui_1d_allocate(2*Max2Sz+1); - for(ib=0; ib<2*Max2Sz+1; ib++){ - HilbertNumToSz[ib]=0; + /* NConserved -> Normal */ + if(X->Def.iFlgCalcSpec == CALCSPEC_NOT){ + if(X->Def.iCalcModel == HubbardNConserved){ + X->Def.iCalcModel = Hubbard; } - - for(ib =0; ibDef.SiteToBit, X->Def.Tpow); - } - list_2_1_Sz[ib]=i2Sz; - HilbertNumToSz[i2Sz+Max2Sz]++; + if(X->Def.iCalcModel == KondoNConserved){ + X->Def.iCalcModel = Kondo; } - jb = 0; - long unsigned int ilftdim=(X->Def.Tpow[X->Def.Nsite-1]*X->Def.SiteToBit[X->Def.Nsite-1])/ihfbit; - for(ib=0;ibDef.SiteToBit, X->Def.Tpow); - } - list_2_2_Sz[ib]=i2Sz; - if((X->Def.Total2Sz- i2Sz +(int)Max2Sz)>=0 && (X->Def.Total2Sz- i2Sz) <= (int)Max2Sz){ - jb += HilbertNumToSz[X->Def.Total2Sz- i2Sz +Max2Sz]; - } - } - - TimeKeeper(X, cFileNameSzTimeKeep, cOMPSzMid, "a"); - TimeKeeper(X, cFileNameTimeKeep, cOMPSzMid, "a"); - - icnt = 0; -#pragma omp parallel for default(none) reduction(+:icnt) private(ib) firstprivate(ilftdim, ihfbit, X) shared(list_1_, list_2_1_, list_2_2_, list_2_1_Sz, list_2_2_Sz,list_jb) - for(ib=0;ibDef.iCalcModel == tJNConserved){ + X->Def.iCalcModel = tJ; } - - free_lui_1d_allocate(HilbertNumToSz); - } - - break; - default: - exitMPI(-1); - - } - i_max=icnt; - //fprintf(stdoutMPI, "Debug: Xicnt=%ld \n",icnt); - TimeKeeper(X, cFileNameSzTimeKeep, cOMPSzFinish, "a"); - TimeKeeper(X, cFileNameTimeKeep, cOMPSzFinish, "a"); - - } - - if(X->Def.iFlgCalcSpec == CALCSPEC_NOT){ - if(X->Def.iCalcModel==HubbardNConserved){ - X->Def.iCalcModel=Hubbard; } - } - - //Error message - //i_max=i_max+1; - if(i_max!=X->Check.idim_max){ - fprintf(stderr, "%s", cErrSz); - fprintf(stderr, cErrSz_ShowDim, i_max, X->Check.idim_max); - strcpy(sdt_err,cFileNameErrorSz); - if(childfopenMPI(sdt_err,"a",&fp_err)!=0){ - exitMPI(-1); - } - fprintf(fp_err, "%s",cErrSz_OutFile); - fclose(fp_err); - exitMPI(-1); - } - - free_li_2d_allocate(comb); - } - fprintf(stdoutMPI, "%s", cProEndCalcSz); - free(list_jb); if(X->Def.iFlgGeneralSpin==TRUE){ free(list_2_1_Sz); free(list_2_2_Sz); } - return 0; + BarrierMPI(); + return 0; } /** @@ -850,35 +495,174 @@ long int Binomial(int n,int k,long int **comb,int Nsite){ int tmp_i,tmp_j; if(n==0 && k==0){ - return 1; - } - else if(n<0 || k<0 || nDef.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 = div_up*div_down; + //printf("Ns %d i %ld j %ld div_up %ld div_down %ld \n",X->Def.Nsite,i,j,div_up,div_down); + if (check_doublon==1){ + break; + } + num_up += div_up; + num_down += div_down; + } + + if(check_doublon==0){ + 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 = 0; + 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 = div_up*div_down; + if (check_doublon==1){ + break; + } + num_up += div_up; + num_down += div_down; + } + if(check_doublon==0){ + 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; + 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 = div_up*div_down; + if (check_doublon==1){ + break; + } + num_up += div_up; + num_down += div_down; + } + if(check_doublon==0){ + 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; + 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 = div_up*div_down; + if (check_doublon==1){ + break; + } + num_up += div_up; + num_down += div_down; + } + if(check_doublon==0){ + 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 * @@ -977,7 +761,7 @@ int omp_sz( } /** * @brief efficient version of calculating restricted Hilbert space for Hubbard systems using snoob - * details of snoob is found in S.H. Warren, Hacker’s Delight, second ed., Addison-Wesley, ISBN: 0321842685, 2012. + * details of snoob is found in S.H. Warren, Hacker's delight(Bs Delight, second ed., Addison-Wesley, ISBN: 0321842685, 2012. * * @param[in] ib upper half bit of i * @param[in] ihfbit 2^(Ns/2) @@ -1064,35 +848,150 @@ int omp_sz_hacker(long unsigned int ib, } ia = snoob(ia); } - } - } - } - } - else if(X->Def.iCalcModel==HubbardNConserved){ - if(tmp_num_up+tmp_num_down <= X->Def.Ne){ //do not exceed Ne - ia = X->Def.Tpow[X->Def.Ne-tmp_num_up-tmp_num_down]-1; - if(ia < X->Check.sdim){ + } + } + } + } + else if(X->Def.iCalcModel==HubbardNConserved){ + if(tmp_num_up+tmp_num_down <= X->Def.Ne){ //do not exceed Ne + ia = X->Def.Tpow[X->Def.Ne-tmp_num_up-tmp_num_down]-1; + if(ia < X->Check.sdim){ + list_1_[ja+jb]=ia+ib*ihfbit; + list_2_1_[ia]=ja+1; + list_2_2_[ib]=jb+1; + ja+=1; + if(ia!=0){ + ia = snoob(ia); + while(ia < X->Check.sdim){ + list_1_[ja+jb]=ia+ib*ihfbit; + list_2_1_[ia]=ja+1; + list_2_2_[ib]=jb+1; + ja+=1; + ia = snoob(ia); + } + } + } + } + } + ja=ja-1; + return ja; +} + +/** + * @brief calculating restricted Hilbert space for KondoNConserved 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_KondoNConserved( + 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 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;iaCheck.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+num_down == X->Def.Ne && icheck_loc==1){ list_1_[ja+jb]=ia+ib*ihfbit; + /* + list_2_1_[ia]=ja; + list_2_2_[ib]=jb; + */ 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; - if(ia!=0){ - ia = snoob(ia); - while(ia < X->Check.sdim){ - list_1_[ja+jb]=ia+ib*ihfbit; - list_2_1_[ia]=ja+1; - list_2_2_[ib]=jb+1; - ja+=1; - ia = snoob(ia); - } - } - } + } } } + //printf("XXX ib %lu ja %lu\n",ib,ja); ja=ja-1; return ja; } + + /** * @brief calculating restricted Hilbert space for Kondo systems * @@ -1478,7 +1377,7 @@ int omp_sz_spin( /** * @brief efficient version of calculating restricted Hilbert space for spin-1/2 systems - * details of snoob is found in S.H. Warren, Hacker’s Delight, second ed., Addison-Wesley, ISBN: 0321842685, 2012. + * details of snoob is found in S.H. Warren, Hacker's delight(Bs Delight, second ed., Addison-Wesley, ISBN: 0321842685, 2012. * * @param[in] ib upper half bit of i * @param[in] ihfbit 2^(Ns/2) @@ -1638,6 +1537,9 @@ int Read_sz case Kondo: sprintf(sdt,"ListForKondo_Ns%d_Ncond%d.dat",X->Def.Nsite,X->Def.Ne); break; + case KondoNConserved: + sprintf(sdt,"ListForKondo_Ns%d_Ncond%d.dat",X->Def.Nsite,X->Def.Ne); + break; } if(childfopenMPI(sdt,"r", &fp)!=0){ exitMPI(-1); @@ -1683,3 +1585,744 @@ int Read_sz return 0; } + +int count_localized_spins(struct BindStruct *X){ + int num_loc = 0; + for (int j = X->Def.Nsite / 2; j < X->Def.Nsite; j++) { /*counting # of localized spins*/ + if(X->Def.LocSpn[j] != ITINERANT) {/*ITINERANT ==0 -> itinerant*/ + num_loc += 1; + } + } + return num_loc; +} + +void calculate_jb_GeneralSpin(struct BindStruct *X, long unsigned int *list_jb, long int *list_2_1_Sz,long int *list_2_2_Sz, long unsigned int ihfbit,long unsigned int ilftdim,unsigned int N){ + unsigned int Max2Sz=0; + unsigned int irghtsite=1; + long unsigned int itmpSize=1; + int i2Sz=0; + long unsigned int *HilbertNumToSz; + long unsigned int ib,jb,j; + + for(j=0; jDef.Nsite; j++){ + itmpSize *= X->Def.SiteToBit[j]; + if(itmpSize==ihfbit){ + break; + } + irghtsite++; + } + for(j=0; jDef.Nsite; j++){ + Max2Sz += X->Def.LocSpn[j]; + } + + HilbertNumToSz = lui_1d_allocate(2*Max2Sz+1); + for(ib=0; ib<2*Max2Sz+1; ib++){ + HilbertNumToSz[ib]=0; + } + + for(ib =0; ibDef.SiteToBit, X->Def.Tpow); + } + list_2_1_Sz[ib]=i2Sz; + HilbertNumToSz[i2Sz+Max2Sz]++; + } + jb = 0; + for(ib=0;ibDef.SiteToBit, X->Def.Tpow); + } + list_2_2_Sz[ib]=i2Sz; + if((X->Def.Total2Sz- i2Sz +(int)Max2Sz)>=0 && (X->Def.Total2Sz- i2Sz) <= (int)Max2Sz){ + jb += HilbertNumToSz[X->Def.Total2Sz- i2Sz +Max2Sz]; + } + } + free_lui_1d_allocate(HilbertNumToSz); +} + +void calculate_jb_Spin_m1(struct BindStruct *X, long unsigned int *list_jb, long unsigned int *list_1_, long unsigned int *list_2_1_,long unsigned int *list_2_2_,\ +long unsigned int ihfbit,long unsigned int irght,long unsigned int ilft,long unsigned int ibpatn, unsigned int N){ + long unsigned int ia,ja,ib,jb,div_up,i,j,tmp_1; + long unsigned int tmp_pow,tmp_i,tmp_j,icnt,max_tmp_i; + int num_up,all_up; + + icnt = 1; + tmp_pow = 1; + tmp_i = 0; + jb = 0; + ja = 0; + while(tmp_pow < X->Def.Tpow[X->Def.Ne]){ + tmp_i += tmp_pow; + tmp_pow = tmp_pow*2; + } + if(X->Def.Nsite%2==0){ + max_tmp_i = X->Check.sdim*X->Check.sdim; + }else{ + max_tmp_i = X->Check.sdim*X->Check.sdim*2-1; + } + while(tmp_iDef.Nsite+1,X->Def.Nsite+1); + jb = 0; + for(ib=0;ibCheck.sdim;ib++){ + list_jb[ib] = jb; + i = ib*ihfbit; + num_up = 0; + for(j=0;jDef.Tpow[j]; + div_up = div_up/X->Def.Tpow[j]; + num_up += div_up; + } + all_up = (X->Def.Nsite+1)/2; + tmp_1 = Binomial(all_up,X->Def.Ne-num_up,comb,all_up); + jb += tmp_1; + } + free_li_2d_allocate(comb); +} + +void calculate_jb_Spin_Old(struct BindStruct *X, long unsigned int *list_jb, long unsigned int ihfbit,unsigned int N){ + long unsigned int ib,jb,div_up,i,j,tmp_1; + int num_up,all_up; + long int **comb; + comb = li_2d_allocate(X->Def.Nsite+1,X->Def.Nsite+1); + for(ib=0;ibCheck.sdim;ib++){ + list_jb[ib] = jb; + i = ib*ihfbit; + num_up=0; + for(j=0;jDef.Tpow[j]; + div_up = div_up/X->Def.Tpow[j]; + num_up += div_up; + } + all_up = (X->Def.Nsite+1)/2; + tmp_1 = Binomial(all_up,X->Def.Ne-num_up,comb,all_up); + jb += tmp_1; + } + free_li_2d_allocate(comb); +} + +void calculate_jb_Kondo(struct BindStruct *X, long unsigned int *list_jb, long unsigned int ihfbit){ + long unsigned int jb = 0,i,ib,j; + long unsigned int div_up, div_down,ihfSpinDown; + int N_all_up, N_all_down; + 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; + + 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); + fprintf(stdoutMPI, cStateNupNdown, N_all_up,N_all_down); + + jb = 0; + num_loc = 0; + for(j=X->Def.Nsite/2; j< X->Def.Nsite ;j++){// counting localized # of spins + if(X->Def.LocSpn[j] != ITINERANT){ + num_loc += 1; + } + } + //printf("DEBUG num_loc %d \n",num_loc); + for(ib=0;ibCheck.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); + + 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]; + 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{ + /*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_up = (X->Def.Nsite+tmp_res)/2-all_loc; + all_down = (X->Def.Nsite-tmp_res)/2-all_loc; + if(X->Def.Nsite%2==1 && X->Def.LocSpn[X->Def.Nsite/2] != ITINERANT){ + 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); + + for(num_loc_up=0; num_loc_up <= all_loc; num_loc_up++){ + tmp_1 = Binomial(all_loc, num_loc_up, comb, all_loc); + if( X->Def.Nsite%2==1 && X->Def.LocSpn[X->Def.Nsite/2] != ITINERANT){ + if(ihfSpinDown !=0){ + tmp_2 = Binomial(all_up, X->Def.Nup-num_up-num_loc_up, comb, all_up); + tmp_3 = Binomial(all_down, X->Def.Ndown-num_down-(all_loc-num_loc_up), comb, all_down); + }else{ + tmp_2 = Binomial(all_up, X->Def.Nup-num_up-num_loc_up, comb, all_up); + tmp_3 = Binomial(all_down, X->Def.Ndown-num_down-(all_loc-num_loc_up), comb, all_down); + } + }else{ + tmp_2 = Binomial(all_up, X->Def.Nup-num_up-num_loc_up, comb, all_up); + tmp_3 = Binomial(all_down, X->Def.Ndown-num_down-(all_loc-num_loc_up), comb, all_down); + } + jb += tmp_1*tmp_2*tmp_3; + } + } + } + free_li_2d_allocate(comb); +} + +void calculate_jb_KondoNConserved(struct BindStruct *X, long unsigned int *list_jb, long unsigned int ihfbit){ + long unsigned int jb = 0,i,ib,j; + long unsigned int div_up, div_down,ihfSpinDown; + int Ne; + 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; + int all_cond,all_cond_up,all_cond_down,num_cond; + + Ne = X->Def.Ne; + + jb = 0; + num_loc = 0; + for(j=X->Def.Nsite/2; j< X->Def.Nsite ;j++){// counting localized # of spins + if(X->Def.LocSpn[j] != ITINERANT){ + num_loc += 1; + } + } + //printf("DEBUG num_loc %d \n",num_loc); + all_loc = X->Def.NLocSpn; + for(ib=0;ibCheck.sdim;ib++){ //sdim = 2^(N/2) + list_jb[ib] = jb; + i = ib*ihfbit; // ihfbit=pow(2,((Nsite+1)/2)) + 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{ + /*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_up = (X->Def.Nsite+tmp_res)/2-all_loc; + all_down = (X->Def.Nsite-tmp_res)/2-all_loc; + if(X->Def.Nsite%2==1 && X->Def.LocSpn[X->Def.Nsite/2] != ITINERANT){ + 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 Ne %d Nsite %d\n",all_loc,all_up,all_down,X->Def.Ne,X->Def.Nsite); + if (num_up+num_down==X->Def.Ne-all_loc){ + jb += X->Def.Tpow[all_loc]; + } + } + } +/* + list_jb[ib] = jb; + i = ib*ihfbit; // ihfbit=pow(2,((Nsite+1)/2)) + num_up = 0; + num_down = 0; + 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]; + 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; + } + } + if (num_up+num_down==X->Def.Ne-all_loc){ + jb += X->Def.Tpow[all_loc]; + } +*/ +} + +void calculate_jb_KondoGC(struct BindStruct *X, int num_loc, long unsigned int *list_jb, long unsigned int ihfbit){ + long unsigned int jb = 0; + for(long unsigned int ib=0;ib < X->Check.sdim;ib++){ + list_jb[ib] = jb; + long unsigned int i = ib*ihfbit; + int icheck_loc = 1; + for(long unsigned int j=(X->Def.Nsite+1)/2; j< X->Def.Nsite ;j++){ + long unsigned int div_up = i & X->Def.Tpow[2*j]; + div_up = div_up/X->Def.Tpow[2*j]; + long unsigned int 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){ + 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){ + if(X->Def.Nsite%2==1 && X->Def.LocSpn[X->Def.Nsite/2] != ITINERANT){ + jb += X->Def.Tpow[X->Def.Nsite-1-(X->Def.NLocSpn-num_loc)]; + }else{ + jb += X->Def.Tpow[X->Def.Nsite-(X->Def.NLocSpn-num_loc)]; + } + } + } +} + + +void calculate_jb_Hubbard(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,i,j,tmp_1,tmp_2; + long int **comb; + int num_up,num_down; + 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; + //[s] counting # of up and down electrons + num_up = 0; + for(j=0;j<=N2-2;j+=2){ // even -> up spin + div = i & X->Def.Tpow[j]; + div = div/X->Def.Tpow[j]; + num_up += div; + } + num_down=0; + for(j=1;j<=N2-1;j+=2){ // odd -> down spin + div=i & X->Def.Tpow[j]; + div=div/X->Def.Tpow[j]; + num_down+=div; + } + //[e] counting # of up and down electrons + 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.Ndown-num_down,comb,all_down); + jb += tmp_1*tmp_2; + } + free_li_2d_allocate(comb); + /*[e] this part can not be parallelized*/ +} + +void calculate_jb_Hubbard_Hacker(struct BindStruct *X,long unsigned int *list_jb, long unsigned int ihfbit, unsigned int N2){ + long unsigned int sdim_div,sdim_rest,ib_start,ib_end; + long unsigned int i,ib,j,jb,div,tmp_res,tmp_1,tmp_2; + int mythread,num_up,num_down,all_up,all_down; + long int **comb2; + long unsigned int *jbthread; + + jbthread = lui_1d_allocate(nthreads); + #pragma omp parallel default(none) \ + shared(X,list_jb,ihfbit,N2,nthreads,jbthread) \ + private(ib,i,j,num_up,num_down,div,tmp_res,tmp_1,tmp_2,jb,all_up,all_down, \ + comb2,mythread,sdim_div,sdim_rest,ib_start,ib_end) + { + jb = 0; + #ifdef _OPENMP + mythread = omp_get_thread_num(); + #else + mythread = 0; + #endif + comb2 = li_2d_allocate(X->Def.Nsite+1,X->Def.Nsite+1); + /* explict loop decomposition is nessesary to fix the asignment to each thread*/ + sdim_div = X->Check.sdim / nthreads; + sdim_rest = X->Check.sdim % nthreads; + if(mythread < sdim_rest){ + ib_start = sdim_div*mythread + mythread; + ib_end = ib_start + sdim_div + 1; + }else{ + ib_start = sdim_div*mythread + sdim_rest; + ib_end = ib_start + sdim_div; + } + for(ib=ib_start;ibDef.Tpow[j]; + div = div/X->Def.Tpow[j]; + num_up += div; + } + num_down=0; + for(j=1;j<=N2-1;j+=2){ + div=i & X->Def.Tpow[j]; + div=div/X->Def.Tpow[j]; + num_down+=div; + } + + 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,comb2,all_up); + tmp_2 = Binomial(all_down,X->Def.Ndown-num_down,comb2,all_down); + jb += tmp_1*tmp_2; + } + free_li_2d_allocate(comb2); + if(mythread != nthreads-1) jbthread[mythread+1] = jb; + #pragma omp barrier + #pragma omp single + { + jbthread[0] = 0; + for(j=1;jDef.Nsite+1,X->Def.Nsite+1); + // this part can not be parallelized + jb = 0; + iSpnup = 0; + 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(ib=0;ibCheck.sdim;ib++){ + list_jb[ib] = jb; + i = ib*ihfbit; + num_up = 0; + for(j=0;j<=N2-2;j+=2){ + div = i & X->Def.Tpow[j]; + div = div/X->Def.Tpow[j]; + num_up += div; + } + num_down=0; + for(j=1;j<=N2-1;j+=2){ + div = i & X->Def.Tpow[j]; + div = div/X->Def.Tpow[j]; + num_down += div; + } + 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, X->Def.Ne-iSpnup-num_down,comb,all_down); + jb += tmp_1*tmp_2; + } + } + free_li_2d_allocate(comb); +} + +void calculate_jb_HubbardNCoserved_Hacker(struct BindStruct *X,long unsigned int *list_jb, long unsigned int ihfbit, unsigned int N2){ + long unsigned int sdim_div,sdim_rest,ib_start,ib_end; + long unsigned int i,ib,j,jb,div,tmp_res,tmp_1,tmp_2; + int mythread,num_up,num_down,all_up,all_down; + long int **comb2; + long unsigned int *jbthread; + int iSpnup, iMinup,iAllup; + + iMinup = 0; + iAllup = X->Def.Ne; + if(X->Def.Ne > X->Def.Nsite){ + iMinup = X->Def.Ne-X->Def.Nsite; + iAllup = X->Def.Nsite; + } + jbthread = lui_1d_allocate(nthreads); + #pragma omp parallel default(none) \ + shared(X,iMinup,iAllup,list_jb,ihfbit,N2,nthreads,jbthread) \ + private(ib,i,j,num_up,num_down,div,tmp_res,tmp_1,tmp_2,iSpnup,jb,all_up,all_down,comb2, \ + mythread,sdim_rest,sdim_div,ib_start,ib_end) + { + jb = 0; + iSpnup=0; + #ifdef _OPENMP + mythread = omp_get_thread_num(); + #else + mythread = 0; + #endif + comb2 = li_2d_allocate(X->Def.Nsite+1,X->Def.Nsite+1); + /* explict loop decomposition is nessesary to fix the asignment to each thread*/ + sdim_div = X->Check.sdim / nthreads; + sdim_rest = X->Check.sdim % nthreads; + if(mythread < sdim_rest){ + ib_start = sdim_div*mythread + mythread; + ib_end = ib_start + sdim_div + 1; + }else{ + ib_start = sdim_div*mythread + sdim_rest; + ib_end = ib_start + sdim_div; + } + for(ib=ib_start;ibDef.Tpow[j]; + div = div/X->Def.Tpow[j]; + num_up += div; + } + num_down=0; + for(j=1;j<=N2-1;j+=2){ + div=i & X->Def.Tpow[j]; + div=div/X->Def.Tpow[j]; + num_down+=div; + } + 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,comb2,all_up); + tmp_2 = Binomial(all_down, X->Def.Ne-iSpnup-num_down,comb2,all_down); + jb += tmp_1*tmp_2; + } + } + free_li_2d_allocate(comb2); + if(mythread != nthreads-1) jbthread[mythread+1] = jb; + #pragma omp barrier + #pragma omp single + { + jbthread[0] = 0; + for(j=1;jDef.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 = 0; + //[s] counting # of up and down electrons + num_up = 0; + num_down = 0; + for(j=0;j<(N2/2);j++){ + 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==1){ + 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==0){ + 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);/* tJ all_down-(X->Def.Nup-num_up)*/ + jb += tmp_1*tmp_2; + //printf("DBB jb=%ld %ld %ld: num_up %d num_down %d : %d %d\n",jb,tmp_1,tmp_2,num_up,num_down,(X->Def.Nup-num_up),X->Def.Ndown-num_down); + } + } + 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 = 0; + //[s] counting # of up and down electrons + num_up = 0; + num_down = 0; + for(j=0;j<(N2/2);j++){ + 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==1){ + 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==0){ + 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);/* 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 = 0; + //[s] counting # of up and down electrons + num_up = 0; + num_down = 0; + for(j=0;j<(N2/2);j++){ + 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==1){ + 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==0){ + 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*/ +} + diff --git a/src/time.c b/src/time.c index f40e28b4..1bb2b3a1 100644 --- a/src/time.c +++ b/src/time.c @@ -194,6 +194,8 @@ void OutputTimer(struct BindStruct *X) { break; case Hubbard: + case tJ: + case tJGC: StampTime(fp," Hubbard", 300); StampTime(fp," trans in Hubbard", 310); StampTime(fp," double", 311); diff --git a/src/xsetmem.c b/src/xsetmem.c index 6ea4bc57..6b54ca58 100644 --- a/src/xsetmem.c +++ b/src/xsetmem.c @@ -294,7 +294,11 @@ int setmem_large case Hubbard: case HubbardNConserved: case Kondo: + case KondoNConserved: case KondoGC: + case tJ: + case tJNConserved: + case tJGC: if (X->Def.iFlgGeneralSpin == FALSE) { if (X->Def.iCalcModel == Spin && X->Def.Nsite % 2 == 1) { X->Large.SizeOflist_2_1 = X->Check.sdim * 2 + 2;