diff --git a/src/expec_cisajscktaltdc.c b/src/expec_cisajscktaltdc.c index 4cf310eb..5cfee092 100644 --- a/src/expec_cisajscktaltdc.c +++ b/src/expec_cisajscktaltdc.c @@ -782,7 +782,7 @@ int expec_cisajscktalt_SpinHalf(struct BindStruct *X,double complex *vec, FILE * * */ -int expec_ThreeBody_SpinGeneral(struct BindStruct *X,double complex *vec, FILE **_fp){ +int expec_Threebody_SpinGeneral(struct BindStruct *X,double complex *vec, FILE **_fp){ long unsigned int i,j; long unsigned int org_isite1,org_isite2,org_isite3,org_isite4,org_isite5,org_isite6; long unsigned int org_sigma1,org_sigma2,org_sigma3,org_sigma4,org_sigma5,org_sigma6; @@ -804,20 +804,26 @@ int expec_ThreeBody_SpinGeneral(struct BindStruct *X,double complex *vec, FILE * for(i=0;iDef.NTBody;i++){ vec_pr = cd_1d_allocate(i_max + 1); - tmp_org_isite1 = X->Def.CisAjtCkuAlvDC[i][0]+1; - tmp_org_sigma1 = X->Def.CisAjtCkuAlvDC[i][1]; - tmp_org_isite2 = X->Def.CisAjtCkuAlvDC[i][2]+1; - tmp_org_sigma2 = X->Def.CisAjtCkuAlvDC[i][3]; - tmp_org_isite3 = X->Def.CisAjtCkuAlvDC[i][4]+1; - tmp_org_sigma3 = X->Def.CisAjtCkuAlvDC[i][5]; - tmp_org_isite4 = X->Def.CisAjtCkuAlvDC[i][6]+1; - tmp_org_sigma4 = X->Def.CisAjtCkuAlvDC[i][7]; + tmp_org_isite1 = X->Def.TBody[i][0]+1; + tmp_org_sigma1 = X->Def.TBody[i][1]; + tmp_org_isite2 = X->Def.TBody[i][2]+1; + tmp_org_sigma2 = X->Def.TBody[i][3]; + tmp_org_isite3 = X->Def.TBody[i][4]+1; + tmp_org_sigma3 = X->Def.TBody[i][5]; + tmp_org_isite4 = X->Def.TBody[i][6]+1; + tmp_org_sigma4 = X->Def.TBody[i][7]; /*[s]For three body*/ - org_isite5 = X->Def.TBody[i][8]+1; - org_sigma5 = X->Def.TBody[i][9]; - org_isite6 = X->Def.TBody[i][10]+1; - org_sigma6 = X->Def.TBody[i][11]; + tmp_org_isite5 = X->Def.TBody[i][8]+1; + tmp_org_sigma5 = X->Def.TBody[i][9]; + tmp_org_isite6 = X->Def.TBody[i][10]+1; + tmp_org_sigma6 = X->Def.TBody[i][11]; + /**/ + org_isite5 = X->Def.TBody[i][8]+1; + org_sigma5 = X->Def.TBody[i][9]; + org_isite6 = X->Def.TBody[i][10]+1; + org_sigma6 = X->Def.TBody[i][11]; + dam_pr = 0.0; X->Large.mode = M_MLTPLY; /* |vec_pr>= c5a6|phi>*/ @@ -825,86 +831,95 @@ int expec_ThreeBody_SpinGeneral(struct BindStruct *X,double complex *vec, FILE * X->Large.mode = M_CORR; /*[e]For three body*/ - if(Rearray_Interactions(i, &org_isite1, &org_isite2, &org_isite3, &org_isite4, &org_sigma1, &org_sigma2, &org_sigma3, &org_sigma4, &tmp_V, X,2)!=0){ + if(Rearray_Interactions(i, &org_isite1, &org_isite2, &org_isite3, &org_isite4, &org_sigma1, &org_sigma2, &org_sigma3, &org_sigma4, &tmp_V, X,3)!=0){ fprintf(*_fp," %4ld %4ld %4ld %4ld %4ld %4ld %4ld %4ld %.10lf %.10lf \n",tmp_org_isite1-1, tmp_org_sigma1, tmp_org_isite2-1, tmp_org_sigma2, tmp_org_isite3-1,tmp_org_sigma3, tmp_org_isite4-1, tmp_org_sigma4,0.0,0.0); continue; } - tmp_Sz=0; - - for(j=0;j<2; j++) { - tmp_org = X->Def.CisAjtCkuAlvDC[i][4*j+1]*X->Def.Tpow[X->Def.CisAjtCkuAlvDC[i][4 * j]]; - tmp_Sz += GetLocal2Sz(X->Def.CisAjtCkuAlvDC[i][4 * j] + 1, tmp_org, X->Def.SiteToBit, X->Def.Tpow); - tmp_org = X->Def.CisAjtCkuAlvDC[i][4*j+3]*X->Def.Tpow[X->Def.CisAjtCkuAlvDC[i][4 * j+2]]; - tmp_Sz -= GetLocal2Sz(X->Def.CisAjtCkuAlvDC[i][4 * j+2] + 1, tmp_org, X->Def.SiteToBit, X->Def.Tpow); - } - if(tmp_Sz !=0){ // not Sz conserved - fprintf(*_fp," %4ld %4ld %4ld %4ld %4ld %4ld %4ld %4ld %.10lf %.10lf \n",tmp_org_isite1-1, tmp_org_sigma1, tmp_org_isite2-1, tmp_org_sigma2, tmp_org_isite3-1,tmp_org_sigma3, tmp_org_isite4-1, tmp_org_sigma4,0.0,0.0); - continue; - } - - dam_pr = 0.0; - if(org_isite1 >X->Def.Nsite && org_isite3>X->Def.Nsite){ + if(org_isite1 > X->Def.Nsite && org_isite3 > X->Def.Nsite){ if(org_sigma1==org_sigma2 && org_sigma3==org_sigma4 ){ //diagonal - dam_pr=child_CisAisCjuAju_GeneralSpin_MPIdouble(org_isite1-1, org_sigma1, org_isite3-1, org_sigma3, tmp_V, X, vec, vec); + dam_pr=child_GC_CisAisCjuAju_GeneralSpin_MPIdouble(org_isite1-1, org_sigma1, org_isite3-1, org_sigma3, tmp_V, X, vec, vec_pr); } - else if(org_sigma1 != org_sigma2 && org_sigma3 != org_sigma4){ - dam_pr=child_CisAitCjuAjv_GeneralSpin_MPIdouble(org_isite1-1, org_sigma1, org_sigma2, org_isite3-1, org_sigma3, org_sigma4,tmp_V, X, vec, vec); + else if(org_sigma1 == org_sigma2 && org_sigma3 != org_sigma4){ + dam_pr=child_GC_CisAisCjuAjv_GeneralSpin_MPIdouble(org_isite1-1, org_sigma1, org_isite3-1, org_sigma3, org_sigma4, tmp_V, X, vec, vec_pr); } - else{ - dam_pr=0.0; + else if(org_sigma1 != org_sigma2 && org_sigma3 == org_sigma4){ + dam_pr=child_GC_CisAitCjuAju_GeneralSpin_MPIdouble(org_isite1-1, org_sigma1, org_sigma2, org_isite3-1, org_sigma3, tmp_V, X, vec, vec_pr); + } + else if(org_sigma1 != org_sigma2 && org_sigma3 != org_sigma4){ + dam_pr=child_GC_CisAitCjuAjv_GeneralSpin_MPIdouble(org_isite1-1, org_sigma1, org_sigma2, org_isite3-1, org_sigma3, org_sigma4,tmp_V, X, vec, vec_pr); } } else if(org_isite3 > X->Def.Nsite || org_isite1 > X->Def.Nsite){ if(org_sigma1==org_sigma2 && org_sigma3==org_sigma4 ){ //diagonal - dam_pr=child_CisAisCjuAju_GeneralSpin_MPIsingle(org_isite1-1, org_sigma1, org_isite3-1, org_sigma3, tmp_V, X, vec, vec); + dam_pr=child_GC_CisAisCjuAju_GeneralSpin_MPIsingle(org_isite1-1, org_sigma1, org_isite3-1, org_sigma3, tmp_V, X, vec, vec_pr); } - else if(org_sigma1 != org_sigma2 && org_sigma3 != org_sigma4){ - dam_pr=child_CisAitCjuAjv_GeneralSpin_MPIsingle(org_isite1-1, org_sigma1, org_sigma2, org_isite3-1, org_sigma3, org_sigma4,tmp_V, X, vec, vec); + else if(org_sigma1 == org_sigma2 && org_sigma3 != org_sigma4){ + dam_pr=child_GC_CisAisCjuAjv_GeneralSpin_MPIsingle(org_isite1-1, org_sigma1, org_isite3-1, org_sigma3, org_sigma4, tmp_V, X, vec, vec_pr); } - else{ - dam_pr=0.0; + else if(org_sigma1 != org_sigma2 && org_sigma3 == org_sigma4){ + dam_pr=child_GC_CisAitCjuAju_GeneralSpin_MPIsingle(org_isite1-1, org_sigma1, org_sigma2, org_isite3-1, org_sigma3, tmp_V, X, vec, vec_pr); + } + else if(org_sigma1 != org_sigma2 && org_sigma3 != org_sigma4){ + dam_pr=child_GC_CisAitCjuAjv_GeneralSpin_MPIsingle(org_isite1-1, org_sigma1, org_sigma2, org_isite3-1, org_sigma3, org_sigma4,tmp_V, X, vec, vec_pr); } } else{ if(org_sigma1==org_sigma2 && org_sigma3==org_sigma4 ){ //diagonal -#pragma omp parallel for default(none) reduction(+:dam_pr) private(j, num1) firstprivate(i_max,X,org_isite1, org_sigma1,org_isite3, org_sigma3, tmp_V) shared(vec,list_1) + #pragma omp parallel for default(none) reduction(+:dam_pr) private(j, num1) firstprivate(i_max,X,org_isite1, org_sigma1,org_isite3, org_sigma3, tmp_V) shared(vec,vec_pr) for(j=1;j<=i_max;j++){ - num1=BitCheckGeneral(list_1[j], org_isite1, org_sigma1, X->Def.SiteToBit, X->Def.Tpow); + num1=BitCheckGeneral(j-1, org_isite1, org_sigma1, X->Def.SiteToBit, X->Def.Tpow); if(num1 != FALSE){ - num1=BitCheckGeneral(list_1[j], org_isite3, org_sigma3, X->Def.SiteToBit, X->Def.Tpow); + num1=BitCheckGeneral(j-1, org_isite3, org_sigma3, X->Def.SiteToBit, X->Def.Tpow); if(num1 != FALSE){ - dam_pr += tmp_V*conj(vec[j])*vec[j]; + dam_pr += tmp_V*conj(vec[j])*vec_pr[j]; } } } - } - else if(org_sigma1 != org_sigma2 && org_sigma3 != org_sigma4){ -#pragma omp parallel for default(none) reduction(+:dam_pr) private(j, num1) firstprivate(i_max,X, org_isite1, org_isite3, org_sigma1, org_sigma2, org_sigma3, org_sigma4, tmp_off, tmp_off_2, list1_off, myrank, tmp_V) shared(vec, list_1) + }else if(org_sigma1 == org_sigma2 && org_sigma3 != org_sigma4){ + #pragma omp parallel for default(none) reduction(+:dam_pr) private(j, num1) firstprivate(i_max,X, org_isite1, org_isite3, org_sigma1,org_sigma3,org_sigma4, tmp_off, tmp_V) shared(vec,vec_pr) for(j=1;j<=i_max;j++){ - num1 = GetOffCompGeneralSpin(list_1[j], org_isite3, org_sigma4, org_sigma3, &tmp_off, X->Def.SiteToBit, X->Def.Tpow); - if(num1 != FALSE) { - num1 = GetOffCompGeneralSpin(tmp_off, org_isite1, org_sigma2, org_sigma1, &tmp_off_2, - X->Def.SiteToBit, X->Def.Tpow); - if (num1 != FALSE) { - ConvertToList1GeneralSpin(tmp_off_2, X->Check.sdim, &list1_off); - dam_pr += tmp_V * conj(vec[list1_off]) * vec[j]; + num1 = GetOffCompGeneralSpin(j-1, org_isite3, org_sigma4, org_sigma3, &tmp_off, X->Def.SiteToBit, X->Def.Tpow); + if(num1 != FALSE){ + num1=BitCheckGeneral(tmp_off, org_isite1, org_sigma1, X->Def.SiteToBit, X->Def.Tpow); + if(num1 != FALSE){ + dam_pr += tmp_V*conj(vec[tmp_off+1])*vec_pr[j]; + } } - } } - //printf("DEBUG: rank=%d, dam_pr=%lf\n", myrank, creal(dam_pr)); - } - else{ - dam_pr=0.0; + }else if(org_sigma1 != org_sigma2 && org_sigma3 == org_sigma4){ + #pragma omp parallel for default(none) reduction(+:dam_pr) private(j, num1) firstprivate(i_max,X, org_isite1, org_isite3, org_sigma1,org_sigma2, org_sigma3, tmp_off, tmp_V) shared(vec,vec_pr) + for(j=1;j<=i_max;j++){ + num1 = BitCheckGeneral(j-1, org_isite3, org_sigma3, X->Def.SiteToBit, X->Def.Tpow); + if(num1 != FALSE){ + num1 = GetOffCompGeneralSpin(j-1, org_isite1, org_sigma2, org_sigma1, &tmp_off, X->Def.SiteToBit, X->Def.Tpow); + if(num1 != FALSE){ + dam_pr += tmp_V*conj(vec[tmp_off+1])*vec_pr[j]; + } + } + } + }else if(org_sigma1 != org_sigma2 && org_sigma3 != org_sigma4){ + #pragma omp parallel for default(none) reduction(+:dam_pr) private(j, num1) firstprivate(i_max,X, org_isite1, org_isite3, org_sigma1, org_sigma2, org_sigma3, org_sigma4, tmp_off, tmp_off_2, tmp_V) shared(vec,vec_pr) + for(j=1;j<=i_max;j++){ + num1 = GetOffCompGeneralSpin(j-1, org_isite3, org_sigma4, org_sigma3, &tmp_off, X->Def.SiteToBit, X->Def.Tpow); + if(num1 != FALSE){ + num1 = GetOffCompGeneralSpin(tmp_off, org_isite1, org_sigma2, org_sigma1, &tmp_off_2, X->Def.SiteToBit, X->Def.Tpow); + if(num1 != FALSE){ + dam_pr += tmp_V*conj(vec[tmp_off_2+1])*vec_pr[j]; + } + } + + } } } dam_pr = SumMPI_dc(dam_pr); - //fprintf(*_fp," %4ld %4ld %4ld %4ld %4ld %4ld %4ld %4ld %.10lf %.10lf \n",tmp_org_isite1-1, tmp_org_sigma1, tmp_org_isite2-1, tmp_org_sigma2, tmp_org_isite3-1, tmp_org_sigma3, tmp_org_isite4-1, tmp_org_sigma4, creal(dam_pr),cimag(dam_pr)); - fprintf(*_fp," %4ld %4ld %4ld %4ld %4ld %4ld %4ld %4ld %4ld %4ld %4ld %4ld %.10lf %.10lf \n", + fprintf(*_fp, " %4ld %4ld %4ld %4ld " + " %4ld %4ld %4ld %4ld " + " %4ld %4ld %4ld %4ld " + " %.10lf %.10lf \n", tmp_org_isite1-1, tmp_org_sigma1, tmp_org_isite2-1, tmp_org_sigma2, - tmp_org_isite3-1, tmp_org_sigma3, tmp_org_isite4-1, tmp_org_sigma4, - tmp_org_isite5-1, tmp_org_sigma5, tmp_org_isite6-1, tmp_org_sigma6, + tmp_org_isite3-1, tmp_org_sigma3, tmp_org_isite4-1, tmp_org_sigma4, + tmp_org_isite5-1, tmp_org_sigma5, tmp_org_isite6-1, tmp_org_sigma6, creal(dam_pr),cimag(dam_pr)); - free_cd_1d_allocate(vec_pr); } return 0; } @@ -1246,6 +1261,9 @@ int expec_cisajscktalt_SpinGC(struct BindStruct *X,double complex *vec, FILE **_ } } else { info=expec_cisajscktalt_SpinGCGeneral(X,vec, _fp); + if(X->Def.NTBody>0){ + info = expec_Threebody_SpinGeneral(X,vec,_fp_2); + } } return info; } @@ -1688,8 +1706,7 @@ int expec_cisajscktalt_SpinGCGeneral(struct BindStruct *X,double complex *vec, F i_max=X->Check.idim_max; X->Large.mode=M_CORR; - - for(i=0;iDef.NCisAjtCkuAlvDC;i++){ + for(i=0;iDef.NCisAjtCkuAlvDC;i++){ tmp_org_isite1 = X->Def.CisAjtCkuAlvDC[i][0]+1; tmp_org_sigma1 = X->Def.CisAjtCkuAlvDC[i][1]; tmp_org_isite2 = X->Def.CisAjtCkuAlvDC[i][2]+1; @@ -1700,13 +1717,13 @@ int expec_cisajscktalt_SpinGCGeneral(struct BindStruct *X,double complex *vec, F tmp_org_sigma4 = X->Def.CisAjtCkuAlvDC[i][7]; dam_pr = 0.0; - if(Rearray_Interactions(i, &org_isite1, &org_isite2, &org_isite3, &org_isite4, &org_sigma1, &org_sigma2, &org_sigma3, &org_sigma4, &tmp_V, X,2)!=0){ - //error message will be added - fprintf(*_fp," %4ld %4ld %4ld %4ld %4ld %4ld %4ld %4ld %.10lf %.10lf \n",tmp_org_isite1-1, tmp_org_sigma1, tmp_org_isite2-1, tmp_org_sigma2, tmp_org_isite3-1,tmp_org_sigma3, tmp_org_isite4-1, tmp_org_sigma4,0.0,0.0); - continue; - } + if(Rearray_Interactions(i, &org_isite1, &org_isite2, &org_isite3, &org_isite4, &org_sigma1, &org_sigma2, &org_sigma3, &org_sigma4, &tmp_V, X,2)!=0){ + //error message will be added + fprintf(*_fp," %4ld %4ld %4ld %4ld %4ld %4ld %4ld %4ld %.10lf %.10lf \n",tmp_org_isite1-1, tmp_org_sigma1, tmp_org_isite2-1, tmp_org_sigma2, tmp_org_isite3-1,tmp_org_sigma3, tmp_org_isite4-1, tmp_org_sigma4,0.0,0.0); + continue; + } - if(org_isite1 > X->Def.Nsite && org_isite3 > X->Def.Nsite){ + if(org_isite1 > X->Def.Nsite && org_isite3 > X->Def.Nsite){ if(org_sigma1==org_sigma2 && org_sigma3==org_sigma4 ){ //diagonal dam_pr=child_GC_CisAisCjuAju_GeneralSpin_MPIdouble(org_isite1-1, org_sigma1, org_isite3-1, org_sigma3, tmp_V, X, vec, vec); } @@ -1736,7 +1753,7 @@ int expec_cisajscktalt_SpinGCGeneral(struct BindStruct *X,double complex *vec, F } else{ if(org_sigma1==org_sigma2 && org_sigma3==org_sigma4 ){ //diagonal -#pragma omp parallel for default(none) reduction(+:dam_pr) private(j, num1) firstprivate(i_max,X,org_isite1, org_sigma1,org_isite3, org_sigma3, tmp_V) shared(vec) + #pragma omp parallel for default(none) reduction(+:dam_pr) private(j, num1) firstprivate(i_max,X,org_isite1, org_sigma1,org_isite3, org_sigma3, tmp_V) shared(vec) for(j=1;j<=i_max;j++){ num1=BitCheckGeneral(j-1, org_isite1, org_sigma1, X->Def.SiteToBit, X->Def.Tpow); if(num1 != FALSE){ @@ -1747,7 +1764,7 @@ int expec_cisajscktalt_SpinGCGeneral(struct BindStruct *X,double complex *vec, F } } }else if(org_sigma1 == org_sigma2 && org_sigma3 != org_sigma4){ -#pragma omp parallel for default(none) reduction(+:dam_pr) private(j, num1) firstprivate(i_max,X, org_isite1, org_isite3, org_sigma1,org_sigma3,org_sigma4, tmp_off, tmp_V) shared(vec) + #pragma omp parallel for default(none) reduction(+:dam_pr) private(j, num1) firstprivate(i_max,X, org_isite1, org_isite3, org_sigma1,org_sigma3,org_sigma4, tmp_off, tmp_V) shared(vec) for(j=1;j<=i_max;j++){ num1 = GetOffCompGeneralSpin(j-1, org_isite3, org_sigma4, org_sigma3, &tmp_off, X->Def.SiteToBit, X->Def.Tpow); if(num1 != FALSE){ @@ -1758,7 +1775,7 @@ int expec_cisajscktalt_SpinGCGeneral(struct BindStruct *X,double complex *vec, F } } }else if(org_sigma1 != org_sigma2 && org_sigma3 == org_sigma4){ -#pragma omp parallel for default(none) reduction(+:dam_pr) private(j, num1) firstprivate(i_max,X, org_isite1, org_isite3, org_sigma1,org_sigma2, org_sigma3, tmp_off, tmp_V) shared(vec) + #pragma omp parallel for default(none) reduction(+:dam_pr) private(j, num1) firstprivate(i_max,X, org_isite1, org_isite3, org_sigma1,org_sigma2, org_sigma3, tmp_off, tmp_V) shared(vec) for(j=1;j<=i_max;j++){ num1 = BitCheckGeneral(j-1, org_isite3, org_sigma3, X->Def.SiteToBit, X->Def.Tpow); if(num1 != FALSE){ @@ -1769,7 +1786,7 @@ int expec_cisajscktalt_SpinGCGeneral(struct BindStruct *X,double complex *vec, F } } }else if(org_sigma1 != org_sigma2 && org_sigma3 != org_sigma4){ -#pragma omp parallel for default(none) reduction(+:dam_pr) private(j, num1) firstprivate(i_max,X, org_isite1, org_isite3, org_sigma1, org_sigma2, org_sigma3, org_sigma4, tmp_off, tmp_off_2, tmp_V) shared(vec) + #pragma omp parallel for default(none) reduction(+:dam_pr) private(j, num1) firstprivate(i_max,X, org_isite1, org_isite3, org_sigma1, org_sigma2, org_sigma3, org_sigma4, tmp_off, tmp_off_2, tmp_V) shared(vec) for(j=1;j<=i_max;j++){ num1 = GetOffCompGeneralSpin(j-1, org_isite3, org_sigma4, org_sigma3, &tmp_off, X->Def.SiteToBit, X->Def.Tpow); if(num1 != FALSE){ diff --git a/src/mltplySpin.c b/src/mltplySpin.c index d5442c16..65f018f8 100644 --- a/src/mltplySpin.c +++ b/src/mltplySpin.c @@ -466,7 +466,6 @@ void mltplyGeneralSpinGC_mini( isite2 = site_j + 1; sigma1 = spin_i; sigma2 = spin_j; - tmp_trans = 1.0; // transverse magnetic field //dam_pr = 0.0; #pragma omp parallel for default(none) reduction(+:dam_pr) \ @@ -484,7 +483,17 @@ void mltplyGeneralSpinGC_mini( } }// sigma1 != sigma2 else{ // sigma1 = sigma2 - fprintf(stderr, "Error: Transverse_Diagonal component must be absorbed !"); + if (isite1 > X->Def.Nsite) { + dam_pr = child_GC_CisAis_GeneralSpin_MPIdouble(isite1 - 1, sigma1, tmp_trans, X, tmp_v0, tmp_v1); + }else{ + // longitudinal magnetic field + #pragma omp parallel for default(none) private(j, num1) firstprivate(i_max, isite1, sigma1, X, tmp_trans) shared(tmp_v0, tmp_v1) + for (j = 1; j <= i_max; j++) { + num1 = BitCheckGeneral(j - 1, isite1, sigma1, X->Def.SiteToBit, X->Def.Tpow); + tmp_v0[j] += tmp_trans * tmp_v1[j] * num1; + } + } + //fprintf(stderr, "Error: Transverse_Diagonal component must be absorbed !"); } }//isite1 = isite2 //else { // isite1 != isite2