mltplyHubbardCore.cpp

Go to the documentation of this file.

/* HPhi  -  Quantum Lattice Model Simulator */
/* Copyright (C) 2015 The University of Tokyo */
/* This program is free software: you can redistribute it and/or modify */
/* it under the terms of the GNU General Public License as published by */
/* the Free Software Foundation, either version 3 of the License, or */
/* (at your option) any later version. */

/* This program is distributed in the hope that it will be useful, */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the */
/* GNU General Public License for more details. */

/* You should have received a copy of the GNU General Public License */
/* along with this program.  If not, see <http://www.gnu.org/licenses/>. */
#include <bitcalc.hpp>
#include "xsetmem.hpp"
#include "wrapperMPI.hpp"
#include "mltplyCommon.hpp"
#include "mltplyHubbardCore.hpp"

/******************************************************************************/
//[s] GetInfo functions
/******************************************************************************/

int child_general_hopp_GetInfo(
  struct BindStruct *X,
  long int isite1,
  long int isite2,
  long int sigma1,
  long int sigma2
) {
  X->Large.is1_spin = X->Def.Tpow[2 * isite1 - 2 + sigma1];
  X->Large.is2_spin = X->Def.Tpow[2 * isite2 - 2 + sigma2];
  if (isite1 > isite2) {
    X->Large.A_spin = (X->Def.Tpow[2 * isite1 - 2 + sigma1] - X->Def.Tpow[2 * isite2 - 1 + sigma2]);
  }
  else if (isite1 < isite2) {
    X->Large.A_spin = (X->Def.Tpow[2 * isite2 - 2 + sigma2] - X->Def.Tpow[2 * isite1 - 1 + sigma1]);
  }
  else {
    if (sigma1 > sigma2) {
      X->Large.A_spin = (X->Def.Tpow[2 * isite1 - 2 + sigma1] - X->Def.Tpow[2 * isite2 - 1 + sigma2]);
    }
    else {
      X->Large.A_spin = (X->Def.Tpow[2 * isite2 - 2 + sigma2] - X->Def.Tpow[2 * isite1 - 1 + sigma1]);
    }
  }
  X->Large.isA_spin = X->Large.is1_spin + X->Large.is2_spin;
  return 0;
}/*int child_general_hopp_GetInfo*/
int child_general_int_GetInfo(
  struct BindStruct *X,
  long int isite1,
  long int isite2,
  long int isite3,
  long int isite4,
  long int sigma1,
  long int sigma2,
  long int sigma3,
  long int sigma4,
  std::complex<double> tmp_V
) {
  long int is1_spin, is2_spin, is3_spin, is4_spin;
  long int A_spin, B_spin;
  long int isA_spin, isB_spin;
  is1_spin = X->Def.Tpow[2 * isite1 - 2 + sigma1];
  is2_spin = X->Def.Tpow[2 * isite2 - 2 + sigma2];
  if (isite1 > isite2) {
    A_spin = (X->Def.Tpow[2 * isite1 - 2 + sigma1] - X->Def.Tpow[2 * isite2 - 1 + sigma2]);
  }
  else if (isite2 > isite1) {
    A_spin = (X->Def.Tpow[2 * isite2 - 2 + sigma2] - X->Def.Tpow[2 * isite1 - 1 + sigma1]);
  }
  else {//isite1=isite2
    if (sigma1 > sigma2) {
      A_spin = (X->Def.Tpow[2 * isite1 - 2 + sigma1] - X->Def.Tpow[2 * isite2 - 1 + sigma2]);
    }
    else {
      A_spin = (X->Def.Tpow[2 * isite2 - 2 + sigma2] - X->Def.Tpow[2 * isite1 - 1 + sigma1]);
    }
  }
  is3_spin = X->Def.Tpow[2 * isite3 - 2 + sigma3];
  is4_spin = X->Def.Tpow[2 * isite4 - 2 + sigma4];
  if (isite3 > isite4) {
    B_spin = (X->Def.Tpow[2 * isite3 - 2 + sigma3] - X->Def.Tpow[2 * isite4 - 1 + sigma4]);
  }
  else if (isite3 < isite4) {
    B_spin = (X->Def.Tpow[2 * isite4 - 2 + sigma4] - X->Def.Tpow[2 * isite3 - 1 + sigma3]);
  }
  else {//isite3=isite4
    if (sigma3 > sigma4) {
      B_spin = (X->Def.Tpow[2 * isite3 - 2 + sigma3] - X->Def.Tpow[2 * isite4 - 1 + sigma4]);
    }
    else {
      B_spin = (X->Def.Tpow[2 * isite4 - 2 + sigma4] - X->Def.Tpow[2 * isite3 - 1 + sigma3]);
    }
  }
  isA_spin = is1_spin + is2_spin;
  isB_spin = is3_spin + is4_spin;

  X->Large.is1_spin = is1_spin;
  X->Large.is2_spin = is2_spin;
  X->Large.is3_spin = is3_spin;
  X->Large.is4_spin = is4_spin;
  X->Large.isA_spin = isA_spin;
  X->Large.isB_spin = isB_spin;
  X->Large.A_spin = A_spin;
  X->Large.B_spin = B_spin;
  X->Large.tmp_V = tmp_V;
  X->Large.isite1 = isite1;
  X->Large.isite2 = isite2;
  X->Large.isite3 = isite3;
  X->Large.isite4 = isite4;

  return 0;
}/*int child_general_int_GetInfo*/
int child_pairhopp_GetInfo(
  int iPairHopp,
  struct BindStruct *X
) {
  int isite1 = X->Def.PairHopping[iPairHopp][0] + 1;
  int isite2 = X->Def.PairHopping[iPairHopp][1] + 1;
  X->Large.tmp_J = X->Def.ParaPairHopping[iPairHopp];
  X->Large.is1_up = X->Def.Tpow[2 * isite1 - 2];
  X->Large.is1_down = X->Def.Tpow[2 * isite1 - 1];
  X->Large.is2_up = X->Def.Tpow[2 * isite2 - 2];
  X->Large.is2_down = X->Def.Tpow[2 * isite2 - 1];

  return 0;
}/*int child_pairhopp_GetInfo*/
int child_exchange_GetInfo(
  int iExchange,
  struct BindStruct *X
) {
  int isite1 = X->Def.ExchangeCoupling[iExchange][0] + 1;
  int isite2 = X->Def.ExchangeCoupling[iExchange][1] + 1;
  X->Large.tmp_J = -X->Def.ParaExchangeCoupling[iExchange];
  X->Large.is1_up = X->Def.Tpow[2 * isite1 - 2];
  X->Large.is1_down = X->Def.Tpow[2 * isite1 - 1];
  X->Large.is2_up = X->Def.Tpow[2 * isite2 - 2];
  X->Large.is2_down = X->Def.Tpow[2 * isite2 - 1];

  return 0;
}/*int child_exchange_GetInfo*/

/******************************************************************************/
//[e] GetInfo functions
/******************************************************************************/

/******************************************************************************/
//[s] core routines
/******************************************************************************/
void GC_CisAis(
  long int j,
  int nstate, 
  std::complex<double> **tmp_v0,
  std::complex<double> **tmp_v1,
  long int is1_spin,
  std::complex<double> tmp_trans
) {
  long int A_ibit_tmp;
  long int list_1_j;
  std::complex<double> dmv;
  int one = 1;

  list_1_j = j - 1;
  A_ibit_tmp = (list_1_j & is1_spin) / is1_spin;
  dmv = tmp_trans * (std::complex<double>)A_ibit_tmp;
  zaxpy_(&nstate, &dmv, &tmp_v1[j][0], &one, &tmp_v0[j][0], &one);
}/*std::complex<double> GC_CisAis*/
void GC_AisCis(
  long int j,
  int nstate, std::complex<double> **tmp_v0,
  std::complex<double> **tmp_v1,
  long int is1_spin,
  std::complex<double> tmp_trans
) {
  long int A_ibit_tmp;
  long int list_1_j;
  std::complex<double> dmv;
  int one = 1;

  list_1_j = j - 1;
  A_ibit_tmp = (list_1_j & is1_spin) / is1_spin;
  dmv = tmp_trans * (std::complex<double>)(1 - A_ibit_tmp);
  zaxpy_(&nstate, &dmv, &tmp_v1[j][0], &one, &tmp_v0[j][0], &one);
}/*std::complex<double> GC_AisCis*/
int X_CisAis(
  long int list_1_j,
  long int is1_spin
) {
  int A_ibit_tmp;

  // off = j
  A_ibit_tmp = (list_1_j & is1_spin) / is1_spin;
  return A_ibit_tmp;
}
void CisAjt(
  long int j,
  int nstate, std::complex<double> **tmp_v0,
  std::complex<double> **tmp_v1,
  struct BindStruct *X,
  long int is1_spin,
  long int is2_spin,
  long int sum_spin,
  long int diff_spin,
  std::complex<double> tmp_V
) {
  long int ibit_tmp_1, ibit_tmp_2;
  long int bit, iexchg, off;
  int sgn;
  std::complex<double> dmv;
  int one = 1;

  ibit_tmp_1 = (list_1[j] & is1_spin);
  ibit_tmp_2 = (list_1[j] & is2_spin);
  if (ibit_tmp_1 == 0 && ibit_tmp_2 != 0) {
    bit = list_1[j] & diff_spin;
    SgnBit(bit, &sgn); // Fermion sign
    iexchg = list_1[j] ^ sum_spin;

    if(GetOffComp(list_2_1, list_2_2, iexchg, X->Large.irght, X->Large.ilft, X->Large.ihfbit, &off)==FALSE){
      return;
    }
    dmv = (std::complex<double>)sgn * tmp_V;
    zaxpy_(&nstate, &dmv, &tmp_v1[j][0], &one, &tmp_v0[off][0], &one);
  }
  else {
    return;
  }
}
void GC_CisAjt(
  long int j,
  int nstate, std::complex<double> **tmp_v0,
  std::complex<double> **tmp_v1,
  long int is1_spin,
  long int is2_spin,
  long int sum_spin,
  long int diff_spin,
  std::complex<double> tmp_V,
  long int *tmp_off
) {
  long int list_1_j, list_1_off;
  long int ibit_tmp_1, ibit_tmp_2;
  long int bit;
  int sgn;
  std::complex<double> dmv;

  list_1_j = j - 1;
  ibit_tmp_1 = (list_1_j & is1_spin);
  ibit_tmp_2 = (list_1_j & is2_spin);
  *tmp_off = 0;
  int one = 1;

  if (ibit_tmp_1 == 0 && ibit_tmp_2 != 0) {
    bit = list_1_j & diff_spin;
    SgnBit(bit, &sgn); // Fermion sign
    list_1_off = list_1_j ^ sum_spin;
    *tmp_off = list_1_off;
    dmv = (std::complex<double>)sgn * tmp_V;
    zaxpy_(&nstate, &dmv, &tmp_v1[j][0], &one, &tmp_v0[list_1_off + 1][0], &one);
  }
  else {
    return;
  }
}/*std::complex<double> GC_CisAjt*/
int X_CisAjt(
  long int list_1_j,
  struct BindStruct *X,
  long int is1_spin,
  long int is2_spin,
  long int sum_spin,
  long int diff_spin,
  long int *tmp_off
) {
  long int off;
  int sgn = 1;

  sgn = X_GC_CisAjt(list_1_j, is1_spin, is2_spin, sum_spin, diff_spin, tmp_off);
  if (sgn != 0) {
    if(GetOffComp(list_2_1, list_2_2, *tmp_off, X->Large.irght, X->Large.ilft, X->Large.ihfbit, &off)!=TRUE){
      *tmp_off = 0;
      return 0;
    }
    *tmp_off = off;
    return sgn;
  }
  else {
    *tmp_off = 0;
    return 0;
  }
}/*int X_CisAjt*/
int X_GC_CisAjt(
  long int list_1_j,
  long int is1_spin,
  long int is2_spin,
  long int sum_spin,
  long int diff_spin,
  long int *tmp_off
) {
  long int ibit_tmp_1, ibit_tmp_2;
  long int bit, off;
  int sgn = 1;

  ibit_tmp_1 = (list_1_j & is1_spin);
  ibit_tmp_2 = (list_1_j & is2_spin);

  if (ibit_tmp_1 == 0 && ibit_tmp_2 != 0) {
    bit = list_1_j & diff_spin;
    SgnBit(bit, &sgn); // Fermion sign
    off = list_1_j ^ sum_spin;
    *tmp_off = off;
    return sgn; // pm 1
  }
  else {
    *tmp_off = 0;
    return 0;
  }
}
/******************************************************************************/
//[e] core routines
/******************************************************************************/

/******************************************************************************/
//[s] child element functions
/******************************************************************************/
void child_exchange_element(
  long int j,
  int nstate, std::complex<double> **tmp_v0,
  std::complex<double> **tmp_v1,
  struct BindStruct *X,
  long int *tmp_off
) {
  long int off;
  long int ibit1_up, ibit2_up, ibit1_down, ibit2_down;
  std::complex<double> dmv;
  long int iexchg;
  long int is1_up = X->Large.is1_up;
  long int is2_up = X->Large.is2_up;
  long int is1_down = X->Large.is1_down;
  long int is2_down = X->Large.is2_down;
  long int irght = X->Large.irght;
  long int ilft = X->Large.ilft;
  long int ihfbit = X->Large.ihfbit;
  std::complex<double> tmp_J = X->Large.tmp_J;
  int one = 1;

  ibit1_up = list_1[j] & is1_up;
  ibit2_up = list_1[j] & is2_up;
  ibit1_down = list_1[j] & is1_down;
  ibit2_down = list_1[j] & is2_down;

  if (ibit1_up == 0 && ibit1_down != 0 && ibit2_up != 0 && ibit2_down == 0) {
    iexchg = list_1[j] - (is1_down + is2_up);
    iexchg += (is1_up + is2_down);
    if(GetOffComp(list_2_1, list_2_2, iexchg, irght, ilft, ihfbit, &off)!=TRUE){
      return;
    }
    *tmp_off = off;
    dmv = tmp_J;
    zaxpy_(&nstate, &dmv, &tmp_v1[j][0], &one, &tmp_v0[off][0], &one);
  }
  else if (ibit1_up != 0 && ibit1_down == 0 && ibit2_up == 0 && ibit2_down != 0) {
    iexchg = list_1[j] - (is1_up + is2_down);
    iexchg += (is1_down + is2_up);
    if(GetOffComp(list_2_1, list_2_2, iexchg, irght, ilft, ihfbit, &off)!=TRUE){
      return;
    }
    *tmp_off = off;
    dmv = tmp_J;
    zaxpy_(&nstate, &dmv, &tmp_v1[j][0], &one, &tmp_v0[off][0], &one);
  }
}/*std::complex<double> child_exchange_element*/
void child_pairhopp_element(
  long int j,
  int nstate, std::complex<double> **tmp_v0,
  std::complex<double> **tmp_v1,
  struct BindStruct *X,
  long int *tmp_off
) {
  long int off;
  long int ibit1_up, ibit2_up, ibit1_down, ibit2_down;
  std::complex<double> dmv;
  long int iexchg;
  long int is1_up = X->Large.is1_up;
  long int is2_up = X->Large.is2_up;
  long int is1_down = X->Large.is1_down;
  long int is2_down = X->Large.is2_down;
  long int irght = X->Large.irght;
  long int ilft = X->Large.ilft;
  long int ihfbit = X->Large.ihfbit;
  std::complex<double> tmp_J = X->Large.tmp_J;
  int one = 1;

  ibit1_up = list_1[j] & is1_up;
  ibit2_up = list_1[j] & is2_up;
  ibit1_down = list_1[j] & is1_down;
  ibit2_down = list_1[j] & is2_down;

  if (ibit1_up == 0 && ibit1_down == 0 && ibit2_up != 0 && ibit2_down != 0) {
    iexchg = list_1[j] - (is2_up + is2_down);
    iexchg += (is1_up + is1_down);

    if(GetOffComp(list_2_1, list_2_2, iexchg, irght, ilft, ihfbit, &off)!=TRUE){
      return;
    }
    *tmp_off = off;
    dmv = tmp_J;
    zaxpy_(&nstate, &dmv, &tmp_v1[j][0], &one, &tmp_v0[off][0], &one);
  }
}/*std::complex<double> child_pairhopp_element*/
void GC_child_exchange_element(
  long int j,
  int nstate, std::complex<double> **tmp_v0,
  std::complex<double> **tmp_v1,
  struct BindStruct *X,
  long int *tmp_off
) {
  long int ibit1_up, ibit2_up, ibit1_down, ibit2_down;
  std::complex<double> dmv;
  long int iexchg;
  long int is1_up = X->Large.is1_up;
  long int is2_up = X->Large.is2_up;
  long int is1_down = X->Large.is1_down;
  long int is2_down = X->Large.is2_down;
  long int list_1_j, list_1_off;
  std::complex<double> tmp_J = X->Large.tmp_J;
  int one = 1;

  list_1_j = j - 1;
  ibit1_up = list_1_j & is1_up;
  ibit2_up = list_1_j & is2_up;
  ibit1_down = list_1_j & is1_down;
  ibit2_down = list_1_j & is2_down;

  if (ibit1_up == 0 && ibit1_down != 0 && ibit2_up != 0 && ibit2_down == 0) {

    iexchg = list_1_j - (is1_down + is2_up);
    iexchg += (is1_up + is2_down);
    list_1_off = iexchg;
    *tmp_off = list_1_off;

    dmv = tmp_J;
    zaxpy_(&nstate, &dmv, &tmp_v1[j][0], &one, &tmp_v0[list_1_off + 1][0], &one);
  }
  else if (ibit1_up != 0 && ibit1_down == 0 && ibit2_up == 0 && ibit2_down != 0) {
    iexchg = list_1_j - (is1_up + is2_down);
    iexchg += (is1_down + is2_up);
    list_1_off = iexchg;
    *tmp_off = list_1_off;

    dmv = tmp_J;
    zaxpy_(&nstate, &dmv, &tmp_v1[j][0], &one, &tmp_v0[list_1_off + 1][0], &one);
  }
}/*std::complex<double> GC_child_exchange_element*/
void GC_child_pairhopp_element(
  long int j,
  int nstate, std::complex<double> **tmp_v0,
  std::complex<double> **tmp_v1,
  struct BindStruct *X,
  long int *tmp_off
) {
  long int ibit1_up, ibit2_up, ibit1_down, ibit2_down;
  std::complex<double> dmv;
  long int iexchg;
  long int is1_up = X->Large.is1_up;
  long int is2_up = X->Large.is2_up;
  long int is1_down = X->Large.is1_down;
  long int is2_down = X->Large.is2_down;
  long int list_1_j, list_1_off;
  std::complex<double> tmp_J = X->Large.tmp_J;
  int one = 1;

  list_1_j = j - 1;

  ibit1_up = list_1_j & is1_up;

  ibit2_up = list_1_j & is2_up;

  ibit1_down = list_1_j & is1_down;

  ibit2_down = list_1_j & is2_down;

  if (ibit1_up == 0 && ibit1_down == 0 && ibit2_up != 0 && ibit2_down != 0) {
    iexchg = list_1_j - (is2_up + is2_down);
    iexchg += (is1_up + is1_down);
    list_1_off = iexchg;
    *tmp_off = list_1_off;
    dmv = tmp_J;
    zaxpy_(&nstate, &dmv, &tmp_v1[j][0], &one, &tmp_v0[list_1_off + 1][0], &one);
  }
}
void child_CisAisCisAis_element(
  long int j,
  long int isite1,
  long int isite3,
  std::complex<double> tmp_V,
  int nstate, std::complex<double> **tmp_v0,
  std::complex<double> **tmp_v1
) {
  int tmp_sgn;
  std::complex<double> dmv;
  int one = 1;
  tmp_sgn = X_CisAis(list_1[j], isite3);
  tmp_sgn *= X_CisAis(list_1[j], isite1);
  dmv = tmp_V * (std::complex<double>)tmp_sgn;
  zaxpy_(&nstate, &dmv, &tmp_v1[j][0], &one, &tmp_v0[j][0], &one);
}/*std::complex<double> child_CisAisCisAis_element*/
void child_CisAisCjtAku_element(
  long int j,
  long int isite1,
  long int isite3,
  long int isite4,
  long int Bsum,
  long int Bdiff,
  std::complex<double> tmp_V,
  int nstate, std::complex<double> **tmp_v0,
  std::complex<double> **tmp_v1,
  struct BindStruct *X,
  long int *tmp_off
) {
  int tmp_sgn;
  std::complex<double> dmv;
  int one = 1;
  tmp_sgn = X_CisAjt(list_1[j], X, isite3, isite4, Bsum, Bdiff, tmp_off);
  if (tmp_sgn != 0) {
    tmp_sgn *= X_CisAis(list_1[*tmp_off], isite1);
    if (tmp_sgn != 0) {
      dmv = tmp_V * (std::complex<double>)tmp_sgn;
      zaxpy_(&nstate, &dmv, &tmp_v1[j][0], &one, &tmp_v0[*tmp_off][0], &one);
    }
  }
}/*std::complex<double> child_CisAisCjtAku_element*/
void child_CisAjtCkuAku_element(
  long int j,
  long int isite1,
  long int isite2,
  long int isite3,
  long int Asum,
  long int Adiff,
  std::complex<double> tmp_V,
  int nstate, std::complex<double> **tmp_v0,
  std::complex<double> **tmp_v1,
  struct BindStruct *X,
  long int *tmp_off
) {
  int tmp_sgn;
  std::complex<double> dmv;
  int one = 1;
  tmp_sgn = X_CisAis(list_1[j], isite3);
  if (tmp_sgn != 0) {
    tmp_sgn *= X_CisAjt(list_1[j], X, isite1, isite2, Asum, Adiff, tmp_off);
    if (tmp_sgn != 0) {
      dmv = tmp_V * (std::complex<double>)tmp_sgn;
      zaxpy_(&nstate, &dmv, &tmp_v1[j][0], &one, &tmp_v0[*tmp_off][0], &one);
    }
  }
}/*std::complex<double> child_CisAjtCkuAku_element*/
void child_CisAjtCkuAlv_element(
  long int j,
  long int isite1,
  long int isite2,
  long int isite3,
  long int isite4,
  long int Asum,
  long int Adiff,
  long int Bsum,
  long int Bdiff,
  std::complex<double> tmp_V,
  int nstate, std::complex<double> **tmp_v0,
  std::complex<double> **tmp_v1,
  struct BindStruct *X,
  long int *tmp_off_2
) {
  int tmp_sgn;
  long int tmp_off_1;
  int one = 1;

  std::complex<double> dmv;
  tmp_sgn = X_GC_CisAjt(list_1[j], isite3, isite4, Bsum, Bdiff, &tmp_off_1);

  if (tmp_sgn != 0) {
    tmp_sgn *= X_CisAjt(tmp_off_1, X, isite1, isite2, Asum, Adiff, tmp_off_2);
    if (tmp_sgn != 0) {
      dmv = tmp_V * (std::complex<double>)tmp_sgn;
      zaxpy_(&nstate, &dmv, &tmp_v1[j][0], &one, &tmp_v0[*tmp_off_2][0], &one);
    }
  }
}/*std::complex<double> child_CisAjtCkuAlv_element*/
//[s] Grand Canonical
void GC_child_CisAisCisAis_element(
  long int j,
  long int isite1,
  long int isite3,
  std::complex<double> tmp_V,
  int nstate, std::complex<double> **tmp_v0,
  std::complex<double> **tmp_v1
) {
  int tmp_sgn;
  std::complex<double> dmv;
  int one = 1;
  tmp_sgn = X_CisAis(j - 1, isite3);
  tmp_sgn *= X_CisAis(j - 1, isite1);
  if (tmp_sgn != 0) {
    dmv = tmp_V * (std::complex<double>)tmp_sgn;
    zaxpy_(&nstate, &dmv, &tmp_v1[j][0], &one, &tmp_v0[j][0], &one);
  }
}/*std::complex<double> GC_child_CisAisCisAis_element*/
void GC_child_CisAisCjtAku_element(
  long int j,
  long int isite1,
  long int isite3,
  long int isite4,
  long int Bsum,
  long int Bdiff,
  std::complex<double> tmp_V,
  int nstate, std::complex<double> **tmp_v0,
  std::complex<double> **tmp_v1,
  long int *tmp_off
) {
  int tmp_sgn;
  std::complex<double> dmv;
  int one = 1;
  tmp_sgn = X_GC_CisAjt((j - 1), isite3, isite4, Bsum, Bdiff, tmp_off);
  if (tmp_sgn != 0) {
    tmp_sgn *= X_CisAis(*tmp_off, isite1);
    if (tmp_sgn != 0) {
      dmv = tmp_V * (std::complex<double>)tmp_sgn;
      zaxpy_(&nstate, &dmv, &tmp_v1[j][0], &one, &tmp_v0[*tmp_off + 1][0], &one);
    }
  }
}/*std::complex<double> GC_child_CisAisCjtAku_element*/
void GC_child_CisAjtCkuAku_element(
  long int j,
  long int isite1,
  long int isite2,
  long int isite3,
  long int Asum,
  long int Adiff,
  std::complex<double> tmp_V,
  int nstate, std::complex<double> **tmp_v0,
  std::complex<double> **tmp_v1,
  long int *tmp_off
) {
  int tmp_sgn;
  std::complex<double> dmv;
  int one = 1;
  tmp_sgn = X_CisAis(j - 1, isite3);
  if (tmp_sgn != 0) {
    tmp_sgn *= X_GC_CisAjt(j - 1, isite1, isite2, Asum, Adiff, tmp_off);
    if (tmp_sgn != 0) {
      dmv = tmp_V * (std::complex<double>)tmp_sgn;
      zaxpy_(&nstate, &dmv, &tmp_v1[j][0], &one, &tmp_v0[*tmp_off + 1][0], &one);
    }/*if (tmp_sgn != 0)*/
  }/*if (tmp_sgn != 0)*/
}/*std::complex<double> GC_child_CisAjtCkuAku_element*/
void GC_child_CisAjtCkuAlv_element(
  long int j,
  long int isite1,
  long int isite2,
  long int isite3,
  long int isite4,
  long int Asum,
  long int Adiff,
  long int Bsum,
  long int Bdiff,
  std::complex<double> tmp_V,
  int nstate, std::complex<double> **tmp_v0,
  std::complex<double> **tmp_v1,
  long int *tmp_off_2
) {
  int tmp_sgn;
  long int tmp_off_1;
  std::complex<double> dmv;
  int one = 1;

  tmp_sgn = X_GC_CisAjt(j - 1, isite3, isite4, Bsum, Bdiff, &tmp_off_1);
  if (tmp_sgn != 0) {
    tmp_sgn *= X_GC_CisAjt(tmp_off_1, isite1, isite2, Asum, Adiff, tmp_off_2);
    if (tmp_sgn != 0) {
      dmv = tmp_V * (std::complex<double>)tmp_sgn;
      zaxpy_(&nstate, &dmv, &tmp_v1[j][0], &one, &tmp_v0[*tmp_off_2 + 1][0], &one);
    }
  }
}/*std::complex<double> GC_child_CisAjtCkuAlv_element*/
//[e] Grand Canonical
void GC_Cis(
  long int j,
  int nstate, std::complex<double> **tmp_v0,
  std::complex<double> **tmp_v1,
  long int is1_spin,
  std::complex<double> tmp_V,
  long int *tmp_off
) {
  long int list_1_j, list_1_off;
  long int ibit_tmp_1;
  long int bit;
  int sgn, ipsgn;
  std::complex<double> dmv;
  int one = 1;

  list_1_j = j - 1;

  ibit_tmp_1 = (list_1_j & is1_spin);
  // is1_spin >= 1
  // is1_spin = Tpow[2*isite + ispin]

  *tmp_off = 0;

  if (ibit_tmp_1 == 0) {
    // able to create an electron at the is1_spin state
    bit = list_1_j - (list_1_j & (2 * is1_spin - 1));
    SgnBit(bit, &sgn); // Fermion sign
    ipsgn = 1;
#ifdef __MPI
    SgnBit(myrank, &ipsgn); // Fermion sign
#endif
    list_1_off = list_1_j | is1_spin; // OR
    *tmp_off = list_1_off;
    dmv = (std::complex<double>)(ipsgn * sgn) * tmp_V;
    zaxpy_(&nstate, &dmv, &tmp_v1[j][0], &one, &tmp_v0[list_1_off + 1][0], &one);
  }
  else {
    return;
  }
}/*std::complex<double> GC_Cis*/
void GC_Ajt(
  long int j,
  int nstate, std::complex<double> **tmp_v0,
  std::complex<double> **tmp_v1,
  long int is1_spin,
  std::complex<double> tmp_V,
  long int *tmp_off
) {
  long int list_1_j, list_1_off;
  long int ibit_tmp_1;
  long int bit;
  int sgn, ipsgn;
  std::complex<double> dmv;
  int one = 1;

  list_1_j = j - 1;

  ibit_tmp_1 = (list_1_j & is1_spin);
  // is1_spin >= 1

  *tmp_off = 0;

  if (ibit_tmp_1 == is1_spin) {
    // able to create an electron at the is1_spin state
    bit = list_1_j - (list_1_j & (2 * is1_spin - 1));
    SgnBit(bit, &sgn); // Fermion sign
    ipsgn = 1;
#ifdef __MPI
    SgnBit(myrank, &ipsgn); // Fermion sign
#endif
    list_1_off = list_1_j ^ is1_spin;
    *tmp_off = list_1_off;
    dmv = (std::complex<double>)(ipsgn * sgn) * tmp_V;
    zaxpy_(&nstate, &dmv, &tmp_v1[j][0], &one, &tmp_v0[list_1_off + 1][0], &one);
  }
  else {
    return;
  }
}/*std::complex<double> GC_Ajt*/
int X_Cis(
  long int j,
  long int is1_spin,
  long int *tmp_off,
  long int *list_1_org,
  long int *list_2_1_target,
  long int *list_2_2_target,
  long int _irght,
  long int _ilft,
  long int _ihfbit
) {
  long int list_1_j, list_1_off;
  long int ibit_tmp_1;
  long int bit;
  int sgn, ipsgn;

  list_1_j = list_1_org[j];

  ibit_tmp_1 = (list_1_j & is1_spin);
  // is1_spin >= 1
  // is1_spin = Tpow[2*isite + ispin]

  *tmp_off = 0;

  if (ibit_tmp_1 == 0) {
    // able to create an electron at the is1_spin state
    bit = list_1_j - (list_1_j & (2 * is1_spin - 1));
    SgnBit(bit, &sgn); // Fermion sign
    ipsgn = 1;
#ifdef __MPI
    SgnBit(myrank, &ipsgn); // Fermion sign
#endif
    list_1_off = list_1_j | is1_spin; // OR

    if(GetOffComp(list_2_1_target, list_2_2_target, list_1_off, _irght, _ilft, _ihfbit, tmp_off)!=TRUE){
      *tmp_off=0;
      return 0;
    }
    sgn *= ipsgn;
    return (sgn);
  }
  else {
    *tmp_off = 0;
    return 0;
  }
}/*int X_Cis*/
int X_Ajt(
  long int j,
  long int is1_spin,
  long int *tmp_off,
  long int *list_1_org,
  long int *list_2_1_target,
  long int *list_2_2_target,
  long int _irght,
  long int _ilft,
  long int _ihfbit
) {
  long int list_1_j, list_1_off;
  long int ibit_tmp_1;
  long int bit;
  int sgn, ipsgn;

  list_1_j = list_1_org[j];

  ibit_tmp_1 = (list_1_j & is1_spin);
// is1_spin >= 1
// is1_spin = Tpow[2*isite + ispin]

  *tmp_off = 0;
  if (ibit_tmp_1 != 0) {
    // able to delete an electron at the is1_spin state
    bit = list_1_j - (list_1_j & (2 * is1_spin - 1));
    SgnBit(bit, &sgn); // Fermion sign
    ipsgn = 1;
#ifdef __MPI
    SgnBit(myrank, &ipsgn); // Fermion sign
#endif
    list_1_off = list_1_j ^ is1_spin;
    if(GetOffComp(list_2_1_target, list_2_2_target, list_1_off, _irght, _ilft, _ihfbit, tmp_off)!=TRUE){
      *tmp_off=0;
      return 0;
    }
    sgn *= ipsgn;
    return(sgn);
  }
  else {
    *tmp_off = 0;
    return 0;
  }
}/*std::complex<double> X_Ajt*/

/******************************************************************************/
//[e] child element functions
/******************************************************************************/