PBC.c File Reference

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

Function Documentation

getETtime()

double getETtime

(

)

                   {
  struct timeval tv;
  gettimeofday(&tv, NULL);
  return tv.tv_sec +(double)tv.tv_usec *1e-6;
}

main()

int main	(	int	argc,
		char **	argv
	)

                                {
  if(argc != 7) {
    fprintf(stderr, "Usage: 1.This 2.(Number of elements) 3.(Number of blocks) 4.(Max elements on a site)  5.|U|*100 6.|V|*100 7.OutDir\n");
    exit(EX_USAGE);
  }
  const MKL_INT N = atoi(argv[1]);
  const MKL_INT L = atoi(argv[2]);
  const MKL_INT MaxOnSite = atoi(argv[3]);
  const double U = -atoi(argv[4])/100.0;
  const double V = -atoi(argv[5])/100.0;
 
  MKL_INT dim_tot, dim_odd, dim_even, NumPairs, **basis=NULL, *order=NULL;
  double *eigenEnergy;
  MKL_Complex16 **Hamiltonian_Odd, **Hamiltonian_Even;
 
  const MKL_INT BUFSIZE = 1000;
  char dstr[BUFSIZE], str[BUFSIZE], stemp[BUFSIZE];
  char *filename = str;
  double start, end;
 
  //******************** Check for the directory structure ********************
  strcpy(dstr, argv[6]);
  strcat(dstr, "/BICData"); mkdir(dstr, 0744);
  strcat(dstr, "/PBC");     mkdir(dstr, 0744);
  if(U > 0) {
    sprintf(stemp, "/N%lld_U%lld_nV%lld_Max%lld_L%lld", N, (MKL_INT)(100*U), (MKL_INT)(-100*V), MaxOnSite, L);
  } else {
    sprintf(stemp, "/N%lld_nU%lld_nV%lld_Max%lld_L%lld", N, (MKL_INT)(-100*U), (MKL_INT)(-100*V), MaxOnSite, L);
  }
  strcat(dstr, stemp);  mkdir(dstr, 0744);
  //******************** (END)Check for the directory structure ********************
 
  //******************** Allocation ********************//
  fprintf(stdout, "# Generating basis vectors...\n");
  dim_tot  = alloc_BasisVectors(N, L, MaxOnSite, &basis, &order, &NumPairs);
  dim_odd  = NumPairs;
  dim_even = dim_tot -dim_odd;
  fprintf(stdout, "# dim_tot = %lld, dim_even = %lld, dim_odd = %lld\n", dim_tot, dim_even, dim_odd);
  eigenEnergy = alloc_dvector(dim_tot);
  Hamiltonian_Odd  = alloc_zmatrix(dim_odd, dim_odd);
  Hamiltonian_Even = alloc_zmatrix(dim_even, dim_even);
  //******************** (END)Allocation ********************//
 
  //****************************************//
  //****************************************//
  fprintf(stdout, "# Constructing Hamiltonian...\n");
  Construct_Hamiltonian(Hamiltonian_Odd, dim_odd, Hamiltonian_Even, dim_even, L, basis, U, V);
 
  start = getETtime();
    fprintf(stdout, "# Diagonalizing Hamiltonian...\n"); fflush(stdout);
    Diagonalize_zmatrix(dim_even, Hamiltonian_Even, &(eigenEnergy[0]));
    Diagonalize_zmatrix(dim_odd,  Hamiltonian_Odd,  &(eigenEnergy[dim_even]));
  end = getETtime();
  fprintf(stderr, "# (%lld %lld), Real time(sec): %.6lf sec\n", N, L, end-start);
  fflush(stdout); fflush(stderr);
 
  //************************************************************//
  //************************************************************//
  //************************************************************//
  start = getETtime();
  FILE *Outfp;
  double **Observable, *KineticEnergy, *ImpurityEnergy;
  KineticEnergy     = alloc_dvector(dim_tot);
  ImpurityEnergy    = alloc_dvector(dim_tot);
  for(int i = 0;i < dim_tot; ++i) {
    KineticEnergy[i]     = 0;
    ImpurityEnergy[i]    = 0;
  }
 
  fprintf(stdout, "# Calculating probability distribution...\n"); fflush(stdout);
  Observable = Probability(Hamiltonian_Even, Hamiltonian_Odd, N, L, basis, dim_tot, NumPairs);
  strcpy(str, dstr);
  sprintf(stemp, "/Probability.txt");
  strcat(str, stemp);
  Outfp = fopen(filename, "w");
    if(Outfp == NULL) {
      printf("Can't open a file %s.", filename);
      exit(EX_CANTCREAT);
    }
  fprintf(Outfp, "# N=%lld, L=%lld, MaxOnSite=%lld\n", N, L, MaxOnSite);
  fprintf(Outfp, "# U=%+lf, V=%+lf\n", U, V);
  for(int alpha = 0;alpha < dim_tot; ++alpha) {
    double Cumulative = 0;
    for(int j = 0;j < L/2; ++j) {
      Cumulative += Observable[alpha][j];
      fprintf(Outfp, "%d %+lld %le %d %le %le\n", alpha, 1-2*(alpha/dim_even), eigenEnergy[alpha], j, Observable[alpha][j], 1-Cumulative);
    }
    fprintf(Outfp, "\n");
  }
  free_dmatrix(Observable); Observable = NULL;
 
 
  fprintf(stdout, "# Calculating particle density distribution...\n"); fflush(stdout);
  Observable = ParticleDensity(Hamiltonian_Even, Hamiltonian_Odd, N, L, basis, dim_tot, NumPairs);
  strcpy(str, dstr);
  sprintf(stemp, "/ParticleDensity.txt");
  strcat(str, stemp);
  Outfp = fopen(filename, "w");
    if(Outfp == NULL) {
      printf("Can't open a file %s.", filename);
      exit(EX_CANTCREAT);
    }
  fprintf(Outfp, "# N=%lld, L=%lld, MaxOnSite=%lld\n", N, L, MaxOnSite);
  fprintf(Outfp, "# U=%+lf, V=%+lf\n", U, V);
  for(int alpha = 0;alpha < dim_tot; ++alpha) {
    for(int j = 0;j < L; ++j) fprintf(Outfp, "%d %+lld %le %d %le\n", alpha, 1-2*(alpha/dim_even), eigenEnergy[alpha], j, Observable[alpha][j]);
    fprintf(Outfp, "\n");
    ImpurityEnergy[alpha] = V*N*Observable[alpha][L/2];
  }
  free_dmatrix(Observable); Observable = NULL;
 
  fprintf(stdout, "# Calculating pseudo-momentum distribution...\n"); fflush(stdout);
  Observable = MomentumDistribution(Hamiltonian_Even, Hamiltonian_Odd, N, L, basis, order, dim_tot, NumPairs);
  strcpy(str, dstr);
  sprintf(stemp, "/MomentumDistribution.txt");
  strcat(str, stemp);
  Outfp = fopen(filename, "w");
    if(Outfp == NULL) {
      printf("Can't open a file %s.", filename);
      exit(EX_CANTCREAT);
    }
  fprintf(Outfp, "# N=%lld, L=%lld, MaxOnSite=%lld\n", N, L, MaxOnSite);
  fprintf(Outfp, "# U=%+lf, V=%+lf\n", U, V);
  for(int alpha = 0;alpha < dim_tot; ++alpha) {
    for(int j = 0;j < L; ++j) {
      fprintf(Outfp, "%d %+lld %le %d %le\n", alpha, 1-2*(alpha/dim_even), eigenEnergy[alpha], j, Observable[alpha][j]);
      KineticEnergy[alpha] += -2*cos(2*M_PI*j/L)*Observable[alpha][j]*N/(double)L;
    }
    fprintf(Outfp, "\n");
  }
 
 
  strcpy(str, dstr);
  sprintf(stemp, "/Energy.txt");
  strcat(str, stemp);
  Outfp = fopen(filename, "w");
    if(Outfp == NULL) {
      printf("Can't open a file %s.", filename);
      exit(EX_CANTCREAT);
    }
  fprintf(Outfp, "# N=%lld, L=%lld, MaxOnSite=%lld\n", N, L, MaxOnSite);
  fprintf(Outfp, "# U=%+lf, V=%+lf\n", U, V);
  fprintf(Outfp, "# 1.No 2.Parity 3.EigenEnergy 4.KineticEnergy 5.InteractionEnergy 6.ImpurityEnergy\n");
  for(int alpha = 0;alpha < dim_tot; ++alpha) {
    double InteractionEnergy = eigenEnergy[alpha] -(KineticEnergy[alpha] +ImpurityEnergy[alpha]);
    fprintf(Outfp, "%d %+lld %le %le %le %le\n", alpha, 1-2*(alpha/dim_even), eigenEnergy[alpha], KineticEnergy[alpha], InteractionEnergy, ImpurityEnergy[alpha]);
  }
 
  //******************** Free ********************//
  free_dvector(eigenEnergy);
  free_zmatrix(Hamiltonian_Odd);
  free_zmatrix(Hamiltonian_Even);
  free_BasisVectors(basis, order);
  //******************** (END)Free ********************//
 
  return 0;
}