From d8baf5fb0e80c5dc7f5d652ebeb4c3f68e7f5eb4 Mon Sep 17 00:00:00 2001 From: tmisawa Date: Thu, 2 May 2024 16:25:48 +0900 Subject: [PATCH] add three body G for general spin --- src/expec_cisajscktaltdc.c | 142 +++++++++++++++++++++++++++++ src/include/expec_cisajscktaltdc.h | 1 + src/include/mltplySpin.h | 10 ++ src/mltplySpin.c | 87 ++++++++++++++++++ 4 files changed, 240 insertions(+) diff --git a/src/expec_cisajscktaltdc.c b/src/expec_cisajscktaltdc.c index 2b7603494..4cf310eb8 100644 --- a/src/expec_cisajscktaltdc.c +++ b/src/expec_cisajscktaltdc.c @@ -771,6 +771,148 @@ int expec_cisajscktalt_SpinHalf(struct BindStruct *X,double complex *vec, FILE * return 0; } +/** + * @brief Child function to calculate three-body green functions for general Spin GC model + * + * @param X [in] data list for calculation + * @param vec [in] eigenvectors + * @param _fp [in] output file name + * @retval 0 normally finished + * @retval -1 abnormally finished + * + */ + +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; + long unsigned int tmp_org_isite1,tmp_org_isite2,tmp_org_isite3,tmp_org_isite4,tmp_org_isite5,tmp_org_isite6; + long unsigned int tmp_org_sigma1,tmp_org_sigma2,tmp_org_sigma3,tmp_org_sigma4,tmp_org_sigma5,tmp_org_sigma6; + long unsigned int tmp_off=0; + long unsigned int tmp_off_2=0; + long unsigned int list1_off=0; + int num1; + double complex tmp_V; + double complex dam_pr; + long int i_max; + int tmp_Sz; + long unsigned int tmp_org=0; + double complex *vec_pr; + vec[0]=0; + i_max=X->Check.idim_max; + X->Large.mode=M_CORR; + + 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]; + + /*[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]; + + X->Large.mode = M_MLTPLY; + /* |vec_pr>= c5a6|phi>*/ + mltplyGeneralSpinGC_mini(X,tmp_org_isite5-1,tmp_org_sigma5,tmp_org_isite6-1,tmp_org_sigma6,vec_pr,vec); + 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){ + 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_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); + } + 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{ + dam_pr=0.0; + } + } + 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); + } + 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{ + dam_pr=0.0; + } + } + 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) + for(j=1;j<=i_max;j++){ + num1=BitCheckGeneral(list_1[j], 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); + if(num1 != FALSE){ + dam_pr += tmp_V*conj(vec[j])*vec[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) + 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]; + } + } + } + //printf("DEBUG: rank=%d, dam_pr=%lf\n", myrank, creal(dam_pr)); + } + else{ + dam_pr=0.0; + } + } + 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", + 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, + creal(dam_pr),cimag(dam_pr)); + free_cd_1d_allocate(vec_pr); + } + return 0; +} + + + + + /** * @brief Child function to calculate two-body green's functions for General Spin model * diff --git a/src/include/expec_cisajscktaltdc.h b/src/include/expec_cisajscktaltdc.h index 96f939d8b..e80c7718b 100644 --- a/src/include/expec_cisajscktaltdc.h +++ b/src/include/expec_cisajscktaltdc.h @@ -69,5 +69,6 @@ void expec_cisajscktaltdc_alldiag_spin( double complex *vec ); +int expec_Threebody_SpinGeneral(struct BindStruct *X,double complex *vec, FILE **_fp); int expec_Threebody_SpinGCHalf(struct BindStruct *X,double complex *vec, FILE **_fp); int expec_Fourbody_SpinGCHalf(struct BindStruct *X,double complex *vec, FILE **_fp); diff --git a/src/include/mltplySpin.h b/src/include/mltplySpin.h index 3937fbd68..ef3ef69eb 100644 --- a/src/include/mltplySpin.h +++ b/src/include/mltplySpin.h @@ -35,6 +35,16 @@ void mltplyHalfSpinGC_mini( double complex *tmp_v1//!<[in] Input producted vector ); +void mltplyGeneralSpinGC_mini( + struct BindStruct *X,//!<[inout] + int site_i, + int spin_i, + int site_j, + int spin_j, + double complex *tmp_v0,//!<[inout] Result vector + double complex *tmp_v1//!<[in] Input producted vector +); + int mltplySpinGC(struct BindStruct *X, double complex *tmp_v0,double complex *tmp_v1); int mltplyHalfSpinGC(struct BindStruct *X, double complex *tmp_v0,double complex *tmp_v1); diff --git a/src/mltplySpin.c b/src/mltplySpin.c index e8b52e106..d5442c164 100644 --- a/src/mltplySpin.c +++ b/src/mltplySpin.c @@ -408,6 +408,93 @@ int mltplySpinGC( return iret; }/*int mltplySpinGC*/ + +/** +@brief Driver function for General Spin Hamiltonian (grandcanonical) +@return error code +@author Takahiro Misawa (The University of Tokyo) +*/ +void mltplyGeneralSpinGC_mini( + struct BindStruct *X,//!<[inout] + int site_i, + int spin_i, + int site_j, + int spin_j, + double complex *tmp_v0,//!<[inout] Result vector + double complex *tmp_v1//!<[in] Input producted vector +) { + long unsigned int j; + long unsigned int i; + long unsigned int off = 0; + long unsigned int tmp_off = 0; + long unsigned int isite1, isite2, sigma1, sigma2; + long unsigned int sigma3, sigma4; + double complex dam_pr; + double complex tmp_trans; + long int tmp_sgn; + double num1 = 0; + /*[s] For InterAll */ + double complex tmp_V; + double complex dmv=0; + /*[e] For InterAll */ + + long unsigned int i_max; + i_max = X->Check.idim_max; + + int ihermite=0; + int idx=0; + + StartTimer(510); + //for (i = 0; i < X->Def.EDNTransfer; i += 2) { + isite1 = site_i + 1; + isite2 = site_j + 1; + sigma1 = spin_i; + sigma2 = spin_j; + tmp_trans = 1.0; + dam_pr = 0.0; + if (isite1 == isite2) { + if (sigma1 != sigma2) { + if (isite1 > X->Def.Nsite) { + dam_pr = child_GC_CisAit_GeneralSpin_MPIdouble(isite1 - 1, sigma1, sigma2, tmp_trans, X, tmp_v0, tmp_v1); + //X->Large.prdct += dam_pr; + }/*if (isite1 > X->Def.Nsite)*/ + else { + //for (ihermite = 0; ihermite<2; ihermite++) { + idx = i + ihermite; + + isite1 = site_i + 1; + 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) \ + private(j, tmp_sgn, num1) firstprivate(i_max, isite1, sigma1, sigma2, X, off, tmp_trans) \ + shared(tmp_v0, tmp_v1) + for (j = 1; j <= i_max; j++) { + num1 = GetOffCompGeneralSpin(j - 1, isite1, sigma2, sigma1, &off, X->Def.SiteToBit, X->Def.Tpow); + if (num1 != 0) { // for multply + tmp_v0[off + 1] += tmp_v1[j] * tmp_trans; + //dam_pr += conj(tmp_v1[off + 1]) * tmp_v1[j] * tmp_trans; + }/*if (num1 != 0)*/ + }/*for (j = 1; j <= i_max; j++)*/ + //X->Large.prdct += dam_pr; + //}/*for (ihermite = 0; ihermite<2; ihermite++)*/ + } + }// sigma1 != sigma2 + else{ // sigma1 = sigma2 + fprintf(stderr, "Error: Transverse_Diagonal component must be absorbed !"); + } + }//isite1 = isite2 + //else { // isite1 != isite2 + // hopping is not allowed in localized spin system + // return -1; + //} + //}/*for (i = 0; i < X->Def.EDNTransfer; i += 2)*/ + StopTimer(510); +} + /** @brief Driver function for Spin 1/2 Hamiltonian (grandcanonical) @return error code