Perturbation.c File Reference

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

Function Documentation

main()

int main	(	int	argc,
		char **	argv
	)

                                {
  if(argc < 5) {
    printf("Usage: 1.This 2.n_max 3.num_h 4.num_op 5.Output dir\n");
    exit(1);
  }
 
  int i,j;
  int dim_tot, dim_sub, s; //全ヒルベルト空間の次元
  int dim_loc = 2; // サイト上のヒルベルト空間の次元
  int n,n_min;
  int n_max   = atoi(argv[1]); //  全スピンの数
  int num_h   = atoi(argv[2]); // 局所ハミルトニアンに含まれるサイト数
  int num_op  = atoi(argv[3]); // 局所演算子に含まれるサイト数
  int dloc_h  = pow(dim_loc, num_h); // 局所ハミルトニアンの空間の次元
  int dloc_op = pow(dim_loc, num_op); // 局所演算子の空間の次元
  int count, ndata; // カウンタ
  const int seed    = 0;
  int div;
  const int divMax = 100;
 
  int MinShell, MaxShell;
  double complex **h, **h_int, **loc;   // 局所ハミルトニアン、局所物理量
  double complex **h_tot, **loc_tot;   // 運動量セクターにおけるハミルトニアンと物理量
  double lambda, res, MicroAverage;
  double Stddev_inShell, Max_inShell;
  const double delta = 0.02; // ミクロカノニカル平均を計算するエネルギーシェルの幅
 
  double dE, Erel, range;
  double *eigenvalue, *expvalue;
  double complex expectation, *intermediate;
 
  //********** For LAPACK and BLASS
  const char e_value = 'V';
  const char upper = 'U';
  double complex *work;
  double *rwork;
  int info;
  const char trans = 'T', ntrans = 'N', ctrans = 'C';
  const double complex alpha = 1;
  const double complex beta  = 0;
  const int incx = 1;
  const int incy = 1;
  //********** (END) For LAPACK and BLASS
 
  const int BUFSIZE = 1000;
  char d_str[BUFSIZE], str[BUFSIZE], stemp[BUFSIZE];
  char *filename = str;
  FILE* FpOperators = NULL;
 
  clock_t start, end, temp_t;
 
  //******************** Check for the directory structure ********************
  strcpy(d_str, argv[4]);
  sprintf(stemp,"/Perturbation_Local%d_Seed%d_Max%d", num_h, seed, n_max);
  strcat(d_str, stemp);
  mkdir(d_str, 0744);
  //******************** (END) Check for the directory structure ********************
 
  //******************** Allocation & Initialization ********************
  if(num_op > num_h) n_min = num_op; else n_min = num_h;
  dim_tot = pow(dim_loc, 2*n_max);
  s       = 3*dim_tot;
 
  expvalue      = alloc_dvector(dim_tot);
  eigenvalue    = alloc_dvector(dim_tot);
  intermediate  = alloc_zvector(dim_tot);
 
  printf("dloc_h=%d, dloc_op=%d\n", dloc_h, dloc_op);
  h     = alloc_zmatrix(dloc_h, dloc_h);
  h_int = alloc_zmatrix(dloc_h, dloc_h);
  loc   = alloc_zmatrix(dloc_op, dloc_op);
  work  = alloc_zvector(s);
  rwork = alloc_dvector(s);
 
  init_genrand(seed);
  fflush(stdout);
  //******************** (END) Allocation & Initialization ********************
 
  start = clock();
  end   = start;
  //****************************** 局所演算子に関する性質を保存するファイルを用意 ******************************//
  strcpy(str, d_str);
  sprintf(stemp, "/Operators_%d.txt", seed);
  strcat(str, stemp);
  FpOperators = fopen(filename, "w");
  if(FpOperators == NULL) {
    printf("Can't open a file %s.", filename);
    exit(1);
  }
  fprintf(FpOperators, "#dim_loc=%d, num_h=%d, num_op=%d, delta=%lf\n", dim_loc, num_h, num_op, delta);
  fprintf(FpOperators, "#Fixed disorder strength\n");
  fprintf(FpOperators, "#1.lambda, 2.System size, 3.Energy density, 4.Relative Pos, 5.Exp value, 6.Micro average, 7.f(E,0) in shell, 8.MaxFluc in shell, 9.Dim of shell\n\n");
  //****************************** (END) 局所演算子に関する性質を保存するファイルを用意 ******************************//
 
  //****************************** 局所ハミルトニアンと局所物理量をランダムにとる ******************************//
  //printf("\rCreating Random matrices. dloc_h=%d, dloc_op=%d", dloc_h, dloc_op);
  //fflush(stdout);
  prepare_h(h, dloc_h , -1);
  prepare_h(loc, dloc_op, -1);
 
  lambda = 0;
  for(div = 0;div <= divMax; ++div) {
    lambda = (double)div/divMax;
    printf("lambda=%.2f\n", lambda);
    for(n = n_max;n >= n_min; --n) {
      dim_sub = pow(dim_loc, 2*n);
      h_tot   = alloc_zmatrix(dim_sub, dim_sub);
      loc_tot = alloc_zmatrix(dim_sub, dim_sub);
      for(i = 0;i < dim_sub; ++i) for(j = 0;j < dim_sub; ++j) {
        h_tot[i][j] = 0;
        loc_tot[i][j] = 0;
      }
 
      for(i = 0;i < dloc_h; ++i) for(j = 0;j < dloc_h; ++j) h_int[i][j] = lambda *h[i][j];
      //****************************** nにおける全系のハミルトニアンと演算子を構成する ******************************//
      // printf("n=%d\n", n);
      for(i = 0;i <= n-num_h; ++i) {
        imbed_locOp(i, h, h_tot, num_h, n, dim_sub, dim_loc);
        // printf("%d ", i);
      }
      for(i = n-num_h+1;i < n; ++i) {
        imbed_locOp(i, h_int, h_tot, num_h, n, dim_sub, dim_loc);
        // printf("%d ", i);
      }
      for(i = n;i <= 2*n-num_h; ++i) {
        imbed_locOp(i, h, h_tot, num_h, n, dim_sub, dim_loc);
        // printf("%d ", i);
      }
      // printf("\n");
      local_to_global(loc, loc_tot, num_op, n, dim_sub, dim_loc, "open", 0.0);
 
 
      //****************************** (END) nにおける全系のハミルトニアンと演算子を構成する ******************************//
 
      //****************************** 全系のハミルトニアンを対角化する ******************************//
      zheev_(&e_value, &upper, &dim_sub, &h_tot[0][0], &dim_sub, &eigenvalue[0], &work[0], &s, &rwork[0], &info);
      if(info != 0) {
        fprintf(stderr, "Error: LAPACK::zheev failed\n");
        exit(1);
      }
      for(i = 0;i < dim_sub; ++i) eigenvalue[i] /= (double)n;
      range = eigenvalue[dim_sub-1] -eigenvalue[0];
      dE    = delta*range;
      //****************************** (END)全系のハミルトニアンを対角化する ******************************//
 
      //****************************** 全系の物理量のエネルギー固有状態に関する期待値を計算 ******************************//
      for(i = 0;i < dim_sub; ++i) {
        for(j = 0;j < dim_sub; ++j) intermediate[j] = 0;
        expectation = 0;
        zgemv_(&ntrans, &dim_sub, &dim_sub, &alpha, &loc_tot[0][0], &dim_sub, &h_tot[i][0], &incx, &beta, &intermediate[0], &incy);
        // SOMETHING IS WRONG IN BLAS zdotc_ // zdotc_(&expectation, &dim_sub, &h_tot[i][0], &incx, intermediate, &incy);
        for(int p=0;p < dim_sub; ++p) expectation += conj(h_tot[i][p])*intermediate[p];
        expvalue[i] = creal(expectation)/(double)n;
      }
      //****************************** (END) 全系の物理量のエネルギー固有状態に関する期待値を計算 ******************************//
 
      //****************************** ミクロカノニカル平均まわりのゆらぎを計算 ******************************//
      ndata = 0;
      for(i = 0;i < dim_sub; ++i) {
        //********** ミクロカノニカル平均を計算
        MicroAverage = 0;
        count = 0;
        MinShell = dim_sub; MaxShell = 0;
        for(j = 0;j<dim_sub && eigenvalue[j]<eigenvalue[i]+dE; ++j) {
          if( eigenvalue[j] < eigenvalue[i]-dE ) continue;
          if (j < MinShell) MinShell = j;
          if (j > MaxShell) MaxShell = j;
          MicroAverage += expvalue[j];
          count += 1;
        }
        if (count!= MaxShell-MinShell+1) {
          printf("(n=%d, dim_sub=%d) count=%d, (MinShell, MaxShell) = (%d,%d), delta=%d\n", n, dim_sub, count, MinShell, MaxShell, MaxShell-MinShell+1);
          exit(1);
        }
        MicroAverage /= (double)count;
        //********** (END) ミクロカノニカル平均を計算
        Erel = (eigenvalue[i]-eigenvalue[0])/range;
        //********** シェル内のゆらぎの大きさを計算
        res = 0;
        Stddev_inShell = 0;
           Max_inShell = 0;
        for(j = MinShell;j <= MaxShell; ++j) {
          res = expvalue[j] -MicroAverage;
          Stddev_inShell += res*res;
          if( fabs(res) > Max_inShell ) Max_inShell = fabs(res);
        }
        Stddev_inShell /= (MaxShell-MinShell);
        Stddev_inShell  = sqrt(Stddev_inShell);
        //********** (END) シェル内のゆらぎの大きさを計算
        fprintf(FpOperators, "%le %d %le %le %le %le %le %le %d\n", lambda, n, eigenvalue[i], Erel, expvalue[i], MicroAverage, Stddev_inShell*sqrt((double)count), Max_inShell, count);
      }
      //****************************** (END) ミクロカノニカル平均まわりのゆらぎを計算 ******************************//
      fprintf(FpOperators, "\n");
      //******************** Free in for(n) ********************//
      free_zmatrix(h_tot);
      free_zmatrix(loc_tot);
      //******************** (END) Free in for(n) ********************//
    }
  }
 
  fclose(FpOperators);
  fflush(stdout);
 
  //******************** Free ********************//
  free_dvector(expvalue);
  free_dvector(eigenvalue);
  free_zvector(intermediate);
 
  free_zmatrix(h);
  free_zmatrix(loc);
  free_zvector(work);
  free_dvector(rwork);
  //******************** (END) Free ********************//
 
  return 0;
}