From 5f802f86b54f2acb9f74bdcebcafcf907b9abdf6 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Tue, 15 Dec 2015 15:59:01 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14366 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- tools/phonon/Makefile | 25 +- tools/phonon/README | 73 +++-- tools/phonon/dynmat.cpp | 217 +++++++-------- tools/phonon/dynmat.h | 2 +- tools/phonon/green.cpp | 47 ++-- tools/phonon/interpolate.cpp | 53 ++-- tools/phonon/phonon.cpp | 505 ++++++++++++----------------------- tools/phonon/phonon.h | 2 + tools/phonon/timer.cpp | 16 +- tools/phonon/version.h | 2 +- 10 files changed, 399 insertions(+), 543 deletions(-) diff --git a/tools/phonon/Makefile b/tools/phonon/Makefile index 6c39f55cea..0aacb1e086 100644 --- a/tools/phonon/Makefile +++ b/tools/phonon/Makefile @@ -9,23 +9,20 @@ 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 @@ -39,7 +36,7 @@ SRC = $(wildcard *.cpp) OBJ = $(SRC:.cpp=.o) #==================================================================== -all: ${EXE} +all: ver ${EXE} ${EXE}: $(OBJ) $(LINK) $(OFLAGS) $(OBJ) $(LIB) -o $@ @@ -51,7 +48,7 @@ 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: diff --git a/tools/phonon/README b/tools/phonon/README index 61e458e5fe..ae6383b6bd 100644 --- a/tools/phonon/README +++ b/tools/phonon/README @@ -1,31 +1,48 @@ -phana = post-processing program for the fix phonon command +#------------------------------------------------------------------------------- + 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". -This program reads the binary file created by the fix phonon command -and helps to analyse the phonon related info. +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$. -The clapack library is needed to solve the eigen problems, which could -be downloaded from: http://www.netlib.org/clapack/ +4. Updates + For updates of phana, please check: + http://nes.sjtu.edu.cn/english/research/software.htm + or + https://github.com/lingtikong/phana.git -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 +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 1886901647..eb91647ed2 100644 --- a/tools/phonon/dynmat.cpp +++ b/tools/phonon/dynmat.cpp @@ -1,8 +1,7 @@ #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) @@ -14,6 +13,8 @@ DynMat::DynMat(int narg, char **arg) DM_q = DM_all = NULL; binfile = funit = dmfile = NULL; + attyp = NULL; + basis = NULL; flag_reset_gamma = flag_skip = 0; // analyze the command line options @@ -74,23 +75,23 @@ DynMat::DynMat(int narg, char **arg) 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 { @@ -100,9 +101,9 @@ DynMat::DynMat(int narg, char **arg) } // 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)){ @@ -112,48 +113,36 @@ DynMat::DynMat(int narg, char **arg) } // 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. + // Enforcing Austic Sum Rule + EnforceASR(); - 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; idimcreate(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= 0.) w[i] = sqrt(w[i]); else w[i] = -sqrt(-w[i]); @@ -298,24 +287,17 @@ 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= 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); } } @@ -501,7 +487,7 @@ void DynMat::GaussJordan(int n, double *Mat) } ipiv[icol] += 1; if (irow != icol){ - for (l=0; 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 #include "green.h" #include - -#define MAXLINE 256 +#include "global.h" /******************************************************************************* * The class of Green is designed to evaluate the LDOS via the Green's Function @@ -45,7 +44,7 @@ Green::Green(const int ntm, const int sdim, const int niter, const double min, c 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; @@ -58,9 +57,9 @@ Green::Green(const int ntm, const int sdim, const int niter, const double min, c // 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(); @@ -100,35 +99,35 @@ void Green::Lanczos() int ipos = iatom*sysdim; // Loop over dimension - for (int idim=0; idim(w*w, epson); - for (int idim=0; idim (sr, si); - for (int j=0; j(w*w, epson); - for (int idim=0; idim(0.,0.); - for (int j=0; j(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 @@ -15,7 +12,7 @@ Interpolate::Interpolate(int nx, int ny, int nz, int ndm, doublecomplex **DM) Nz = nz; Npt = Nx*Ny*Nz; ndim = ndm; - memory = new Memory; + memory = new Memory(); which = UseGamma = 0; @@ -33,13 +30,13 @@ 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; } @@ -47,9 +44,9 @@ void Interpolate::tricubic_init() // get the derivatives int n=0; const double half = 0.5, one4 = 0.25, one8 = 0.125; - for (int ii=0; ii= 1.) q[i] -= 1.; } @@ -150,8 +147,8 @@ void Interpolate::tricubic(double *qin, doublecomplex *DMq) 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= 1.) q[i] -= 1.; } @@ -220,7 +219,7 @@ void Interpolate::trilinear(double *qin, doublecomplex *DMq) 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); @@ -233,10 +232,10 @@ void Interpolate::trilinear(double *qin, doublecomplex *DMq) fac[7] = x*y*z; // now to do the interpolation - for (int idim=0; idim(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. @@ -42,7 +39,10 @@ Phonon::Phonon(DynMat *dm) // 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"); @@ -52,14 +52,15 @@ Phonon::Phonon(DynMat *dm) 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(); @@ -71,6 +72,7 @@ Phonon::Phonon(DynMat *dm) 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; } @@ -113,11 +115,10 @@ void Phonon::pdos() char str[MAXLINE]; fmin = fmax = eigs[0][0]; - for (int iq=0; iqdestroy(dos); - dos = memory->create(dos, ndos, "pdos:dos"); + memory->create(dos, ndos, "pdos:dos"); for (int i=0; i 0.){ - for (int j=0; j=0 && idxeml2f*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; iDM_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= 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); @@ -332,7 +341,7 @@ void Phonon::ldos_rsgf() 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; @@ -366,104 +375,6 @@ void Phonon::ldos_rsgf() 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 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; iigetDMq(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 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; igeteigen(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; igeteigen(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; inucell; 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; idimnx,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)); @@ -602,16 +512,17 @@ void Phonon::DMdisp() 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; iigetDMq(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; } @@ -625,26 +536,26 @@ void Phonon::smooth(double *array, const int npt) 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; idestroy(tmp); memory->destroy(table); @@ -682,9 +593,9 @@ void Phonon::therm() 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 0.); fclose(fp); @@ -737,15 +649,15 @@ void Phonon::local_therm() 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){ @@ -755,11 +667,11 @@ void Phonon::local_therm() // 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; inx == 1) nx = 1; if (dynmat->ny == 1) ny = 1; if (dynmat->nz == 1) nz = 1; @@ -855,10 +767,10 @@ void Phonon::QMesh() // 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); @@ -869,149 +781,52 @@ void Phonon::QMesh() #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; icreate(atpos, dynmat->nucell,3,"QMesh:atpos"); - attyp = memory->create(attyp, dynmat->nucell, "QMesh:attyp"); + } else { + if (atpos == NULL) memory->create(atpos, dynmat->nucell,3,"QMesh:atpos"); + if (attyp == NULL) memory->create(attyp, dynmat->nucell, "QMesh:attyp"); - for (int i=0; inucell; i++) - for (int idim=0; idim<3; idim++) atpos[i][idim] = 0.; - for (int i=0; i<3; i++) latvec[i][i] = 1.; + 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.; - int flag_lat_info_read = dynmat->flag_latinfo; + // get atomic type info + for (int i = 0; i < num_atom; ++i) attyp[i] = dynmat->attyp[i]; - 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; inucell; i++) - for (int idim=0; idim<3; idim++) atpos[i][idim] = 0.; - for (int i=0; i<3; i++) latvec[i][i] = 1.; + // 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 atomic type info - for (int i=0; iattyp[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; ibasis[i][idim]; + // 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 NUMATOM) printf(" ... (%d atoms omitted.)\n", num_atom - NUMATOM); - 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 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= 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; idestroy(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"); @@ -1091,7 +905,7 @@ void Phonon::ldos_egv() if (nlocal < 1) return; printf("Local PDOS for atom(s):"); - for (int i=0; idestroy(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 10) nprint = nq/10; @@ -1137,7 +951,7 @@ void Phonon::ldos_egv() doublecomplex **egvec = dynmat->DM_q; printf("\nNow to compute the phonons and DOSs "); fflush(stdout); - for (int iq=0; iqgetDMq(qpts[iq], &wt[iq]); @@ -1145,14 +959,14 @@ void Phonon::ldos_egv() dynmat->geteigen(&egval[0], 1); - for (int idim=0; idim= 0 && hit 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. @@ -1188,24 +1029,24 @@ void Phonon::Normalize() double odd, even, sum; if (dos){ odd = even = 0.; - for (int i=1; idestroy(eigs); - eigs = memory->create(eigs, nq,ndim,"QMesh_eigs"); + memory->create(eigs, nq,ndim,"QMesh_eigs"); - for (int iq=0; iqgetDMq(qpts[iq], &wt[iq]); @@ -1246,7 +1087,7 @@ 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; diff --git a/tools/phonon/phonon.h b/tools/phonon/phonon.h index 02552fc409..80d52429be 100644 --- a/tools/phonon/phonon.h +++ b/tools/phonon/phonon.h @@ -43,6 +43,8 @@ private: void DMdisp(); void vecanyq(); + void ShowCell(); + void smooth(double *, const int); void writeDOS(); void writeLDOS(); diff --git a/tools/phonon/timer.cpp b/tools/phonon/timer.cpp index 0b66a9572e..af1d9950d9 100644 --- a/tools/phonon/timer.cpp +++ b/tools/phonon/timer.cpp @@ -1,5 +1,7 @@ #include "timer.h" - +/* ----------------------------------------------------------------------------- + * Initialization of time + * -------------------------------------------------------------------------- */ Timer::Timer() { flag = 0; @@ -7,6 +9,9 @@ Timer::Timer() return; } +/* ----------------------------------------------------------------------------- + * public function, start the timer + * -------------------------------------------------------------------------- */ void Timer::start() { t1 = clock(); @@ -15,6 +20,9 @@ void Timer::start() return; } +/* ----------------------------------------------------------------------------- + * public function, stop the timer + * -------------------------------------------------------------------------- */ void Timer::stop() { if ( flag&1 ) { @@ -24,6 +32,9 @@ void Timer::stop() return; } +/* ----------------------------------------------------------------------------- + * public function, print the total time used after timer stops + * -------------------------------------------------------------------------- */ void Timer::print() { if ( (flag&3) != 3) return; @@ -34,6 +45,9 @@ void Timer::print() return; } +/* ----------------------------------------------------------------------------- + * public function, return the total time used up to now, in seconds + * -------------------------------------------------------------------------- */ double Timer::elapse() { if ( (flag&3) != 3) return 0.; diff --git a/tools/phonon/version.h b/tools/phonon/version.h index 69817dc824..8ed0e80aa7 100644 --- a/tools/phonon/version.h +++ b/tools/phonon/version.h @@ -1 +1 @@ -#define VERSION 75 +#define VERSION 7