Skip to content

Commit

Permalink
fixed bugs
Browse files Browse the repository at this point in the history
  • Loading branch information
三澤貴宏 authored and 三澤貴宏 committed May 13, 2024
1 parent d8baf5f commit ee29091
Show file tree
Hide file tree
Showing 2 changed files with 103 additions and 77 deletions.
167 changes: 92 additions & 75 deletions src/expec_cisajscktaltdc.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -804,107 +804,122 @@ int expec_ThreeBody_SpinGeneral(struct BindStruct *X,double complex *vec, FILE *

for(i=0;i<X->Def.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>*/
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){
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;
}
Expand Down Expand Up @@ -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;
}
Expand Down Expand Up @@ -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;i<X->Def.NCisAjtCkuAlvDC;i++){
for(i=0;i<X->Def.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;
Expand All @@ -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);
}
Expand Down Expand Up @@ -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){
Expand All @@ -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){
Expand All @@ -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){
Expand All @@ -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){
Expand Down
13 changes: 11 additions & 2 deletions src/mltplySpin.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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) \
Expand All @@ -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
Expand Down

0 comments on commit ee29091

Please sign in to comment.