PairExSpin.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 "PairExSpin.hpp"
#include "bitcalc.hpp"
#include "wrapperMPI.hpp"
#include "mltplyMPISpinCore.hpp"
#include "mltplySpinCore.hpp"
#include "mltplyCommon.hpp"
#ifdef __MPI
#include "common/setmemory.hpp"
#endif
int GetPairExcitedStateHalfSpinGC(
  struct BindStruct *X,
  int nstate,
  std::complex<double> **tmp_v0, 
  std::complex<double> **tmp_v1, 
  int iEx
) {
  long int i, j;
  long int isite1;
  long int org_isite1, org_isite2, org_sigma1, org_sigma2;
  long int tmp_off = 0;

  std::complex<double> tmp_trans = 0, dmv;
  long int i_max;
  int tmp_sgn, one = 1;

  i_max = X->Check.idim_maxOrg;

  for (i = 0; i < X->Def.NPairExcitationOperator[iEx]; i++) {
    org_isite1 = X->Def.PairExcitationOperator[iEx][i][0] + 1;
    org_isite2 = X->Def.PairExcitationOperator[iEx][i][2] + 1;
    org_sigma1 = X->Def.PairExcitationOperator[iEx][i][1];
    org_sigma2 = X->Def.PairExcitationOperator[iEx][i][3];
    tmp_trans = X->Def.ParaPairExcitationOperator[iEx][i];
    if (org_isite1 == org_isite2) {
      if (org_isite1 > X->Def.Nsite) {
        if (org_sigma1 == org_sigma2) {  // longitudinal magnetic field
          if (X->Def.PairExcitationOperator[iEx][i][4] == 0) {
            X_GC_child_AisCis_spin_MPIdouble(org_isite1 - 1, org_sigma1, -tmp_trans, X, nstate, tmp_v0, tmp_v1);
          }
          else {
            X_GC_child_CisAis_spin_MPIdouble(org_isite1 - 1, org_sigma1, tmp_trans, X, nstate, tmp_v0, tmp_v1);
          }
        }
        else {  // transverse magnetic field
            //fprintf(stdoutMPI, "Debug: test, org_isite1=%d, org_sigma1=%d, orgsima_2=%d\n", org_isite1, org_sigma1, org_sigma2);
          X_GC_child_CisAit_spin_MPIdouble(org_isite1 - 1, org_sigma1, org_sigma2, tmp_trans, X, nstate, tmp_v0, tmp_v1);
        }
      }
      else {
        isite1 = X->Def.Tpow[org_isite1 - 1];
        if (org_sigma1 == org_sigma2) {
          if (X->Def.PairExcitationOperator[iEx][i][4] == 0) {
            // longitudinal magnetic field
#pragma omp parallel for default(none) private(j, tmp_sgn,dmv) \
  firstprivate(i_max, isite1, org_sigma1, X,tmp_trans) shared(one,nstate,tmp_v0, tmp_v1)
            for (j = 1; j <= i_max; j++) {
              dmv = (1.0 - X_SpinGC_CisAis(j, isite1, org_sigma1))* (-tmp_trans);
              zaxpy_(&nstate, &dmv, tmp_v1[j], &one, tmp_v0[j], &one);
            }
          }
          else {
            // longitudinal magnetic field
#pragma omp parallel for default(none) private(j, tmp_sgn,dmv)             \
  firstprivate(i_max, isite1, org_sigma1, X,tmp_trans) shared(tmp_v0, tmp_v1,one,nstate)
            for (j = 1; j <= i_max; j++) {
              dmv = (std::complex<double>)X_SpinGC_CisAis(j, isite1, org_sigma1)* tmp_trans;
              zaxpy_(&nstate, &dmv, tmp_v1[j], &one, tmp_v0[j], &one);
            }
          }
        }
        else {
          // transverse magnetic field
          // fprintf(stdoutMPI, "Debug: isite1=%d, org_sigma2=%d\n", isite1, org_sigma2);
#pragma omp parallel for default(none) private(j, tmp_sgn, tmp_off,dmv)    \
  firstprivate(i_max, isite1, org_sigma2, X, tmp_trans) shared(tmp_v0, tmp_v1,one,nstate)
          for (j = 1; j <= i_max; j++) {
            tmp_sgn = X_SpinGC_CisAit(j, isite1, org_sigma2, &tmp_off);
            if (tmp_sgn != 0) {
              dmv = (std::complex<double>)tmp_sgn * tmp_trans;
              zaxpy_(&nstate, &dmv, tmp_v1[j], &one, tmp_v0[tmp_off + 1], &one);
            }
          }
        }
      }
    }
    else {
      fprintf(stdoutMPI, "ERROR: hopping is not allowed in localized spin system\n");
      return FALSE;
    }
  }
  return TRUE;
}
int GetPairExcitedStateGeneralSpinGC(
  struct BindStruct *X,
  int nstate,
  std::complex<double> **tmp_v0, 
  std::complex<double> **tmp_v1, 
  int iEx
) {
  long int i, j;
  int num1;
  long int org_isite1, org_isite2, org_sigma1, org_sigma2;
  long int tmp_off = 0;
  int one = 1;
  std::complex<double> tmp_trans = 0, dmv;
  long int i_max;
  i_max = X->Check.idim_maxOrg;

  for (i = 0; i < X->Def.NPairExcitationOperator[iEx]; i++) {
    org_isite1 = X->Def.PairExcitationOperator[iEx][i][0] + 1;
    org_isite2 = X->Def.PairExcitationOperator[iEx][i][2] + 1;
    org_sigma1 = X->Def.PairExcitationOperator[iEx][i][1];
    org_sigma2 = X->Def.PairExcitationOperator[iEx][i][3];
    tmp_trans = X->Def.ParaPairExcitationOperator[iEx][i];
    if (org_isite1 == org_isite2) {
      if (org_isite1 > X->Def.Nsite) {
        if (org_sigma1 == org_sigma2) {
          if (X->Def.PairExcitationOperator[iEx][i][4] == 0) {
            // longitudinal magnetic field
            X_GC_child_AisCis_GeneralSpin_MPIdouble(org_isite1 - 1, org_sigma1, -tmp_trans, X, nstate, tmp_v0, tmp_v1);
          }
          else {
            X_GC_child_CisAis_GeneralSpin_MPIdouble(org_isite1 - 1, org_sigma1, tmp_trans, X, nstate, tmp_v0, tmp_v1);
          }
        }
        else {
          // transverse magnetic field
          X_GC_child_CisAit_GeneralSpin_MPIdouble(org_isite1 - 1, org_sigma1, org_sigma2, tmp_trans, X, nstate, tmp_v0, tmp_v1);
        }
      }
      else {//org_isite1 <= X->Def.Nsite
        if (org_sigma1 == org_sigma2) {
          if (X->Def.PairExcitationOperator[iEx][i][4] == 0) {
            // longitudinal magnetic field
#pragma omp parallel for default(none) private(j,num1,dmv) \
firstprivate(i_max,org_isite1,org_sigma1,X,tmp_trans) shared(tmp_v0,tmp_v1,one,nstate)
            for (j = 1; j <= i_max; j++) {
              num1 = BitCheckGeneral(j - 1, org_isite1, org_sigma1, X->Def.SiteToBit, X->Def.Tpow);
              dmv = -tmp_trans * (1.0 - num1);
              zaxpy_(&nstate, &dmv, tmp_v1[j], &one, tmp_v0[j], &one);
            }
          }
          else {
            // longitudinal magnetic field
#pragma omp parallel for default(none) private(j,num1,dmv) \
  firstprivate(i_max,org_isite1,org_sigma1,X,tmp_trans) shared(tmp_v0,tmp_v1,one,nstate)
            for (j = 1; j <= i_max; j++) {
              num1 = BitCheckGeneral(j - 1, org_isite1, org_sigma1, X->Def.SiteToBit, X->Def.Tpow);
              dmv = tmp_trans * (std::complex<double>)num1;
              zaxpy_(&nstate, &dmv, tmp_v1[j], &one, tmp_v0[j], &one);
            }
          }
        }
        else {
          // transverse magnetic field
#pragma omp parallel for default(none) private(j,num1,dmv) \
firstprivate(i_max,org_isite1,org_sigma1,org_sigma2,X,tmp_off,tmp_trans) \
shared(tmp_v0,tmp_v1,one,nstate)
          for (j = 1; j <= i_max; j++) {
            num1 = GetOffCompGeneralSpin(j - 1, org_isite1, org_sigma2, org_sigma1, &tmp_off, X->Def.SiteToBit, X->Def.Tpow);
            if (num1 != 0) {
              dmv = tmp_trans * (std::complex<double>)num1;
              zaxpy_(&nstate, &dmv, tmp_v1[j], &one, tmp_v0[tmp_off + 1], &one);
            }
          }
        }
      }
    }
    else {
      fprintf(stdoutMPI, "ERROR: hopping is not allowed in localized spin system\n");
      return FALSE;
    }
  }
  return TRUE;
}
int GetPairExcitedStateSpinGC(
  struct BindStruct *X,
  int nstate, std::complex<double> **tmp_v0, 
  std::complex<double> **tmp_v1, 
  int iEx
) {

  int iret = 0;
  if (X->Def.iFlgGeneralSpin == FALSE) {
    iret = GetPairExcitedStateHalfSpinGC(X, nstate, tmp_v0, tmp_v1, iEx);
  }
  else {
    iret = GetPairExcitedStateGeneralSpinGC(X, nstate, tmp_v0, tmp_v1, iEx);
  }
  return iret;
}
int GetPairExcitedStateHalfSpin(
  struct BindStruct *X,
  int nstate, 
  std::complex<double> **tmp_v0, 
  std::complex<double> **tmp_v1, 
  int iEx
)
{
  long int i, j;
  long int isite1;
  long int org_isite1, org_isite2, org_sigma1, org_sigma2;
  long int tmp_off = 0;
  std::complex<double> tmp_trans = 0, dmv;
  long int i_max;
  int num1, one = 1;
  long int ibit1;
  long int is1_up;

  i_max = X->Check.idim_maxOrg;

  for (i = 0; i < X->Def.NPairExcitationOperator[iEx]; i++) {
    org_isite1 = X->Def.PairExcitationOperator[iEx][i][0] + 1;
    org_isite2 = X->Def.PairExcitationOperator[iEx][i][2] + 1;
    org_sigma1 = X->Def.PairExcitationOperator[iEx][i][1];
    org_sigma2 = X->Def.PairExcitationOperator[iEx][i][3];
    tmp_trans = X->Def.ParaPairExcitationOperator[iEx][i];
    if (org_sigma1 == org_sigma2) {
      if (org_isite1 == org_isite2) {
        if (org_isite1 > X->Def.Nsite) {
          is1_up = X->Def.Tpow[org_isite1 - 1];
          ibit1 = X_SpinGC_CisAis((long int) myrank + 1, is1_up, org_sigma1);
          if (X->Def.PairExcitationOperator[iEx][i][4] == 0) {
            if (ibit1 == 0) {
              dmv = -tmp_trans;
#pragma omp parallel for default(none) shared(tmp_v0, tmp_v1,one,nstate,dmv) \
  firstprivate(i_max, tmp_trans) private(j)
              for (j = 1; j <= i_max; j++) {
                zaxpy_(&nstate, &dmv, tmp_v1[j], &one, tmp_v0[j], &one);
              }
            }
          }
          else {
            if (ibit1 != 0) {
#pragma omp parallel for default(none) shared(tmp_v0, tmp_v1,one,nstate)    \
  firstprivate(i_max, tmp_trans) private(j)
              for (j = 1; j <= i_max; j++) 
                zaxpy_(&nstate, &tmp_trans, tmp_v1[j], &one, tmp_v0[j], &one);
            }
          }
        }// org_isite1 > X->Def.Nsite
        else {
          isite1 = X->Def.Tpow[org_isite1 - 1];
          if (org_isite1 == org_isite2 && org_sigma1 == org_sigma2 &&
            X->Def.PairExcitationOperator[iEx][i][4] == 0) {
#pragma omp parallel for default(none) private(j,dmv) \
firstprivate(i_max,isite1,org_sigma1,X,tmp_trans) shared(tmp_v0,tmp_v1,one,nstate)
            for (j = 1; j <= i_max; j++) {
              dmv = (1.0 - X_Spin_CisAis(j, isite1, org_sigma1)) * (-tmp_trans);
              zaxpy_(&nstate, &dmv, tmp_v1[j], &one, tmp_v0[j], &one);
            }
          }
          else {
#pragma omp parallel for default(none) private(j,dmv) \
firstprivate(i_max,isite1,org_sigma1,X,tmp_trans) shared(tmp_v0,tmp_v1,one,nstate)
            for (j = 1; j <= i_max; j++) {
              dmv = (std::complex<double>)X_Spin_CisAis(j, isite1, org_sigma1) * tmp_trans;
              zaxpy_(&nstate, &dmv, tmp_v1[j], &one, tmp_v0[j], &one);
            }
          }
        }
      }
      else {
        fprintf(stdoutMPI, "Error: isite1 must be equal to isite2 for Spin system. \n");
        return FALSE;
      }
    }
    else { //org_sigma1 != org_sigma2             // for the canonical case
      if (org_isite1 > X->Def.Nsite) {//For MPI
        X_child_CisAit_spin_MPIdouble(org_isite1 - 1, org_sigma2, tmp_trans, 
          X, nstate, tmp_v0, tmp_v1, i_max);
      }
      else {
        isite1 = X->Def.Tpow[org_isite1 - 1];
#pragma omp parallel for default(none) private(j,tmp_off,num1,dmv) \
firstprivate(i_max,isite1,org_sigma2,X,tmp_trans,list_1_org,list_1,list_2_1,list_2_2) \
shared(tmp_v0,tmp_v1,one,nstate)
        for (j = 1; j <= i_max; j++) {
          num1 = X_Spin_CisAit(j, X, isite1, org_sigma2, &tmp_off);
          if (num1 != 0) {
            dmv = tmp_trans*(double)num1;
            zaxpy_(&nstate, &dmv, tmp_v1[j], &one, tmp_v0[tmp_off], &one);
          }
        }
      }
    }
  }
  return TRUE;
}
int GetPairExcitedStateGeneralSpin(
  struct BindStruct *X,
  int nstate,
  std::complex<double> **tmp_v0, 
  std::complex<double> **tmp_v1, 
  int iEx
)
{
  long int i, j;
  long int org_isite1, org_isite2, org_sigma1, org_sigma2;
  long int tmp_off = 0;
  long int off = 0;
  std::complex<double> tmp_trans = 0, dmv;
  long int i_max;
  int tmp_sgn, num1, one = 1;
  i_max = X->Check.idim_maxOrg;

  for (i = 0; i < X->Def.NPairExcitationOperator[iEx]; i++) {
    org_isite1 = X->Def.PairExcitationOperator[iEx][i][0] + 1;
    org_isite2 = X->Def.PairExcitationOperator[iEx][i][2] + 1;
    org_sigma1 = X->Def.PairExcitationOperator[iEx][i][1];
    org_sigma2 = X->Def.PairExcitationOperator[iEx][i][3];
    tmp_trans = X->Def.ParaPairExcitationOperator[iEx][i];
    if (org_isite1 == org_isite2) {
      if (org_isite1 > X->Def.Nsite) {
        if (org_sigma1 == org_sigma2) {
          // longitudinal magnetic field
          num1 = BitCheckGeneral((long int) myrank,
            org_isite1, org_sigma1, X->Def.SiteToBit, X->Def.Tpow);
          if (X->Def.PairExcitationOperator[iEx][i][4] == 0) {
            if (num1 == 0) {
#pragma omp parallel for default(none) private(j,dmv)  \
  firstprivate(i_max, tmp_trans)                  \
  shared(tmp_v0, tmp_v1,one,nstate)
              for (j = 1; j <= i_max; j++) {
                dmv = -tmp_trans;
                zaxpy_(&nstate, &dmv, tmp_v1[j], &one, tmp_v0[j], &one);
              }
            }
          }
          else {
            if (num1 != 0) {
#pragma omp parallel for default(none) private(j) \
  firstprivate(i_max, tmp_trans) shared(tmp_v0, tmp_v1,one,nstate)
              for (j = 1; j <= i_max; j++) {
                zaxpy_(&nstate, &tmp_trans, tmp_v1[j], &one, tmp_v0[j], &one);
              }
            }
          }
        }//org_sigma1=org_sigma2
        else {//org_sigma1 != org_sigma2
          X_child_CisAit_GeneralSpin_MPIdouble(org_isite1 - 1, org_sigma1, org_sigma2, 
            tmp_trans, X, nstate, tmp_v0, tmp_v1, i_max);
        }
      }
      else {//org_isite1 <= X->Def.Nsite
        if (org_sigma1 == org_sigma2) {
          // longitudinal magnetic field
          if (X->Def.PairExcitationOperator[iEx][i][4] == 0) {
#pragma omp parallel for default(none) private(j, num1,dmv) \
  firstprivate(i_max, org_isite1, org_sigma1, X, tmp_trans) \
  shared(tmp_v0, tmp_v1, list_1,one,nstate)
            for (j = 1; j <= i_max; j++) {
              num1 = BitCheckGeneral(list_1[j], org_isite1, org_sigma1, X->Def.SiteToBit, X->Def.Tpow);
              dmv = -tmp_trans * (1.0 - num1);
              zaxpy_(&nstate, &dmv, tmp_v1[j], &one, tmp_v0[j], &one);
            }
          }
          else {
#pragma omp parallel for default(none) private(j, num1,dmv)    \
  firstprivate(i_max, org_isite1, org_sigma1, X, tmp_trans) \
  shared(tmp_v0, tmp_v1, list_1,one,nstate)
            for (j = 1; j <= i_max; j++) {
              num1 = BitCheckGeneral(list_1[j], org_isite1, org_sigma1, X->Def.SiteToBit, X->Def.Tpow);
              dmv = tmp_trans * (std::complex<double>)num1;
              zaxpy_(&nstate, &dmv, tmp_v1[j], &one, tmp_v0[j], &one);
            }
          }
        }//org_sigma1=org_sigma2
        else {//org_sigma1 != org_sigma2
#pragma omp parallel for default(none) private(j, tmp_sgn, tmp_off)      \
  firstprivate(i_max, org_isite1, org_sigma1, org_sigma2, X, off, tmp_trans, myrank) \
  shared(tmp_v0, tmp_v1, list_1_org, list_1,one,nstate)
          for (j = 1; j <= i_max; j++) {
            tmp_sgn = GetOffCompGeneralSpin(list_1_org[j], org_isite1, org_sigma2, org_sigma1, &off,
              X->Def.SiteToBit, X->Def.Tpow);
            if (tmp_sgn != FALSE) {
              ConvertToList1GeneralSpin(off, X->Large.ihfbit, &tmp_off);
              zaxpy_(&nstate, &tmp_trans, tmp_v1[j], &one, tmp_v0[tmp_off], &one);
            }
          }
        }
      }
    }
    else {
      fprintf(stdoutMPI, "ERROR: hopping is not allowed in localized spin system\n");
      return FALSE;
    }//org_isite1 != org_isite2
  }

  return TRUE;
}
int GetPairExcitedStateSpin(
  struct BindStruct *X,
  int nstate,
  std::complex<double> **tmp_v0, 
  std::complex<double> **tmp_v1, 
  int iEx
) {
  int iret = 0;
  if (X->Def.iFlgGeneralSpin == FALSE) {
    iret = GetPairExcitedStateHalfSpin(X, nstate, tmp_v0, tmp_v1, iEx);
  }
  else {
    iret = GetPairExcitedStateGeneralSpin(X, nstate, tmp_v0, tmp_v1, iEx);
  }
  return iret;
}