diff --git a/src/expec_cisajscktaltdc.c b/src/expec_cisajscktaltdc.c
index 2b760349..4cf310eb 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;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];
+
+        /*[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 96f939d8..e80c7718 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 3937fbd6..ef3ef69e 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 e8b52e10..d5442c16 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