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