Fluc_randHOp.cpp File Reference

Functions
int	main (int argc, char **argv)

Function Documentation

main()

int main	(	int	argc,
		char **	argv
	)

                                {
  if(argc != N_ARGS) {
    std::cerr << USAGE << std::endl;
    std::exit(EX_USAGE);
  }
  Integer_t const MAX = 100000; // 1つのファイルに書き込むデータ数の上限
  Integer_t const dim_loc = 2; // サイト上のヒルベルト空間の次元
  Integer_t const div     = (Integer_t)std::atoll(argv[5]);
  Integer_t const num_h   = (Integer_t)std::atoll(argv[3]); // 局所ハミルトニアンに含まれるサイト数
  Integer_t const num_op  = (Integer_t)std::atoll(argv[4]); // 局所演算子に含まれるサイト数
  Integer_t const n_max   = (Integer_t)std::atoll(argv[1]); //  全スピンの数
  Integer_t const n_min   = (Integer_t)std::atoll(argv[2]);
  Integer_t const dloc_h  = pow(dim_loc, num_h); // 局所ハミルトニアンの空間の次元
  Integer_t const dloc_op = pow(dim_loc, num_op); // 局所演算子の空間の次元
  Integer_t const seed    = 0;
  Integer_t const repMin  = MAX*div;
  Integer_t const repMax  = repMin +MAX-1;
 
  Integer_t dim_sub; //全ヒルベルト空間の次元
 
  Integer_t ndata, info, failed=0; // カウンタ
  double gE, EnergyRange, OpRange, OpMin;
  double Fluc, MaxFluc, MaxPos; // ETHの指標
  double const cutoff = 0.45;
  double const dE     = 0.02; // ミクロカノニカル平均を計算するエネルギーシェルの幅
  double const Emin   = 0 +cutoff;
  double const Emax   = 1 -cutoff;
 
  std::stringstream buff;
 
  std::cout << "dloc_h=" << dloc_h << ", dloc_op=" << dloc_op << "\n";
  if(n_min > n_max) {
      std::cerr << "Error: n_min(" << n_min << ") > n_max(" << n_max << ")." << std::endl;
      std::exit(EX_USAGE);
  }
  if( !Initialize(argc, argv) ) {
    std::cerr << "Error: Initialization failed." << std::endl;
    std::exit(EX_USAGE);
  }
  #ifndef NDEBUG
    std::cerr << "# Successfully initialized." << std::endl;
  #endif
 
  //******************** Check for the directory structure ********************
  #ifndef NDEBUG
    std::cerr << "# Checking for the directory structure." << std::endl;
  #endif
  std::string d_str1(argv[6]); d_str1 += "/Probability/"; d_str1 += QUOTE(ENSEMBLE);
  std::string d_str2(d_str1);
  buff.str("");
  buff.clear(std::stringstream::goodbit);
  buff << "/" << TYPE << "randHOp_" << "Local" << num_h << "_Cut" << (Integer_t)(100*cutoff) << "_Seed" << seed << "_Min" << n_min << "_Max" << n_max;
  d_str1 += buff.str();
  std::filesystem::create_directories(d_str1);
 
  if(div == 0) {
    buff.str("");
    buff.clear(std::stringstream::goodbit);
    buff << "/" << TYPE << "Structure_" << "Local" << num_h << "_Seed" << seed << "_Min" << n_min << "_Max" << n_max;
    d_str2 += buff.str();
    std::filesystem::create_directories(d_str2);
  }
 
  std::ofstream OutMeasure[n_max+1];
  for(Integer_t n = n_min;n <= n_max; ++n) {
    buff.str("");
    buff.clear(std::stringstream::goodbit);
    buff << "/MeasureOfETH_N" << n << "_" << div << ".txt";
    std::string filename(d_str1); filename += buff.str();
    OutMeasure[n].open(filename);
    if( OutMeasure[n].fail() ) {
      std::cerr << "Can't open a file" << filename << "." << std::endl;
      std::exit(EX_CANTCREAT);
    }
    OutMeasure[n] << std::right;
  }
 
  std::ofstream OutSample;
  //******************** (END) Check for the directory structure ********************
 
  for(Integer_t n = n_max;n >= n_min; --n) {
    dim_sub = pow(dim_loc, n);
    std::cout << "(n=" << n << ") dim = " << dim_sub << "\n";
 
    OutMeasure[n] << "# (n_max=" << n_max << ") dim=" << dim_sub << "\n"
                  << "# 1.(Sample No.) 2.(Fluc_2) 3.(Fluc_infty) 4.(Pos_max_flux) 5.(ndata)\n"
                  << std::endl << std::scientific;
  }
 
  //***************************************************************************
  //******************** Allocation & Initialization **************************
  //***************************************************************************
  #ifndef NDEBUG
    std::cerr << "# Allocating CPU memories." << std::endl;
  #endif
  Integer_t const dim_max = pow(dim_loc, n_max);
  std::vector<Integer_t> NdataInShell(dim_max);
  std::vector<double>     eigenEnergy(dim_max);
  std::vector<double>        EXPvalue(dim_max);
  std::vector<double>           MCAverage(dim_max);
  matrix<MKL_Complex16>        h(dloc_h , dloc_h );
  matrix<MKL_Complex16>      loc(dloc_op, dloc_op);
  matrix<MKL_Complex16>    h_tot(dim_max, dim_max);
  matrix<MKL_Complex16>  loc_tot(dim_max, dim_max);
 
  double start, t_int, end, temp_t;
  double T_diag=0, T_post=0, T_pre=0;
  init_genrand(seed);
  start = getETtime();
  for(Integer_t repetition = 0;repetition < repMin; ++repetition) {
    generateLocal_h( h,   dloc_h,  -1);
    generateLocal_op(loc, dloc_op, -1);
  }
  end = getETtime();
  std::cout << "(init_genrand): time=" << std::fixed << (end-start) << std::endl;
  //***************************************************************************
  //******************** (END) Allocation & Initialization ********************
  //***************************************************************************
 
  start = getETtime();
  end   = start;
  for(Integer_t repetition = repMin;repetition <= repMax; ++repetition) {
    // 局所ハミルトニアンと局所物理量をランダムにとる ******************************//
    generateLocal_h(  h,  dloc_h,  -1);
    generateLocal_op(loc, dloc_op, -1);
 
    if(repetition < 10000) {
      //****************************** 局所演算子に関する性質を保存するファイルを用意 ******************************//
      buff.str("");
      buff.clear(std::stringstream::goodbit);
      buff << "/Sample_No" << repetition << ".txt";
      std::string filename(d_str2); filename += "/RawData";
      std::filesystem::create_directory(filename);
      filename += buff.str();
      OutSample.open(filename);
        if( OutSample.fail() ) {
          std::cerr << "Can't open a file" << filename;
          std::exit(EX_CANTCREAT);
        }
      OutSample << "# dim_loc=" << dim_loc << ", num_h=" << num_h << ", num_op=" << num_op << "\n"
                // << "# wave_number=" << wave_number << "\n"
                << "# 1.(System size) 2.(Energy density) 3.(Relative Pos) 4.(Exp. value)\n" << std::endl
                << std::scientific;
      //****************************** (END) 局所演算子に関する性質を保存するファイルを用意 ******************************//
    }
    for(Integer_t n = n_max;n >= n_min; --n) {
      dim_sub = pow(dim_loc, n);
 
      #ifndef NDEBUG
        std::cout << "# (rep,n)=(" << repetition << "," << n << ")" << std::endl;
      #endif
      temp_t  = getETtime();
      #pragma omp parallel
      {
        std::cout << "parallel(rep,n)=(" << repetition << "," << n << ")" << std::endl;
        #pragma omp sections
        {
          #pragma omp section
          {
            #pragma omp critical
            std::cout << "section1(rep,n)=(" << repetition << "," << n << ")" << std::endl;
            constructGlobal_h(   h_tot, n, dim_loc, num_h,    h);
          }
          #pragma omp section
          {
            #pragma omp critical
            std::cout << "section2(rep,n)=(" << repetition << "," << n << ")" << std::endl;
            constructGlobal_op(loc_tot, n, dim_loc, num_op, loc);
          }
        }
      }
      T_pre += getETtime() -temp_t;
 
      #ifndef NDEBUG
        std::cerr << "# Calculating eigenstate expectation values." << std::endl;
      #endif
      temp_t = getETtime();
        info = EigenExpValue(EXPvalue, eigenEnergy, dim_sub, h_tot, loc_tot);
          if(info != 0) {
            failed += 1;
            buff.str("");
            buff.clear(std::stringstream::goodbit);
            buff << "/Failed_N" << n << "_No" << repetition << ".txt";
            std::string filename(d_str1); filename += buff.str();
            std::ofstream ErrFp;
            ErrFp.open(filename);
            if( !ErrFp.is_open() ) {
              std::cerr << "# Warning: Can't open a file " << filename << ". Continue..." << std::endl;
            } else {
              ErrFp << "# info=" << info << "\n\n"
                       "# seed=" << seed << ", div=" << div << ", N=" << n << ", repetition=" << repetition << "\n"
                       "# Local Hamiltonian\n";
              h.print(ErrFp, dloc_h, dloc_h); ErrFp << "\n";
              ErrFp << "# Local operator\n";
              loc.print(ErrFp, dloc_op, dloc_op); ErrFp << "\n";
              ErrFp.close();
            }
            continue;
          }
        OpRange = SpectralRange(OpMin, dim_sub, loc_tot);
      T_diag += getETtime() -temp_t;
 
      #ifndef NDEBUG
        std::cerr << "# Calculating measures of the ETH." << std::endl;
      #endif
      temp_t  = getETtime();
        for(Integer_t i = 0;i < dim_sub; ++i) EXPvalue[i] /= OpRange;
        EnergyRange = eigenEnergy[dim_sub-1] -eigenEnergy[0];
        gE = eigenEnergy[0];
        for(Integer_t i = 0;i < dim_sub; ++i) eigenEnergy[i] = (eigenEnergy[i] -gE)/EnergyRange;
        // Calculate Microcanonical averages with energy corresponding to each eigenstate.
        MCAverage.resize(dim); NdataInShell.resize(dim);
        std::fill(   MCAverage.begin(),    MCAverage.end(), 0.0);
        std::fill(NdataInShell.begin(), NdataInShell.end(), 0.0);
        for(Integer_t i = 0;i < dim; ++i) {
          MCAverage[i] = MicroCanonicalAverage(NdataInShell[i], eigenEnergy[i], dE, dim_sub, eigenEnergy, EXPvalue, i);
        }
 
        ndata = MeasureOfETH(Fluc, MaxFluc, MaxPos, Emin, Emax, dim_sub, eigenEnergy, EXPvalue, MCAverage);
      T_post += getETtime() -temp_t;
 
      OutMeasure[n] << std::setw(6) << repetition << " " << std::showpos
                    << std::setw(13) << Fluc << " "
                    << std::setw(13) << MaxFluc << " "
                    << std::setw(13) << MaxPos << " " << std::noshowpos << ndata << "\n";
 
      if(repetition < 10000) {
        for(Integer_t i = 0;i < dim_sub; ++i) {
          OutSample << std::setw(2) << n << " " << std::showpos
                    << std::setw(13) << eigenEnergy[i]*(double)EnergyRange +(double)gE << " "
                    << std::setw(13) << eigenEnergy[i] << " "
                    << std::setw(13) << EXPvalue[i]    << "\n" << std::noshowpos;
        }
      }
    }
    if(repetition < 10000) OutSample.close();
    if(repetition%10 == 9) {
      t_int = end; end = getETtime();
      std::cerr << "(total=" << repetition+1
                << "): timeINT="    << std::setprecision(6) << (end-t_int)
                << ", timeTOT="     << std::setprecision(6) << (end-start)
                << ", T_construct=" << std::setprecision(6) << T_pre  << "(" << std::setprecision(1) << 100*T_pre /(end-start) << "%)"
                << ", T_diag="      << std::setprecision(6) << T_diag << "(" << std::setprecision(1) << 100*T_diag/(end-start) << "%)"
                << ", T_process="   << std::setprecision(6) << T_post << "(" << std::setprecision(1) << 100*T_post/(end-start) << "%)"
                << std::endl;
      for(Integer_t n = n_min;n <= n_max; ++n) OutMeasure[n].flush();
    }
  }
 
  Finalize(argc, argv);
  return 0;
}