diff --git a/tools/phonon/Makefile b/tools/phonon/Makefile
index 6c39f55ce..0aacb1e08 100644
--- a/tools/phonon/Makefile
+++ b/tools/phonon/Makefile
@@ -1,64 +1,61 @@
 .SUFFIXES : .o .cpp
 # compiler and flags
 CC     = g++ -Wno-unused-result
 LINK   = $(CC) -static
 CFLAGS = -O3 $(DEBUG) $(UFLAG)
 #
 OFLAGS = -O3 $(DEBUG)
 INC    = $(LPKINC) $(TCINC) $(SPGINC)
 LIB    = $(LPKLIB) $(TCLIB) $(SPGLIB)
 #
 # cLapack library needed
-LPKINC = -I/opt/clapack/3.2.1/include
-LPKLIB = -L/opt/clapack/3.2.1/lib -lclapack -lblas -lf2c #-lm
+LPKINC = -I/opt/libs/clapack/3.2.1/include
+LPKLIB = -L/opt/libs/clapack/3.2.1/lib -lclapack -lblas -lf2c #-lm
 #
 # Tricubic library needed
-TCINC = -I/opt/tricubic/1.0/include
-TCLIB = -L/opt/tricubic/1.0/lib -ltricubic
+TCINC = -I/opt/libs/tricubic/1.0/include
+TCLIB = -L/opt/libs/tricubic/1.0/lib -ltricubic
 #
-# spglib 0.7.1, used to get the irreducible q-points
+# spglib 1.8.2, used to get the irreducible q-points
 # if UFLAG is not set, spglib won't be used.
 UFLAG  = -DUseSPG
-SPGINC = -I/opt/spglib/0.7.1/include
-SPGLIB = -L/opt/spglib/0.7.1/lib -lsymspg
-# if spglib > 0.7.1 is used, please 
-# 1) modify file phonon.cpp, instruction can be found by searching 0.7.1
-# 2) uncomment the following two lines
-#SPGINC = -I/opt/spglib/1.1.2/include
-#SPGLIB = -L/opt/spglib/1.1.2/lib -lsymspg
+SPGINC = -I/opt/libs/spglib/1.8.2/include
+SPGLIB = -L/opt/libs/spglib/1.8.2/lib -lsymspg
+# if spglib other than version 1.8.2 is used, please 
+# modify file phonon.cpp, instruction can be found by searching 1.8.2
 
 # Debug flags
 #DEBUG = -g -DDEBUG
 #====================================================================
 ROOT   = phana
 # executable name
 EXE    = $(ROOT)
 #====================================================================
 # source and rules
 SRC = $(wildcard *.cpp)
 OBJ = $(SRC:.cpp=.o)
 
 #====================================================================
-all:  ${EXE}
+all:  ver ${EXE}
 
 ${EXE}: $(OBJ)
 	$(LINK) $(OFLAGS) $(OBJ) $(LIB) -o $@
 
 clean: 
 	rm -f *.o *~ *.mod ${EXE}
 
 tar:
 	rm -f ${ROOT}.tar; tar -czvf ${ROOT}.tar.gz *.cpp  *.h Makefile README
 
 ver:
-	@echo "#define VERSION `svn info|grep '^Revision'|cut -d: -f2`" > version.h
+	@echo "#define VERSION `git log|grep '^commit'|wc -l`" > version.h
 
 #====================================================================
 .f.o:
 	$(FC) $(FFLAGS) $(FREE) $(MPI) ${INC} -c $<
 .f90.o:
 	$(FC) $(FFLAGS) $(FREE) $(MPI) ${INC} -c $<
 .c.o:
 	$(CC) $(CFLAGS) -c $<
 .cpp.o:
 	$(CC) $(CFLAGS) $(INC) -c $<
diff --git a/tools/phonon/README b/tools/phonon/README
index 61e458e5f..ae6383b6b 100644
--- a/tools/phonon/README
+++ b/tools/phonon/README
@@ -1,31 +1,48 @@
-phana = post-processing program for the fix phonon command
-
-This program reads the binary file created by the fix phonon command
-and helps to analyse the phonon related info.
-
-The clapack library is needed to solve the eigen problems, which could
-be downloaded from: http://www.netlib.org/clapack/
-
-The tricubic library is also needed to to tricubic interpolations,
-which could be obtained from:
-http://orca.princeton.edu/francois/software/tricubic/ or
-http://code.google.com/p/fix-phonon/downloads/list
-
-The spglib (version 0.7.1) is optionally needed, enabling one to
-evaluate the phonon density of states or vibrational thermal
-properties using only the irreducible q-points in the first Brillouin
-zone.
-
-To compile the code, one needs therefore to install the above
-libraries and set the paths correctly in the Makefile.
-
-The units of the output frequencies by this code is THz for LAMMPS
-units "real", "si", "metal", and "cgs"; in these cases, the
-frequencies are $\nu$ instead of $\omega$.
-
-One is encouraged to visit http://code.google.com/p/fix-phonon/ to
-check out the latest revision on fix-phonon and the post-processing
-code.
-
-Author: Ling-Ti Kong (konglt at sjtu.edu.cn)
-Feb 2013
+#-------------------------------------------------------------------------------
+   phana
+#
+   This program reads the binary file created by fix_phonon and helps to
+   analyse the phonon related information.
+#-------------------------------------------------------------------------------
+1. Dependencies
+   The clapack library is needed to solve the eigen problems,
+   which could be downloaded from:
+   http://www.netlib.org/clapack/
+   
+   The tricubic library is also needed to do tricubic interpolations,
+   which could be obtained from:
+      http://orca.princeton.edu/francois/software/tricubic/
+   or
+      http://1drv.ms/1J2WFYk
+   
+   The spglib is optionally needed, enabling one to evaluate the
+   phonon density of states or vibrational thermal properties
+   using only the irreducible q-points in the first Brillouin zone,
+   as well as to evaluate the phonon dispersion curvers with the
+   automatic mode. Currently, the 1.8.3 version of spglib is used.
+   It can be obtained from:
+   http://spglib.sourceforge.net/
+   
+2. Compilation
+   To compile the code, one needs therefore to install the above
+   libraries and set the paths correctly in the Makefile.
+   Once this is done, by typing
+     make
+   will yield the executable "phana".
+
+3. Unit system
+   The units of the output frequencies by this code is THz for
+   LAMMPS units "real", "si", "metal", and "cgs"; in these cases,
+   the frequencies are $\nu$ instead of $\omega$.
+
+4. Updates
+   For updates of phana, please check:
+     http://nes.sjtu.edu.cn/english/research/software.htm
+   or
+     https://github.com/lingtikong/phana.git
+
+5. Bug report
+   If any bug found, please drop a line to: konglt(at)sjtu.edu.cn
+#-------------------------------------------------------------------------------
+Author: Ling-Ti Kong, konglt(at)sjtu.edu.cn
+Oct 2015
diff --git a/tools/phonon/dynmat.cpp b/tools/phonon/dynmat.cpp
index 188690164..eb91647ed 100644
--- a/tools/phonon/dynmat.cpp
+++ b/tools/phonon/dynmat.cpp
@@ -1,602 +1,589 @@
 #include "dynmat.h"
 #include "math.h"
 #include "version.h"
-
-#define MAXLINE 256
+#include "global.h"
 
 // to intialize the class
 DynMat::DynMat(int narg, char **arg)
 {
   attyp = NULL;
   memory = NULL;
   M_inv_sqrt = NULL;
   interpolate = NULL;
   DM_q = DM_all = NULL;
   binfile = funit = dmfile = NULL;
 
+  attyp = NULL;
+  basis = NULL;
   flag_reset_gamma = flag_skip = 0;
 
   // analyze the command line options
   int iarg = 1;
   while (narg > iarg){
     if (strcmp(arg[iarg], "-s") == 0){
       flag_reset_gamma = flag_skip = 1;
 
     } else if (strcmp(arg[iarg], "-r") == 0){
       flag_reset_gamma = 1;
 
     } else if (strcmp(arg[iarg], "-h") == 0){
       help();
 
     } else {
       if (binfile) delete []binfile;
       int n = strlen(arg[iarg]) + 1;
       binfile = new char[n];
       strcpy(binfile, arg[iarg]); 
     }
 
     iarg++;
   }
 
   ShowVersion();
   // get the binary file name from user input if not found in command line
   char str[MAXLINE];
   if (binfile == NULL) {
     char *ptr = NULL;
     printf("\n");
     do {
       printf("Please input the binary file name from fix_phonon: ");
       fgets(str,MAXLINE,stdin);
       ptr = strtok(str, " \n\t\r\f");
     } while (ptr == NULL);
 
     int n = strlen(ptr) + 1;
     binfile = new char[n];
     strcpy(binfile, ptr);
   }
 
   // open the binary file
   FILE *fp = fopen(binfile, "rb");
   if (fp == NULL) {
     printf("\nFile %s not found! Programe terminated.\n", binfile);
     help();
   }
 
   // read header info from the binary file
   if ( fread(&sysdim, sizeof(int),    1, fp) != 1) {printf("\nError while reading sysdim from file: %s\n", binfile); fclose(fp); exit(2);}
   if ( fread(&nx,     sizeof(int),    1, fp) != 1) {printf("\nError while reading nx from file: %s\n", binfile); fclose(fp); exit(2);}
   if ( fread(&ny,     sizeof(int),    1, fp) != 1) {printf("\nError while reading ny from file: %s\n", binfile); fclose(fp); exit(2);}
   if ( fread(&nz,     sizeof(int),    1, fp) != 1) {printf("\nError while reading nz from file: %s\n", binfile); fclose(fp); exit(2);}
   if ( fread(&nucell, sizeof(int),    1, fp) != 1) {printf("\nError while reading nucell from file: %s\n", binfile); fclose(fp); exit(2);}
   if ( fread(&boltz,  sizeof(double), 1, fp) != 1) {printf("\nError while reading boltz from file: %s\n", binfile); fclose(fp); exit(2);}
 
   fftdim = sysdim*nucell; fftdim2 = fftdim*fftdim;
   npt = nx*ny*nz;
 
   // display info related to the read file
-  printf("\n"); for (int i=0; i<80; i++) printf("="); printf("\n");
+  printf("\n"); for (int i = 0; i < 80; ++i) printf("="); printf("\n");
   printf("Dynamical matrix is read from file: %s\n", binfile);
   printf("The system size in three dimension: %d x %d x %d\n", nx, ny, nz);
   printf("Number of atoms per unit cell     : %d\n", nucell);
   printf("System dimension                  : %d\n", sysdim);
   printf("Boltzmann constant in used units  : %g\n", boltz);
-  for (int i=0; i<80; i++) printf("="); printf("\n");
-  if (sysdim<1||sysdim>3||nx<1||ny<1||nz<1||nucell<1){
+  for (int i = 0; i < 80; ++i) printf("="); printf("\n");
+  if (sysdim < 1||sysdim > 3||nx < 1||ny < 1||nz < 1||nucell < 1){
     printf("Wrong values read from header of file: %s, please check the binary file!\n", binfile);
     fclose(fp); exit(3);
   }
 
   funit = new char[4];
   strcpy(funit, "THz");
-  if (boltz == 1.){eml2f = 1.; delete funit; funit=new char[22]; strcpy(funit,"sqrt(epsilon/(m.sigma^2))");}
-  else if (boltz == 0.0019872067) eml2f = 3.256576161;
-  else if (boltz == 8.617343e-5)  eml2f = 15.63312493;
+  if (boltz == 1.){eml2f = 1.; delete funit; funit = new char[27]; strcpy(funit,"sqrt(epsilon/(m.sigma^2))");}
+  else if (boltz == 0.0019872067)  eml2f = 3.256576161;
+  else if (boltz == 8.617343e-5)   eml2f = 15.63312493;
   else if (boltz == 1.3806504e-23) eml2f = 1.;
   else if (boltz == 1.3806504e-16) eml2f = 1.591549431e-14;
   else {
     printf("WARNING: Because of float precision, I cannot get the factor to convert sqrt(E/ML^2)\n");
     printf("into THz, instead, I set it to be 1; you should check the unit used by LAMMPS.\n");
     eml2f = 1.;
   }
 
   // now to allocate memory for DM
-  memory = new Memory;
-  DM_all = memory->create(DM_all, npt, fftdim2, "DynMat:DM_all");
-  DM_q   = memory->create(DM_q, fftdim,fftdim,"DynMat:DM_q");
+  memory = new Memory();
+  memory->create(DM_all, npt, fftdim2, "DynMat:DM_all");
+  memory->create(DM_q, fftdim,fftdim,"DynMat:DM_q");
 
   // read all dynamical matrix info into DM_all
   if ( fread(DM_all[0], sizeof(doublecomplex), npt*fftdim2, fp) != size_t(npt*fftdim2)){
     printf("\nError while reading the DM from file: %s\n", binfile);
     fclose(fp);
     exit(1);
   }
 
   // now try to read unit cell info from the binary file
-  flag_latinfo = 0;
-  basis = memory->create(basis,nucell,sysdim,"DynMat:basis");
-  attyp = memory->create(attyp,nucell, "DynMat:attyp");
-  M_inv_sqrt = memory->create(M_inv_sqrt, nucell, "DynMat:M_inv_sqrt");
-  int flag_mass_read = 0;
+  memory->create(basis, nucell, sysdim, "DynMat:basis");
+  memory->create(attyp, nucell,         "DynMat:attyp");
+  memory->create(M_inv_sqrt, nucell,    "DynMat:M_inv_sqrt");
   
-  if ( fread(&Tmeasure,   sizeof(double), 1, fp) == 1) flag_latinfo |= 1;
-  if ( fread(&basevec[0], sizeof(double), 9, fp) == 9) flag_latinfo |= 2;
-  if ( fread(basis[0],    sizeof(double), fftdim, fp) == fftdim) flag_latinfo |= 4;
-  if ( fread(&attyp[0],   sizeof(int),    nucell, fp) == nucell) flag_latinfo |= 8;
-  if ( fread(&M_inv_sqrt[0], sizeof(double), nucell, fp) == nucell) flag_mass_read = 1;
+  if ( fread(&Tmeasure,      sizeof(double), 1,      fp) != 1     ){printf("\nError while reading temperature from file: %s\n",   binfile); fclose(fp); exit(3);}
+  if ( fread(&basevec[0],    sizeof(double), 9,      fp) != 9     ){printf("\nError while reading lattice info from file: %s\n",  binfile); fclose(fp); exit(3);}
+  if ( fread(basis[0],       sizeof(double), fftdim, fp) != fftdim){printf("\nError while reading basis info from file: %s\n",    binfile); fclose(fp); exit(3);}
+  if ( fread(&attyp[0],      sizeof(int),    nucell, fp) != nucell){printf("\nError while reading atom types from file: %s\n",    binfile); fclose(fp); exit(3);}
+  if ( fread(&M_inv_sqrt[0], sizeof(double), nucell, fp) != nucell){printf("\nError while reading atomic masses from file: %s\n", binfile); fclose(fp); exit(3);}
   fclose(fp);
 
-  if ((flag_latinfo&15) == 15){
-    car2dir(flag_mass_read);
-    real2rec();
-
-    flag_latinfo = 1;
-  } else {
-    Tmeasure = 0.;
-    flag_latinfo = 0;
-  }
+  car2dir();
+  real2rec();
 
   // initialize interpolation
   interpolate = new Interpolate(nx,ny,nz,fftdim2,DM_all);
   if (flag_reset_gamma) interpolate->reset_gamma();
 
-  if ( flag_mass_read ){ // M_inv_sqrt info read, the data stored are force constant matrix instead of dynamical matrix.
-
-    EnforceASR();
-
-    // get the dynamical matrix from force constant matrix: D = 1/M x Phi
-    for (int idq=0; idq< npt; idq++){
-      int ndim =0;
-      for (int idim=0; idim<fftdim; idim++){
-        for (int jdim=0; jdim<fftdim; jdim++){
-          double inv_mass = M_inv_sqrt[idim/sysdim]*M_inv_sqrt[jdim/sysdim];
-          DM_all[idq][ndim].r *= inv_mass;
-          DM_all[idq][ndim].i *= inv_mass;
-          ndim++;
-        }
-      }
+  // Enforcing Austic Sum Rule
+  EnforceASR();
+
+  // get the dynamical matrix from force constant matrix: D = 1/M x Phi
+  for (int idq = 0; idq < npt; ++idq){
+    int ndim =0;
+    for (int idim = 0; idim < fftdim; ++idim)
+    for (int jdim = 0; jdim < fftdim; ++jdim){
+      double inv_mass = M_inv_sqrt[idim/sysdim]*M_inv_sqrt[jdim/sysdim];
+      DM_all[idq][ndim].r *= inv_mass;
+      DM_all[idq][ndim].i *= inv_mass;
+      ndim++;
     }
   }
 
   // ask for the interpolation method
   interpolate->set_method();
 
   return;
 }
 
 // to destroy the class
 DynMat::~DynMat()
 {
  // destroy all memory allocated
  if (funit) delete []funit;
  if (dmfile) delete []dmfile;
  if (binfile) delete []binfile;
  if (interpolate) delete interpolate;
 
  memory->destroy(DM_q);
  memory->destroy(attyp);
  memory->destroy(basis);
  memory->destroy(DM_all);
  memory->destroy(M_inv_sqrt);
  if (memory) delete memory;
 }
 
 /* ----------------------------------------------------------------------------
  * method to write DM_q to file, single point
  * ---------------------------------------------------------------------------- */
 void DynMat::writeDMq(double *q)
 {
   FILE *fp;
   // only ask for file name for the first time
   // other calls will append the result to the file.
   if (dmfile == NULL){
     char str[MAXLINE], *ptr;
     printf("\n");
-    while (1){
+    while ( 1 ){
       printf("Please input the filename to output the DM at selected q: ");
       fgets(str,MAXLINE,stdin);
       ptr = strtok(str, " \r\t\n\f");
       if (ptr) break;
     }
 
     int n = strlen(ptr) + 1;
     dmfile = new char[n];
     strcpy(dmfile, ptr);
     fp = fopen(dmfile,"w");
 
   } else {
     fp = fopen(dmfile,"a");
   }
   fprintf(fp,"# q = [%lg %lg %lg]\n", q[0], q[1], q[2]);
 
-  for (int i=0; i<fftdim; i++){
-    for (int j=0; j<fftdim; j++) fprintf(fp,"%lg %lg\t", DM_q[i][j].r, DM_q[i][j].i);
+  for (int i = 0; i < fftdim; ++i){
+    for (int j = 0; j < fftdim; ++j) fprintf(fp,"%lg %lg\t", DM_q[i][j].r, DM_q[i][j].i);
     fprintf(fp,"\n");
   }
   fprintf(fp,"\n");
   fclose(fp);
 return;
 }
 
 /* ----------------------------------------------------------------------------
  * method to write DM_q to file, dispersion-like
  * ---------------------------------------------------------------------------- */
 void DynMat::writeDMq(double *q, const double qr, FILE *fp)
 {
 
   fprintf(fp, "%lg %lg %lg %lg ", q[0], q[1], q[2], qr);
 
-  for (int i=0; i<fftdim; i++){
-    for (int j=0; j<fftdim; j++) fprintf(fp,"%lg %lg\t", DM_q[i][j].r, DM_q[i][j].i);
-  }
+  for (int i = 0; i < fftdim; ++i)
+  for (int j = 0; j < fftdim; ++j) fprintf(fp,"%lg %lg\t", DM_q[i][j].r, DM_q[i][j].i);
+
   fprintf(fp,"\n");
 return;
 }
 
 /* ----------------------------------------------------------------------------
  * method to evaluate the eigenvalues of current q-point;
  * return the eigenvalues in egv.
  * cLapack subroutine zheevd is employed.
  * ---------------------------------------------------------------------------- */
 int DynMat::geteigen(double *egv, int flag)
 {
   char jobz, uplo;
   integer n, lda, lwork, lrwork, *iwork, liwork, info;
   doublecomplex *work;
   doublereal *w = &egv[0], *rwork;
 
   n     = fftdim;
   if (flag) jobz = 'V';
   else jobz = 'N';
 
   uplo = 'U';
   lwork = (n+2)*n;
   lrwork = 1 + (5+n+n)*n;
   liwork = 3 + 5*n;
   lda    = n;
 
-  work  = memory->create(work,  lwork,  "geteigen:work");
-  rwork = memory->create(rwork, lrwork, "geteigen:rwork");
-  iwork = memory->create(iwork, liwork, "geteigen:iwork");
+  memory->create(work,  lwork,  "geteigen:work");
+  memory->create(rwork, lrwork, "geteigen:rwork");
+  memory->create(iwork, liwork, "geteigen:iwork");
 
   zheevd_(&jobz, &uplo, &n, DM_q[0], &lda, w, work, &lwork, rwork, &lrwork, iwork, &liwork, &info);
  
   // to get w instead of w^2; and convert w into v (THz hopefully)
-  for (int i=0; i<n; i++){
+  for (int i = 0; i < n; ++i){
     if (w[i]>= 0.) w[i] = sqrt(w[i]);
     else w[i] = -sqrt(-w[i]);
 
     w[i] *= eml2f;
   }
 
   memory->destroy(work);
   memory->destroy(rwork);
   memory->destroy(iwork);
 
 return info;
 }
 
 /* ----------------------------------------------------------------------------
  * method to get the Dynamical Matrix at q
  * ---------------------------------------------------------------------------- */
 void DynMat::getDMq(double *q)
 {
   interpolate->execute(q, DM_q[0]);
 return;
 }
 
 /* ----------------------------------------------------------------------------
  * method to get the Dynamical Matrix at q
  * ---------------------------------------------------------------------------- */
 void DynMat::getDMq(double *q, double *wt)
 {
   interpolate->execute(q, DM_q[0]);
 
   if (flag_skip && interpolate->UseGamma ) wt[0] = 0.;
 return;
 }
 
 /* ----------------------------------------------------------------------------
  * private method to convert the cartisan coordinate of basis into fractional
  * ---------------------------------------------------------------------------- */
-void DynMat::car2dir(int flag)
+void DynMat::car2dir()
 {
-  if (!flag){ // in newer version, this is done in fix-phonon
-    for (int i=0; i<3; i++){
-      basevec[i]   /= double(nx);
-      basevec[i+3] /= double(ny);
-      basevec[i+6] /= double(nz);
-    }
-  }
   double mat[9];
-  for (int idim=0; idim<9; idim++) mat[idim] = basevec[idim];
+  for (int idim = 0; idim < 9; ++idim) mat[idim] = basevec[idim];
   GaussJordan(3, mat);
 
-  for (int i=0; i<nucell; i++){
+  for (int i = 0; i < nucell; ++i){
     double x[3];
     x[0] = x[1] = x[2] = 0.;
-    for (int idim=0; idim<sysdim; idim++) x[idim] = basis[i][idim];
-    for (int idim=0; idim<sysdim; idim++)
+    for (int idim = 0; idim < sysdim; idim++) x[idim] = basis[i][idim];
+    for (int idim = 0; idim < sysdim; idim++)
       basis[i][idim] = x[0]*mat[idim] + x[1]*mat[3+idim] + x[2]*mat[6+idim];
   }
 
 return;
 }
 
 /* ----------------------------------------------------------------------------
  * private method to enforce the acoustic sum rule on force constant matrix at G
  * ---------------------------------------------------------------------------- */
 void DynMat::EnforceASR()
 {
   char str[MAXLINE];
   int nasr = 20;
   if (nucell <= 1) nasr = 1;
-  printf("\n"); for (int i=0; i<80; i++) printf("=");
+  printf("\n"); for (int i = 0; i < 80; ++i) printf("=");
 
   // compute and display eigenvalues of Phi at gamma before ASR
   if (nucell > 100){
     printf("\nYour unit cell is rather large, eigenvalue evaluation takes some time...");
     fflush(stdout);
   }
 
   double egvs[fftdim];
-  for (int i=0; i<fftdim; i++)
-  for (int j=0; j<fftdim; j++) DM_q[i][j] = DM_all[0][i*fftdim+j];
+  for (int i = 0; i < fftdim; ++i)
+  for (int j = 0; j < fftdim; ++j) DM_q[i][j] = DM_all[0][i*fftdim+j];
   geteigen(egvs, 0);
   printf("\nEigenvalues of Phi at gamma before enforcing ASR:\n");
-  for (int i=0; i<fftdim; i++){
+  for (int i = 0; i < fftdim; ++i){
     printf("%lg ", egvs[i]);
     if (i%10 == 9) printf("\n");
     if (i == 99){ printf("...... (%d more skipped)\n", fftdim-100); break;}
   }
   printf("\n\n");
 
   // ask for iterations to enforce ASR
   printf("Please input the # of iterations to enforce ASR [%d]: ", nasr);
   fgets(str,MAXLINE,stdin);
   char *ptr = strtok(str," \t\n\r\f");
   if (ptr) nasr = atoi(ptr);
-  if (nasr < 1){return; for (int i=0; i<80; i++) printf("="); printf("\n");}
+  if (nasr < 1){
+    for (int i=0; i<80; i++) printf("="); printf("\n");
+    return;
+  }
 
-  for (int iit=0; iit<nasr; iit++){
+  for (int iit = 0; iit < nasr; ++iit){
     // simple ASR; the resultant matrix might not be symmetric
-    for (int a=0; a<sysdim; a++)
-    for (int b=0; b<sysdim; b++){
-      for (int k=0; k<nucell; k++){
+    for (int a = 0; a < sysdim; ++a)
+    for (int b = 0; b < sysdim; ++b){
+      for (int k = 0; k < nucell; ++k){
         double sum = 0.;
-        for (int kp=0; kp<nucell; kp++){
+        for (int kp = 0; kp < nucell; ++kp){
           int idx = (k*sysdim+a)*fftdim+kp*sysdim+b;
           sum += DM_all[0][idx].r;
         }
         sum /= double(nucell);
-        for (int kp=0; kp<nucell; kp++){
+        for (int kp = 0; kp < nucell; ++kp){
           int idx = (k*sysdim+a)*fftdim+kp*sysdim+b;
           DM_all[0][idx].r -= sum;
         }
       }
     }
    
     // symmetrize
-    for (int k=0; k<nucell; k++)
-    for (int kp=k; kp<nucell; kp++){
+    for (int k  = 0; k  < nucell; ++k)
+    for (int kp = k; kp < nucell; ++kp){
       double csum = 0.;
-      for (int a=0; a<sysdim; a++)
-      for (int b=0; b<sysdim; b++){
+      for (int a = 0; a < sysdim; ++a)
+      for (int b = 0; b < sysdim; ++b){
         int idx = (k*sysdim+a)*fftdim+kp*sysdim+b;
         int jdx = (kp*sysdim+b)*fftdim+k*sysdim+a;
         csum = (DM_all[0][idx].r + DM_all[0][jdx].r )*0.5;
         DM_all[0][idx].r = DM_all[0][jdx].r = csum;
       }
     }
   }
 
   // symmetric ASR
-  for (int a=0; a<sysdim; a++)
-  for (int b=0; b<sysdim; b++){
-    for (int k=0; k<nucell; k++){
+  for (int a = 0; a < sysdim; ++a)
+  for (int b = 0; b < sysdim; ++b){
+    for (int k = 0; k < nucell; ++k){
       double sum = 0.;
-      for (int kp=0; kp<nucell; kp++){
+      for (int kp = 0; kp < nucell; ++kp){
         int idx = (k*sysdim+a)*fftdim+kp*sysdim+b;
         sum += DM_all[0][idx].r;
       }
       sum /= double(nucell-k);
-      for (int kp=k; kp<nucell; kp++){
+      for (int kp = k; kp < nucell; ++kp){
         int idx = (k*sysdim+a)*fftdim+kp*sysdim+b;
         int jdx = (kp*sysdim+b)*fftdim+k*sysdim+a;
         DM_all[0][idx].r -= sum;
         DM_all[0][jdx].r  = DM_all[0][idx].r;
       }
     }
   }
 
   // compute and display eigenvalues of Phi at gamma after ASR
-  for (int i=0; i<fftdim; i++)
-  for (int j=0; j<fftdim; j++) DM_q[i][j] = DM_all[0][i*fftdim+j];
+  for (int i = 0; i < fftdim; ++i)
+  for (int j = 0; j < fftdim; ++j) DM_q[i][j] = DM_all[0][i*fftdim+j];
   geteigen(egvs, 0);
   printf("Eigenvalues of Phi at gamma after enforcing ASR:\n");
-  for (int i=0; i<fftdim; i++){
+  for (int i = 0; i < fftdim; ++i){
     printf("%lg ", egvs[i]);
     if (i%10 == 9) printf("\n");
     if (i == 99){ printf("...... (%d more skiped)", fftdim-100); break;}
   }
-  printf("\n"); for (int i=0; i<80; i++) printf("="); printf("\n\n");
+  printf("\n");
+  for (int i = 0; i < 80; ++i) printf("="); printf("\n\n");
 
 return;
 }
 
 /* ----------------------------------------------------------------------------
  * private method to get the reciprocal lattice vectors from the real space ones
  * ---------------------------------------------------------------------------- */
 void DynMat::real2rec()
 {
   ibasevec[0] = basevec[4]*basevec[8] - basevec[5]*basevec[7];
   ibasevec[1] = basevec[5]*basevec[6] - basevec[3]*basevec[8];
   ibasevec[2] = basevec[3]*basevec[7] - basevec[4]*basevec[6];
 
   ibasevec[3] = basevec[7]*basevec[2] - basevec[8]*basevec[1];
   ibasevec[4] = basevec[8]*basevec[0] - basevec[6]*basevec[2];
   ibasevec[5] = basevec[6]*basevec[1] - basevec[7]*basevec[0];
 
   ibasevec[6] = basevec[1]*basevec[5] - basevec[2]*basevec[4];
   ibasevec[7] = basevec[2]*basevec[3] - basevec[0]*basevec[5];
   ibasevec[8] = basevec[0]*basevec[4] - basevec[1]*basevec[3];
 
   double vol = 0.;
-  for (int i=0; i<sysdim; i++) vol += ibasevec[i] * basevec[i];
+  for (int i = 0; i < sysdim; ++i) vol += ibasevec[i] * basevec[i];
   vol = 8.*atan(1.)/vol;
 
-  for (int i=0; i<9; i++) ibasevec[i] *= vol;
+  for (int i = 0; i < 9; ++i) ibasevec[i] *= vol;
 
-  printf("\n"); for (int i=0; i<80; i++) printf("=");
+  printf("\n"); for (int i = 0; i < 80; ++i) printf("=");
   printf("\nBasis vectors of the unit cell in real space:");
-  for (int i=0; i<sysdim; i++){
+  for (int i = 0; i < sysdim; ++i){
     printf("\n     A%d: ", i+1);
-    for (int j=0; j<sysdim; j++) printf("%8.4f ", basevec[i*3+j]);
+    for (int j = 0; j < sysdim; ++j) printf("%8.4f ", basevec[i*3+j]);
   }
   printf("\nBasis vectors of the corresponding reciprocal cell:");
-  for (int i=0; i<sysdim; i++){
+  for (int i = 0; i < sysdim; ++i){
     printf("\n     B%d: ", i+1);
-    for (int j=0; j<sysdim; j++) printf("%8.4f ", ibasevec[i*3+j]);
+    for (int j = 0; j < sysdim; ++j) printf("%8.4f ", ibasevec[i*3+j]);
   }
-  printf("\n"); for (int i=0; i<80; i++) printf("="); printf("\n");
+  printf("\n"); for (int i = 0; i < 80; ++i) printf("="); printf("\n");
 
 return;
 }
 
 /* ----------------------------------------------------------------------
  * private method, to get the inverse of a double matrix by means of
  * Gaussian-Jordan Elimination with full pivoting; square matrix required.
  *
  * Adapted from the Numerical Recipes in Fortran.
  * --------------------------------------------------------------------*/
 void DynMat::GaussJordan(int n, double *Mat)
 {
   int i,icol,irow,j,k,l,ll,idr,idc;
   int *indxc,*indxr,*ipiv;
   double big, nmjk;
   double dum, pivinv;
 
   indxc = new int[n];
   indxr = new int[n];
   ipiv  = new int[n];
 
-  for (i=0; i<n; i++) ipiv[i] = 0;
-  for (i=0; i<n; i++){
+  for (i = 0; i < n; ++i) ipiv[i] = 0;
+  for (i = 0; i < n; ++i){
     big = 0.;
-    for (j=0; j<n; j++){
+    for (j = 0; j < n; ++j){
       if (ipiv[j] != 1){
-        for (k=0; k<n; k++){
+        for (k = 0; k < n; ++k){
           if (ipiv[k] == 0){
-            idr = j*n+k;
+            idr = j * n + k;
             nmjk = abs(Mat[idr]);
             if (nmjk >= big){
               big  = nmjk;
               irow = j;
               icol = k;
             }
-          }else if (ipiv[k]>1){
+          } else if (ipiv[k] > 1){
             printf("DynMat: Singular matrix in double GaussJordan!\n"); exit(1);
           }
         }
       }
     }
     ipiv[icol] += 1;
     if (irow != icol){
-      for (l=0; l<n; l++){
+      for (l = 0; l < n; ++l){
         idr  = irow*n+l;
         idc  = icol*n+l;
         dum  = Mat[idr];
         Mat[idr] = Mat[idc];
         Mat[idc] = dum;
       }
     }
     indxr[i] = irow;
     indxc[i] = icol;
-    idr = icol*n+icol;
+    idr = icol * n + icol;
     if (Mat[idr] == 0.){
       printf("DynMat: Singular matrix in double GaussJordan!");
       exit(1);
     }
     
     pivinv = 1./ Mat[idr];
     Mat[idr] = 1.;
     idr = icol*n;
-    for (l=0; l<n; l++) Mat[idr+l] *= pivinv;
-    for (ll=0; ll<n; ll++){
+    for (l = 0; l < n; ++l) Mat[idr+l] *= pivinv;
+    for (ll = 0; ll < n; ++ll){
       if (ll != icol){
-        idc = ll*n+icol;
+        idc = ll * n + icol;
         dum = Mat[idc];
         Mat[idc] = 0.;
         idc -= icol;
-        for (l=0; l<n; l++) Mat[idc+l] -= Mat[idr+l]*dum;
+        for (l = 0; l < n; ++l) Mat[idc+l] -= Mat[idr+l]*dum;
       }
     }
   }
-  for (l=n-1; l>=0; l--){
+  for (l = n-1; l >= 0; --l){
     int rl = indxr[l];
     int cl = indxc[l];
     if (rl != cl){
-      for (k=0; k<n; k++){
-        idr = k*n+rl;
-        idc = k*n+cl;
+      for (k = 0; k < n; ++k){
+        idr = k * n + rl;
+        idc = k * n + cl;
         dum = Mat[idr];
         Mat[idr] = Mat[idc];
         Mat[idc] = dum;
       }
     }
   }
   delete []indxr;
   delete []indxc;
   delete []ipiv;
+
 return;
 }
 
 /* ----------------------------------------------------------------------------
  * Public method to reset the interpolation method
  * ---------------------------------------------------------------------------- */
 void DynMat::reset_interp_method()
 {
   interpolate->set_method();
 
 return;
 }
 
 /* ----------------------------------------------------------------------------
  * Private method to display help info
  * ---------------------------------------------------------------------------- */
 void DynMat::help()
 {
   ShowVersion();
   printf("\nUsage:\n  phana [options] [file]\n\n");
   printf("Available options:\n");
   printf("  -r          To reset the dynamical matrix at the gamma point by a 4th order\n");
   printf("              polynomial interpolation along the [100] direction; this might be\n");
   printf("              useful for systems with charges. As for charged system, the dynamical\n");
   printf("              matrix at Gamma is far from accurate because of the long range nature\n");
   printf("              of Coulombic interaction. By reset it by interpolation, will partially\n");
   printf("              elliminate the unwanted behavior, but the result is still inaccurate.\n");
   printf("              By default, this is not set; and not expected for uncharged systems.\n\n");
   printf("  -s          This will reset the dynamical matrix at the gamma point, too, but it\n");
   printf("              will also inform the code to skip all q-points that is in the vicinity\n");
   printf("              of the gamma point when evaluating phonon DOS and/or phonon dispersion.\n\n");
   printf("              By default, this is not set; and not expected for uncharged systems.\n\n");
   printf("  -h          To print out this help info.\n\n");
   printf("  file        To define the filename that carries the binary dynamical matrice generated\n");
   printf("              by fix-phonon. If not provided, the code will ask for it.\n");
   printf("\n\n");
   exit(0);
 }
 
 /* ----------------------------------------------------------------------------
  * Private method to display the version info
  * ---------------------------------------------------------------------------- */
 void DynMat::ShowVersion()
 {
   printf("                ____  _   _    __    _  _    __   \n");
   printf("               (  _ \\( )_( )  /__\\  ( \\( )  /__\\  \n");
   printf("                )___/ ) _ (  /(__)\\  )  (  /(__)\\ \n");
   printf("               (__)  (_) (_)(__)(__)(_)\\_)(__)(__)\n");
-  printf("\nPHonon ANAlyzer for Fix-Phonon, version 2.%d, compiled on %s.\n", VERSION, __DATE__);
+  printf("\nPHonon ANAlyzer for Fix-Phonon, version 2.%02d, compiled on %s.\n", VERSION, __DATE__);
 
 return;
 }
 /* --------------------------------------------------------------------*/
diff --git a/tools/phonon/dynmat.h b/tools/phonon/dynmat.h
index 95cb2d7a3..1d6e71658 100644
--- a/tools/phonon/dynmat.h
+++ b/tools/phonon/dynmat.h
@@ -1,66 +1,66 @@
 #ifndef DYNMAT_H
 #define DYNMAT_H
 
 #include "stdio.h"
 #include "stdlib.h"
 #include "string.h"
 #include "memory.h"
 #include "interpolate.h"
 
 extern "C"{
 #include "f2c.h"
 #include "clapack.h"
 }
 
 using namespace std;
 
 class DynMat {
 public:
 
   DynMat(int, char**);
   ~DynMat();
 
   int nx, ny, nz, nucell;
   int sysdim, fftdim;
   double eml2f;
   char *funit;
 
   void getDMq(double *);
   void getDMq(double *, double *);
   void writeDMq(double *);
   void writeDMq(double *, const double, FILE *fp);
   int geteigen(double *, int);
   void reset_interp_method();
 
   doublecomplex **DM_q;
 
   int flag_latinfo;
   double Tmeasure, basevec[9], ibasevec[9];
   double **basis;
   int *attyp;
 
 private:
 
   int flag_skip, flag_reset_gamma;
   Interpolate *interpolate;
   
   Memory *memory;
   int npt, fftdim2;
 
   int nasr;
   void EnforceASR();
 
   char *binfile, *dmfile;
   double boltz, q[3];
   double *M_inv_sqrt;
 
   doublecomplex **DM_all;
 
-  void car2dir(int); // to convert basis from cartisian coordinate into factional.
+  void car2dir();      // to convert basis from cartisian coordinate into factional.
   void real2rec();
   void GaussJordan(int, double *);
 
   void help();
   void ShowVersion();
 };
 #endif
diff --git a/tools/phonon/green.cpp b/tools/phonon/green.cpp
index 5136ecc65..8f8946dc4 100644
--- a/tools/phonon/green.cpp
+++ b/tools/phonon/green.cpp
@@ -1,250 +1,249 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
 #include <math.h>
 #include "green.h"
 #include <complex>
-
-#define MAXLINE 256
+#include "global.h"
 
 /*******************************************************************************
  * The class of Green is designed to evaluate the LDOS via the Green's Function
  * method. The meaning of input/output parameters are as follows:
  *
  *   ntm     (input, value)  total number of atoms in system
  *   sdim    (input, value)  dimension of the system; usually 3
  *   niter   (input, value)  maximum iterations during Lanczos diagonalization
  *   min     (input, value)  minimum value for the angular frequency
  *   max     (input, value)  maximum value for the angular frequency
  *   ndos    (input, value)  total number of points in LDOS
  *   eps     (input, value)  epson that govens the width of delta-function
  *   Hessian (input, pointer of pointer) mass-weighted force constant matrix, of
  *                           dimension [natom*sysdim][natm*sysdim]; it is actually
  *                           the dynamical matrix at gamma point
  *   itm     (input, value)  index of the atom to evaluate local phonon DOS, from 0
  *   lpdos   (output, array) double array of size (ndos, sdim)
  *******************************************************************************
  * References:
  *  1. Z. Tang and N. R. Aluru, Phys. Rev. B 74, 235441 (2006).
  *  2. C. Hudon, R. Meyer, and L.J. Lewis, Phys. Rev. B 76, 045409 (2007).
  *  3. L.T. Kong and L.J. Lewis, Phys. Rev. B 77, 165422 (2008).
  *
  * NOTE: The real-space Green's function method is not expected to work accurately
  *       for small systems, say, system (unit cell) less than 500 atoms.
  *******************************************************************************/
 
 /*------------------------------------------------------------------------------
  * Constructor is used as the main driver
  *----------------------------------------------------------------------------*/
 Green::Green(const int ntm, const int sdim, const int niter, const double min, const double max,
              const int ndos, const double eps, double **Hessian, const int itm, double **lpdos)
 {
   const double tpi = 8.*atan(1.);
   natom = ntm; sysdim = sdim; nit = niter; epson = eps;
   wmin = min*tpi; wmax = max*tpi; nw = ndos + (ndos+1)%2;
   H = Hessian; iatom = itm;
   ldos = lpdos;
 
-  memory = new Memory;
+  memory = new Memory();
   if (natom < 1 || iatom < 0 || iatom >= natom){
     printf("\nError: Wrong number of total atoms or wrong index of interested atom!\n");
     return;
   }
   ndim = natom * sysdim;
 
   if (nit < 1){printf("\nError: Wrong input of maximum iterations!\n"); return;}
   if (nit > ndim){printf("\nError: # Lanczos iterations is not expected to exceed the degree of freedom!\n"); return;}
   if (nw  < 1){printf("\nError: Wrong input of points in LDOS!\n"); return;}
 
   // initialize variables and allocate local memories
   dw = (wmax - wmin)/double(nw-1);
-  alpha = memory->create(alpha, sysdim,nit,  "Green_Green:alpha");
-  beta  = memory->create(beta,  sysdim,nit+1,"Green_Green:beta");
-  //ldos  = memory->create(ldos,  nw,sysdim, "Green_Green:ldos");
+  memory->create(alpha, sysdim,nit,  "Green_Green:alpha");
+  memory->create(beta,  sysdim,nit+1,"Green_Green:beta");
+  //memory->create(ldos,  nw,sysdim, "Green_Green:ldos");
 
   // use Lanczos algorithm to diagonalize the Hessian
   Lanczos();
 
   // Get the inverser of the treated hessian by continued fractional method
   Recursion();
 
 return;
 }
 
 /*------------------------------------------------------------------------------
  * Deconstructor is used to free memory
  *----------------------------------------------------------------------------*/
 Green::~Green()
 {
   H = NULL;
   ldos = NULL;
   memory->destroy(alpha);
   memory->destroy(beta);
 
   delete memory;
 
 return;
 }
       
 /*------------------------------------------------------------------------------
  * Private method to diagonalize a matrix by the Lanczos algorithm
  *----------------------------------------------------------------------------*/
 void Green::Lanczos()
 {
   double *vp, *v, *w, *ptr;
 
   vp = new double [ndim];
   v  = new double [ndim];
   w  = new double [ndim];
   
   int ipos = iatom*sysdim;
 
   // Loop over dimension
-  for (int idim=0; idim<sysdim; idim++){
+  for (int idim = 0; idim < sysdim; ++idim){
     beta[idim][0] = 0.;
-    for (int i=0; i<ndim; i++){vp[i] = v[i] = 0.;}
+    for (int i = 0; i < ndim; ++i) vp[i] = v[i] = 0.;
     v[ipos+idim] = 1.;
 
     // Loop on fraction levels
-    for (int i=0; i<nit; i++){
+    for (int i = 0; i < nit; ++i){
       double sum_a = 0.;
-      for (int j=0; j<ndim; j++){
+      for (int j = 0; j < ndim; ++j){
         double sumHv = 0.;
-        for (int k=0; k<ndim; k++) sumHv += H[j][k]*v[k];
+        for (int k = 0; k < ndim; ++k) sumHv += H[j][k]*v[k];
         w[j] = sumHv - beta[idim][i]*vp[j];
         sum_a += w[j]*v[j];
       }
       alpha[idim][i] = sum_a;
 
-      for (int k=0; k<ndim; k++) w[k] -= alpha[idim][i]*v[k];
+      for (int k = 0; k < ndim; ++k) w[k] -= alpha[idim][i]*v[k];
 
       double gamma = 0.;
-      for (int k=0; k<ndim; k++) gamma += w[k]*v[k];
-      for (int k=0; k<ndim; k++) w[k] -= gamma*v[k];
+      for (int k = 0; k < ndim; ++k) gamma += w[k]*v[k];
+      for (int k = 0; k < ndim; ++k) w[k] -= gamma*v[k];
 
       double sum_b = 0.;
-      for (int k=0; k<ndim; k++) sum_b += w[k]*w[k];
+      for (int k = 0; k < ndim; ++k) sum_b += w[k]*w[k];
       beta[idim][i+1] = sqrt(sum_b);
 
       ptr = vp; vp = v; v = ptr;
       double tmp = 1./beta[idim][i+1];    
-      for (int k=0; k<ndim; k++) v[k] = w[k]*tmp;
+      for (int k = 0; k < ndim; ++k) v[k] = w[k] * tmp;
     }
   }
 
   ptr = NULL;
   delete []vp;
   delete []v;
   delete []w;
 
 return;
 }
 
 /*------------------------------------------------------------------------------
  * Private method to compute the LDOS via the recusive method for system with
  * many atoms
  *----------------------------------------------------------------------------*/
 void Green::Recursion()
 {
   // local variables
   double *alpha_inf, *beta_inf, *xmin, *xmax;
   alpha_inf = new double [sysdim];
   beta_inf  = new double [sysdim];
   xmin  = new double [sysdim];
   xmax  = new double [sysdim];
 
   int nave = nit/4;
-  for (int idim=0; idim<sysdim; idim++){
+  for (int idim = 0; idim < sysdim; ++idim){
     alpha_inf[idim] = beta_inf[idim] = 0.;
 
-    for (int i= nit-nave; i<nit; i++){
+    for (int i = nit-nave; i < nit; ++i){
       alpha_inf[idim] += alpha[idim][i];
       beta_inf[idim] += beta[idim][i+1];
     }
 
     alpha_inf[idim] /= double(nave);
     beta_inf[idim]  /= double(nave);
 
     xmin[idim] = alpha_inf[idim] - 2.*beta_inf[idim];
     xmax[idim] = alpha_inf[idim] + 2.*beta_inf[idim];
   }
 
   std::complex<double> Z, z_m_a, r_x, rec_x, rec_x_inv;
   double sr, si;
 
   double w = wmin;
-  for (int i=0; i<nw; i++){
+  for (int i = 0; i < nw; ++i){
     double a = w*w, ax, bx;
     Z = std::complex<double>(w*w, epson);
 
-    for (int idim=0; idim<sysdim; idim++){
+    for (int idim = 0; idim < sysdim; ++idim){
       double two_b = 2.*beta_inf[idim]*beta_inf[idim];
       double rtwob = 1./two_b;
 
       z_m_a = Z - alpha_inf[idim]*alpha_inf[idim];
 
       if ( a < xmin[idim] ){
         r_x = sqrt(-2.*two_b + z_m_a);
         ax = std::real(r_x) * rtwob;
         bx = std::imag(r_x) * rtwob;
       } else if (a > xmax[idim]) {
         r_x = sqrt(-2.*two_b + z_m_a);
         ax = -std::real(r_x) * rtwob;
         bx = -std::imag(r_x) * rtwob;
       } else {
         r_x = sqrt(2.*two_b - z_m_a);
         ax =  std::imag(r_x) * rtwob;
         bx = -std::real(r_x) * rtwob;
       }
 
       sr = (a - alpha_inf[idim])*rtwob + ax;
       si = epson * rtwob + bx;
       rec_x = std::complex<double> (sr, si);
 
-      for (int j=0; j<nit; j++){
+      for (int j = 0; j < nit; ++j){
         rec_x_inv = Z - alpha[idim][nit-j-1] - beta[idim][nit-j]*beta[idim][nit-j]*rec_x;
         rec_x = 1./rec_x_inv;
       }
       ldos[i][idim] = std::imag(rec_x)*w;
     }
     w += dw;
   }
   delete []alpha_inf;
   delete []beta_inf;
   delete []xmin;
   delete []xmax;
 
 return;
 }
 
 /*------------------------------------------------------------------------------
  * Private method to compute the LDOS via the recusive method for system with
  * a few atoms (less than NMAX)
  *----------------------------------------------------------------------------*/
 void Green::recursion()
 {
   // local variables
   std::complex<double> Z, rec_x, rec_x_inv;
   std::complex<double> cunit = std::complex<double>(0.,1.);
 
   double w = wmin;
 
-  for (int i=0; i<nw; i++){
+  for (int i = 0; i < nw; ++i){
     Z = std::complex<double>(w*w, epson);
 
-    for (int idim=0; idim<sysdim; idim++){
+    for (int idim = 0; idim < sysdim; ++idim){
       rec_x = std::complex<double>(0.,0.);
 
-      for (int j=0; j<nit; j++){
+      for (int j = 0; j < nit; ++j){
         rec_x_inv = Z - alpha[idim][nit-j-1] - beta[idim][nit-j]*beta[idim][nit-j]*rec_x;
         rec_x = 1./rec_x_inv;
       }
       ldos[i][idim] = std::imag(rec_x)*w;
     }
     w += dw;
   }
 
 return;
 }
 
 /*------------------------------------------------------------------------------*/
diff --git a/tools/phonon/interpolate.cpp b/tools/phonon/interpolate.cpp
index 35b015e50..8c0cbde1c 100644
--- a/tools/phonon/interpolate.cpp
+++ b/tools/phonon/interpolate.cpp
@@ -1,310 +1,309 @@
 #include "interpolate.h"
 #include "math.h"
-
-#define MAXLINE 256
-#define MIN(a,b) ((a)>(b)?(b):(a))
-#define MAX(a,b) ((a)>(b)?(a):(b))
+#include "global.h"
 
 /* ----------------------------------------------------------------------------
  * Constructor used to get info from caller, and prepare other necessary data
  * ---------------------------------------------------------------------------- */
 Interpolate::Interpolate(int nx, int ny, int nz, int ndm, doublecomplex **DM)
 {
   Nx = nx;
   Ny = ny;
   Nz = nz;
   Npt = Nx*Ny*Nz;
   ndim = ndm;
-  memory = new Memory;
+  memory = new Memory();
 
   which = UseGamma = 0;
 
   data = DM;
   Dfdx = Dfdy = Dfdz = D2fdxdy = D2fdxdz = D2fdydz = D3fdxdydz = NULL;
   flag_reset_gamma = flag_allocated_dfs = 0;
 
 return;
 }
 
 /* ----------------------------------------------------------------------------
  * Private method to initialize tricubic interpolations
  * ---------------------------------------------------------------------------- */
 void Interpolate::tricubic_init()
 {
   // prepare necessary data for tricubic
   if (flag_allocated_dfs == 0){
-    Dfdx = memory->create(Dfdx, Npt, ndim, "Interpolate_Interpolate:Dfdx");
-    Dfdy = memory->create(Dfdy, Npt, ndim, "Interpolate_Interpolate:Dfdy");
-    Dfdz = memory->create(Dfdz, Npt, ndim, "Interpolate_Interpolate:Dfdz");
-    D2fdxdy = memory->create(D2fdxdy, Npt, ndim, "Interpolate_Interpolate:D2fdxdy");
-    D2fdxdz = memory->create(D2fdxdz, Npt, ndim, "Interpolate_Interpolate:D2fdxdz");
-    D2fdydz = memory->create(D2fdydz, Npt, ndim, "Interpolate_Interpolate:D2fdydz");
-    D3fdxdydz = memory->create(D3fdxdydz, Npt, ndim, "Interpolate_Interpolate:D2fdxdydz");
+    memory->create(Dfdx, Npt, ndim, "Interpolate_Interpolate:Dfdx");
+    memory->create(Dfdy, Npt, ndim, "Interpolate_Interpolate:Dfdy");
+    memory->create(Dfdz, Npt, ndim, "Interpolate_Interpolate:Dfdz");
+    memory->create(D2fdxdy, Npt, ndim, "Interpolate_Interpolate:D2fdxdy");
+    memory->create(D2fdxdz, Npt, ndim, "Interpolate_Interpolate:D2fdxdz");
+    memory->create(D2fdydz, Npt, ndim, "Interpolate_Interpolate:D2fdydz");
+    memory->create(D3fdxdydz, Npt, ndim, "Interpolate_Interpolate:D2fdxdydz");
 
     flag_allocated_dfs = 1;
   }
 
   // get the derivatives
   int n=0;
   const double half = 0.5, one4 = 0.25, one8 = 0.125;
-  for (int ii=0; ii<Nx; ii++)
-  for (int jj=0; jj<Ny; jj++)
-  for (int kk=0; kk<Nz; kk++){
+  for (int ii = 0; ii < Nx; ++ii)
+  for (int jj = 0; jj < Ny; ++jj)
+  for (int kk = 0; kk < Nz; ++kk){
 
     int ip = (ii+1)%Nx, jp = (jj+1)%Ny, kp = (kk+1)%Nz;
     int im = (ii-1+Nx)%Nx, jm = (jj-1+Ny)%Ny, km = (kk-1+Nz)%Nz;
 
     int p100 = (ip*Ny+jj)*Nz+kk;
     int p010 = (ii*Ny+jp)*Nz+kk;
     int p001 = (ii*Ny+jj)*Nz+kp;
     int p110 = (ip*Ny+jp)*Nz+kk;
     int p101 = (ip*Ny+jj)*Nz+kp;
     int p011 = (ii*Ny+jp)*Nz+kp;
     int pm00 = (im*Ny+jj)*Nz+kk;
     int p0m0 = (ii*Ny+jm)*Nz+kk;
     int p00m = (ii*Ny+jj)*Nz+km;
     int pmm0 = (im*Ny+jm)*Nz+kk;
     int pm0m = (im*Ny+jj)*Nz+km;
     int p0mm = (ii*Ny+jm)*Nz+km;
     int p1m0 = (ip*Ny+jm)*Nz+kk;
     int p10m = (ip*Ny+jj)*Nz+km;
     int p01m = (ii*Ny+jp)*Nz+km;
     int pm10 = (im*Ny+jp)*Nz+kk;
     int pm01 = (im*Ny+jj)*Nz+kp;
     int p0m1 = (ii*Ny+jm)*Nz+kp;
     int p111 = (ip*Ny+jp)*Nz+kp;
     int pm11 = (im*Ny+jp)*Nz+kp;
     int p1m1 = (ip*Ny+jm)*Nz+kp;
     int p11m = (ip*Ny+jp)*Nz+km;
     int pm1m = (im*Ny+jp)*Nz+km;
     int p1mm = (ip*Ny+jm)*Nz+km;
     int pmm1 = (im*Ny+jm)*Nz+kp;
     int pmmm = (im*Ny+jm)*Nz+km;
 
     for (int idim=0; idim<ndim; idim++){
       Dfdx[n][idim].r = (data[p100][idim].r - data[pm00][idim].r) * half;
       Dfdx[n][idim].i = (data[p100][idim].i - data[pm00][idim].i) * half;
       Dfdy[n][idim].r = (data[p010][idim].r - data[p0m0][idim].r) * half;
       Dfdy[n][idim].i = (data[p010][idim].i - data[p0m0][idim].i) * half;
       Dfdz[n][idim].r = (data[p001][idim].r - data[p00m][idim].r) * half;
       Dfdz[n][idim].i = (data[p001][idim].i - data[p00m][idim].i) * half;
       D2fdxdy[n][idim].r = (data[p110][idim].r - data[p1m0][idim].r - data[pm10][idim].r + data[pmm0][idim].r) * one4;
       D2fdxdy[n][idim].i = (data[p110][idim].i - data[p1m0][idim].i - data[pm10][idim].i + data[pmm0][idim].i) * one4;
       D2fdxdz[n][idim].r = (data[p101][idim].r - data[p10m][idim].r - data[pm01][idim].r + data[pm0m][idim].r) * one4;
       D2fdxdz[n][idim].i = (data[p101][idim].i - data[p10m][idim].i - data[pm01][idim].i + data[pm0m][idim].i) * one4;
       D2fdydz[n][idim].r = (data[p011][idim].r - data[p01m][idim].r - data[p0m1][idim].r + data[p0mm][idim].r) * one4;
       D2fdydz[n][idim].i = (data[p011][idim].i - data[p01m][idim].i - data[p0m1][idim].i + data[p0mm][idim].i) * one4;
       D3fdxdydz[n][idim].r = (data[p111][idim].r-data[pm11][idim].r - data[p1m1][idim].r - data[p11m][idim].r +
                               data[p1mm][idim].r+data[pm1m][idim].r + data[pmm1][idim].r - data[pmmm][idim].r) * one8;
       D3fdxdydz[n][idim].i = (data[p111][idim].i-data[pm11][idim].i - data[p1m1][idim].i - data[p11m][idim].i +
                               data[p1mm][idim].i+data[pm1m][idim].i + data[pmm1][idim].i - data[pmmm][idim].i) * one8;
     }
     n++;
   }
 return;
 }
 
 /* ----------------------------------------------------------------------------
  * Deconstructor used to free memory
  * ---------------------------------------------------------------------------- */
 Interpolate::~Interpolate()
 {
   data = NULL;
   memory->destroy(Dfdx);
   memory->destroy(Dfdy);
   memory->destroy(Dfdz);
   memory->destroy(D2fdxdy);
   memory->destroy(D2fdxdz);
   memory->destroy(D2fdydz);
   memory->destroy(D3fdxdydz);
   delete memory;
 }
 
 /* ----------------------------------------------------------------------------
  * Tricubic interpolation, by calling the tricubic library
  * ---------------------------------------------------------------------------- */
 void Interpolate::tricubic(double *qin, doublecomplex *DMq)
 {
   // qin should be in unit of 2*pi/L
   double q[3];
-  for (int i=0; i<3; i++) q[i] = qin[i];
-  for (int i=0; i<3; i++){
+  for (int i = 0; i < 3; ++i) q[i] = qin[i];
+  for (int i = 0; i < 3; ++i){
     while (q[i] < 0.)  q[i] += 1.;
     while (q[i] >= 1.) q[i] -= 1.;
   }
   
   int ix = int(q[0]*double(Nx));
   int iy = int(q[1]*double(Ny));
   int iz = int(q[2]*double(Nz));
   double x = q[0]*double(Nx)-double(ix);
   double y = q[1]*double(Ny)-double(iy);
   double z = q[2]*double(Nz)-double(iz);
   int ixp = (ix+1)%Nx, iyp = (iy+1)%Ny, izp = (iz+1)%Nz;
   vidx[0] = (ix*Ny+iy)*Nz+iz;
   vidx[1] = (ixp*Ny+iy)*Nz+iz;
   vidx[2] = (ix*Ny+iyp)*Nz+iz;
   vidx[3] = (ixp*Ny+iyp)*Nz+iz;
   vidx[4] = (ix*Ny+iy)*Nz+izp;
   vidx[5] = (ixp*Ny+iy)*Nz+izp;
   vidx[6] = (ix*Ny+iyp)*Nz+izp;
   vidx[7] = (ixp*Ny+iyp)*Nz+izp;
   for (int i=0; i<8; i++) if (vidx[i] == 0) UseGamma = 1;
 
-  for (int idim=0; idim<ndim; idim++){
-    for (int i=0; i<8; i++){
+  for (int idim = 0; idim < ndim; ++idim){
+    for (int i = 0; i < 8; ++i){
       f[i] = data[vidx[i]][idim].r;
       dfdx[i] = Dfdx[vidx[i]][idim].r;
       dfdy[i] = Dfdy[vidx[i]][idim].r;
       dfdz[i] = Dfdz[vidx[i]][idim].r;
       d2fdxdy[i] = D2fdxdy[vidx[i]][idim].r;
       d2fdxdz[i] = D2fdxdz[vidx[i]][idim].r;
       d2fdydz[i] = D2fdydz[vidx[i]][idim].r;
       d3fdxdydz[i] = D3fdxdydz[vidx[i]][idim].r;
     }
     tricubic_get_coeff(&a[0],&f[0],&dfdx[0],&dfdy[0],&dfdz[0],&d2fdxdy[0],&d2fdxdz[0],&d2fdydz[0],&d3fdxdydz[0]); 
     DMq[idim].r = tricubic_eval(&a[0],x,y,z);
     
-    for (int i=0; i<8; i++){
+    for (int i = 0; i < 8; ++i){
       f[i] = data[vidx[i]][idim].i;
       dfdx[i] = Dfdx[vidx[i]][idim].i;
       dfdy[i] = Dfdy[vidx[i]][idim].i;
       dfdz[i] = Dfdz[vidx[i]][idim].i;
       d2fdxdy[i] = D2fdxdy[vidx[i]][idim].i;
       d2fdxdz[i] = D2fdxdz[vidx[i]][idim].i;
       d2fdydz[i] = D2fdydz[vidx[i]][idim].i;
       d3fdxdydz[i] = D3fdxdydz[vidx[i]][idim].i;
     }
     tricubic_get_coeff(&a[0],&f[0],&dfdx[0],&dfdy[0],&dfdz[0],&d2fdxdy[0],&d2fdxdz[0],&d2fdydz[0],&d3fdxdydz[0]); 
     DMq[idim].i = tricubic_eval(&a[0],x,y,z);
   }
+
+return;
 }
 
 /* ----------------------------------------------------------------------------
  * method to interpolate the DM at an arbitrary q point;
  * the input q should be a vector in unit of (2pi/a 2pi/b 2pi/c).
  * All q components will be rescaled into [0 1).
  * ---------------------------------------------------------------------------- */
 void Interpolate::trilinear(double *qin, doublecomplex *DMq)
 {
   // rescale q[i] into [0 1)
   double q[3];
-  for (int i=0; i<3; i++) q[i] = qin[i];
-  for (int i=0; i<3; i++){
+  for (int i = 0; i < 3; ++i) q[i] = qin[i];
+  for (int i = 0; i < 3; ++i){
     while (q[i] < 0.)  q[i] += 1.;
     while (q[i] >= 1.) q[i] -= 1.;
   }
 
   // find the index of the eight vertice
   int ix, iy, iz, ixp, iyp, izp;
   double x, y, z;
   q[0] *= double(Nx);
   q[1] *= double(Ny);
   q[2] *= double(Nz);
 
   ix = int(q[0])%Nx;
   iy = int(q[1])%Ny;
   iz = int(q[2])%Nz;
   ixp = (ix+1)%Nx;
   iyp = (iy+1)%Ny;
   izp = (iz+1)%Nz;
   x = q[0] - double(ix);
   y = q[1] - double(iy);
   z = q[2] - double(iz);
 
 //--------------------------------------
   vidx[0] = ((ix*Ny)+iy)*Nz + iz;
   vidx[1] = ((ixp*Ny)+iy)*Nz + iz;
   vidx[2] = ((ix*Ny)+iyp)*Nz + iz;
   vidx[3] = ((ix*Ny)+iy)*Nz + izp;
   vidx[4] = ((ixp*Ny)+iy)*Nz + izp;
   vidx[5] = ((ix*Ny)+iyp)*Nz + izp;
   vidx[6] = ((ixp*Ny)+iyp)*Nz + iz;
   vidx[7] = ((ixp*Ny)+iyp)*Nz + izp;
-  for (int i=0; i<8; i++) if (vidx[i] == 0) UseGamma = 1;
+  for (int i = 0; i < 8; ++i) if (vidx[i] == 0) UseGamma = 1;
 
   double fac[8];
   fac[0] = (1.-x)*(1.-y)*(1.-z);
   fac[1] = x*(1.-y)*(1.-z);
   fac[2] = (1.-x)*y*(1.-z);
   fac[3] = (1.-x)*(1.-y)*z;
   fac[4] = x*(1.-y)*z;
   fac[5] = (1.-x)*y*z;
   fac[6] = x*y*(1.-z);
   fac[7] = x*y*z;
   
   // now to do the interpolation
-  for (int idim=0; idim<ndim; idim++){
+  for (int idim = 0; idim < ndim; ++idim){
     DMq[idim].r = 0.;
     DMq[idim].i = 0.;
-    for (int i=0; i<8; i++){
+    for (int i = 0; i < 8; ++i){
       DMq[idim].r += data[vidx[i]][idim].r*fac[i];
       DMq[idim].i += data[vidx[i]][idim].i*fac[i];
     }
   }
 
 return;
 }
 
 /* ----------------------------------------------------------------------------
  * To invoke the interpolation
  * ---------------------------------------------------------------------------- */
 void Interpolate::execute(double *qin, doublecomplex *DMq)
 {
   UseGamma = 0;
   if (which == 1) // 1: tricubic
     tricubic(qin, DMq);
   else       // otherwise: trilinear
     trilinear(qin, DMq);
 return;
 }
 
 /* ----------------------------------------------------------------------------
  * Public method, to set/reset the interpolation method
  * ---------------------------------------------------------------------------- */
 void Interpolate::set_method()
 {
   char str[MAXLINE];
   int im = 1;
   printf("\n");for(int i=0; i<80; i++) printf("=");
   printf("\nWhich interpolation method would you like to use?\n");
   printf("  1. Tricubic;\n  2. Trilinear;\n");
-  printf("Your choice[1]: ");
+  printf("Your choice [1]: ");
   fgets(str,MAXLINE,stdin);
   char *ptr = strtok(str," \t\n\r\f");
   if (ptr) im = atoi(ptr);
 
   which =2-im%2;
-  printf("Your selection: %d\n", which);
+  printf("Your  selection: %d\n", which);
   for(int i=0; i<80; i++) printf("="); printf("\n\n");
 
   if (which == 1) tricubic_init();
 
 return;
 }
 
 /* ----------------------------------------------------------------------------
  * Public method, to reset gamma point data; in this case, the gamma point data
  * will be meaningless. should only be called once.
  * ---------------------------------------------------------------------------- */
 void Interpolate::reset_gamma()
 {
   if (flag_reset_gamma) return;
   flag_reset_gamma = 1;
 
   int p1 = 1%Nx, p2 = 2%Nx;
   int m1 = (Nx-1), m2 = (Nx-2+Nx)%Nx;
 
   int ip1 = p1*Ny*Nz, ip2 = p2*Ny*Nz;
   int im1 = m1*Ny*Nz, im2 = m2*Ny*Nz;
 
   double const one6 = -1./6., two3 = 2./3.;
 
   for (int idim=0; idim<ndim; idim++){
     data[0][idim].i = 0.;
     data[0][idim].r = (data[im2][idim].r + data[ip2][idim].r) * one6
                     + (data[im1][idim].r + data[ip1][idim].r) * two3;
   }
 
 return;
 }
 /* ---------------------------------------------------------------------------- */
diff --git a/tools/phonon/phonon.cpp b/tools/phonon/phonon.cpp
index 5ee34a6ce..43bea111b 100644
--- a/tools/phonon/phonon.cpp
+++ b/tools/phonon/phonon.cpp
@@ -1,1266 +1,1107 @@
 #include <vector>
 #include "string.h"
 #include "phonon.h"
 #include "green.h"
 #include "timer.h"
+#include "global.h"
 
 #ifdef UseSPG
 extern "C"{
 #include "spglib.h"
 }
 #endif
 
-#define MAXLINE 256
-#define MIN(a,b) ((a)>(b)?(b):(a))
-#define MAX(a,b) ((a)>(b)?(a):(b))
-
 /* ----------------------------------------------------------------------------
  * Class Phonon is the main driver to calculate phonon DOS, phonon
  * dispersion curve and some other things.
  * ---------------------------------------------------------------------------- */
 Phonon::Phonon(DynMat *dm)
 {
   // create memory 
   memory = new Memory();
 
   // pass the class from main
   dynmat = dm;
   sysdim = dynmat->sysdim;
   ndim   = dynmat->fftdim;
   dos  = NULL;
   ldos = NULL;
   qpts = NULL;
   wt   = NULL;
   eigs = NULL;
   locals = NULL;
 
 #ifdef UseSPG
   attyp = NULL;
   atpos = NULL;
 #endif
 
   // display the menu
   char str[MAXLINE];
   while ( 1 ){
-    printf("\n"); for (int i=0; i<37;i++) printf("="); printf(" Menu "); for (int i=0; i<37;i++) printf("="); printf("\n");
+    printf("\n");
+    for (int i = 0; i < 37; ++i) printf("=");
+    printf(" Menu ");
+    for (int i = 0; i < 37; ++i) printf("="); printf("\n");
     printf("  1. Phonon DOS evaluation;\n");
     printf("  2. Phonon dispersion curves;\n");
     printf("  3. Dynamical matrix at arbitrary q;\n");
     printf("  4. Vibrational frequencies at arbitrary q;\n");
     printf("  5. Dispersion-like curve for dynamical matrix;\n");
     printf("  6. Vibrational thermodynamical properties;\n");
     printf("  7. Local phonon DOS from eigenvectors;\n");
     printf("  8. Local phonon DOS by RSGF method;\n");
     printf("  9. Freqs and eigenvectors at arbitrary q;\n");
+    printf(" 10. Show information related to the unit cell;\n");
     printf(" -1. Reset the interpolation method;\n");
     printf("  0. Exit.\n");
     // read user choice
     int job = 0;
-    printf("Your choice[0]: ");
+    printf("Your choice [0]: ");
     if (count_words(fgets(str,MAXLINE,stdin)) > 0) job = atoi(strtok(str," \t\n\r\f"));
-    printf("\nYour selection: %d\n", job);
-    for (int i=0; i<80;i++) printf("=");printf("\n\n");
+    printf("\nYour  selection: %d\n", job);
+    for (int i = 0; i < 80; ++i) printf("=");printf("\n\n");
 
     // now to do the job according to user's choice
     if      (job == 1) pdos();
     else if (job == 2) pdisp();
     else if (job == 3) dmanyq(); 
     else if (job == 4) vfanyq(); 
     else if (job == 5) DMdisp(); 
     else if (job == 6) therm(); 
     else if (job == 7) ldos_egv(); 
     else if (job == 8) ldos_rsgf(); 
     else if (job == 9) vecanyq(); 
+    else if (job ==10) ShowCell();
     else if (job ==-1) dynmat->reset_interp_method();
     else break;
   }
 return;
 }
 
 /* ----------------------------------------------------------------------------
  * Deconstructor to free memory
  * ---------------------------------------------------------------------------- */
 Phonon::~Phonon()
 {
   dynmat = NULL;
 
   memory->destroy(wt);
   memory->destroy(qpts);
   memory->destroy(eigs);
 
   memory->destroy(locals);
 
   memory->destroy(dos);
   memory->destroy(ldos);
 
 #ifdef UseSPG
   memory->destroy(attyp);
   memory->destroy(atpos);
 #endif
   delete memory;
 }
 
 /* ----------------------------------------------------------------------------
  * Private method to calculate the phonon DOS
  * ---------------------------------------------------------------------------- */
 void Phonon::pdos()
 {
   // get frequencies on a q-mesh
   QMesh();       // generate q-points, hopefully irreducible
   ComputeAll();  // get all eigen values ==> frequencies
 
   // now to get the frequency range
   char str[MAXLINE];
 
   fmin = fmax = eigs[0][0];
-  for (int iq=0; iq<nq; iq++){
-    for (int j=0; j<ndim; j++){
-      fmin = MIN(fmin, eigs[iq][j]);
-      fmax = MAX(fmax, eigs[iq][j]);
-    }
+  for (int iq = 0; iq < nq; ++iq)
+  for (int j  = 0; j  < ndim; ++j){
+    fmin = MIN(fmin, eigs[iq][j]);
+    fmax = MAX(fmax, eigs[iq][j]);
   }
 
   // Now to ask for the output frequency range
   printf("\nThe frequency range of all q-points are: [%g %g]\n", fmin, fmax);
   printf("Please input the desired range to get DOS [%g %g]: ", fmin, fmax);
   if (count_words(fgets(str,MAXLINE,stdin)) >= 2){
     fmin = atof(strtok(str," \t\n\r\f"));
     fmax = atof(strtok(NULL," \t\n\r\f"));
   }
   if (fmin > fmax){double swap = fmin; fmin = fmax; fmax = swap;}
   printf("The fequency range for your phonon DOS is [%g %g].\n", fmin, fmax);
 
   ndos = 201;
   printf("Please input the number of intervals [%d]: ", ndos);
   if (count_words(fgets(str,MAXLINE,stdin)) > 0) ndos = atoi(strtok(str," \t\n\r\f"));
   ndos += (ndos+1)%2; ndos = MAX(2,ndos);
 
   
   df  = (fmax-fmin)/double(ndos-1);
   rdf = 1./df;
   memory->destroy(dos);
-  dos = memory->create(dos, ndos, "pdos:dos");
+  memory->create(dos, ndos, "pdos:dos");
   for (int i=0; i<ndos; i++) dos[i] = 0.;
 
   // now to calculate the DOS
   double offset = fmin-0.5*df;
-  for (int iq=0; iq<nq; iq++){
+  for (int iq = 0; iq < nq; ++iq){
     if (wt[iq] > 0.){
-      for (int j=0; j<ndim; j++){
+      for (int j = 0; j < ndim; ++j){
         int idx = int((eigs[iq][j]-offset)*rdf);
         if (idx>=0 && idx<ndos) dos[idx] += wt[iq];
       }
     }
   }
 
   // smooth dos ?
   printf("Would you like to smooth the phonon dos? (y/n)[n]: ");
   if (count_words(fgets(str,MAXLINE,stdin)) > 0){
     char *flag = strtok(str," \t\n\r\f");
     if (strcmp(flag,"y") == 0 || strcmp(flag,"Y") == 0){
       smooth(dos, ndos);
     }
   }
 
   // normalize dos to 1
   Normalize();
 
   // output DOS
   writeDOS();
 
 return;
 }
 
 /* ----------------------------------------------------------------------------
  * Private method to write the phonon DOS to file
  * ---------------------------------------------------------------------------- */
 void Phonon::writeDOS()
 {
   if (dos == NULL) return;
 
   char str[MAXLINE];
   // now to output the phonon DOS
   printf("\nPlease input the filename to write DOS [pdos.dat]: ");
   if (count_words(fgets(str,MAXLINE,stdin)) < 1) strcpy(str, "pdos.dat");
   char *fname = strtok(str," \t\n\r\f");
 
   printf("The total phonon DOS will be written to file: %s\n", fname);
 
   FILE *fp = fopen(fname, "w");
   fprintf(fp,"# frequency  DOS\n");
   fprintf(fp,"#%s  number\n", dynmat->funit);
   double freq = fmin;
   for (int i=0; i<ndos; i++){
     fprintf(fp,"%lg %lg\n", freq, dos[i]);
     freq += df;
   }
   fclose(fp);
 
   // also write the gnuplot script to generate the figure
   fp = fopen("pdos.gnuplot", "w");
   fprintf(fp,"set term post enha colo 20\nset out %cpdos.eps%c\n\n",char(34),char(34));
   fprintf(fp,"set xlabel %cfrequency (THz)%c\n",char(34),char(34));
   fprintf(fp,"set ylabel %cPhonon DOS%c\n",char(34),char(34));
   fprintf(fp,"unset key\n");
   fprintf(fp,"plot %c%s%c u 1:2 w l\n",char(34),fname,char(34));
   fclose(fp);
 
   fname = NULL;
   
 return;
 }
 
 /* ----------------------------------------------------------------------------
  * Private method to write the local DOS to files.
  * ---------------------------------------------------------------------------- */
 void Phonon::writeLDOS()
 {
   if (ldos == NULL) return;
 
   printf("The phonon LDOSs will be written to file(s) : pldos_?.dat\n\n");
   const double one3 = 1./double(sysdim);
   char str[MAXLINE];
-  for (int ilocal=0; ilocal<nlocal; ilocal++){
+  for (int ilocal = 0; ilocal < nlocal; ++ilocal){
     sprintf(str,"pldos_%d.dat", locals[ilocal]);
     char *fname = strtok(str," \t\n\r\f");
 
     FILE *fp = fopen(fname, "w"); fname = NULL;
     fprintf(fp,"#freq xDOS yDOS zDOS total\n");
     double freq = fmin;
-    for (int i=0; i<ndos; i++){
+    for (int i = 0; i < ndos; ++i){
       fprintf(fp,"%lg", freq);
       double total = 0.;
-      for (int idim=0; idim<sysdim; idim++) {fprintf(fp," %lg",ldos[ilocal][i][idim]); total += ldos[ilocal][i][idim];}
+      for (int idim = 0; idim < sysdim; ++idim){
+        fprintf(fp," %lg",ldos[ilocal][i][idim]);
+        total += ldos[ilocal][i][idim];
+      }
       fprintf(fp," %lg\n", total*one3);
       freq += df;
     }
     fclose(fp);
   }
 
 return;
 }
 
 /* ----------------------------------------------------------------------------
  * Private method to calculate the local phonon DOS via the real space Green's
  * function method
  * ---------------------------------------------------------------------------- */
 void Phonon::ldos_rsgf()
 {
   char str[MAXLINE];
   const double tpi = 8.*atan(1.);
   double **Hessian, scale;
   scale = dynmat->eml2f*tpi; scale *= scale;
-  Hessian = memory->create(Hessian, ndim, ndim, "phonon_ldos:Hessian");
+  memory->create(Hessian, ndim, ndim, "phonon_ldos:Hessian");
 
   double q0[3];
   q0[0] = q0[1] = q0[2] = 0.;
 
   dynmat->getDMq(q0);
-  for (int i=0; i<ndim; i++)
-  for (int j=0; j<ndim; j++) Hessian[i][j] = dynmat->DM_q[i][j].r*scale;
+  for (int i = 0; i < ndim; ++i)
+  for (int j = 0; j < ndim; ++j) Hessian[i][j] = dynmat->DM_q[i][j].r*scale;
 
   if (ndim < 300){
     double *egvs = new double [ndim];
     dynmat->geteigen(egvs, 0);
 
     fmin = fmax = egvs[0];
-    for (int i=1; i<ndim; i++){fmin = MIN(fmin, egvs[i]); fmax = MAX(fmax, egvs[i]);}
+    for (int i = 1; i < ndim; ++i){fmin = MIN(fmin, egvs[i]); fmax = MAX(fmax, egvs[i]);}
     delete []egvs;
   } else {
     fmin = 0.; fmax = 20.;
   }
 
   ndos = 201;
   int ik = 0, nit = MAX(ndim*0.1, MIN(ndim,50));
   double eps = 12.; // for Cu with 1000+ atoms, 12 is enough; for small system, eps should be large.
 
   while (1) {
     int istr, iend, iinc;
     // ask for relevant info
     printf("\nThere are %d atoms in each unit cell of your lattice.\n", dynmat->nucell);
     printf("Please input the index/index range/index range and increment of atom(s)\n");
     printf("in the unit cell to evaluate LDOS, q to exit [%d]: ", ik);
     int nr = count_words( fgets(str,MAXLINE,stdin) );
     if (nr < 1){
       istr = iend = ik;
       iinc = 1;
+
     } else if (nr == 1) {
       char *ptr = strtok(str," \t\n\r\f");
       if (strcmp(ptr,"q") == 0) break;
 
       ik = atoi(ptr);
-      if (ik < 0 || ik >= dynmat->nucell) break;
+      ik = MAX(0, MIN(ik, dynmat->nucell-1));
       istr = iend = ik;
       iinc = 1;
+
     } else if (nr == 2) {
       istr = atoi(strtok(str," \t\n\r\f"));
       iend = atoi(strtok(NULL," \t\n\r\f"));
       iinc = 1;
-      if (istr < 0||iend >= dynmat->nucell||istr > iend) break;
+      istr = MAX(0, MIN(istr, dynmat->nucell-1));
+      iend = MAX(0, MIN(iend, dynmat->nucell-1));
+
     } else if (nr >= 3) {
       istr = atoi(strtok(str," \t\n\r\f"));
       iend = atoi(strtok(NULL," \t\n\r\f"));
       iinc = atoi(strtok(NULL," \t\n\r\f"));
-      if (istr<0 || iend >= dynmat->nucell || istr > iend || iinc<1) break;
+      istr = MAX(0, MIN(istr, dynmat->nucell-1));
+      iend = MAX(0, MIN(iend, dynmat->nucell-1));
     }
 
     printf("Please input the frequency range to evaluate LDOS [%g %g]: ", fmin, fmax);
     if (count_words(fgets(str,MAXLINE,stdin)) >= 2){
       fmin = atof(strtok(str," \t\n\r\f"));
       fmax = atof(strtok(NULL," \t\n\r\f"));
     }
     if (fmax < fmin) break;
     printf("The frequency range for your LDOS is [%g %g].\n", fmin, fmax);
 
     printf("Please input the desired number of points in LDOS [%d]: ", ndos);
     if (count_words(fgets(str,MAXLINE,stdin)) > 0) ndos = atoi(strtok(str," \t\n\r\f"));
     if (ndos < 2) break;
     ndos += (ndos+1)%2;
 
     printf("Please input the maximum # of Lanczos iterations  [%d]: ", nit);
     if (count_words(fgets(str,MAXLINE,stdin)) > 0) nit = atoi(strtok(str," \t\n\r\f"));
     if (nit < 1) break;
 
     printf("Please input the value of epsilon for delta-function [%g]: ", eps);
     if (count_words(fgets(str,MAXLINE,stdin)) > 0) eps = atof(strtok(str," \t\n\r\f"));
     if (eps <= 0.) break;
     
     // prepare array for local pdos
     nlocal = 0;
     for (ik = istr; ik <= iend; ik += iinc) nlocal++;
     memory->destroy(ldos);
     ldos = memory->create(ldos,nlocal,ndos,dynmat->sysdim,"ldos_rsgf:ldos");
 
     memory->destroy(locals);
-    locals = memory->create(locals, nlocal, "ldos_rsgf:locals");
+    memory->create(locals, nlocal, "ldos_rsgf:locals");
 
     df  = (fmax-fmin)/double(ndos-1);
     rdf = 1./df;
 
     // to measure the LDOS via real space Green's function method
     int ilocal = 0;
     for (ik = istr; ik <= iend; ik += iinc){
       locals[ilocal] = ik;
 
       // time info
       Timer *time = new Timer();
       printf("\nNow to compute the LDOS for atom %d by Real Space Greens function method ...\n", ik);
       fflush(stdout);
   
       // run real space green's function calculation
       Green *green = new Green(dynmat->nucell, dynmat->sysdim, nit, fmin, fmax, ndos, eps, Hessian, ik, ldos[ilocal++]);
       delete green;
   
       time->stop(); time->print(); delete time;
     }
 
     Normalize();
     writeLDOS();
 
     // evaluate the local vibrational thermal properties optionally
     local_therm();
 
   }
   memory->destroy(Hessian);
 
 return;
 }
 
-/*------------------------------------------------------------------------------
- * Private method to evaluate the phonon dispersion curves
- *----------------------------------------------------------------------------*/
-void Phonon::pdisp()
-{
-  // ask the output file name and write the header.
-  char str[MAXLINE];
-  printf("Please input the filename to output the dispersion data [pdisp.dat]:");
-  if (count_words(fgets(str,MAXLINE,stdin)) < 1) strcpy(str, "pdisp.dat");
-  char *ptr = strtok(str," \t\n\r\f");
-  char *fname = new char[strlen(ptr)+1];
-  strcpy(fname,ptr);
-
-  FILE *fp = fopen(fname, "w");
-  fprintf(fp,"# q     qr    freq\n");
-  fprintf(fp,"# 2pi/L  2pi/L %s\n", dynmat->funit);
-
-  // to store the nodes of the dispersion curve
-  std::vector<double> nodes; nodes.clear();
-
-  // now the calculate the dispersion curve
-  double qstr[3], qend[3], q[3], qinc[3], qr=0., dq;
-  int nq = MAX(MAX(dynmat->nx,dynmat->ny),dynmat->nz)/2+1;
-  qend[0] = qend[1] = qend[2] = 0.;
-
-  double *egvs = new double [ndim];
-  while (1){
-    for (int i=0; i<3; i++) qstr[i] = qend[i];
-
-    int quit = 0;
-    printf("\nPlease input the start q-point in unit of B1->B3, q to exit [%g %g %g]: ", qstr[0], qstr[1], qstr[2]);
-    int n = count_words(fgets(str,MAXLINE,stdin));
-    ptr = strtok(str," \t\n\r\f");
-    if ((n == 1) && (strcmp(ptr,"q") == 0)) break;
-    else if (n >= 3){
-      qstr[0] = atof(ptr);
-      qstr[1] = atof(strtok(NULL," \t\n\r\f"));
-      qstr[2] = atof(strtok(NULL," \t\n\r\f"));
-    }
-
-    do printf("Please input the end q-point in unit of B1->B3: ");
-    while (count_words(fgets(str,MAXLINE,stdin)) < 3);
-    qend[0] = atof(strtok(str," \t\n\r\f"));
-    qend[1] = atof(strtok(NULL," \t\n\r\f"));
-    qend[2] = atof(strtok(NULL," \t\n\r\f"));
-
-    printf("Please input the # of points along the line [%d]: ", nq);
-    if (count_words(fgets(str,MAXLINE,stdin)) > 0) nq = atoi(strtok(str," \t\n\r\f"));
-    nq = MAX(nq,2);
-
-    for (int i=0; i<3; i++) qinc[i] = (qend[i]-qstr[i])/double(nq-1);
-    dq = sqrt(qinc[0]*qinc[0]+qinc[1]*qinc[1]+qinc[2]*qinc[2]);
-
-    nodes.push_back(qr);
-    for (int i=0; i<3; i++) q[i] = qstr[i];
-    for (int ii=0; ii<nq; ii++){
-      double wii = 1.;
-      dynmat->getDMq(q, &wii);
-      if (wii > 0.){
-        dynmat->geteigen(egvs, 0);
-        fprintf(fp,"%lg %lg %lg %lg ", q[0], q[1], q[2], qr);
-        for (int i=0; i<ndim; i++) fprintf(fp," %lg", egvs[i]);
-      }
-      fprintf(fp,"\n");
-
-      for (int i=0; i<3; i++) q[i] += qinc[i];
-      qr += dq;
-    }
-    qr -= dq;
-  }
-  if (qr > 0.) nodes.push_back(qr);
-  fclose(fp);
-  delete []egvs;
-
-  // write the gnuplot script which helps to visualize the result
-  int nnd = nodes.size();
-  if (nnd > 1){
-    fp = fopen("pdisp.gnuplot", "w");
-    fprintf(fp,"set term post enha colo 20\nset out %cpdisp.eps%c\n\n",char(34),char(34));
-    fprintf(fp,"set xlabel %cq%c\n",char(34),char(34));
-    fprintf(fp,"set ylabel %cfrequency (THz)%c\n\n",char(34),char(34));
-    fprintf(fp,"set xrange [0:%lg]\nset yrange [0:*]\n\n", nodes[nnd-1]);
-    fprintf(fp,"set grid xtics\n");
-    fprintf(fp,"# {/Symbol G} will give you letter gamma in the label\nset xtics (");
-    for (int i=0; i<nnd-1; i++) fprintf(fp,"%c%c %lg, ", char(34),char(34),nodes[i]);
-    fprintf(fp,"%c%c %lg)\n\n",char(34),char(34),nodes[nnd-1]);
-    fprintf(fp,"unset key\n\n");
-    fprintf(fp,"plot %c%s%c u 4:5 w l lt 1",char(34),fname,char(34));
-    for (int i=1; i<ndim; i++) fprintf(fp,",\\\n%c%c u 4:%d w l lt 1",char(34),char(34),i+5);
-    fclose(fp);
-  }
-
-  delete []fname;
-  nodes.clear();
-
-return;
-}
-
 /* ----------------------------------------------------------------------------
  * Private method to write out the dynamical matrix at selected q
  * ---------------------------------------------------------------------------- */
 void Phonon::dmanyq()
 {
   char str[MAXLINE];
   double q[3];
   do printf("Please input the q-point to output the dynamical matrix:");
   while (count_words(fgets(str,MAXLINE,stdin)) < 3);
   q[0] = atof(strtok(str," \t\n\r\f"));
   q[1] = atof(strtok(NULL," \t\n\r\f"));
   q[2] = atof(strtok(NULL," \t\n\r\f"));
 
   dynmat->getDMq(q);
   dynmat->writeDMq(q);
 
 return;
 }
 
 /* ----------------------------------------------------------------------------
  * Private method to get the vibrational frequencies at selected q
  * ---------------------------------------------------------------------------- */
 void Phonon::vfanyq()
 {
   char str[MAXLINE];
   double q[3], egvs[ndim];
   
-  while (1){
+  while ( 1 ){
     printf("Please input the q-point to compute the frequencies, q to exit: ");
     if (count_words(fgets(str,MAXLINE,stdin)) < 3) break;
 
     q[0] = atof(strtok(str, " \t\n\r\f"));
     q[1] = atof(strtok(NULL," \t\n\r\f"));
     q[2] = atof(strtok(NULL," \t\n\r\f"));
   
     dynmat->getDMq(q);
     dynmat->geteigen(egvs, 0);
     printf("q-point: [%lg %lg %lg], ", q[0], q[1], q[2]);
     printf("vibrational frequencies at this q-point:\n");
-    for (int i=0; i<ndim; i++) printf("%lg ", egvs[i]); printf("\n\n");
+    for (int i = 0; i < ndim; ++i) printf("%lg ", egvs[i]); printf("\n\n");
   }
 
 return;
 }
 
 /* ----------------------------------------------------------------------------
  * Private method to get the vibrational frequencies and eigenvectors at selected q
  * ---------------------------------------------------------------------------- */
 void Phonon::vecanyq()
 {
   char str[MAXLINE];
   double q[3], egvs[ndim];
   doublecomplex **eigvec = dynmat->DM_q;
   printf("Please input the filename to output the result [eigvec.dat]: ");
   if (count_words(fgets(str,MAXLINE,stdin)) < 1) strcpy(str,"eigvec.dat");
   FILE *fp = fopen(strtok(str," \t\n\r\f"), "w");
 
-  while (1){
+  while ( 1 ){
     printf("Please input the q-point to compute the frequencies, q to exit: ");
     if (count_words(fgets(str,MAXLINE,stdin)) < 3) break;
 
     q[0] = atof(strtok(str, " \t\n\r\f"));
     q[1] = atof(strtok(NULL," \t\n\r\f"));
     q[2] = atof(strtok(NULL," \t\n\r\f"));
   
     dynmat->getDMq(q);
     dynmat->geteigen(egvs, 1);
     fprintf(fp,"# q-point: [%lg %lg %lg], sysdim: %d, # of atoms per cell: %d\n",
       q[0],q[1],q[2], sysdim, dynmat->nucell);
-    for (int i=0; i<ndim; i++){
+    for (int i = 0; i < ndim; ++i){
       fprintf(fp,"# frequency %d at [%lg %lg %lg]: %lg\n",i+1,q[0],q[1],q[2],egvs[i]);
       fprintf(fp,"# atom eigenvector : |e|\n");
-      for (int j=0; j<dynmat->nucell; j++){
+      for (int j = 0; j < dynmat->nucell; ++j){
         int ipos = j * sysdim;
         double sum = 0.;
         fprintf(fp,"%d", j+1);
-        for (int idim=0; idim<sysdim; idim++){
+        for (int idim = 0; idim < sysdim; ++idim){
           fprintf(fp,"  %lg %lg", eigvec[i][ipos+idim].r, eigvec[i][ipos+idim].i);
           sum += eigvec[i][ipos+idim].r * eigvec[i][ipos+idim].r + eigvec[i][ipos+idim].i * eigvec[i][ipos+idim].i;
         }
         fprintf(fp,"  : %lg\n", sqrt(sum));
       }
       fprintf(fp,"\n");
     }
     fprintf(fp,"\n");
   }
   fclose(fp);
   eigvec = NULL;
 return;
 }
 
 /* ----------------------------------------------------------------------------
  * Private method to get the dispersion-like data for dynamical matrix
  * ---------------------------------------------------------------------------- */
 void Phonon::DMdisp()
 {
   // ask the output file name and write the header.
   char str[MAXLINE];
 
   printf("Please input the filename to output the DM data [DMDisp.dat]: ");
   if (count_words(fgets(str,MAXLINE,stdin)) < 1) strcpy(str, "DMDisp.dat");
   char *fname = strtok(str," \t\n\r\f");
 
   FILE *fp = fopen(fname, "w"); fname = NULL;
   fprintf(fp,"# q     qr    D\n");
 
   // now the calculate the dispersion-like curve
   double qstr[3], qend[3], q[3], qinc[3], qr=0., dq;
   int nq = MAX(MAX(dynmat->nx,dynmat->ny),dynmat->nz)/2;
   qend[0] = qend[1] = qend[2] = 0.;
 
-  while (1){
-
-    for (int i=0; i<3; i++) qstr[i] = qend[i];
+  while ( 1 ){
+    for (int i = 0; i < 3; ++i) qstr[i] = qend[i];
 
     printf("\nPlease input the start q-point in unit of B1->B3, q to exit [%g %g %g]: ", qstr[0], qstr[1], qstr[2]);
     int n = count_words(fgets(str,MAXLINE,stdin));
     char *ptr = strtok(str," \t\n\r\f");
     if ((n == 1) && (strcmp(ptr,"q") == 0)) break;
     else if (n >= 3){
       qstr[0] = atof(ptr);
       qstr[1] = atof(strtok(NULL," \t\n\r\f"));
       qstr[2] = atof(strtok(NULL," \t\n\r\f"));
     }
 
     do printf("Please input the end q-point in unit of B1->B3: ");
     while (count_words(fgets(str,MAXLINE,stdin)) < 3);
     qend[0] = atof(strtok(str," \t\n\r\f"));
     qend[1] = atof(strtok(NULL," \t\n\r\f"));
     qend[2] = atof(strtok(NULL," \t\n\r\f"));
 
     printf("Please input the # of points along the line [%d]: ", nq);
     if (count_words(fgets(str,MAXLINE,stdin)) > 0) nq = atoi(strtok(str," \t\n\r\f"));
     nq = MAX(nq,2);
 
     for (int i=0; i<3; i++) qinc[i] = (qend[i]-qstr[i])/double(nq-1);
     dq = sqrt(qinc[0]*qinc[0]+qinc[1]*qinc[1]+qinc[2]*qinc[2]);
 
-    for (int i=0; i<3; i++) q[i] = qstr[i];
-    for (int ii=0; ii<nq; ii++){
+    for (int i = 0; i < 3; ++i) q[i] = qstr[i];
+    for (int ii = 0; ii < nq; ++ii){
       dynmat->getDMq(q);
       dynmat->writeDMq(q, qr, fp);
-      for (int i=0; i<3; i++) q[i] += qinc[i];
+      for (int i = 0; i < 3; ++i) q[i] += qinc[i];
       qr += dq;
     }
     qr -= dq;
   }
   fclose(fp);
+
 return;
 }
 
 /* ----------------------------------------------------------------------------
  * Private method to smooth the dos
  * ---------------------------------------------------------------------------- */
 void Phonon::smooth(double *array, const int npt)
 {
   if (npt < 4) return;
 
   int nlag = npt/4;
 
   double *tmp, *table;
-  tmp   = memory->create(tmp, npt, "smooth:tmp");
-  table = memory->create(table, nlag+1, "smooth:table");
+  memory->create(tmp, npt, "smooth:tmp");
+  memory->create(table, nlag+1, "smooth:table");
   
   double fnorm = -1.;
   double sigma = 4., fac = 1./(sigma*sigma);
-  for (int jj=0; jj<= nlag; jj++){
+  for (int jj = 0; jj <= nlag; ++jj){
     table[jj] = exp(-double(jj*jj)*fac);
     fnorm += table[jj];
   }
   fnorm = 1./fnorm;
 
-  for (int i=0; i<npt; i++){
+  for (int i = 0; i < npt; ++i){
     tmp[i] = 0.;
-    for (int jj=-nlag; jj<= nlag; jj++){
+    for (int jj = -nlag; jj <= nlag; ++jj){
       int j = (i+jj+npt)%npt; // assume periodical data
 
       tmp [i] += array[j]*table[abs(jj)];
     }
   }
-  for (int i=0; i<npt; i++) array[i] = tmp[i]*fnorm;
+  for (int i = 0; i < npt; ++i) array[i] = tmp[i]*fnorm;
 
   memory->destroy(tmp);
   memory->destroy(table);
 
 return;
 }
 
 /* ----------------------------------------------------------------------------
  * Private method to calculate the thermal properties
  * ---------------------------------------------------------------------------- */
 void Phonon::therm()
 {
   // get frequencies on a q-mesh
   QMesh();
   ComputeAll();
  
   // get the filename to output thermal properties
   char str[MAXLINE];
 
   printf("\nPlease input the filename to output thermal properties [therm.dat]:");
   if (count_words(fgets(str,MAXLINE,stdin)) < 1) strcpy(str, "therm.dat");
   char *fname = strtok(str," \t\n\r\f");
   FILE *fp = fopen(fname, "a"); fname = NULL;
   // header line 
   fprintf(fp,"#Temp   Uvib    Svib     Fvib    ZPE      Cvib\n");
   fprintf(fp,"# K      eV      Kb       eV      eV       Kb\n");
 
   // constants          J.s             J/K                J
   const double h = 6.62606896e-34, Kb = 1.380658e-23, eV = 1.60217733e-19;
 
   // first temperature
   double T = dynmat->Tmeasure;
   do {
     // constants under the same temperature; assuming angular frequency in THz
     double h_o_KbT = h/(Kb*T)*1.e12, KbT_in_eV = Kb*T/eV;
 
     double Uvib = 0., Svib = 0., Fvib = 0., Cvib = 0., ZPE = 0.;
-    for (int iq=0; iq<nq; iq++){
+    for (int iq = 0; iq < nq; ++iq){
       double Utmp = 0., Stmp = 0., Ftmp = 0., Ztmp = 0., Ctmp = 0.;
-      for (int i=0; i<ndim; i++){
+      for (int i = 0; i < ndim; ++i){
         if (eigs[iq][i] <= 0.) continue;
         double x = eigs[iq][i] * h_o_KbT;
         double expterm = 1./(exp(x)-1.);
         Stmp += x*expterm - log(1.-exp(-x));
         Utmp += (0.5+expterm)*x;
         Ftmp += log(2.*sinh(0.5*x));
         Ctmp += x*x*exp(x)*expterm*expterm;
         Ztmp += 0.5*h*eigs[iq][i];
       }
       Svib += wt[iq]*Stmp;
       Uvib += wt[iq]*Utmp;
       Fvib += wt[iq]*Ftmp;
       Cvib += wt[iq]*Ctmp;
       ZPE  += wt[iq]*Ztmp;
     }
     Uvib *= KbT_in_eV;
     Fvib *= KbT_in_eV;
     ZPE  /= eV*1.e-12;
     // output result under current temperature
     fprintf(fp,"%lg %lg %lg %lg %lg %lg\n", T, Uvib, Svib, Fvib, ZPE, Cvib);
 
     printf("Please input the desired temperature (K), enter to exit: ");
     if (count_words(fgets(str,MAXLINE,stdin)) < 1) break;
     T = atof(strtok(str," \t\n\r\f"));
+
   } while (T > 0.);
   fclose(fp);
 
 return;
 }
 
 /* ----------------------------------------------------------------------------
  * Private method to calculate the local thermal properties
  * ---------------------------------------------------------------------------- */
 void Phonon::local_therm()
 {
   char str[MAXLINE];
   printf("\nWould you like to compute the local thermal properties (y/n)[n]: ");
   if (count_words(fgets(str,MAXLINE,stdin)) < 1) return;
   char *ptr = strtok(str," \t\n\r\f");
   if (strcmp(ptr,"y") != 0 && strcmp(ptr, "Y") != 0 && strcmp(ptr, "yes") != 0) return;
 
   printf("Please input the filename to output vibrational thermal info [localtherm.dat]: ");
   if (count_words(fgets(str,MAXLINE,stdin)) < 1) strcpy(str, "localtherm.dat");
 
   FILE *fp = fopen(strtok(str," \t\n\r\f"), "w");
   fprintf(fp,"# atom Temp  U_vib (eV)    S_vib (kB)    F_vib (eV)    C_vib (kB)     ZPE (eV)\n");
   fprintf(fp,"#           ------------  ------------  -----------   -----------   ------------\n");
   fprintf(fp,"#            Ux Uy Uz Ut   Sx Sy Sz St   Fx Fy Fz Ft   Cx Cy Cz Ct   Zx Zy Zz Zt\n");
   fprintf(fp,"#  1   2     3  4  5  6    7  8  9  10   11 12 13 14   15 16 17 18   19 20 21 22\n");
   fprintf(fp,"#-------------------------------------------------------------------------------\n");
 
   double **Uvib, **Svib, **Fvib, **Cvib, **ZPE;
-  Uvib = memory->create(Uvib,nlocal,sysdim,"local_therm:Uvib");
-  Svib = memory->create(Svib,nlocal,sysdim,"local_therm:Svib");
-  Fvib = memory->create(Fvib,nlocal,sysdim,"local_therm:Fvib");
-  Cvib = memory->create(Cvib,nlocal,sysdim,"local_therm:Cvib");
-  ZPE  = memory->create(ZPE ,nlocal,sysdim,"local_therm:ZPE");
+  memory->create(Uvib,nlocal,sysdim,"local_therm:Uvib");
+  memory->create(Svib,nlocal,sysdim,"local_therm:Svib");
+  memory->create(Fvib,nlocal,sysdim,"local_therm:Fvib");
+  memory->create(Cvib,nlocal,sysdim,"local_therm:Cvib");
+  memory->create(ZPE ,nlocal,sysdim,"local_therm:ZPE");
   // constants          J.s             J/K                J
   const double h = 6.62606896e-34, Kb = 1.380658e-23, eV = 1.60217733e-19;
   double T = dynmat->Tmeasure;
-  while (1){
+  while ( 1 ){
     printf("\nPlease input the temperature at which to evaluate the local vibrational\n");
     printf("thermal properties, non-positive number to exit [%g]: ", T);
     if (count_words(fgets(str,MAXLINE,stdin)) > 0){
       T = atoi(strtok(str," \t\n\r\f"));
       if (T <= 0.) break;
     }
     // constants under the same temperature; assuming angular frequency in THz
     double h_o_KbT = h/(Kb*T)*1.e12, KbT_in_eV = Kb*T/eV;
 
-    for (int i=0; i<nlocal; i++)
-    for (int j=0; j<sysdim; j++) Uvib[i][j] = Svib[i][j] = Fvib[i][j] = Cvib[i][j] = ZPE[i][j] = 0.;
+    for (int i = 0; i <nlocal; ++i)
+    for (int j = 0; j <sysdim; ++j) Uvib[i][j] = Svib[i][j] = Fvib[i][j] = Cvib[i][j] = ZPE[i][j] = 0.;
   
     double freq = fmin-df;
-    for (int i=0; i<ndos; i++){
+    for (int i = 0; i < ndos; ++i){
       freq += df;
       if (freq <= 0.) continue;
   
       double x = freq * h_o_KbT;
       double expterm = 1./(exp(x)-1.);
   
       double Stmp = x*expterm - log(1.-exp(-x));
       double Utmp = (0.5+expterm)*x;
       double Ftmp = log(2.*sinh(0.5*x));
       double Ctmp = x*x*exp(x)*expterm*expterm;
       double Ztmp = 0.5*h*freq;
   
-      for (int il=0; il<nlocal; il++)
-      for (int idim=0; idim<sysdim; idim++){
+      for (int il = 0; il < nlocal; ++il)
+      for (int idim = 0; idim < sysdim; ++idim){
         Uvib[il][idim] += ldos[il][i][idim]*Utmp;
         Svib[il][idim] += ldos[il][i][idim]*Stmp;
         Fvib[il][idim] += ldos[il][i][idim]*Ftmp;
         Cvib[il][idim] += ldos[il][i][idim]*Ctmp;
         ZPE [il][idim] += ldos[il][i][idim]*Ztmp;
       }
     }
-    for (int il=0; il<nlocal; il++)
-    for (int idim=0; idim<sysdim; idim++){
+    for (int il = 0; il < nlocal; ++il)
+    for (int idim = 0; idim < sysdim; ++idim){
       Uvib[il][idim] *= KbT_in_eV*df;
       Svib[il][idim] *= df;
       Fvib[il][idim] *= KbT_in_eV*df;
       Cvib[il][idim] *= df;
       ZPE [il][idim] /= eV*1.e-12*rdf;
     }
 
     // output result under current temperature
-    for (int il=0; il<nlocal; il++){
+    for (int il = 0; il < nlocal; ++il){
       fprintf(fp,"%d %g ", locals[il], T);
       double total = 0.;
-      for (int idim=0; idim<sysdim; idim++){
+      for (int idim = 0; idim < sysdim; ++idim){
         fprintf(fp,"%g ", Uvib[il][idim]);
         total += Uvib[il][idim];
       }
       fprintf(fp,"%g ", total); total = 0.;
 
-      for (int idim=0; idim<sysdim; idim++){
+      for (int idim = 0; idim < sysdim; ++idim){
         fprintf(fp,"%g ", Svib[il][idim]);
         total += Svib[il][idim];
       }
       fprintf(fp,"%g ", total); total = 0.;
 
-      for (int idim=0; idim<sysdim; idim++){
+      for (int idim = 0; idim < sysdim; ++idim){
         fprintf(fp,"%g ", Fvib[il][idim]);
         total += Fvib[il][idim];
       }
       fprintf(fp,"%g ", total); total = 0.;
 
-      for (int idim=0; idim<sysdim; idim++){
+      for (int idim = 0; idim < sysdim; ++idim){
         fprintf(fp,"%g ", Cvib[il][idim]);
         total += Cvib[il][idim];
       }
       fprintf(fp,"%g ", total); total = 0.;
 
-      for (int idim=0; idim<sysdim; idim++){
+      for (int idim = 0; idim < sysdim; ++idim){
         fprintf(fp,"%g ", ZPE[il][idim]);
         total += ZPE[il][idim];
       }
       fprintf(fp,"%g\n", total);
     }
   }
   fclose(fp);
 
 return;
 }
 
 /* ----------------------------------------------------------------------------
  * Private method to generate the q-points from a uniform q-mesh
  * ---------------------------------------------------------------------------- */
 void Phonon::QMesh()
 {
   // ask for mesh info
   char str[MAXLINE];
   int nx = dynmat->nx, ny = dynmat->ny, nz = dynmat->nz;
   printf("\nThe q-mesh size from the read dynamical matrix is: %d x %d x %d\n", nx, ny, nz);
   printf("A denser mesh can be interpolated, but NOTE a too dense mesh can cause segmentation fault.\n");
   printf("Please input your desired q-mesh size [%d %d %d]: ", nx, ny, nz);
   if (count_words(fgets(str,MAXLINE,stdin)) >= 3){
     nx = atoi(strtok(str," \t\n\r\f"));
     ny = atoi(strtok(NULL," \t\n\r\f"));
     nz = atoi(strtok(NULL," \t\n\r\f"));
   }
-  if (nx<1||ny<1||nz<1) return;
+  if (nx < 1||ny < 1||nz < 1) return;
   if (dynmat->nx == 1) nx = 1;
   if (dynmat->ny == 1) ny = 1;
   if (dynmat->nz == 1) nz = 1;
 
 #ifdef UseSPG
   // ask method to generate q-points
   int method = 2;
   printf("Please select your method to generate the q-points:\n");
-  printf("  1. uniform;\n  2. Monkhost-Pack mesh;\nYour choice[2]: ");
+  printf("  1. uniform;\n  2. Monkhost-Pack mesh;\nYour choice [2]: ");
   if (count_words(fgets(str,MAXLINE,stdin)) > 0) method = atoi(strtok(str," \t\n\r\f"));
-  method = 2-method%2;
-  printf("Your selection: %d\n", method);
+  method = 2 - method%2;
+  printf("Your  selection: %d\n", method);
 #endif
  
   memory->destroy(wt);
   memory->destroy(qpts);
 
 #ifdef UseSPG
   if (method == 1){
 #endif
     nq = nx*ny*nz;
     double w = 1./double(nq);
-    wt   = memory->create(wt,   nq, "QMesh:wt");
-    qpts = memory->create(qpts, nq, 3, "QMesh:qpts");
+    memory->create(wt,   nq, "QMesh:wt");
+    memory->create(qpts, nq, 3, "QMesh:qpts");
 
     int iq = 0;
-    for (int i=0; i<nx; i++)
-    for (int j=0; j<ny; j++)
-    for (int k=0; k<nz; k++){
+    for (int i = 0; i < nx; ++i)
+    for (int j = 0; j < ny; ++j)
+    for (int k = 0; k < nz; ++k){
       qpts[iq][0] = double(i)/double(nx);
       qpts[iq][1] = double(j)/double(ny);
       qpts[iq][2] = double(k)/double(nz);
       wt[iq++] = w;
     }
 #ifdef UseSPG
-  }
-  if ((method == 2) && (atpos == NULL)){
-    atpos = memory->create(atpos, dynmat->nucell,3,"QMesh:atpos");
-    attyp = memory->create(attyp, dynmat->nucell,  "QMesh:attyp");
-
-    for (int i=0; i<dynmat->nucell; i++)
-    for (int idim=0; idim<3; idim++) atpos[i][idim] = 0.;
-    for (int i=0; i<3; i++) latvec[i][i] = 1.;
-
-    int flag_lat_info_read = dynmat->flag_latinfo;
-
-    if ( flag_lat_info_read ){ // get unit cell info from binary file; done by dynmat
-      num_atom = dynmat->nucell;
-      // set default, in case system dimension under study is not 3.
-      for (int i=0; i<dynmat->nucell; i++)
-      for (int idim=0; idim<3; idim++) atpos[i][idim] = 0.;
-      for (int i=0; i<3; i++) latvec[i][i] = 1.;
-
-      // get atomic type info
-      for (int i=0; i<num_atom; i++) attyp[i] = dynmat->attyp[i];
-      // get unit cell vector info
-      int ndim = 0;
-      for (int idim=0; idim<3; idim++)
-      for (int jdim=0; jdim<3; jdim++) latvec[jdim][idim] = dynmat->basevec[ndim++];
-      // get atom position in unit cell; fractional
-      for (int i=0; i<num_atom; i++)
-      for (int idim=0; idim<sysdim; idim++) atpos[i][idim] = dynmat->basis[i][idim];
-
-      // display the unit cell info read
-      printf("\n");for (int ii=0; ii<80; ii++) printf("="); printf("\n");
-      printf("The basis vectors of the unit cell:\n");
-      for (int idim=0; idim<3; idim++) printf("  A%d = %lg %lg %lg\n", idim+1, latvec[0][idim], latvec[1][idim], latvec[2][idim]);
-      printf("Atom(s) in the unit cell:\n");
-      printf("  No.  type  sx  sy sz\n");
-      for (int i=0; i<num_atom; i++) printf("  %d %d %lg %lg %lg\n", i+1, attyp[i], atpos[i][0], atpos[i][1], atpos[i][2]);
-      printf("\nIs the above info correct? (y/n)[y]: ");
-      fgets(str,MAXLINE,stdin);
-      char *ptr = strtok(str," \t\n\r\f");
-      if ( (ptr) && ( (strcmp(ptr,"y") != 0) && (strcmp(ptr,"Y") != 0)) ) flag_lat_info_read = 0;
-    }
-
-    if (flag_lat_info_read == 0) { // get unit cell info from file or user input
-      int latsrc = 1;
-      printf("\nPlease select the way to provide the unit cell info:\n");
-      printf("  1. By file;\n  2. Read in.\nYour choice [1]: ");
-      if (count_words(fgets(str,MAXLINE,stdin)) > 0) latsrc = atoi(strtok(str," \t\n\r\f"));
-      latsrc = 2-latsrc%2;
-      /*----------------------------------------------------------------
-       * Ask for lattice info from the user; the format of the file is:
-       * A1_x A1_y A1_z
-       * A2_x A2_y A2_z
-       * A3_x A3_y A3_z
-       * natom
-       * Type_1 sx_1 sy_1 sz_1
-       * ...
-       * Type_n sx_n sy_n sz_n
-       *----------------------------------------------------------------*/
-      if (latsrc == 1){ // to read unit cell info from file; get file name first
-        do printf("Please input the file name containing the unit cell info: ");
-        while (count_words(fgets(str,MAXLINE,stdin)) < 1);
-        char *fname = strtok(str," \t\n\r\f");
-        FILE *fp = fopen(fname,"r"); fname = NULL;
-  
-        if (fp == NULL) latsrc = 2;
-        else {
-          for (int i=0; i<3; i++){ // read unit cell vector info; # of atoms per unit cell
-            if (count_words(fgets(str,MAXLINE,fp)) < 3){
-              latsrc = 2;
-              break;
-            }
-            latvec[0][i] = atof(strtok(str, " \t\n\r\f"));
-            latvec[1][i] = atof(strtok(NULL," \t\n\r\f"));
-            latvec[2][i] = atof(strtok(NULL," \t\n\r\f"));
-          }
-          if (count_words(fgets(str,MAXLINE,fp)) < 1) latsrc = 2;
-          else {
-            num_atom = atoi(strtok(str," \t\n\r\f"));
-            if (num_atom > dynmat->nucell){
-              printf("\nError: # of atoms read from file (%d) is bigger than that given by the dynamical matrix (%d)!\n", num_atom, dynmat->nucell);
-              return;
-            }
-      
-            for (int i=0; i<num_atom; i++){ // read atomic type and fractional positions
-              if (count_words(fgets(str,MAXLINE,fp)) < 4){
-                latsrc = 2;
-                break;
-
-              } else {
-                attyp[i] = atoi(strtok(str," \t\n\r\f"));
-
-                atpos[i][0] = atof(strtok(NULL," \t\n\r\f"));
-                atpos[i][1] = atof(strtok(NULL," \t\n\r\f"));
-                atpos[i][2] = atof(strtok(NULL," \t\n\r\f"));
-              }
-            }
-          }
-        }
-        fclose(fp);
-      }
-      if (latsrc == 2){
-        for (int i=0; i<3; i++){
-          do printf("Please input the vector A%d: ", i+1);
-          while (count_words(fgets(str,MAXLINE,stdin)) < 3);
-          latvec[0][i] = atof(strtok(str," \t\n\r\f"));
-          latvec[1][i] = atof(strtok(NULL," \t\n\r\f"));
-          latvec[2][i] = atof(strtok(NULL," \t\n\r\f"));
-        }
-  
-        do printf("please input the number of atoms per unit cell: ");
-        while (count_words(fgets(str,MAXLINE,stdin)) < 1);
-        num_atom = atoi(strtok(str," \t\n\r\f"));
-        if (num_atom > dynmat->nucell){
-          printf("\nError: # of atoms input (%d) is bigger than that given by the dynamical matrix (%d)!\n", num_atom, dynmat->nucell);
-          return;
-        }
-  
-        for (int i=0; i<num_atom; i++){
-          do printf("Please input the type, and fractional coordinate of atom No.%d: ", i+1);
-          while (count_words(fgets(str,MAXLINE,stdin)) < 4);
-          attyp[i] = atoi(strtok(str," \t\n\r\f"));
-  
-          atpos[i][0] = atof(strtok(NULL," \t\n\r\f"));
-          atpos[i][1] = atof(strtok(NULL," \t\n\r\f"));
-          atpos[i][2] = atof(strtok(NULL," \t\n\r\f"));
-        }
-      }
-    } // end of read from file or input
-  } // end of if (method == 2 && ...
+  } else {
+    if (atpos == NULL) memory->create(atpos, dynmat->nucell,3,"QMesh:atpos");
+    if (attyp == NULL) memory->create(attyp, dynmat->nucell,  "QMesh:attyp");
+
+    num_atom = dynmat->nucell;
+    // set default, in case system dimension under study is not 3.
+    for (int i = 0; i < dynmat->nucell; ++i)
+    for (int idim = 0; idim < 3; ++idim) atpos[i][idim] = 0.;
+    for (int i = 0; i < 3; ++i)
+    for (int j = 0; j < 3; ++j) latvec[i][j] = 0.;
+    for (int i = 0; i < 3; ++i) latvec[i][i] = 1.;
+
+    // get atomic type info
+    for (int i = 0; i < num_atom; ++i) attyp[i] = dynmat->attyp[i];
+
+    // get unit cell vector info
+    int ndim = 0;
+    for (int idim = 0; idim < 3; ++idim)
+    for (int jdim = 0; jdim < 3; ++jdim) latvec[jdim][idim] = dynmat->basevec[ndim++];
+
+    // get atom position in unit cell; fractional
+    for (int i = 0; i < num_atom; ++i)
+    for (int idim = 0; idim < sysdim; ++idim) atpos[i][idim] = dynmat->basis[i][idim];
+
+    // display the unit cell info read
+    printf("\n");for (int ii = 0; ii < 80; ++ii) printf("="); printf("\n");
+    printf("The basis vectors of the unit cell:\n");
+    for (int idim = 0; idim < 3; ++idim) printf("  A%d = %lg %lg %lg\n", idim+1, latvec[0][idim], latvec[1][idim], latvec[2][idim]);
+    printf("Atom(s) in the unit cell:\n");
+    printf("  No.  type  sx  sy sz\n");
+    for (int i = 0; i < MIN(num_atom, NUMATOM); ++i) printf("  %d %d %lg %lg %lg\n", i+1, attyp[i], atpos[i][0], atpos[i][1], atpos[i][2]);
+    if (num_atom > NUMATOM) printf("  ... (%d atoms omitted.)\n", num_atom - NUMATOM);
 
-  if (method == 2){
     int mesh[3], shift[3], is_time_reversal = 0;
     mesh[0] = nx; mesh[1] = ny; mesh[2] = nz;
     shift[0] = shift[1] = shift[2] = 0;
     int num_grid = mesh[0]*mesh[1]*mesh[2];
     int grid_point[num_grid][3], map[num_grid];
     double symprec = 1.e-4, pos[num_atom][3];
 
-    for (int i=0; i<num_atom; i++)
-    for (int j=0; j<3; j++) pos[i][j] = atpos[i][j];
+    for (int i = 0; i < num_atom; ++i)
+    for (int j = 0; j < 3; ++j) pos[i][j] = atpos[i][j];
 
-    //spg_show_symmetry(latvec, pos, attyp,  num_atom, symprec);
 
     // if spglib 0.7.1 is used
-    nq = spg_get_ir_reciprocal_mesh(grid_point, map, num_grid, mesh, shift, is_time_reversal, latvec, pos, attyp, num_atom, symprec);
+    // nq = spg_get_ir_reciprocal_mesh(grid_point, map, num_grid, mesh, shift, is_time_reversal, latvec, pos, attyp, num_atom, symprec);
 
     // if spglib >= 1.0.3 is used
-    //nq = spg_get_ir_reciprocal_mesh(grid_point, map, mesh, shift, is_time_reversal, latvec, pos, attyp, num_atom, symprec);
+    nq = spg_get_ir_reciprocal_mesh(grid_point, map, mesh, shift, is_time_reversal, latvec, pos, attyp, num_atom, symprec);
 
-    wt   = memory->create(wt,   nq, "QMesh:wt");
-    qpts = memory->create(qpts, nq,3,"QMesh:qpts");
+    memory->create(wt,   nq, "QMesh:wt");
+    memory->create(qpts, nq,3,"QMesh:qpts");
 
     int *iq2idx = new int[num_grid];
     int numq = 0;
-    for (int i=0; i<num_grid; i++){
+    for (int i = 0; i < num_grid; ++i){
       int iq = map[i];
       if (iq == i) iq2idx[iq] = numq++;
     }
-    for (int iq=0; iq<nq; iq++) wt[iq] = 0.;
+    for (int iq = 0; iq < nq; ++iq) wt[iq] = 0.;
     numq = 0;
-    for (int i=0; i<num_grid; i++){
+    for (int i = 0; i < num_grid; ++i){
       int iq = map[i];
       if (iq == i){
         qpts[numq][0] = double(grid_point[i][0])/double(mesh[0]);
         qpts[numq][1] = double(grid_point[i][1])/double(mesh[1]);
         qpts[numq][2] = double(grid_point[i][2])/double(mesh[2]);
         numq++;
       }
       wt[iq2idx[iq]] += 1.;
     }
     delete []iq2idx;
 
     double wsum = 0.;
-    for (int iq=0; iq<nq; iq++) wsum += wt[iq];
-    for (int iq=0; iq<nq; iq++) wt[iq] /= wsum;
+    for (int iq = 0; iq < nq; ++iq) wsum += wt[iq];
+    for (int iq = 0; iq < nq; ++iq) wt[iq] /= wsum;
     
   }
 #endif
   printf("Your new q-mesh size would be: %d x %d x %d => %d points\n", nx,ny,nz,nq);
 
 return;
 }
 
 /* ----------------------------------------------------------------------------
  * Private method to calculate the local phonon DOS and total phonon DOS based
  * on the eigenvectors
  * ---------------------------------------------------------------------------- */
 void Phonon::ldos_egv()
 {
   // get local position info
   char str[MAXLINE], *ptr;
   printf("\nThe # of atoms per cell is: %d, please input the atom IDs to compute\n", dynmat->nucell);
   printf("local PDOS, IDs begin with 0: ");
   int nmax = count_words(fgets(str,MAXLINE,stdin));
   if (nmax < 1) return;
 
   memory->destroy(locals);
-  locals = memory->create(locals, nmax, "ldos_egv:locals");
+  memory->create(locals, nmax, "ldos_egv:locals");
 
   nlocal = 0;
   ptr = strtok(str," \t\n\r\f");
   while (ptr != NULL){
     int id = atoi(ptr);
     if (id >= 0 && id < dynmat->nucell) locals[nlocal++] = id;
 
     ptr = strtok(NULL," \t\n\r\f");
   }
   if (nlocal < 1) return;
 
   printf("Local PDOS for atom(s):");
-  for (int i=0; i<nlocal; i++) printf(" %d", locals[i]);
+  for (int i = 0; i < nlocal; ++i) printf(" %d", locals[i]);
   printf("  will be computed.\n");
 
   fmin = 0.; fmax = 10.;
   printf("Please input the freqency (nv, THz) range to compute PDOS [%g %g]: ", fmin, fmax);
   if (count_words(fgets(str,MAXLINE,stdin)) >= 2) {
     fmin = atof(strtok(str," \t\n\r\f"));
     fmax = atof(strtok(NULL," \t\n\r\f"));
   }
   if (fmax < 0. || fmax < fmin) return;
 
   ndos = 201;
   printf("Please input your desired # of points in PDOS [%d]: ", ndos);
   if (count_words(fgets(str,MAXLINE,stdin)) > 0) ndos = atoi(strtok(str," \t\n\r\f"));
   if (ndos < 2) return;
   ndos += (ndos+1)%2;
 
   df = (fmax-fmin)/double(ndos-1);
   rdf = 1./df;
 
   // get the q-points
   QMesh();
 
   // allocate memory for DOS and LDOS
   memory->destroy(dos);
   memory->destroy(ldos);
 
-  dos  = memory->create(dos, ndos,"ldos_egv:dos");
-  ldos = memory->create(ldos,nlocal,ndos,sysdim,"ldos_egv:ldos");
+  memory->create(dos, ndos,"ldos_egv:dos");
+  memory->create(ldos,nlocal,ndos,sysdim,"ldos_egv:ldos");
 
-  for (int i=0; i<ndos; i++) dos[i] = 0.;
+  for (int i = 0; i < ndos; ++i) dos[i] = 0.;
 
-  for (int ilocal=0; ilocal<nlocal; ilocal++)
-  for (int i=0; i<ndos; i++)
-  for (int idim=0; idim<sysdim; idim++) ldos[ilocal][i][idim] = 0.;
+  for (int ilocal = 0; ilocal < nlocal; ++ilocal)
+  for (int i = 0; i < ndos; ++i)
+  for (int idim = 0; idim < sysdim; ++idim) ldos[ilocal][i][idim] = 0.;
 
   int nprint;
   if (nq > 10) nprint = nq/10;
   else nprint = 1;
   Timer *time = new Timer();
 
   // memory and pointer for eigenvalues and eigenvectors
   double egval[ndim], offset=fmin-0.5*df;
   doublecomplex **egvec = dynmat->DM_q;
 
   printf("\nNow to compute the phonons and DOSs "); fflush(stdout);
-  for (int iq=0; iq<nq; iq++){
+  for (int iq = 0; iq < nq; ++iq){
     if ((iq+1)%nprint == 0) {printf("."); fflush(stdout);}
 
     dynmat->getDMq(qpts[iq], &wt[iq]);
     if (wt[iq] <= 0.) continue;
 
     dynmat->geteigen(&egval[0], 1);
 
-    for (int idim=0; idim<ndim; idim++){
+    for (int idim = 0; idim < ndim; ++idim){
       int hit = int((egval[idim] - offset)*rdf);
       if (hit >= 0 && hit <ndos){
         dos[hit] += wt[iq];
 
-        for (int ilocal=0; ilocal<nlocal; ilocal++){
+        for (int ilocal = 0; ilocal < nlocal; ++ilocal){
           int ipos = locals[ilocal]*sysdim;
-          for (int jdim=0; jdim<sysdim; jdim++){
+          for (int jdim = 0; jdim < sysdim; ++jdim){
             double dr = egvec[idim][ipos+jdim].r, di = egvec[idim][ipos+jdim].i;
             double norm = dr * dr + di * di;
             ldos[ilocal][hit][jdim] += wt[iq] * norm;
           }
         }
       }
     }
   }
   egvec = NULL;
   printf("Done!\nNow to normalize the DOSs ..."); fflush(stdout);
 
   // normalize the measure DOS and LDOS
   Normalize();
   printf("Done! ");
   time->stop(); time->print(); delete time;
 
   // to write the DOSes
   writeDOS();
   writeLDOS();
 
   // evaluate the local vibrational thermal properties optionally
   local_therm();
 
 return;
 }
 
+/* ----------------------------------------------------------------------------
+ * Private method to show the unit cell info
+ * ---------------------------------------------------------------------------- */
+void Phonon::ShowCell()
+{
+  printf("\n");
+  for (int i = 0; i < 30; ++i) printf("=");
+  printf("   Unit Cell Info   ");
+  for (int i = 0; i < 30; ++i) printf("="); printf("\n");
+  printf("Number of atoms in the unit cell: %d\n", dynmat->nucell);
+  printf("Basis  vectors  of the unit cell:\n");
+  printf("  %15.8f  %15.8f  %15.8f\n", dynmat->basevec[0],  dynmat->basevec[1],  dynmat->basevec[2]);
+  printf("  %15.8f  %15.8f  %15.8f\n", dynmat->basevec[3],  dynmat->basevec[4],  dynmat->basevec[5]);
+  printf("  %15.8f  %15.8f  %15.8f\n", dynmat->basevec[6],  dynmat->basevec[7],  dynmat->basevec[8]);
+  printf("Basis  vectors  of the reciprocal:\n");
+  printf("  %15.8f  %15.8f  %15.8f\n", dynmat->ibasevec[0],  dynmat->ibasevec[1],  dynmat->ibasevec[2]);
+  printf("  %15.8f  %15.8f  %15.8f\n", dynmat->ibasevec[3],  dynmat->ibasevec[4],  dynmat->ibasevec[5]);
+  printf("  %15.8f  %15.8f  %15.8f\n", dynmat->ibasevec[6],  dynmat->ibasevec[7],  dynmat->ibasevec[8]);
+  printf("Atomic type and fractional coordinates:\n");
+  for (int i = 0; i < dynmat->nucell; ++i)
+    printf("%4d %12.8f %12.8f %12.8f\n", dynmat->attyp[i], dynmat->basis[i][0], dynmat->basis[i][1], dynmat->basis[i][2]);
+  for (int i = 0; i < 80; ++i) printf("=");
+  printf("\n");
+
+return;
+}
+
 /* ----------------------------------------------------------------------------
  * Private method to normalize the DOS and/or Local DOS.
  * Simpson's rule is used for the integration.
  * ---------------------------------------------------------------------------- */
 void Phonon::Normalize()
 {
   double odd, even, sum;
   if (dos){
     odd = even = 0.;
-    for (int i=1; i<ndos-1; i +=2) odd  += dos[i];
-    for (int i=2; i<ndos-1; i +=2) even += dos[i];
+    for (int i = 1; i < ndos-1; i +=2) odd  += dos[i];
+    for (int i = 2; i < ndos-1; i +=2) even += dos[i];
     sum = dos[0] + dos[ndos-1];
     sum += 4.*odd + 2.*even;
     sum = 3.*rdf/sum;
-    for (int i=0; i<ndos; i++) dos[i] *= sum;
+    for (int i = 0; i < ndos; ++i) dos[i] *= sum;
   }
 
   if (ldos){
-    for (int ilocal=0; ilocal<nlocal; ilocal++)
-    for (int idim=0; idim<sysdim; idim++){
+    for (int ilocal = 0; ilocal < nlocal; ++ilocal)
+    for (int idim = 0; idim < sysdim; ++idim){
       odd = even = 0.;
-      for (int i=1; i<ndos-1; i +=2) odd  += ldos[ilocal][i][idim];
-      for (int i=2; i<ndos-1; i +=2) even += ldos[ilocal][i][idim];
+      for (int i = 1; i < ndos-1; i += 2) odd  += ldos[ilocal][i][idim];
+      for (int i = 2; i < ndos-1; i += 2) even += ldos[ilocal][i][idim];
       sum = ldos[ilocal][0][idim] + ldos[ilocal][ndos-1][idim];
       sum += 4.*odd + 2.*even;
       sum = 3.*rdf/sum;
-      for (int i=0; i<ndos; i++) ldos[ilocal][i][idim] *= sum;
+      for (int i = 0; i < ndos; ++i) ldos[ilocal][i][idim] *= sum;
     }
   }
 
 return;
 }
 
 /* ----------------------------------------------------------------------------
  * Private method to calculate vibrational frequencies for all q-points
  * ---------------------------------------------------------------------------- */
 void Phonon::ComputeAll()
 {
   int nprint;
   if (nq > 10) nprint = nq/10;
   else nprint = 1;
   Timer *time = new Timer();
 
   printf("\nNow to compute the phonons "); fflush(stdout);
   // now to calculate the frequencies at all q-points
   memory->destroy(eigs);
-  eigs = memory->create(eigs, nq,ndim,"QMesh_eigs");
+  memory->create(eigs, nq,ndim,"QMesh_eigs");
   
-  for (int iq=0; iq<nq; iq++){
+  for (int iq = 0; iq < nq; ++iq){
     if ((iq+1)%nprint == 0) {printf("."); fflush(stdout);}
 
     dynmat->getDMq(qpts[iq], &wt[iq]);
     if (wt[iq] > 0.) dynmat->geteigen(eigs[iq], 0);
   }
   printf("Done!\n");
   time->stop(); time->print(); delete time;
 
 return;
 }
 
 /*------------------------------------------------------------------------------
  * Method to count # of words in a string, without destroying the string
  *----------------------------------------------------------------------------*/
 int Phonon::count_words(const char *line)
 {
   int n = strlen(line) + 1;
   char *copy;
-  copy = memory->create(copy, n, "count_words:copy");
+  memory->create(copy, n, "count_words:copy");
   strcpy(copy,line);
 
   char *ptr;
   if (ptr = strchr(copy,'#')) *ptr = '\0';
 
   if (strtok(copy," \t\n\r\f") == NULL) {
     memory->destroy(copy);
     return 0;
   }
   n = 1;
   while (strtok(NULL," \t\n\r\f")) n++;
 
   memory->destroy(copy);
   return n;
 }
 
 /*----------------------------------------------------------------------------*/
diff --git a/tools/phonon/phonon.h b/tools/phonon/phonon.h
index 02552fc40..80d52429b 100644
--- a/tools/phonon/phonon.h
+++ b/tools/phonon/phonon.h
@@ -1,59 +1,61 @@
 #ifndef PHONON_H
 #define PHONON_H
 
 #include "stdio.h"
 #include "stdlib.h"
 #include <complex>
 #include "dynmat.h"
 #include "memory.h"
 
 using namespace std;
 
 class Phonon{
 public:
   Phonon(DynMat *);
   ~Phonon();
 
   DynMat *dynmat;
 
 private:
   int nq, ndim, sysdim;
   double **qpts, *wt;
   double **eigs;
 
   int ndos, nlocal, *locals;
   double *dos, fmin, fmax, df, rdf;
   double ***ldos;
 
   Memory *memory;
 
   void QMesh();
   void ComputeAll();
 
   void pdos();
   void pdisp();
   void therm();
 
   void ldos_egv();
   void ldos_rsgf();
   void local_therm();
 
   void dmanyq();
   void vfanyq();
   void DMdisp();
   void vecanyq();
 
+  void ShowCell();
+
   void smooth(double *, const int);
   void writeDOS();
   void writeLDOS();
   void Normalize();
 
   int count_words(const char *);
 
 #ifdef UseSPG
   int num_atom, *attyp;
   double latvec[3][3], **atpos;
 #endif
 };
 
 #endif
diff --git a/tools/phonon/timer.cpp b/tools/phonon/timer.cpp
index 0b66a9572..af1d9950d 100644
--- a/tools/phonon/timer.cpp
+++ b/tools/phonon/timer.cpp
@@ -1,41 +1,55 @@
 #include "timer.h"
-
+/* -----------------------------------------------------------------------------
+ * Initialization of time
+ * -------------------------------------------------------------------------- */
 Timer::Timer()
 {
   flag = 0;
   start();
 return;
 }
 
+/* -----------------------------------------------------------------------------
+ * public function, start the timer
+ * -------------------------------------------------------------------------- */
 void Timer::start()
 {
  t1 = clock();
  flag |= 1;
 
 return;
 }
 
+/* -----------------------------------------------------------------------------
+ * public function, stop the timer
+ * -------------------------------------------------------------------------- */
 void Timer::stop()
 {
   if ( flag&1 ) {
     t2 = clock();
     flag |= 2;
   }
 return;
 }
 
+/* -----------------------------------------------------------------------------
+ * public function, print the total time used after timer stops
+ * -------------------------------------------------------------------------- */
 void Timer::print()
 {
   if ( (flag&3) != 3) return;
 
   cpu_time_used = ((double) (t2 - t1)) / CLOCKS_PER_SEC;
   printf("Total CPU time used: %g seconds.\n", cpu_time_used);
 
 return;
 }
 
+/* -----------------------------------------------------------------------------
+ * public function, return the total time used up to now, in seconds
+ * -------------------------------------------------------------------------- */
 double Timer::elapse()
 {
   if ( (flag&3) != 3) return 0.;
   else return ((double) (t2 - t1)) / CLOCKS_PER_SEC;
 }
diff --git a/tools/phonon/version.h b/tools/phonon/version.h
index 69817dc82..8ed0e80aa 100644
--- a/tools/phonon/version.h
+++ b/tools/phonon/version.h
@@ -1 +1 @@
-#define VERSION  75
+#define VERSION 7