phys.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 "phys.hpp"
#include "expec_energy_flct.hpp"
#include "expec_totalspin.hpp"
#include "expec_cisajs.hpp"
#include "expec_cisajscktaltdc.hpp"
#include "wrapperMPI.hpp"
#ifdef _SCALAPACK
#include "matrixscalapack.hpp"
#endif

void phys(struct BindStruct *X, 
          long int neig 
) {
  long int i;
  double tmp_N;
#ifdef _SCALAPACK
  std::complex<double> *vec_tmp;
  int ictxt, ierr, rank;
  long int j, i_max;

  i_max = X->Check.idim_max;

  if(use_scalapack){
  fprintf(stdoutMPI, "In scalapack fulldiag, total spin is not calculated !\n");
  vec_tmp = malloc(i_max*sizeof(std::complex<double>));
  }
  for (i = 0; i < neig; i++) {
    for (j = 0; j < i_max; j++) {
      v0[j + 1] = 0.0;
    }
    if (use_scalapack) {
      MPI_Comm_rank(MPI_COMM_WORLD, &rank);
      GetEigenVector(i, i_max, Z_vec, descZ_vec, vec_tmp);
      if (rank == 0) {
        for (j = 0; j < i_max; j++) {
          v0[j + 1] = vec_tmp[j];
        }
      }
      else {
        for (j = 0; j < i_max; j++) {
          v0[j + 1] = 0.0;
        }
      }
    }
    else {
      if (X->Def.iCalcType == FullDiag) {
        if (myrank == 0) {
          for (j = 0; j < i_max; j++) {
            v0[j + 1] = v1[i][j];
          }
        }
      }
      else {
        for (j = 0; j < i_max; j++) {
          v0[j + 1] = v1[i][j];
        }
      }
    }
  }/*for (i = 0; i < neig; i++)*/
#endif

  if (expec_energy_flct(X, neig, v0, v1) != 0) {
    fprintf(stderr, "Error: calc expec_energy.\n");
    exitMPI(-1);
  }
  if (expec_cisajs(X, neig, v0, v1) != 0) {
    fprintf(stderr, "Error: calc OneBodyG.\n");
    exitMPI(-1);
  }
  if (expec_cisajscktaltdc(X, neig, v0, v1) != 0) {
    fprintf(stderr, "Error: calc TwoBodyG.\n");
    exitMPI(-1);
  }
    
#ifdef _SCALAPACK
  if (use_scalapack) {
    if (X->Def.iCalcType == FullDiag) {
      X->Phys.s2 = 0.0;
      X->Phys.Sz = 0.0;
    }
  }
  else {
    if (X->Def.iCalcType == FullDiag) {
      if (expec_totalspin(X, v1) != 0) {
        fprintf(stderr, "Error: calc TotalSpin.\n");
        exitMPI(-1);
      }
    }
  }
#else
  if (X->Def.iCalcType == FullDiag) {
    if (expec_totalspin(X, neig, v1) != 0) {
      fprintf(stderr, "Error: calc TotalSpin.\n");
      exitMPI(-1);
    }
  }
#endif

  for (i = 0; i < neig; i++) {
    if (X->Def.iCalcModel == Spin || X->Def.iCalcModel == SpinGC) {
      tmp_N = X->Def.NsiteMPI;
    }
    else {
      tmp_N = X->Phys.num_up[i] + X->Phys.num_down[i];
    }
    if (X->Def.iCalcType == FullDiag) {
#ifdef _SCALAPACK
      if (use_scalapack) {
        fprintf(stdoutMPI, "i=%5ld Energy=%10lf N=%10lf Sz=%10lf Doublon=%10lf \n", i, X->Phys.energy, tmp_N,
          X->Phys.Sz, X->Phys.doublon);
      }
      else {
        fprintf(stdoutMPI, "i=%5ld Energy=%10lf N=%10lf Sz=%10lf S2=%10lf Doublon=%10lf \n", i, X->Phys.energy, tmp_N,
          X->Phys.Sz, X->Phys.s2, X->Phys.doublon);
      }
#else
      fprintf(stdoutMPI, "i=%5ld Energy=%10lf N=%10lf Sz=%10lf S2=%10lf Doublon=%10lf \n",
        i, X->Phys.energy[i], tmp_N, X->Phys.Sz[i], X->Phys.s2[i], X->Phys.doublon[i]);
#endif      
    }
    else if (X->Def.iCalcType == CG)
      fprintf(stdoutMPI, "i=%5ld Energy=%10lf N=%10lf Sz=%10lf Doublon=%10lf \n",
        i, X->Phys.energy[i], tmp_N, X->Phys.Sz[i], X->Phys.doublon[i]);
  }
#ifdef _SCALAPACK
  if(use_scalapack) free(vec_tmp);
#endif  
}