expec_totalspin.c File Reference

File for calculating total spin. More...

#include <bitcalc.h>
#include "mltplyCommon.h"
#include "mltplySpinCore.h"
#include "wrapperMPI.h"
#include "mltplyMPISpin.h"
#include "mltplyMPISpinCore.h"
#include "expec_totalspin.h"

Include dependency graph for expec_totalspin.c:

Go to the source code of this file.

Functions
int	expec_totalspin (struct BindStruct *X, double complex *vec)
	Parent function of calculation of total spin. More...
void	totalspin_Hubbard (struct BindStruct *X, double complex *vec)
	function of calculating totalspin for Hubbard model More...
void	totalspin_HubbardGC (struct BindStruct *X, double complex *vec)
	function of calculating totalspin for Hubbard model in grand canonical ensemble More...
void	totalspin_Spin (struct BindStruct *X, double complex *vec)
	function of calculating totalspin for spin model More...
void	totalspin_SpinGC (struct BindStruct *X, double complex *vec)
	function of calculating totalspin for spin model in grand canonical ensemble More...
int	expec_totalSz (struct BindStruct *X, double complex *vec)
void	totalSz_HubbardGC (struct BindStruct *X, double complex *vec)
	function of calculating totalSz for Hubbard model in grand canonical ensemble More...
void	totalSz_SpinGC (struct BindStruct *X, double complex *vec)
	function of calculating totalSz for Spin model in grand canonical ensemble More...

Detailed Description

File for calculating total spin.

Version

0.2

modify to treat the case of general spin

Version

0.1

Author

Takahiro Misawa (The University of Tokyo)

Kazuyoshi Yoshimi (The University of Tokyo)

Definition in file expec_totalspin.c.

Function Documentation

expec_totalspin()

int expec_totalspin	(	struct BindStruct *	X,
		double complex *	vec
	)

Parent function of calculation of total spin.

Parameters

[in,out]	X	data list of calculation parameters
[in]	vec	eigenvector

Author

Takahiro Misawa (The University of Tokyo)

Kazuyoshi Yoshimi (The University of Tokyo)

Return values

0	calculation is normally finished

Definition at line 50 of file expec_totalspin.c.

References BindStruct::Def, DefineList::iCalcModel, BindStruct::Large, LargeList::mode, BindStruct::Phys, PhysList::s2, PhysList::Sz, DefineList::Total2SzMPI, totalspin_Hubbard(), totalspin_HubbardGC(), totalspin_Spin(), and totalspin_SpinGC().

Referenced by phys().

{
  X->Large.mode = M_TOTALS;
  switch(X->Def.iCalcModel){
  case Spin:
     totalspin_Spin(X,vec);
     X->Phys.Sz=X->Def.Total2SzMPI/2.;
     break;
  case SpinGC:
     totalspin_SpinGC(X,vec);
     break;
   case Hubbard:
   case Kondo:
     totalspin_Hubbard(X,vec);
     break;
  case HubbardGC:
  case KondoGC:
     totalspin_HubbardGC(X,vec);
     break;
  default:
    X->Phys.s2=0.0;
    X->Phys.Sz=0.0;
  }
  return 0;
}

Here is the call graph for this function:

Here is the caller graph for this function:

expec_totalSz()

int expec_totalSz	(	struct BindStruct *	X,
		double complex *	vec
	)

Definition at line 649 of file expec_totalspin.c.

References BindStruct::Def, DefineList::iCalcModel, BindStruct::Large, LargeList::mode, BindStruct::Phys, PhysList::Sz, DefineList::Total2SzMPI, totalSz_HubbardGC(), and totalSz_SpinGC().

Referenced by CalcByLanczos().

  {
  X->Large.mode = M_TOTALS;
  switch (X->Def.iCalcModel) {
    case Spin:
      X->Phys.Sz = X->Def.Total2SzMPI / 2.;
      break;
    case SpinGC:
      totalSz_SpinGC(X, vec);
      break;
    case Hubbard:
    case Kondo:
      X->Phys.Sz = X->Def.Total2SzMPI / 2.;

      break;
    case HubbardGC:
    case KondoGC:
      totalSz_HubbardGC(X, vec);
      break;
    default:
      X->Phys.Sz = 0.0;
  }

  return 0;
}

Here is the call graph for this function:

Here is the caller graph for this function:

totalspin_Hubbard()

void totalspin_Hubbard	(	struct BindStruct *	X,
		double complex *	vec
	)

function of calculating totalspin for Hubbard model

Parameters

[in,out]	X	data list of calculation parameters
	vec	eigenvector

Version

0.1

Author

Takahiro Misawa (The University of Tokyo)

Kazuyoshi Yoshimi (The University of Tokyo)

Definition at line 88 of file expec_totalspin.c.

References BindStruct::Check, BindStruct::Def, GetOffComp(), GetSplitBitByModel(), DefineList::iCalcModel, CheckList::idim_max, list_1, list_2_1, list_2_2, DefineList::Nsite, DefineList::NsiteMPI, BindStruct::Phys, PhysList::s2, PhysList::Sz, and DefineList::Tpow.

Referenced by expec_totalspin().

                                                                 {
  long unsigned int j;
  long unsigned int irght, ilft, ihfbit;
  long unsigned int isite1, isite2;
  long unsigned int is1_up, is2_up, is1_down, is2_down;
  long unsigned int iexchg, off;
  int num1_up, num2_up;
  int num1_down, num2_down;
  long unsigned int ibit1_up, ibit2_up, ibit1_down, ibit2_down;
  double complex spn_z, tmp_spn_z;
  double complex spn;
  long unsigned i_max;
  i_max = X->Check.idim_max;

  GetSplitBitByModel(X->Def.Nsite, X->Def.iCalcModel, &irght, &ilft, &ihfbit);
  spn = 0.0;
  spn_z = 0.0;
  for (isite1 = 1; isite1 <= X->Def.NsiteMPI; isite1++) {
    is1_up = X->Def.Tpow[2 * isite1 - 2];
    is1_down = X->Def.Tpow[2 * isite1 - 1];

    for (isite2 = 1; isite2 <= X->Def.NsiteMPI; isite2++) {
      is2_up = X->Def.Tpow[2 * isite2 - 2];
      is2_down = X->Def.Tpow[2 * isite2 - 1];

#pragma omp parallel for reduction(+:spn, spn_z) default(none) firstprivate(i_max, is1_up, is1_down, is2_up, is2_down, irght, ilft, ihfbit, isite1, isite2) private(ibit1_up, num1_up, ibit2_up, num2_up, ibit1_down, num1_down, ibit2_down, num2_down, tmp_spn_z, iexchg, off) shared(vec, list_1, list_2_1, list_2_2)
      for (j = 1; j <= i_max; j++) {

        ibit1_up = list_1[j] & is1_up;
        num1_up = ibit1_up / is1_up;
        ibit2_up = list_1[j] & is2_up;
        num2_up = ibit2_up / is2_up;

        ibit2_down = list_1[j] & is2_down;
        num2_down = ibit2_down / is2_down;
        ibit1_down = list_1[j] & is1_down;
        num1_down = ibit1_down / is1_down;

        tmp_spn_z = (num1_up - num1_down) * (num2_up - num2_down);
        spn += conj(vec[j]) * vec[j] * tmp_spn_z / 4.0;
        if (isite1 == isite2) {
          spn += conj(vec[j]) * vec[j] * (num1_up + num1_down - 2 * num1_up * num1_down) / 2.0;
          spn_z += conj(vec[j]) * vec[j] * (num1_up - num1_down) / 2.0;
        } else {
          if (ibit1_up != 0 && ibit1_down == 0 && ibit2_up == 0 && ibit2_down != 0) {
            iexchg = list_1[j] - (is1_up + is2_down);
            iexchg += (is2_up + is1_down);
            GetOffComp(list_2_1, list_2_2, iexchg, irght, ilft, ihfbit, &off);
            spn += conj(vec[j]) * vec[off] / 2.0;
          } else if (ibit1_up == 0 && ibit1_down != 0 && ibit2_up != 0 && ibit2_down == 0) {
            iexchg = list_1[j] - (is1_down + is2_up);
            iexchg += (is2_down + is1_up);
            GetOffComp(list_2_1, list_2_2, iexchg, irght, ilft, ihfbit, &off);
            spn += conj(vec[j]) * vec[off] / 2.0;
          }
        }
      }
    }
  }
  X->Phys.s2 = creal(spn);
  X->Phys.Sz = creal(spn_z);
}

Here is the call graph for this function:

Here is the caller graph for this function:

totalspin_HubbardGC()

void totalspin_HubbardGC	(	struct BindStruct *	X,
		double complex *	vec
	)

function of calculating totalspin for Hubbard model in grand canonical ensemble

Parameters

[in,out]	X	data list of calculation parameters
	vec	eigenvector

Version

0.1

Author

Takahiro Misawa (The University of Tokyo)

Kazuyoshi Yoshimi (The University of Tokyo)

Definition at line 161 of file expec_totalspin.c.

References BindStruct::Check, BindStruct::Def, CheckList::idim_max, DefineList::NsiteMPI, BindStruct::Phys, PhysList::s2, PhysList::Sz, and DefineList::Tpow.

Referenced by expec_totalspin().

                                                                   {
  long unsigned int j;
  long unsigned int isite1, isite2;
  long unsigned int is1_up, is2_up, is1_down, is2_down;
  long unsigned int iexchg, off;
  int num1_up, num2_up;
  int num1_down, num2_down;
  long unsigned int ibit1_up, ibit2_up, ibit1_down, ibit2_down, list_1_j;
  double complex spn_z, tmp_spn_z;
  double complex spn;
  long unsigned int i_max;

  i_max = X->Check.idim_max;

  spn = 0.0;
  spn_z = 0.0;
  for (isite1 = 1; isite1 <= X->Def.NsiteMPI; isite1++) {
    for (isite2 = 1; isite2 <= X->Def.NsiteMPI; isite2++) {
      is1_up = X->Def.Tpow[2 * isite1 - 2];
      is1_down = X->Def.Tpow[2 * isite1 - 1];
      is2_up = X->Def.Tpow[2 * isite2 - 2];
      is2_down = X->Def.Tpow[2 * isite2 - 1];

#pragma omp parallel for reduction(+:spn, spn_z) default(none) firstprivate(i_max, is1_up, is1_down, is2_up, is2_down, isite1, isite2) private(list_1_j, ibit1_up, num1_up, ibit2_up, num2_up, ibit1_down, num1_down, ibit2_down, num2_down, tmp_spn_z, iexchg, off) shared(vec)
      for (j = 1; j <= i_max; j++) {
        list_1_j = j - 1;
        ibit1_up = list_1_j & is1_up;
        num1_up = ibit1_up / is1_up;
        ibit2_up = list_1_j & is2_up;
        num2_up = ibit2_up / is2_up;

        ibit1_down = list_1_j & is1_down;
        num1_down = ibit1_down / is1_down;
        ibit2_down = list_1_j & is2_down;
        num2_down = ibit2_down / is2_down;

        tmp_spn_z = (num1_up - num1_down) * (num2_up - num2_down);
        spn += conj(vec[j]) * vec[j] * tmp_spn_z / 4.0;
        if (isite1 == isite2) {
          spn += conj(vec[j]) * vec[j] * (num1_up + num1_down - 2 * num1_up * num1_down) / 2.0;
          spn_z += conj(vec[j]) * vec[j] * (num1_up - num1_down) / 2.0;
        } else {
          if (ibit1_up != 0 && ibit1_down == 0 && ibit2_up == 0 && ibit2_down != 0) {
            iexchg = list_1_j - (is1_up + is2_down);
            iexchg += (is2_up + is1_down);
            off = iexchg + 1;
            spn += conj(vec[j]) * vec[off] / 2.0;
          } else if (ibit1_up == 0 && ibit1_down != 0 && ibit2_up != 0 && ibit2_down == 0) {
            iexchg = list_1_j - (is1_down + is2_up);
            iexchg += (is2_down + is1_up);
            off = iexchg + 1;
            spn += conj(vec[j]) * vec[off] / 2.0;
          }
        }
      }
    }
  }
  X->Phys.s2 = creal(spn);
  X->Phys.Sz = creal(spn_z);
}

Here is the caller graph for this function:

totalspin_Spin()

void totalspin_Spin	(	struct BindStruct *	X,
		double complex *	vec
	)

function of calculating totalspin for spin model

Parameters

[in,out]	X	data list of calculation parameters
	vec	eigenvector

Version

0.2

modify for hybrid parallel

Version

0.1

Author

Takahiro Misawa (The University of Tokyo)

Kazuyoshi Yoshimi (The University of Tokyo)

Definition at line 233 of file expec_totalspin.c.

References BindStruct::Check, ConvertToList1GeneralSpin(), BindStruct::Def, FALSE, GetBitGeneral(), GetLocal2Sz(), GetOffComp(), GetOffCompGeneralSpin(), GetSplitBitByModel(), DefineList::iCalcModel, CheckList::idim_max, DefineList::iFlgGeneralSpin, list_1, list_2_1, list_2_2, myrank, DefineList::Nsite, DefineList::NsiteMPI, BindStruct::Phys, PhysList::s2, CheckList::sdim, DefineList::SiteToBit, SumMPI_dc(), PhysList::Sz, DefineList::Tpow, TRUE, X_child_general_int_spin_MPIsingle(), X_child_general_int_spin_TotalS_MPIdouble(), and X_SpinGC_CisAis().

Referenced by expec_totalspin().

                                                              {

  long unsigned int j;
  long unsigned int irght, ilft, ihfbit;
  long unsigned int isite1, isite2;
  long unsigned int tmp_isite1, tmp_isite2;

  long unsigned int is1_up, is2_up;
  long unsigned int iexchg, off, off_2;

  int num1_up, num2_up;
  int num1_down, num2_down;
  int sigma_1, sigma_2;
  long unsigned int ibit1_up, ibit2_up, ibit_tmp, is_up;
  double complex spn_z = 0.0;
  double complex spn_z1, spn_z2, spn_zd;
  double complex spn = 0.0;
  long unsigned int i_max;

  i_max = X->Check.idim_max;
  if (X->Def.iFlgGeneralSpin == FALSE) {
    GetSplitBitByModel(X->Def.Nsite, X->Def.iCalcModel, &irght, &ilft, &ihfbit);
    spn = 0.0;
    spn_z = 0.0;
    spn_zd = 0.0;
    for (isite1 = 1; isite1 <= X->Def.NsiteMPI; isite1++) {
      for (isite2 = 1; isite2 <= X->Def.NsiteMPI; isite2++) {

        if (isite1 > X->Def.Nsite && isite2 > X->Def.Nsite) {
#ifdef MPI
          is1_up = X->Def.Tpow[isite1 - 1];
          is2_up = X->Def.Tpow[isite2 - 1];
          is_up = is1_up + is2_up;
          num1_up = X_SpinGC_CisAis((unsigned long int) myrank + 1, X, is1_up, 1);
          num1_down = 1 - num1_up;
          num2_up = X_SpinGC_CisAis((unsigned long int) myrank + 1, X, is2_up, 1);
          num2_down = 1 - num2_up;
          spn_z = (num1_up - num1_down) * (num2_up - num2_down);

#pragma omp parallel for default(none) reduction (+:spn_zd) shared(vec)  \
  firstprivate(i_max, spn_z) private(j)
          for (j = 1; j <= i_max; j++) {
            spn_zd += conj(vec[j]) * vec[j] * spn_z / 4.0;
          }
          if (isite1 == isite2) {
#pragma omp parallel for default(none) reduction (+:spn_zd) shared(vec)  \
  firstprivate(i_max) private(j)
            for (j = 1; j <= i_max; j++) {
              spn_zd += conj(vec[j]) * vec[j] / 2.0;
            }
          } else {//off diagonal
            spn += X_child_general_int_spin_TotalS_MPIdouble(isite1 - 1, isite2 - 1, X, vec, vec);
          }
#endif
        } else if (isite1 > X->Def.Nsite || isite2 > X->Def.Nsite) {
#ifdef MPI
          if (isite1 < isite2) {
            tmp_isite1 = isite1;
            tmp_isite2 = isite2;
          } else {
            tmp_isite1 = isite2;
            tmp_isite2 = isite1;
          }

          is1_up = X->Def.Tpow[tmp_isite1 - 1];
          is2_up = X->Def.Tpow[tmp_isite2 - 1];
          num2_up = X_SpinGC_CisAis((unsigned long int) myrank + 1, X, is2_up, 1);
          num2_down = 1 - num2_up;

          //diagonal
#pragma omp parallel for reduction(+: spn_zd) default(none) firstprivate(i_max, is1_up, num2_up, num2_down) private(ibit1_up, num1_up, num1_down, spn_z) shared(list_1, vec)
          for (j = 1; j <= i_max; j++) {
            ibit1_up = list_1[j] & is1_up;
            num1_up = ibit1_up / is1_up;
            num1_down = 1 - num1_up;
            spn_z = (num1_up - num1_down) * (num2_up - num2_down);
            spn_zd += conj(vec[j]) * vec[j] * spn_z / 4.0;
          }
          if (isite1 < isite2) {
            spn += X_child_general_int_spin_MPIsingle(isite1 - 1, 0, 1, isite2 - 1, 1, 0, 1.0, X, vec, vec);
          } else {
            spn += conj(X_child_general_int_spin_MPIsingle(isite2 - 1, 1, 0, isite1 - 1, 0, 1, 1.0, X, vec, vec));
          }
#endif
        }//isite1 > Nsite || isite2 > Nsite
        else {
          is1_up = X->Def.Tpow[isite1 - 1];
          is2_up = X->Def.Tpow[isite2 - 1];
          is_up = is1_up + is2_up;

#pragma omp parallel for reduction(+: spn, spn_zd) default(none) firstprivate(i_max, is_up, is1_up, is2_up, irght, ilft, ihfbit, isite1, isite2) private(ibit1_up, num1_up, ibit2_up, num2_up, num1_down, num2_down, spn_z, iexchg, off, ibit_tmp) shared(list_1, list_2_1, list_2_2, vec)
          for (j = 1; j <= i_max; j++) {
            ibit1_up = list_1[j] & is1_up;
            num1_up = ibit1_up / is1_up;
            num1_down = 1 - num1_up;
            ibit2_up = list_1[j] & is2_up;
            num2_up = ibit2_up / is2_up;
            num2_down = 1 - num2_up;

            spn_z = (num1_up - num1_down) * (num2_up - num2_down);
            spn_zd += conj(vec[j]) * vec[j] * spn_z / 4.0;

            if (isite1 == isite2) {
              spn_zd += conj(vec[j]) * vec[j] / 2.0;
            } else {
              ibit_tmp = (num1_up) ^ (num2_up);
              if (ibit_tmp != 0) {
                iexchg = list_1[j] ^ (is_up);
                GetOffComp(list_2_1, list_2_2, iexchg, irght, ilft, ihfbit, &off);
                spn += conj(vec[j]) * vec[off] / 2.0;
              }
            }
          }// j
        }
      }//isite2
    }//isite1
  }//generalspin=FALSE
  else {
    double S1 = 0;
    double S2 = 0;
    for (isite1 = 1; isite1 <= X->Def.NsiteMPI; isite1++) {
      for (isite2 = 1; isite2 <= X->Def.NsiteMPI; isite2++) {
        S1 = 0.5 * (X->Def.SiteToBit[isite1 - 1] - 1);
        S2 = 0.5 * (X->Def.SiteToBit[isite2 - 1] - 1);
        if (isite1 == isite2) {
#pragma omp parallel for reduction(+: spn, spn_z) default(none) firstprivate(i_max, isite1, X, S1) private (spn_z1)shared(vec, list_1)
          for (j = 1; j <= i_max; j++) {
            spn_z1 = 0.5 * GetLocal2Sz(isite1, list_1[j], X->Def.SiteToBit, X->Def.Tpow);
            spn += conj(vec[j]) * vec[j] * S1 * (S1 + 1.0);
            spn_z += conj(vec[j]) * vec[j] * spn_z1;
          }
        } else {
#pragma omp parallel for reduction(+: spn) default(none) firstprivate(i_max, isite1, isite2, X, S1, S2) private(spn_z1, spn_z2, off, off_2, ibit_tmp, sigma_1, sigma_2) shared(vec, list_1)
          for (j = 1; j <= i_max; j++) {
            spn_z1 = 0.5 * GetLocal2Sz(isite1, list_1[j], X->Def.SiteToBit, X->Def.Tpow);
            spn_z2 = 0.5 * GetLocal2Sz(isite2, list_1[j], X->Def.SiteToBit, X->Def.Tpow);
            spn += conj(vec[j]) * vec[j] * spn_z1 * spn_z2;

            sigma_1 = GetBitGeneral(isite1, list_1[j], X->Def.SiteToBit, X->Def.Tpow);
            sigma_2 = GetBitGeneral(isite2, list_1[j], X->Def.SiteToBit, X->Def.Tpow);

            ibit_tmp = GetOffCompGeneralSpin(list_1[j], isite2, sigma_2, sigma_2 + 1, &off, X->Def.SiteToBit,
                                             X->Def.Tpow);
            if (ibit_tmp == TRUE) {
              ibit_tmp = GetOffCompGeneralSpin(off, isite1, sigma_1, sigma_1 - 1, &off_2, X->Def.SiteToBit,
                                               X->Def.Tpow);
              if (ibit_tmp == TRUE) {
                ConvertToList1GeneralSpin(off_2, X->Check.sdim, &off);
                spn += conj(vec[j]) * vec[off] * sqrt(S2 * (S2 + 1) - spn_z2 * (spn_z2 + 1)) *
                       sqrt(S1 * (S1 + 1) - spn_z1 * (spn_z1 - 1)) / 2.0;
              }
            }

            ibit_tmp = GetOffCompGeneralSpin(list_1[j], isite2, sigma_2, sigma_2 - 1, &off, X->Def.SiteToBit,
                                             X->Def.Tpow);
            if (ibit_tmp == TRUE) {
              ibit_tmp = GetOffCompGeneralSpin(off, isite1, sigma_1, sigma_1 + 1, &off_2, X->Def.SiteToBit,
                                               X->Def.Tpow);
              if (ibit_tmp == TRUE) {
                ConvertToList1GeneralSpin(off_2, X->Check.sdim, &off);
                spn += conj(vec[j]) * vec[off] * sqrt(S2 * (S2 + 1) - spn_z2 * (spn_z2 - 1.0)) *
                       sqrt(S1 * (S1 + 1) - spn_z1 * (spn_z1 + 1)) / 2.0;
              }
            }
          }
        }
      }
    }
  }

  spn = SumMPI_dc(spn);
  spn_zd = SumMPI_dc(spn_zd);
  spn_z = SumMPI_dc(spn_z);
  spn += spn_zd;
  X->Phys.s2 = creal(spn);
  X->Phys.Sz = creal(spn_z);
}

Here is the call graph for this function:

Here is the caller graph for this function:

totalspin_SpinGC()

void totalspin_SpinGC	(	struct BindStruct *	X,
		double complex *	vec
	)

function of calculating totalspin for spin model in grand canonical ensemble

Parameters

[in,out]	X	data list of calculation parameters
	vec	eigenvector

Version

0.2

add function to treat a calculation of total spin for general spin

Version

0.1

Author

Takahiro Misawa (The University of Tokyo)

Kazuyoshi Yoshimi (The University of Tokyo)

Definition at line 424 of file expec_totalspin.c.

References BindStruct::Check, BindStruct::Def, FALSE, GetBitGeneral(), GetLocal2Sz(), GetOffCompGeneralSpin(), CheckList::idim_max, DefineList::iFlgGeneralSpin, BindStruct::Large, LargeList::mode, myrank, DefineList::Nsite, DefineList::NsiteMPI, BindStruct::Phys, PhysList::s2, DefineList::SiteToBit, SumMPI_dc(), PhysList::Sz, DefineList::Tpow, X_GC_child_CisAitCiuAiv_spin_MPIdouble(), X_GC_child_CisAitCiuAiv_spin_MPIsingle(), and X_SpinGC_CisAis().

Referenced by expec_totalspin().

                                                               {

  long unsigned int j;
  long unsigned int isite1,isite2, tmp_isite1, tmp_isite2;
  long unsigned int is1_up,is2_up;
  long unsigned int iexchg, off, off_2;
  int num1_up,num2_up;
  int num1_down,num2_down;
  int sigma_1, sigma_2;
  long unsigned int ibit1_up,ibit2_up,ibit_tmp,is_up;
  double complex spn_z;
  double complex spn_z1, spn_z2;
  double complex spn, spn_d;
  long unsigned int list_1_j;
  long unsigned int i_max;
  i_max=X->Check.idim_max;
  X->Large.mode = M_TOTALS;
  spn=0.0;
  spn_d=0.0;
  spn_z=0.0;
  if(X->Def.iFlgGeneralSpin==FALSE){
    for(isite1=1;isite1<=X->Def.NsiteMPI;isite1++){
      if(isite1 > X->Def.Nsite){
  is1_up      = X->Def.Tpow[isite1-1];
  ibit1_up = myrank&is1_up;
  num1_up = ibit1_up/is1_up;
  num1_down =1-num1_up;
#pragma omp parallel for reduction(+: spn_z) default(none) firstprivate(i_max, is1_up,  num1_up, num1_down) shared(vec)
  for(j=1;j<=i_max;j++){
    spn_z  +=  conj(vec[j])*vec[j]*(num1_up-num1_down)/2.0;
  }
      }
      else{
  is1_up      = X->Def.Tpow[isite1-1];
#pragma omp parallel for reduction(+: spn_z) default(none) firstprivate(i_max, is1_up) private(list_1_j, ibit1_up, num1_up, num1_down) shared(vec)
  for(j=1;j<=i_max;j++){
    list_1_j=j-1;
    ibit1_up  = list_1_j&is1_up;
    num1_up   = ibit1_up/is1_up;
    num1_down = 1-num1_up;
    spn_z  +=  conj(vec[j])*vec[j]*(num1_up-num1_down)/2.0;
  }

      }
      for(isite2=1;isite2<=X->Def.NsiteMPI;isite2++){

  if(isite1 > X->Def.Nsite && isite2 > X->Def.Nsite){
    is1_up      = X->Def.Tpow[isite1-1];
    is2_up      = X->Def.Tpow[isite2-1];
    num1_up = X_SpinGC_CisAis((unsigned long int)myrank + 1, X, is1_up, 1);
    num1_down = 1-num1_up;
    num2_up = X_SpinGC_CisAis((unsigned long int)myrank + 1, X, is2_up, 1);
    num2_down = 1-num2_up;
    spn_z2  = (num1_up-num1_down)*(num2_up-num2_down)/4.0;
#pragma omp parallel for default(none) reduction (+:spn_d) shared(vec)  \
  firstprivate(i_max, spn_z2) private(j)
    for (j = 1; j <= i_max; j++) {
      spn_d   += conj(vec[j])*vec[j]*spn_z2;
    }
    if(isite1 == isite2){
#pragma omp parallel for default(none) reduction (+:spn_d) shared(vec)  \
  firstprivate(i_max) private(j)
      for (j = 1; j <= i_max; j++) {
        spn_d      += conj(vec[j])*vec[j]/2.0;
      }
    }//isite1 = isite2
    else{//off diagonal
      spn += X_GC_child_CisAitCiuAiv_spin_MPIdouble(isite1-1, 0, 1, isite2-1, 1, 0, 1.0, X, vec, vec)/2.0;
    }
  }
  else if(isite1 > X->Def.Nsite || isite2 > X->Def.Nsite){
    if(isite1 < isite2){
      tmp_isite1=isite1;
      tmp_isite2=isite2;
    }
    else{
      tmp_isite1 = isite2;
      tmp_isite2 = isite1;
    }
    is1_up = X->Def.Tpow[tmp_isite1 - 1];
    is2_up = X->Def.Tpow[tmp_isite2 - 1];
    num2_up = X_SpinGC_CisAis((unsigned long int)myrank + 1, X, is2_up, 1);
    num2_down =1-num2_up;
    //diagonal
#pragma omp parallel for reduction(+: spn_d) default(none) firstprivate(i_max, is1_up, num2_up, num2_down) private(ibit1_up, num1_up, num1_down, spn_z2, list_1_j) shared(vec)
    for(j=1;j<=i_max;j++){
      list_1_j=j-1;
      ibit1_up  = list_1_j&is1_up;
      num1_up   = ibit1_up/is1_up;
      num1_down = 1-num1_up;
      spn_z2  = (num1_up-num1_down)*(num2_up-num2_down);
      spn_d   += conj(vec[j])*vec[j]*spn_z2/4.0;
    }
    if(isite1 < isite2){
      spn += X_GC_child_CisAitCiuAiv_spin_MPIsingle(isite1-1, 0, 1, isite2-1, 1, 0, 1.0, X, vec, vec)/2.0;
    }
    else{
      spn += conj(X_GC_child_CisAitCiuAiv_spin_MPIsingle(isite2-1, 1, 0, isite1-1, 0, 1, 1.0, X, vec, vec))/2.0;
    }
  }
  else{
    is2_up      = X->Def.Tpow[isite2-1];
    is_up       = is1_up+is2_up;
#pragma omp parallel for reduction(+: spn, spn_d) default(none) firstprivate(i_max, is_up, is1_up, is2_up, isite1, isite2) private(list_1_j, ibit1_up, num1_up, ibit2_up, num2_up, num1_down, num2_down, spn_z2, iexchg, off, ibit_tmp) shared(vec)
    for(j=1;j<=i_max;j++){
      list_1_j=j-1;
      ibit1_up  = list_1_j&is1_up;
      num1_up   = ibit1_up/is1_up;
      num1_down = 1-num1_up;
      ibit2_up  = list_1_j&is2_up;
      num2_up   = ibit2_up/is2_up;
      num2_down = 1-num2_up;

      spn_z2  = (num1_up-num1_down)*(num2_up-num2_down);
      spn_d   += conj(vec[j])*vec[j]*spn_z2/4.0;

      if(isite1==isite2){
        spn_d      += conj(vec[j])*vec[j]/2.0;
      }else{
        ibit_tmp  = (num1_up) ^ (num2_up);
        if(ibit_tmp!=0){
    iexchg  = list_1_j ^ (is_up);
    off    = iexchg+1;
    spn    += conj(vec[j])*vec[off]/2.0;
        }

      }
    }//j
  }//else
      }
    }
  }
  else{//general spin
    double S1=0;
    double S2=0;
    spn =0.0;
    spn_z=0.0;
    for(isite1=1;isite1<=X->Def.NsiteMPI;isite1++){
      S1=0.5*(X->Def.SiteToBit[isite1-1]-1);
      if(isite1 > X->Def.Nsite){
  spn_z1  = 0.5*GetLocal2Sz(isite1, (unsigned long int) myrank, X->Def.SiteToBit, X->Def.Tpow);
#pragma omp parallel for reduction(+: spn, spn_z) default(none) firstprivate(S1, spn_z1,i_max) shared(vec)
  for(j=1;j<=i_max;j++){
    spn   += conj(vec[j])*vec[j]*S1*(S1+1.0);
    spn_z += conj(vec[j])*vec[j]*spn_z1;
  }
      }
      else{
#pragma omp parallel for reduction(+: spn, spn_z) default(none) firstprivate(i_max, isite1, X, S1) private(spn_z1) shared(vec)
  for(j=1;j<=i_max;j++){
    spn_z1  = 0.5*GetLocal2Sz(isite1, j-1, X->Def.SiteToBit, X->Def.Tpow);
    spn    += conj(vec[j])*vec[j]*S1*(S1+1.0);
    spn_z += conj(vec[j])*vec[j]*spn_z1;
  }
      }
      for(isite2=1;isite2<=X->Def.NsiteMPI;isite2++){
  if(isite1==isite2) continue;
  S2=0.5*(X->Def.SiteToBit[isite2-1]-1);

  if(isite1 > X->Def.Nsite && isite2 > X->Def.Nsite){
    /*
    spn_z1  = 0.5*GetLocal2Sz(isite1, (unsigned long int) myrank, X->Def.SiteToBit, X->Def.Tpow);
    spn_z2  = 0.5*GetLocal2Sz(isite2, (unsigned long int) myrank, X->Def.SiteToBit, X->Def.Tpow);
#pragma omp parallel for reduction(+: spn, spn_z) default(none) firstprivate(spn_z1, spn_z2, i_max) shared(vec)
    for(j=1;j<=i_max;j++){
      spn   += conj(vec[j])*vec[j]*spn_z1*spn_z2;
    }
    tmp_V= sqrt(S2*(S2+1) - spn_z2*(spn_z2+1))*sqrt(S1*(S1+1) - spn_z1*(spn_z1-1))/2.0;
    spn += X_GC_child_CisAitCjuAjv_GeneralSpin_MPIdouble(isite1-1, sigma_1-1, sigma_1, isite2-1, sigma_2+1, sigma_2, tmp_V, X,vec, vec);
    tmp_V= sqrt(S2*(S2+1) - spn_z2*(spn_z2-1))*sqrt(S1*(S1+1) - spn_z1*(spn_z1+1))/2.0;
    spn += X_GC_child_CisAitCjuAjv_GeneralSpin_MPIdouble(isite1-1, sigma_1+1, sigma_1, isite2-1, sigma_2-1, sigma_2, tmp_V, X,vec, vec);
    */
  }
  else if(isite1 > X->Def.Nsite || isite2 > X->Def.Nsite){
    /*
   if(isite1 < isite2){
      tmp_isite1=isite1;
      tmp_isite2=isite2;
    }
    else{
      tmp_isite1 = isite2;
      tmp_isite2 = isite1;
    }
   spn_z2  = 0.5*GetLocal2Sz(tmp_isite2, (unsigned long int) myrank, X->Def.SiteToBit, X->Def.Tpow);
   sigma_2 = GetBitGeneral(tmp_isite2, (unsigned long int) myrank, X->Def.SiteToBit, X->Def.Tpow);
    */
  }
  else{ //inner-process
#pragma omp parallel for reduction(+: spn) default(none) firstprivate(i_max, isite1, isite2, X, S1, S2) private(spn_z1, spn_z2, off, off_2, ibit_tmp, sigma_1, sigma_2) shared(vec)
    for(j=1;j<=i_max;j++){
      spn_z1  = 0.5*GetLocal2Sz(isite1, j-1, X->Def.SiteToBit, X->Def.Tpow);
      spn_z2  = 0.5*GetLocal2Sz(isite2, j-1, X->Def.SiteToBit, X->Def.Tpow);
      spn    += conj(vec[j])*vec[j]*spn_z1*spn_z2;

      sigma_1=GetBitGeneral(isite1, j-1, X->Def.SiteToBit, X->Def.Tpow);
      sigma_2=GetBitGeneral(isite2, j-1, X->Def.SiteToBit, X->Def.Tpow);

      ibit_tmp = GetOffCompGeneralSpin(j-1, isite2, sigma_2, sigma_2+1, &off, X->Def.SiteToBit, X->Def.Tpow);
      if(ibit_tmp!=0){
        ibit_tmp = GetOffCompGeneralSpin(off, isite1, sigma_1, sigma_1-1,&off_2, X->Def.SiteToBit, X->Def.Tpow);
        if(ibit_tmp!=0){
    spn += conj(vec[j])*vec[off_2+1]*sqrt(S2*(S2+1) - spn_z2*(spn_z2+1))*sqrt(S1*(S1+1) - spn_z1*(spn_z1-1))/2.0;
        }
      }

      ibit_tmp = GetOffCompGeneralSpin(j-1, isite2, sigma_2, sigma_2-1, &off, X->Def.SiteToBit, X->Def.Tpow);
      if(ibit_tmp !=0){
        ibit_tmp = GetOffCompGeneralSpin(off, isite1, sigma_1, sigma_1+1, &off_2, X->Def.SiteToBit, X->Def.Tpow);
        if(ibit_tmp!=0){
    spn += conj(vec[j])*vec[off_2+1]*sqrt(S2*(S2+1) - spn_z2*(spn_z2-1.0))*sqrt(S1*(S1+1)- spn_z1*(spn_z1+1))/2.0;
        }
      }
    }//j
  }//inner-process
      }//isite2
    }//isite1
  }

  spn = SumMPI_dc(spn);
  spn_d = SumMPI_dc(spn_d);
  spn_z = SumMPI_dc(spn_z);
  X->Phys.s2=creal(spn+spn_d);
  X->Phys.Sz=creal(spn_z);
}

Here is the call graph for this function:

Here is the caller graph for this function:

totalSz_HubbardGC()

void totalSz_HubbardGC	(	struct BindStruct *	X,
		double complex *	vec
	)

function of calculating totalSz for Hubbard model in grand canonical ensemble

Parameters

[in,out]	X	data list of calculation parameters
	vec	eigenvector

Version

0.1

Author

Takahiro Misawa (The University of Tokyo)

Kazuyoshi Yoshimi (The University of Tokyo)

Definition at line 688 of file expec_totalspin.c.

References BindStruct::Check, BindStruct::Def, CheckList::idim_max, myrank, DefineList::Nsite, DefineList::NsiteMPI, BindStruct::Phys, SumMPI_dc(), PhysList::Sz, totalSz_SpinGC(), and DefineList::Tpow.

Referenced by expec_totalSz().

   {

  long unsigned int j;
  long unsigned int isite1;
  long unsigned int is1_up, is1_down;
  int num1_up, num1_down, num1_sz;
  long unsigned int ibit1_up, ibit1_down, list_1_j;
  double complex spn_z;
  long unsigned int i_max;

  i_max = X->Check.idim_max;
  spn_z = 0.0;
  for (isite1 = 1; isite1 <= X->Def.NsiteMPI; isite1++) {
    if (isite1 > X->Def.Nsite) {
#ifdef MPI
      is1_up = X->Def.Tpow[2 * isite1 - 2];
      is1_down = X->Def.Tpow[2 * isite1 - 1];
      ibit1_up = (unsigned long int) myrank & is1_up;
      num1_up = ibit1_up / is1_up;
      ibit1_down = (unsigned long int) myrank & is1_down;
      num1_down = ibit1_down / is1_down;
      num1_sz = num1_up - num1_down;
#pragma omp parallel for reduction(+:spn_z) default(none) firstprivate(i_max) private(j) shared(num1_sz,vec)
      for (j = 1; j <= i_max; j++) {
        spn_z += (num1_sz) / 2.0 * conj(vec[j]) * vec[j];
      }
#endif
    } else {//isite1 > X->Def.Nsite
      is1_up = X->Def.Tpow[2 * isite1 - 2];
      is1_down = X->Def.Tpow[2 * isite1 - 1];
#pragma omp parallel for reduction(+:spn_z) default(none) firstprivate(i_max, is1_up, is1_down, isite1) \
  private(list_1_j, ibit1_up, num1_up, ibit1_down, num1_down) shared(vec)
      for (j = 1; j <= i_max; j++) {
        list_1_j = j - 1;
        ibit1_up = list_1_j & is1_up;
        num1_up = ibit1_up / is1_up;
        ibit1_down = list_1_j & is1_down;
        num1_down = ibit1_down / is1_down;
        spn_z += conj(vec[j]) * vec[j] * (num1_up - num1_down) / 2.0;
      }
    }
  }
  spn_z = SumMPI_dc(spn_z);
  X->Phys.Sz = creal(spn_z);
}

Here is the call graph for this function:

Here is the caller graph for this function:

totalSz_SpinGC()

void totalSz_SpinGC	(	struct BindStruct *	X,
		double complex *	vec
	)

function of calculating totalSz for Spin model in grand canonical ensemble

Parameters

[in,out]	X	data list of calculation parameters
	vec	eigenvector

Version

0.1

Author

Takahiro Misawa (The University of Tokyo)

Kazuyoshi Yoshimi (The University of Tokyo)

Definition at line 747 of file expec_totalspin.c.

References BindStruct::Check, BindStruct::Def, FALSE, GetLocal2Sz(), CheckList::idim_max, DefineList::iFlgGeneralSpin, BindStruct::Large, LargeList::mode, myrank, DefineList::Nsite, DefineList::NsiteMPI, BindStruct::Phys, DefineList::SiteToBit, SumMPI_dc(), PhysList::Sz, and DefineList::Tpow.

Referenced by expec_totalSz(), and totalSz_HubbardGC().

   {
  long unsigned int j, list_1_j;
  long unsigned int isite1;
  long unsigned int is1_up;
  int num1_up;
  int num1_down;
  long unsigned int ibit1_up;
  double complex spn_z, spn_z1;
  long unsigned int i_max;
  i_max = X->Check.idim_max;
  X->Large.mode = M_TOTALS;
  spn_z = 0.0;
  if (X->Def.iFlgGeneralSpin == FALSE) {
    for (isite1 = 1; isite1 <= X->Def.NsiteMPI; isite1++) {
      if (isite1 > X->Def.Nsite) {
#ifdef MPI
        is1_up = X->Def.Tpow[isite1 - 1];
        ibit1_up = myrank & is1_up;
        num1_up = ibit1_up / is1_up;
        num1_down = 1 - num1_up;
#pragma omp parallel for reduction(+: spn_z) default(none) firstprivate(i_max, is1_up, num1_up, num1_down) shared(vec)
        for (j = 1; j <= i_max; j++) {
          spn_z += conj(vec[j]) * vec[j] * (num1_up - num1_down) / 2.0;
        }
#endif
      } else {
        is1_up = X->Def.Tpow[isite1 - 1];
#pragma omp parallel for reduction(+: spn_z) default(none) firstprivate(i_max, is1_up) private(list_1_j, ibit1_up, num1_up, num1_down) shared(vec)
        for (j = 1; j <= i_max; j++) {
          list_1_j = j - 1;
          ibit1_up = list_1_j & is1_up;
          num1_up = ibit1_up / is1_up;
          num1_down = 1 - num1_up;
          spn_z += conj(vec[j]) * vec[j] * (num1_up - num1_down) / 2.0;
        }
      }//else
    }
  } else {//general spin
    spn_z = 0.0;
    for (isite1 = 1; isite1 <= X->Def.NsiteMPI; isite1++) {
      if (isite1 > X->Def.Nsite) {
        spn_z1 = 0.5 * GetLocal2Sz(isite1, (unsigned long int) myrank, X->Def.SiteToBit, X->Def.Tpow);
#pragma omp parallel for reduction(+: spn_z) default(none) firstprivate(spn_z1, i_max) shared(vec)
        for (j = 1; j <= i_max; j++) {
          spn_z += conj(vec[j]) * vec[j] * spn_z1;
        }
      } else {
#pragma omp parallel for reduction(+:  spn_z) default(none) firstprivate(i_max, isite1, X) private(spn_z1) shared(vec)
        for (j = 1; j <= i_max; j++) {
          spn_z1 = 0.5 * GetLocal2Sz(isite1, j - 1, X->Def.SiteToBit, X->Def.Tpow);
          spn_z += conj(vec[j]) * vec[j] * spn_z1;
        }
      }
    }//isite1
  }
  spn_z = SumMPI_dc(spn_z);
  X->Phys.Sz = creal(spn_z);
}

Here is the call graph for this function:

Here is the caller graph for this function: