diff --git a/tools/phonon/Makefile b/tools/phonon/Makefile index 0aacb1e08..67f9b91fd 100644 --- a/tools/phonon/Makefile +++ b/tools/phonon/Makefile @@ -1,61 +1,73 @@ .SUFFIXES : .o .cpp # compiler and flags -CC = g++ -Wno-unused-result -LINK = $(CC) -static +CC = g++ -Wall +LINK = $(CC) CFLAGS = -O3 $(DEBUG) $(UFLAG) # OFLAGS = -O3 $(DEBUG) INC = $(LPKINC) $(TCINC) $(SPGINC) LIB = $(LPKLIB) $(TCLIB) $(SPGLIB) # # cLapack library needed -LPKINC = -I/opt/libs/clapack/3.2.1/include -LPKLIB = -L/opt/libs/clapack/3.2.1/lib -lclapack -lblas -lf2c #-lm +LPKINC = +LPKLIB =-llapack # -# Tricubic library needed -TCINC = -I/opt/libs/tricubic/1.0/include -TCLIB = -L/opt/libs/tricubic/1.0/lib -ltricubic # # 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/libs/spglib/1.8.2/include -SPGLIB = -L/opt/libs/spglib/1.8.2/lib -lsymspg + +# UFLAG = -DUseSPG +# 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: ver ${EXE} +all: ${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 `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 $< + +#==================================================================== +# dependencies +disp.o: disp.cpp phonon.h dynmat.h memory.h interpolate.h green.h timer.h \ + global.h +dynmat.o: dynmat.cpp dynmat.h memory.h interpolate.h version.h global.h +green.o: green.cpp green.h memory.h global.h +interpolate.o: interpolate.cpp interpolate.h memory.h global.h +main.o: main.cpp dynmat.h memory.h interpolate.h phonon.h +memory.o: memory.cpp memory.h +phonon.o: phonon.cpp phonon.h dynmat.h memory.h interpolate.h green.h \ + timer.h global.h +timer.o: timer.cpp timer.h diff --git a/tools/phonon/README b/tools/phonon/README index ae6383b6b..b54d96d8a 100644 --- a/tools/phonon/README +++ b/tools/phonon/README @@ -1,48 +1,42 @@ #------------------------------------------------------------------------------- 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 LAPACK library is needed to solve the eigen problems. + http://www.netlib.org/lapack/ + Intel MKL can be used as well. 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/disp.cpp b/tools/phonon/disp.cpp index 2fa603916..218e01e7f 100644 --- a/tools/phonon/disp.cpp +++ b/tools/phonon/disp.cpp @@ -1,2854 +1,2855 @@ #include #include "string.h" #include "phonon.h" #include "green.h" #include "timer.h" #include "global.h" #ifdef UseSPG extern "C"{ #include "spglib.h" } #endif /*------------------------------------------------------------------------------ * Private method to evaluate the phonon dispersion curves *----------------------------------------------------------------------------*/ void Phonon::pdisp() { // ask the output file name and write the header. char str[MAXLINE]; - for (int ii = 0; ii < 80; ++ii) printf("="); printf("\n"); + for (int ii = 0; ii < 80; ++ii) printf("="); + printf("\n"); #ifdef UseSPG // ask method to generate q-lines int method = 2; printf("Please select your method to generate the phonon dispersion:\n"); printf(" 1. Manual, should always work;\n"); printf(" 2. Automatic, works only for 3D crystals (CMS49-299).\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); #endif printf("\nPlease 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); // to store the nodes of the dispersion curve std::vector nodes; nodes.clear(); std::vector ndstr; ndstr.clear(); std::vector qs, qe; qs.clear(); qe.clear(); std::vector nqbin; nqbin.clear(); // now the calculate the dispersion curve double qstr[3], qend[3]; int nq = MAX(MAX(dynmat->nx,dynmat->ny),dynmat->nz)/2+1; qend[0] = qend[1] = qend[2] = 0.; double *egvs = new double [ndim]; #ifdef UseSPG if (method == 1){ #endif 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); double *qtmp = new double [3]; for (int i = 0; i < 3; ++i) qtmp[i] = qstr[i]; qs.push_back(qtmp); qtmp = new double [3]; for (int i = 0; i < 3; ++i) qtmp[i] = qend[i]; qe.push_back(qtmp); nqbin.push_back(nq); ndstr.push_back(""); } ndstr.push_back(""); #ifdef UseSPG } else { memory->grow(atpos, dynmat->nucell, 3, "pdisp:atpos"); memory->grow(attyp, dynmat->nucell, "pdisp: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); char symbol[11]; 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]; int spgnum = spg_get_international(symbol, latvec, pos, attyp, num_atom, symprec); printf("The space group number of your unit cell is: %d => %s\n", spgnum, symbol); for (int ii = 0; ii < 80; ++ii) printf("-"); printf("\n"); // angles const double la = sqrt(latvec[0][0] * latvec[0][0] + latvec[0][1] * latvec[0][1] + latvec[0][2] * latvec[0][2]); const double lb = sqrt(latvec[1][0] * latvec[1][0] + latvec[1][1] * latvec[1][1] + latvec[1][2] * latvec[1][2]); const double lc = sqrt(latvec[2][0] * latvec[2][0] + latvec[2][1] * latvec[2][1] + latvec[2][2] * latvec[2][2]); const double cosa = sqrt(latvec[1][0] * latvec[2][0] + latvec[1][1] * latvec[2][1] + latvec[1][2] * latvec[2][2])/(lb*lc); const double cosg = sqrt(latvec[0][0] * latvec[1][0] + latvec[0][1] * latvec[1][1] + latvec[0][2] * latvec[1][2])/(la*lb); double ivec[3][3]; ndim = 0; for (int idim = 0; idim < 3; ++idim) for (int jdim = 0; jdim < 3; ++jdim) ivec[jdim][idim] = dynmat->ibasevec[ndim++]; const double ka = sqrt(ivec[0][0] * ivec[0][0] + ivec[0][1] * ivec[0][1] + ivec[0][2] * ivec[0][2]); const double kb = sqrt(ivec[1][0] * ivec[1][0] + ivec[1][1] * ivec[1][1] + ivec[1][2] * ivec[1][2]); const double coskg = sqrt(ivec[0][0] * ivec[1][0] + ivec[0][1] * ivec[1][1] + ivec[0][2] * ivec[1][2])/(ka*kb); double *qtmp; if (spgnum >= 1 && spgnum <= 2){ // Triclinic system if (fabs(coskg) > ZERO){ // A.14, TRI1a and TRI2a ndstr.push_back("X"); // X-G qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("{/Symbol G}"); // G-Y qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.; qtmp[1] = 0.5; qe.push_back(qtmp); ndstr.push_back("Y/L"); // L-G qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("{/Symbol G}"); // G-Z qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("Z/N"); // N-G qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("{/Symbol G}"); // G-M qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.; qtmp[1] = qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("M/R"); // R-G qtmp = new double [3]; qtmp[0] = qtmp[2] = qtmp[1] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("{/Symbol G}"); } else { // A.14, TRI1b and TRI2b ndstr.push_back("X"); // X-G qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.; qtmp[1] = -0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("{/Symbol G}"); // G-Y qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[1] = qtmp[2] = 0.; qtmp[0] = 0.5; qe.push_back(qtmp); ndstr.push_back("Y/L"); // L-G qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = -0.5; qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("{/Symbol G}"); // G-Z qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = -0.5; qtmp[1] = 0.; qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("Z/N"); // N-G qtmp = new double [3]; qtmp[0] = qtmp[2] = -0.5; qtmp[1] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("{/Symbol G}"); // G-M qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("M/R"); // R-G qtmp = new double [3]; qtmp[0] = 0.; qtmp[1] = -0.5; qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("{/Symbol G}"); } } else if (spgnum >= 3 && spgnum <= 15){ // Monoclinic if (symbol[0] == 'P'){ // MCL-P const double eta = (1.-lb*cosa/lc)/(2.*(1.-cosa*cosa)); const double niu = 0.5 - eta * lc * cosa / lb; ndstr.push_back("{/Symbol G}"); // G-Y qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("Y"); // Y-H qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.; qtmp[1] = eta; qtmp[2] = 1.-niu; qe.push_back(qtmp); ndstr.push_back("H"); // H-C qtmp = new double [3]; qtmp[0] = 0.; qtmp[1] = eta; qtmp[2] = 1.-niu; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.; qtmp[1] = qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("C"); // C-E qtmp = new double [3]; qtmp[0] = 0.; qtmp[1] = qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("E"); // E-M1 qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = 1.-eta; qtmp[2] = niu; qe.push_back(qtmp); ndstr.push_back("M_1"); // M1-A qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = 1.-eta; qtmp[2] = niu; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("A"); // A-X qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.; qtmp[1] = 0.5; qe.push_back(qtmp); ndstr.push_back("X"); // X-H1 qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.; qtmp[1] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.; qtmp[2] = 1.-eta; qtmp[1] = niu; qe.push_back(qtmp); ndstr.push_back("H_1/M"); // M-D qtmp = new double [3]; qtmp[0] = 0.5; qtmp[2] = eta; qtmp[1] = 1.-niu; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.; qe.push_back(qtmp); ndstr.push_back("D"); // D-Z qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("Z/Y"); // Y-D qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.; qe.push_back(qtmp); ndstr.push_back("D"); } else { // MCL-C if (coskg < 0.){ // MCLC1 const double xi = (2. - lb * cosa / lc) / (4.*(1.-cosa*cosa)); const double eta = 0.5 + 2.*xi*lc/lb*cosa; const double psi = 0.75 - la * la /(4.*lb*lb*(1.-cosa*cosa)); const double phi = psi + (0.75 - psi) * lb / lc * cosa; ndstr.push_back("{/Symbol G}"); // G-Y qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("Y"); // Y-F qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 1.-xi; qtmp[2] = 1.-eta; qe.push_back(qtmp); ndstr.push_back("F"); // F-L qtmp = new double [3]; qtmp[0] = qtmp[1] = 1.-xi; qtmp[2] = 1.-eta; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("L"); // L-I qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = phi; qtmp[1] = 1.-phi; qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("I/I_1"); // I1-Z qtmp = new double [3]; qtmp[0] = 1.-phi; qtmp[1] = -qtmp[0]; qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("Z"); // Z-F1 qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = xi; qtmp[2] = eta; qe.push_back(qtmp); ndstr.push_back("F_1/Y"); // Y-X1 qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = psi; qtmp[1] = 1.-psi; qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("X_1/X"); // X-G qtmp = new double [3]; qtmp[0] = 1.-psi; qtmp[1] = -qtmp[0]; qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("{/Symbol G}"); // G-N qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("N/M"); // M-G qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.5; qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("{/Symbol G}"); } else if ( fabs(coskg) < ZERO) { // MCLC2 const double xi = (2. - lb * cosa / lc) / (4.*(1.-cosa*cosa)); const double eta = 0.5 + 2.*xi*lc/lb*cosa; const double psi = 0.75 - la * la /(4.*lb*lb*(1.-cosa*cosa)); const double phi = psi + (0.75 - psi) * lb / lc * cosa; ndstr.push_back("{/Symbol G}"); // G-Y qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("Y"); // Y-F qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 1.-xi; qtmp[2] = 1.-eta; qe.push_back(qtmp); ndstr.push_back("F"); // F-L qtmp = new double [3]; qtmp[0] = qtmp[1] = 1.-xi; qtmp[2] = 1.-eta; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("L"); // L-I qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = phi; qtmp[1] = 1.-phi; qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("I/I_1"); // I1-Z qtmp = new double [3]; qtmp[0] = 1.-phi; qtmp[1] = -qtmp[0]; qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("Z"); // Z-F1 qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = xi; qtmp[2] = eta; qe.push_back(qtmp); ndstr.push_back("F_1/N"); // N-G qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("{/Symbol G}"); // G-M qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.5; qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("M"); } else { double flag = lb / lc * cosa + lb * lb / (la * la *(1.-cosa*cosa)); if (fabs(flag) < ZERO){ // MCLC4 const double miu = 0.25*(1. + lb * lb / (la *la)); const double del = lb * lc * cosa / (2.*la*la); const double xi = miu - 0.25 + (1. - lb * cosa / lc)/(4.*(1.-cosa*cosa)); const double eta = 0.5 + 2.*xi*lc/lb*cosa; const double phi = 1. + xi - 2.*miu; const double psi = eta - 2.*del; ndstr.push_back("{/Symbol G}"); // G-Y qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = miu; qtmp[2] = del; qe.push_back(qtmp); ndstr.push_back("Y"); // Y-F qtmp = new double [3]; qtmp[0] = qtmp[1] = miu; qtmp[2] = del; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 1.-phi; qtmp[2] = 1.-psi; qe.push_back(qtmp); ndstr.push_back("F"); // F-H qtmp = new double [3]; qtmp[0] = qtmp[1] = 1.-phi; qtmp[2] = 1.-psi; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = xi; qtmp[2] = eta; qe.push_back(qtmp); ndstr.push_back("H"); // H-Z qtmp = new double [3]; qtmp[0] = qtmp[1] = xi; qtmp[2] = eta; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("Z"); // Z-I qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.5; qtmp[1] = -0.5; qe.push_back(qtmp); ndstr.push_back("I/H_1"); // H1-Y1 qtmp = new double [3]; qtmp[0] = 1.-xi; qtmp[1] = -xi; qtmp[2] = 1.-eta; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 1.-miu; qtmp[1] = -miu; qtmp[2] = -del; qe.push_back(qtmp); ndstr.push_back("Y_1"); // Y1-X qtmp = new double [3]; qtmp[0] = 1.-miu; qtmp[1] = -miu; qtmp[2] = -del; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = -0.5; qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("X"); // X-G qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = -0.5; qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("{/Symbol G}"); // G-N qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("N/M"); // M-G qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("{/Symbol G}"); } else if (flag < 1.){ // MCLC3 const double miu = 0.25*(1. + lb * lb / (la *la)); const double del = lb * lc * cosa / (2.*la*la); const double xi = miu - 0.25 + (1. - lb * cosa / lc)/(4.*(1.-cosa*cosa)); const double eta = 0.5 + 2.*xi*lc/lb*cosa; const double phi = 1. + xi - 2.*miu; const double psi = eta - 2.*del; ndstr.push_back("{/Symbol G}"); // G-Y qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = miu; qtmp[2] = del; qe.push_back(qtmp); ndstr.push_back("Y"); // Y-F qtmp = new double [3]; qtmp[0] = qtmp[1] = miu; qtmp[2] = del; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 1.-phi; qtmp[2] = 1.-psi; qe.push_back(qtmp); ndstr.push_back("F"); // F-H qtmp = new double [3]; qtmp[0] = qtmp[1] = 1.-phi; qtmp[2] = 1.-psi; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = xi; qtmp[2] = eta; qe.push_back(qtmp); ndstr.push_back("H"); // H-Z qtmp = new double [3]; qtmp[0] = qtmp[1] = xi; qtmp[2] = eta; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("Z"); // Z-I qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.5; qtmp[1] = -0.5; qe.push_back(qtmp); ndstr.push_back("I/H_1"); // I-F1 qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = phi; qtmp[2] = phi - 1.; qtmp[1] = psi; qe.push_back(qtmp); ndstr.push_back("F_1/H_1"); // H1-Y1 qtmp = new double [3]; qtmp[0] = 1.-xi; qtmp[1] = -xi; qtmp[2] = 1.-eta; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 1.-miu; qtmp[1] = -miu; qtmp[2] = -del; qe.push_back(qtmp); ndstr.push_back("Y_1"); // Y1-X qtmp = new double [3]; qtmp[0] = 1.-miu; qtmp[1] = -miu; qtmp[2] = -del; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = -0.5; qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("X"); // X-G qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = -0.5; qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("{/Symbol G}"); // G-N qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("N/M"); // M-G qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("{/Symbol G}"); } else { // MCLC5 const double xi = (lb*lb/(la*la) + (1.-lb*cosa/lc)/(1.-cosa*cosa))*0.25; const double eta = 0.5 + 2.*xi*lc*cosa/lb; const double miu = 0.5*eta + lb * lb /(4.*la*la) - lb*lc/(2.*la*la)*cosa; const double niu = 2.*miu - xi; const double omg = (4.*niu - 1. - lb*lb*(1.-cosa*cosa)/(la*la))*lc/(2.*lb*cosa); const double del = xi*lc*cosa/lb + omg*0.5 - 0.25; const double rho = 1.-xi*la*la/(lb*lb); ndstr.push_back("{/Symbol G}"); // G-Y qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = miu; qtmp[2] = del; qe.push_back(qtmp); ndstr.push_back("Y"); // Y-F qtmp = new double [3]; qtmp[0] = qtmp[1] = miu; qtmp[2] = del; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = niu; qtmp[2] = omg; qe.push_back(qtmp); ndstr.push_back("F"); // F-L qtmp = new double [3]; qtmp[0] = qtmp[1] = niu; qtmp[2] = omg; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("Y"); // L-I qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = rho; qtmp[1] = 1.-rho; qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("I/I_1"); // I1-Z qtmp = new double [3]; qtmp[0] = 1.-rho; qtmp[1] = -qtmp[0]; qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("Z"); // Z-H qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = xi; qtmp[2] = eta; qe.push_back(qtmp); ndstr.push_back("H"); // H-F1 qtmp = new double [3]; qtmp[0] = qtmp[1] = xi; qtmp[2] = eta; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 1.-niu; qtmp[2] = 1.-omg; qe.push_back(qtmp); ndstr.push_back("F_1/H_1"); // H1-Y1 qtmp = new double [3]; qtmp[0] = 1.-xi; qtmp[1] = -xi; qtmp[2] = 1.-eta; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 1.-miu; qtmp[1] = -miu; qtmp[2] = -del; qe.push_back(qtmp); ndstr.push_back("Y_1"); // Y1-X qtmp = new double [3]; qtmp[0] = 1.-miu; qtmp[1] = -miu; qtmp[2] = -del; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("Y_1"); // X-G qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("{/Symbol G}"); // G-N qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("N/M"); // M-G qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.; qe.push_back(qtmp); ndstr.push_back("{/Symbol G}"); } } } } else if (spgnum >= 16 && spgnum <= 74){ // Orthorhombic if (symbol[0] == 'P'){ // ORC ndstr.push_back("{/Symbol G}"); // G-X qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("X"); // X-S qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("S"); // S-G qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("{/Symbol G}"); // G-Z qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("Z"); // Z-U qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.; qe.push_back(qtmp); ndstr.push_back("U"); // U-R qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("R"); // R-T qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.; qtmp[1] = qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("T"); // T-Z qtmp = new double [3]; qtmp[0] = 0.; qtmp[1] = qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("Z/Y"); // Y-T qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.; qtmp[1] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.; qtmp[1] = qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("T/U"); // U-X qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("X/S"); // S-R qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("R"); } else if (symbol[0] == 'F'){ // ORCF double flag = 1./(la*la) - 1./(lb*lb) - 1./(lc*lc); if (fabs(flag) < ZERO){ // ORCF3 const double xi = 0.25 * (1. + la*la/(lb*lb) - la*la/(lc*lc)); const double eta = 0.25 * (1. + la*la/(lb*lb) + la*la/(lc*lc)); ndstr.push_back("{/Symbol G}"); // G-Y qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.; qe.push_back(qtmp); ndstr.push_back("Y"); // Y-T qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 1.; qtmp[1] = qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("T"); // T-Z qtmp = new double [3]; qtmp[0] = 1.; qtmp[1] = qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("Z"); // Z-G qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("{/Symbol G}"); // G-X qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.; qtmp[1] = qtmp[2] = eta; qe.push_back(qtmp); ndstr.push_back("X"); // X-A1 qtmp = new double [3]; qtmp[0] = 0.; qtmp[1] = qtmp[2] = eta; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = 0.5-xi; qtmp[2] = 1.-xi; qe.push_back(qtmp); ndstr.push_back("A_1"); // A1-Y qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = 0.5-xi; qtmp[2] = 1.-xi; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.; qe.push_back(qtmp); ndstr.push_back("Y/X"); // X-A qtmp = new double [3]; qtmp[0] = 0.; qtmp[1] = qtmp[2] = eta; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = 0.5+xi; qtmp[2] = xi; qe.push_back(qtmp); ndstr.push_back("A"); // A-Z qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = 0.5+xi; qtmp[2] = xi; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("Z/L"); // L-G qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("{/Symbol G}"); } else if (flag > 0.){ // ORCF1 const double xi = 0.25 * (1. + la*la/(lb*lb) - la*la/(lc*lc)); const double eta = 0.25 * (1. + la*la/(lb*lb) + la*la/(lc*lc)); ndstr.push_back("{/Symbol G}"); // G-Y qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.; qe.push_back(qtmp); ndstr.push_back("Y"); // Y-T qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 1.; qtmp[1] = qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("T"); // T-Z qtmp = new double [3]; qtmp[0] = 1.; qtmp[1] = qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("Z"); // Z-G qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("{/Symbol G}"); // G-X qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.; qtmp[1] = qtmp[2] = eta; qe.push_back(qtmp); ndstr.push_back("X"); // X-A1 qtmp = new double [3]; qtmp[0] = 0.; qtmp[1] = qtmp[2] = eta; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = 0.5-xi; qtmp[2] = 1.-xi; qe.push_back(qtmp); ndstr.push_back("A_1"); // A1-Y qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = 0.5-xi; qtmp[2] = 1.-xi; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.; qe.push_back(qtmp); ndstr.push_back("Y/T"); // T-X1 qtmp = new double [3]; qtmp[0] = 1.; qtmp[1] = qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 1.; qtmp[1] = qtmp[2] = 1.-eta; qe.push_back(qtmp); ndstr.push_back("X_1/X"); // X-A qtmp = new double [3]; qtmp[0] = 0.; qtmp[1] = qtmp[2] = eta; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = 0.5+xi; qtmp[2] = xi; qe.push_back(qtmp); ndstr.push_back("A"); // A-Z qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = 0.5+xi; qtmp[2] = xi; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("Z/L"); // L-G qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("{/Symbol G}"); } else { // ORCF2 const double eta = 0.25 * (1. + la*la/(lb*lb) - la*la/(lc*lc)); const double phi = 0.25 * (1. + lc*lc/(lb*lb) - lc*lc/(la*la)); const double del = 0.25 * (1. + lb*lb/(la*la) - lb*lb/(lc*lc)); ndstr.push_back("{/Symbol G}"); // G-Y qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.; qe.push_back(qtmp); ndstr.push_back("Y"); // Y-C qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.5; qtmp[2] = 0.5-eta; qtmp[1] = 1.-eta; qe.push_back(qtmp); ndstr.push_back("C"); // C-D qtmp = new double [3]; qtmp[0] = 0.5; qtmp[2] = 0.5-eta; qtmp[1] = 1.-eta; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.5-del; qtmp[2] = 0.5; qtmp[1] = 1.-del; qe.push_back(qtmp); ndstr.push_back("D"); // D-X qtmp = new double [3]; qtmp[0] = 0.5-del; qtmp[2] = 0.5; qtmp[1] = 1.-del; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.; qtmp[1] = qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("X"); // X-G qtmp = new double [3]; qtmp[0] = 0.; qtmp[1] = qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("{/Symbol G}"); // G-Z qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("Z"); // Z-D1 qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.5+del; qtmp[1] = 0.5; qtmp[2] = del; qe.push_back(qtmp); ndstr.push_back("D_1"); // D1-H qtmp = new double [3]; qtmp[0] = 0.5+del; qtmp[1] = 0.5; qtmp[2] = del; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 1.-phi; qtmp[1] = 0.5-phi; qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("H"); // H-C qtmp = new double [3]; qtmp[0] = 1.-phi; qtmp[1] = 0.5-phi; qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = 0.5-eta; qtmp[2] = 1.-eta; qe.push_back(qtmp); ndstr.push_back("C/C_1"); // C1-Z qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = 0.5+eta; qtmp[2] = eta; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("Z/X"); // X-H1 qtmp = new double [3]; qtmp[0] = 0.; qtmp[1] = qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = phi; qtmp[1] = 0.5+phi; qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("H_1/H"); // H-Y qtmp = new double [3]; qtmp[0] = 1.-phi; qtmp[1] = 0.5-phi; qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.; qe.push_back(qtmp); ndstr.push_back("Y/L"); // L-G qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("{/Symbol G}"); } } else if (symbol[0] == 'C'){ // ORCC const double xi = 0.25 * (1. + la*la/(lb*lb)); ndstr.push_back("{/Symbol G}"); // G-X qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = xi; qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("X"); // X-S qtmp = new double [3]; qtmp[0] = qtmp[1] = xi; qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.; qtmp[1] = 0.5; qe.push_back(qtmp); ndstr.push_back("S"); // S-R qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.; qtmp[1] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.; qtmp[1] = qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("R"); // R-A qtmp = new double [3]; qtmp[0] = 0.; qtmp[1] = qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = xi; qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("A"); // A-Z qtmp = new double [3]; qtmp[0] = qtmp[1] = xi; qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("Z"); // Z-G qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("{/Symbol G}"); // G-Y qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = -0.5; qtmp[1] = 0.5; qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("Y"); // Y-X1 qtmp = new double [3]; qtmp[0] = -0.5; qtmp[1] = 0.5; qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = -xi; qtmp[1] = 1.-xi; qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("X_1"); // X1-A1 qtmp = new double [3]; qtmp[0] = -xi; qtmp[1] = 1.-xi; qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = -xi; qtmp[1] = 1.-xi; qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("A_1"); // A1-T qtmp = new double [3]; qtmp[0] = -xi; qtmp[1] = 1.-xi; qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = -0.5; qtmp[1] = qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("T"); // T-Y qtmp = new double [3]; qtmp[0] = -0.5; qtmp[1] = qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = -0.5; qtmp[1] = 0.5; qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("Y/Z"); // Z-T qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = -0.5; qtmp[1] = qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("T"); } else { // ORCI const double xi = 0.25 * (1. + la*la/(lc*lc)); const double eta = 0.25 * (1. + lb*lb/(lc*lc)); const double del = (lb*lb-la*la)/(4.*lc*lc); const double miu = (lb*lb+la*la)/(4.*lc*lc); ndstr.push_back("{/Symbol G}"); // G-X qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = -xi; qtmp[1] = qtmp[2] = xi; qe.push_back(qtmp); ndstr.push_back("X"); // X-L qtmp = new double [3]; qtmp[0] = -xi; qtmp[1] = qtmp[2] = xi; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = -miu; qtmp[1] = miu; qtmp[2] = 0.5-del; qe.push_back(qtmp); ndstr.push_back("L"); // L-T qtmp = new double [3]; qtmp[0] = -miu; qtmp[1] = miu; qtmp[2] = 0.5-del; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("T"); // T-W qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.25; qe.push_back(qtmp); ndstr.push_back("W"); // W-R qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.25; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.; qtmp[1] = 0.5; qe.push_back(qtmp); ndstr.push_back("R"); // R-X1 qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.; qtmp[1] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = xi; qtmp[1] = 1.-xi; qtmp[2] = -xi; qe.push_back(qtmp); ndstr.push_back("X_1"); // X1-Z qtmp = new double [3]; qtmp[0] = xi; qtmp[1] = 1.-xi; qtmp[2] = -xi; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.5; qtmp[2] = -0.5; qe.push_back(qtmp); ndstr.push_back("Z"); // Z-G qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.5; qtmp[2] = -0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("{/Symbol G}"); // G-Y qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[2] = eta; qtmp[1] = -eta; qe.push_back(qtmp); ndstr.push_back("Y"); // Y-S qtmp = new double [3]; qtmp[0] = qtmp[2] = eta; qtmp[1] = -eta; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("S"); // S-W qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.25; qe.push_back(qtmp); ndstr.push_back("W/L_1"); // L1-Y qtmp = new double [3]; qtmp[0] = miu; qtmp[1] = -miu; qtmp[2] = 0.5+del; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[2] = eta; qtmp[1] = -eta; qe.push_back(qtmp); ndstr.push_back("Y/Y_1"); // Y1-Z qtmp = new double [3]; qtmp[0] = 1.-eta; qtmp[1] = eta; qtmp[2] = -eta; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.5; qtmp[2] = -0.5; qe.push_back(qtmp); ndstr.push_back("Z"); } } else if (spgnum >= 75 && spgnum <= 142){ // Tetragonal if (symbol[0] == 'P'){ // TET ndstr.push_back("{/Symbol G}"); // G-X qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.; qtmp[1] = 0.5; qe.push_back(qtmp); ndstr.push_back("X"); // X-M qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.; qtmp[1] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("M"); // M-G qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("{/Symbol G}"); // G-Z qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("Z"); // Z-R qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.; qtmp[1] = qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("R"); // R-A qtmp = new double [3]; qtmp[0] = 0.; qtmp[1] = qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("A"); // A-Z qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("Z/X"); // X-R qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.; qtmp[1] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.; qtmp[1] = qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("R/M"); // M-A qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("A"); } else { // BCT if (la > lc){ // BCT1 const double eta = 0.25 * (1. + lc*lc/(la*la)); ndstr.push_back("{/Symbol G}"); // G-X qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("X"); // X-M qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = -0.5; qtmp[1] = qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("M"); // M-G qtmp = new double [3]; qtmp[0] = -0.5; qtmp[1] = qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("{/Symbol G}"); // G-Z qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = eta; qtmp[2] = -eta; qe.push_back(qtmp); ndstr.push_back("Z"); // Z-P qtmp = new double [3]; qtmp[0] = qtmp[1] = eta; qtmp[2] = -eta; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.25; qe.push_back(qtmp); ndstr.push_back("P"); // P-N qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.25; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.; qtmp[1] = 0.5; qe.push_back(qtmp); ndstr.push_back("N"); // N-Z1 qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.; qtmp[1] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = -eta; qtmp[1] = 1.-eta; qtmp[2] = eta; qe.push_back(qtmp); ndstr.push_back("Z_1"); // Z1-M qtmp = new double [3]; qtmp[0] = -eta; qtmp[1] = 1.-eta; qtmp[2] = eta; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = -0.5; qtmp[1] = qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("M"); // X-P qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.25; qe.push_back(qtmp); ndstr.push_back("P"); } else { // BCT2 const double eta = 0.25 * (1. + la*la/(lc*lc)); const double xi = la*la/(2.*lc*lc); ndstr.push_back("{/Symbol G}"); // G-X qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("X"); // X-Y qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = -xi; qtmp[1] = xi; qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("Y"); // Y-Sigma qtmp = new double [3]; qtmp[0] = -xi; qtmp[1] = xi; qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = -eta; qtmp[1] = qtmp[2] = eta; qe.push_back(qtmp); ndstr.push_back("{/Symbol S}"); // Sigma-G qtmp = new double [3]; qtmp[0] = -eta; qtmp[1] = qtmp[2] = eta; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("{/Symbol G}"); // G-Z qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.5; qtmp[2] = -0.5; qe.push_back(qtmp); ndstr.push_back("Z"); // Z-Sigma_1 qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.5; qtmp[2] = -0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = eta; qtmp[1] = 1.-eta; qtmp[2] = -eta; qe.push_back(qtmp); ndstr.push_back("{/Symbol S}_1"); // Sigma_1-N qtmp = new double [3]; qtmp[0] = eta; qtmp[1] = 1.-eta; qtmp[2] = -eta; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.; qtmp[1] = 0.5; qe.push_back(qtmp); ndstr.push_back("N"); // N-P qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.; qtmp[1] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[2] = qtmp[1] = 0.25; qe.push_back(qtmp); ndstr.push_back("P"); // P-Y1 qtmp = new double [3]; qtmp[0] = qtmp[2] = qtmp[1] = 0.25; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.5; qtmp[1] = -xi; qe.push_back(qtmp); ndstr.push_back("Y_1"); // Y1-Z qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.5; qtmp[1] = -xi; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.5; qtmp[1] = -0.5; qe.push_back(qtmp); ndstr.push_back("Z"); // X-P qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[2] = qtmp[1] = 0.25; qe.push_back(qtmp); ndstr.push_back("P"); } } } else if (spgnum >= 143 && spgnum <= 167){ // Trigonal if (cosg > 0.){ // RHL1 const double eta = (1.+4.*cosa)/(2.+4.*cosa); const double niu = 0.75 - 0.5*eta; ndstr.push_back("{/Symbol G}"); // G-L qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("L"); // L-B1 qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = 1.-eta; qtmp[2] = eta - 1.; qe.push_back(qtmp); ndstr.push_back("B_1/B"); // B-Z qtmp = new double [3]; qtmp[0] = eta; qtmp[1] = 0.5; qtmp[2] = 1.-eta; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("Z"); // Z-G qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("{/Symbol G}"); // G-X qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = niu; qtmp[1] = 0.; qtmp[2] = -niu; qe.push_back(qtmp); ndstr.push_back("X/Q"); // Q-F qtmp = new double [3]; qtmp[0] = 1.-niu; qtmp[1] = niu; qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("F"); // F-P1 qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 1.-niu; qtmp[2] = 1.-eta; qe.push_back(qtmp); ndstr.push_back("P_1"); // P1-Z qtmp = new double [3]; qtmp[0] = qtmp[1] = 1.-niu; qtmp[2] = 1.-eta; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("Z/L"); // L-P qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = eta; qtmp[1] = qtmp[2] = niu; qe.push_back(qtmp); ndstr.push_back("P"); } else { // RHL2 const double eta = 0.5*(1.+cosa)/(1.-cosa); const double niu = 0.75 - 0.5*eta; ndstr.push_back("{/Symbol G}"); // G-P qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[2] = 1.-niu; qtmp[1] = -niu; qe.push_back(qtmp); ndstr.push_back("P"); // P-Z qtmp = new double [3]; qtmp[0] = qtmp[2] = 1.-niu; qtmp[1] = -niu; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.5; qtmp[1] = -0.5; qe.push_back(qtmp); ndstr.push_back("Z"); // Z-Q qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.5; qtmp[1] = -0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = eta; qe.push_back(qtmp); ndstr.push_back("Q"); // Q-G qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = eta; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("{/Symbol G}"); // G-F qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = -0.5; qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("F"); // F-P1 qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = -0.5; qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = niu; qtmp[1] = qtmp[2] = niu - 1.; qe.push_back(qtmp); ndstr.push_back("P_1"); // P1-Q1 qtmp = new double [3]; qtmp[0] = niu; qtmp[1] = qtmp[2] = niu - 1.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 1.-eta; qtmp[1] = qtmp[2] = -eta; qe.push_back(qtmp); ndstr.push_back("Q_1"); // Q1-L qtmp = new double [3]; qtmp[0] = 1.-eta; qtmp[1] = qtmp[2] = -eta; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("L"); // L-Z qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.5; qtmp[1] = -0.5; qe.push_back(qtmp); ndstr.push_back("Z"); } } else if (spgnum >= 168 && spgnum <= 194){ // Hexagonal ndstr.push_back("{/Symbol G}"); // G-M qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("M"); // M-K qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 1./3.; qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("K"); // K-G qtmp = new double [3]; qtmp[0] = qtmp[1] = 1./3.; qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("{/Symbol G}"); // G-A qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("A"); // A-L qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.; qe.push_back(qtmp); ndstr.push_back("L"); // L-H qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 1./3.; qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("H"); // H-A qtmp = new double [3]; qtmp[0] = qtmp[1] = 1./3.; qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("A/L"); // L-M qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("M/K"); // K-H qtmp = new double [3]; qtmp[0] = qtmp[1] = 1./3.; qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 1./3.; qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("H"); } else if (spgnum >= 195 && spgnum <= 230){ // Cubic if (symbol[0] == 'P'){ // CUB ndstr.push_back("{/Symbol G}"); // G-X qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.; qtmp[1] = 0.5; qe.push_back(qtmp); ndstr.push_back("X"); // X-M qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.; qtmp[1] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("M"); // M-G qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("{/Symbol G}"); // G-R qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("R"); // R-X qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.; qtmp[1] = 0.5; qe.push_back(qtmp); ndstr.push_back("X/M"); // M-R qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.5; qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("R"); } else if (symbol[0] == 'F'){ // FCC ndstr.push_back("{/Symbol G}"); // G-X qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.; qe.push_back(qtmp); ndstr.push_back("X"); // X-W qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.5; qtmp[1] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = 0.25; qtmp[2] = 0.75; qe.push_back(qtmp); ndstr.push_back("W"); // W-K qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = 0.25; qtmp[2] = 0.75; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.375; qtmp[2] = 0.75; qe.push_back(qtmp); ndstr.push_back("K"); // K-G qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.375; qtmp[2] = 0.75; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("{/Symbol G}"); // G-L qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("L"); // L-U qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.625; qtmp[1] = 0.25; qtmp[2] = 0.625; qe.push_back(qtmp); ndstr.push_back("U"); // U-W qtmp = new double [3]; qtmp[0] = 0.625; qtmp[1] = 0.25; qtmp[2] = 0.625; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = 0.25; qtmp[2] = 0.75; qe.push_back(qtmp); ndstr.push_back("W"); // W-L qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = 0.25; qtmp[2] = 0.75; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("L"); // L-K qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.375; qtmp[2] = 0.75; qe.push_back(qtmp); ndstr.push_back("K/U"); // U-X qtmp = new double [3]; qtmp[0] = 0.625; qtmp[1] = 0.25; qtmp[2] = 0.625; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = 0.5; qtmp[1] = 0.; qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("X"); } else { // BCC ndstr.push_back("{/Symbol G}"); // G-H qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.5; qtmp[1] = -0.5; qe.push_back(qtmp); ndstr.push_back("H"); // H-N qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.5; qtmp[1] = -0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("N"); // N-G qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qe.push_back(qtmp); ndstr.push_back("{/Symbol G}"); // G-P qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.25; qe.push_back(qtmp); ndstr.push_back("P"); // P-H qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.25; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[2] = 0.5; qtmp[1] = -0.5; qe.push_back(qtmp); ndstr.push_back("H/P"); // P-N qtmp = new double [3]; qtmp[0] = qtmp[1] = qtmp[2] = 0.25; qs.push_back(qtmp); qtmp = new double [3]; qtmp[0] = qtmp[1] = 0.; qtmp[2] = 0.5; qe.push_back(qtmp); ndstr.push_back("N"); } } else { printf("\nSorry, failed to identify the crystal system, please use the manual mode.\n"); } // to determine the number of points along each line, with a step size of 0.05 const double qs_inv = 1./QSTEP; int nbin = qs.size(); if (nbin > 0) printf("\nPhonon dispersion will be evaluated along lines:\n\t%s", ndstr[0].c_str()); for (int is = 0; is < nbin; ++is){ double *qstr = qs[is]; double *qend = qe[is]; double ql = 0.; for (int i = 0; i < 3; ++i) ql += (qend[i] - qstr[i])*(qend[i] - qstr[i]); int nqpt = MAX(int(sqrt(ql) * qs_inv + 0.5), 2); nqbin.push_back(nqpt); printf("-%s", ndstr[is+1].c_str()); } if (nbin > 0) printf("\n"); } #endif FILE *fp = fopen(fname, "w"); fprintf(fp,"# q qr freq\n"); fprintf(fp,"# 2pi/L 2pi/L %s\n", dynmat->funit); double qr = 0., dq, q[3], qinc[3]; int nbin = qs.size(); for (int is = 0; is < nbin; ++is){ double *qstr = qs[is]; double *qend = qe[is]; int nbin = nqbin[is]; for (int i = 0; i < 3; ++i) qinc[i] = (qend[i]-qstr[i])/double(nbin-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 < nbin; ++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; delete []qstr; delete []qend; } qs.clear(); qe.clear(); 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){ const char qmk = char(34); // " fp = fopen("pdisp.gnuplot", "w"); fprintf(fp,"set term post enha colo 20\nset out %cpdisp.eps%c\n\n", qmk, qmk); fprintf(fp,"set xlabel %cq%c\n", qmk, qmk); fprintf(fp,"set ylabel %cfrequency (THz)%c\n\n", qmk, qmk); 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%s%c %lg, ", qmk, ndstr[i].c_str(), qmk, nodes[i]); fprintf(fp, "%c%s%c %lg)\n\n", qmk, ndstr[nnd-1].c_str(), qmk, nodes[nnd-1]); fprintf(fp, "unset key\n\n"); fprintf(fp, "plot %c%s%c u 4:5 w l lt 1", qmk, fname, qmk); for (int i = 1; i < ndim; ++i) fprintf(fp,",\\\n%c%c u 4:%d w l lt 1", qmk, qmk, i+5); fclose(fp); printf("\nPhonon dispersion data are written to: %s, you can visualize the results\n", fname); printf("by invoking: `gnuplot pdisp.gnuplot; gv pdisp.eps`\n"); } - for (int ii = 0; ii < 80; ++ii) printf("="); printf("\n"); + for (int ii = 0; ii < 80; ++ii) printf("="); + printf("\n"); delete []fname; nodes.clear(); ndstr.clear(); return; } diff --git a/tools/phonon/dynmat.cpp b/tools/phonon/dynmat.cpp index e82f47313..3b7bfe826 100644 --- a/tools/phonon/dynmat.cpp +++ b/tools/phonon/dynmat.cpp @@ -1,589 +1,597 @@ #include "dynmat.h" #include "math.h" #include "version.h" #include "global.h" +extern "C" void zheevd_(char *, char *, long int *, doublecomplex *, + long int *, double *, doublecomplex *, + long int *, double *, long int *, long int *, + long int *, long int *); + // to initialize 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("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"); + 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[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(); 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 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 ){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);} + if ( (int) fread(&Tmeasure, sizeof(double), 1, fp) != 1 ){printf("\nError while reading temperature from file: %s\n", binfile); fclose(fp); exit(3);} + if ( (int) fread(&basevec[0], sizeof(double), 9, fp) != 9 ){printf("\nError while reading lattice info from file: %s\n", binfile); fclose(fp); exit(3);} + if ( (int) fread(basis[0], sizeof(double), fftdim, fp) != fftdim){printf("\nError while reading basis info from file: %s\n", binfile); fclose(fp); exit(3);} + if ( (int) fread(&attyp[0], sizeof(int), nucell, fp) != nucell){printf("\nError while reading atom types from file: %s\n", binfile); fclose(fp); exit(3);} + if ( (int) 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); car2dir(); real2rec(); // initialize interpolation interpolate = new Interpolate(nx,ny,nz,fftdim2,DM_all); if (flag_reset_gamma) interpolate->reset_gamma(); // 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 ){ 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); 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); 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; + long int n, lda, lwork, lrwork, *iwork, liwork, info; doublecomplex *work; - doublereal *w = &egv[0], *rwork; + double *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; 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){ 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() { double mat[9]; for (int idim = 0; idim < 9; ++idim) mat[idim] = basevec[idim]; GaussJordan(3, mat); 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++) 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("="); // 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]; geteigen(egvs, 0); printf("\nEigenvalues of Phi at gamma before enforcing ASR:\n"); 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){ - for (int i=0; i<80; i++) printf("="); printf("\n"); + for (int i=0; i<80; i++) printf("="); + printf("\n"); return; } 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){ double sum = 0.; 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){ 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){ double csum = 0.; 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){ double sum = 0.; 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){ 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]; geteigen(egvs, 0); printf("Eigenvalues of Phi at gamma after enforcing ASR:\n"); 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"); + 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]; vol = 8.*atan(1.)/vol; for (int i = 0; i < 9; ++i) ibasevec[i] *= vol; 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){ printf("\n A%d: ", i+1); 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){ printf("\n B%d: ", i+1); 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"); 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 i,icol=0,irow=0,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){ big = 0.; for (j = 0; j < n; ++j){ if (ipiv[j] != 1){ for (k = 0; k < n; ++k){ if (ipiv[k] == 0){ idr = j * n + k; nmjk = abs(Mat[idr]); if (nmjk >= big){ big = nmjk; irow = j; icol = k; } } 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){ 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; 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){ if (ll != 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 = 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; 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.%02d, compiled on %s.\n", VERSION, __DATE__); return; } /* --------------------------------------------------------------------*/ diff --git a/tools/phonon/dynmat.h b/tools/phonon/dynmat.h index 1d6e71658..f5bd4010b 100644 --- a/tools/phonon/dynmat.h +++ b/tools/phonon/dynmat.h @@ -1,66 +1,61 @@ #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(); // 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 8f8946dc4..35514c03f 100644 --- a/tools/phonon/green.cpp +++ b/tools/phonon/green.cpp @@ -1,249 +1,248 @@ #include #include #include #include #include "green.h" #include #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(); 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); 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){ beta[idim][0] = 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){ double sum_a = 0.; for (int j = 0; j < ndim; ++j){ double sumHv = 0.; 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]; 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]; double sum_b = 0.; 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; } } 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){ alpha_inf[idim] = beta_inf[idim] = 0.; 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 Z, z_m_a, r_x, rec_x, rec_x_inv; double sr, si; double w = wmin; for (int i = 0; i < nw; ++i){ double a = w*w, ax, bx; Z = std::complex(w*w, epson); 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 (sr, si); 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 Z, rec_x, rec_x_inv; - std::complex cunit = std::complex(0.,1.); double w = wmin; for (int i = 0; i < nw; ++i){ Z = std::complex(w*w, epson); for (int idim = 0; idim < sysdim; ++idim){ rec_x = std::complex(0.,0.); 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 8c0cbde1c..954062d41 100644 --- a/tools/phonon/interpolate.cpp +++ b/tools/phonon/interpolate.cpp @@ -1,309 +1,427 @@ #include "interpolate.h" -#include "math.h" +#include #include "global.h" +/////////////////////// +// tricubic library code +static int A[64][64] = { +{ 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, +{ 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, +{-3, 3, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, +{ 2,-2, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, +{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, +{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, +{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, +{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, +{-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, +{ 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, +{ 9,-9,-9, 9, 0, 0, 0, 0, 6, 3,-6,-3, 0, 0, 0, 0, 6,-6, 3,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 2, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, +{-6, 6, 6,-6, 0, 0, 0, 0,-3,-3, 3, 3, 0, 0, 0, 0,-4, 4,-2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-2,-1,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, +{ 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, +{ 0, 0, 0, 0, 0, 0, 0, 0, 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, +{-6, 6, 6,-6, 0, 0, 0, 0,-4,-2, 4, 2, 0, 0, 0, 0,-3, 3,-3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-1,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, +{ 4,-4,-4, 4, 0, 0, 0, 0, 2, 2,-2,-2, 0, 0, 0, 0, 2,-2, 2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, +{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, +{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, +{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, +{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, +{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, +{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0}, +{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0}, +{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0}, +{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, +{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0, 0, 0, 0, 0}, +{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9,-9,-9, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 3,-6,-3, 0, 0, 0, 0, 6,-6, 3,-3, 0, 0, 0, 0, 4, 2, 2, 1, 0, 0, 0, 0}, +{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-6, 6, 6,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3,-3, 3, 3, 0, 0, 0, 0,-4, 4,-2, 2, 0, 0, 0, 0,-2,-2,-1,-1, 0, 0, 0, 0}, +{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, +{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0}, +{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-6, 6, 6,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4,-2, 4, 2, 0, 0, 0, 0,-3, 3,-3, 3, 0, 0, 0, 0,-2,-1,-2,-1, 0, 0, 0, 0}, +{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4,-4,-4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2,-2,-2, 0, 0, 0, 0, 2,-2, 2,-2, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0}, +{-3, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, +{ 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, +{ 9,-9, 0, 0,-9, 9, 0, 0, 6, 3, 0, 0,-6,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6,-6, 0, 0, 3,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 2, 0, 0, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, +{-6, 6, 0, 0, 6,-6, 0, 0,-3,-3, 0, 0, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4, 4, 0, 0,-2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-2, 0, 0,-1,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, +{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, +{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0, 0, 0,-1, 0, 0, 0}, +{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9,-9, 0, 0,-9, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 3, 0, 0,-6,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6,-6, 0, 0, 3,-3, 0, 0, 4, 2, 0, 0, 2, 1, 0, 0}, +{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-6, 6, 0, 0, 6,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3,-3, 0, 0, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4, 4, 0, 0,-2, 2, 0, 0,-2,-2, 0, 0,-1,-1, 0, 0}, +{ 9, 0,-9, 0,-9, 0, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 3, 0,-6, 0,-3, 0, 6, 0,-6, 0, 3, 0,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0, 2, 0, 2, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0}, +{ 0, 0, 0, 0, 0, 0, 0, 0, 9, 0,-9, 0,-9, 0, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 3, 0,-6, 0,-3, 0, 6, 0,-6, 0, 3, 0,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0, 2, 0, 2, 0, 1, 0}, +{-27,27,27,-27,27,-27,-27,27,-18,-9,18, 9,18, 9,-18,-9,-18,18,-9, 9,18,-18, 9,-9,-18,18,18,-18,-9, 9, 9,-9,-12,-6,-6,-3,12, 6, 6, 3,-12,-6,12, 6,-6,-3, 6, 3,-12,12,-6, 6,-6, 6,-3, 3,-8,-4,-4,-2,-4,-2,-2,-1}, +{18,-18,-18,18,-18,18,18,-18, 9, 9,-9,-9,-9,-9, 9, 9,12,-12, 6,-6,-12,12,-6, 6,12,-12,-12,12, 6,-6,-6, 6, 6, 6, 3, 3,-6,-6,-3,-3, 6, 6,-6,-6, 3, 3,-3,-3, 8,-8, 4,-4, 4,-4, 2,-2, 4, 4, 2, 2, 2, 2, 1, 1}, +{-6, 0, 6, 0, 6, 0,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0,-3, 0, 3, 0, 3, 0,-4, 0, 4, 0,-2, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-2, 0,-1, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0}, +{ 0, 0, 0, 0, 0, 0, 0, 0,-6, 0, 6, 0, 6, 0,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0,-3, 0, 3, 0, 3, 0,-4, 0, 4, 0,-2, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-2, 0,-1, 0,-1, 0}, +{18,-18,-18,18,-18,18,18,-18,12, 6,-12,-6,-12,-6,12, 6, 9,-9, 9,-9,-9, 9,-9, 9,12,-12,-12,12, 6,-6,-6, 6, 6, 3, 6, 3,-6,-3,-6,-3, 8, 4,-8,-4, 4, 2,-4,-2, 6,-6, 6,-6, 3,-3, 3,-3, 4, 2, 4, 2, 2, 1, 2, 1}, +{-12,12,12,-12,12,-12,-12,12,-6,-6, 6, 6, 6, 6,-6,-6,-6, 6,-6, 6, 6,-6, 6,-6,-8, 8, 8,-8,-4, 4, 4,-4,-3,-3,-3,-3, 3, 3, 3, 3,-4,-4, 4, 4,-2,-2, 2, 2,-4, 4,-4, 4,-2, 2,-2, 2,-2,-2,-2,-2,-1,-1,-1,-1}, +{ 2, 0, 0, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, +{ 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, +{-6, 6, 0, 0, 6,-6, 0, 0,-4,-2, 0, 0, 4, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0,-3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, +{ 4,-4, 0, 0,-4, 4, 0, 0, 2, 2, 0, 0,-2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, +{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, +{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0}, +{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-6, 6, 0, 0, 6,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4,-2, 0, 0, 4, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0,-3, 3, 0, 0,-2,-1, 0, 0,-2,-1, 0, 0}, +{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4,-4, 0, 0,-4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 0, 0,-2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 2,-2, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0}, +{-6, 0, 6, 0, 6, 0,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4, 0,-2, 0, 4, 0, 2, 0,-3, 0, 3, 0,-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0,-2, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0}, +{ 0, 0, 0, 0, 0, 0, 0, 0,-6, 0, 6, 0, 6, 0,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4, 0,-2, 0, 4, 0, 2, 0,-3, 0, 3, 0,-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0,-2, 0,-1, 0}, +{18,-18,-18,18,-18,18,18,-18,12, 6,-12,-6,-12,-6,12, 6,12,-12, 6,-6,-12,12,-6, 6, 9,-9,-9, 9, 9,-9,-9, 9, 8, 4, 4, 2,-8,-4,-4,-2, 6, 3,-6,-3, 6, 3,-6,-3, 6,-6, 3,-3, 6,-6, 3,-3, 4, 2, 2, 1, 4, 2, 2, 1}, +{-12,12,12,-12,12,-12,-12,12,-6,-6, 6, 6, 6, 6,-6,-6,-8, 8,-4, 4, 8,-8, 4,-4,-6, 6, 6,-6,-6, 6, 6,-6,-4,-4,-2,-2, 4, 4, 2, 2,-3,-3, 3, 3,-3,-3, 3, 3,-4, 4,-2, 2,-4, 4,-2, 2,-2,-2,-1,-1,-2,-2,-1,-1}, +{ 4, 0,-4, 0,-4, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 2, 0,-2, 0,-2, 0, 2, 0,-2, 0, 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0}, +{ 0, 0, 0, 0, 0, 0, 0, 0, 4, 0,-4, 0,-4, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 2, 0,-2, 0,-2, 0, 2, 0,-2, 0, 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0}, +{-12,12,12,-12,12,-12,-12,12,-8,-4, 8, 4, 8, 4,-8,-4,-6, 6,-6, 6, 6,-6, 6,-6,-6, 6, 6,-6,-6, 6, 6,-6,-4,-2,-4,-2, 4, 2, 4, 2,-4,-2, 4, 2,-4,-2, 4, 2,-3, 3,-3, 3,-3, 3,-3, 3,-2,-1,-2,-1,-2,-1,-2,-1}, +{ 8,-8,-8, 8,-8, 8, 8,-8, 4, 4,-4,-4,-4,-4, 4, 4, 4,-4, 4,-4,-4, 4,-4, 4, 4,-4,-4, 4, 4,-4,-4, 4, 2, 2, 2, 2,-2,-2,-2,-2, 2, 2,-2,-2, 2, 2,-2,-2, 2,-2, 2,-2, 2,-2, 2,-2, 1, 1, 1, 1, 1, 1, 1, 1}}; + +static int ijk2n(int i, int j, int k) { + return(i+4*j+16*k); +} + +/* ---------------------------------------------------------------------------- */ + +static void tricubic_get_coeff_stacked(double a[64], double x[64]) { + int i,j; + for (i=0;i<64;i++) { + a[i]=(double)(0.0); + for (j=0;j<64;j++) { + a[i]+=A[i][j]*x[j]; + } + } +} + +static void tricubic_get_coeff(double a[64], double f[8], double dfdx[8], double dfdy[8], double dfdz[8], double d2fdxdy[8], double d2fdxdz[8], double d2fdydz[8], double d3fdxdydz[8]) { + int i; + double x[64]; + for (i=0;i<8;i++) { + x[0+i]=f[i]; + x[8+i]=dfdx[i]; + x[16+i]=dfdy[i]; + x[24+i]=dfdz[i]; + x[32+i]=d2fdxdy[i]; + x[40+i]=d2fdxdz[i]; + x[48+i]=d2fdydz[i]; + x[56+i]=d3fdxdydz[i]; + } + tricubic_get_coeff_stacked(a,x); +} + +static double tricubic_eval(double a[64], double x, double y, double z) { + int i,j,k; + double ret=(double)(0.0); + /* TRICUBIC EVAL + This is the short version of tricubic_eval. It is used to compute + the value of the function at a given point (x,y,z). To compute + partial derivatives of f, use the full version with the extra args. + */ + for (i=0;i<4;i++) { + for (j=0;j<4;j++) { + for (k=0;k<4;k++) { + ret+=a[ijk2n(i,j,k)]*pow(x,i)*pow(y,j)*pow(z,k); + } + } + } + return(ret); +} + /* ---------------------------------------------------------------------------- * 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(); 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){ 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){ 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; idimdestroy(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){ 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){ 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){ 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){ 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; 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){ DMq[idim].r = 0.; DMq[idim].i = 0.; 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]: "); 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); - for(int i=0; i<80; i++) printf("="); printf("\n\n"); + 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 -extern "C"{ -#include "f2c.h" -#include "clapack.h" -} + +extern "C" typedef struct { double r, i; } doublecomplex; using namespace std; class Interpolate{ public: Interpolate(int, int, int, int, doublecomplex **); ~Interpolate(); void set_method(); void execute(double *, doublecomplex *); void reset_gamma(); int UseGamma; private: void tricubic_init(); void tricubic(double *, doublecomplex *); void trilinear(double *, doublecomplex *); Memory *memory; int which; int Nx, Ny, Nz, Npt, ndim; int flag_reset_gamma, flag_allocated_dfs; doublecomplex **data; doublecomplex **Dfdx, **Dfdy, **Dfdz, **D2fdxdy, **D2fdxdz, **D2fdydz, **D3fdxdydz; double a[64], f[8], dfdx[8], dfdy[8], dfdz[8], d2fdxdy[8], d2fdxdz[8], d2fdydz[8], d3fdxdydz[8]; int vidx[8]; }; #endif diff --git a/tools/phonon/phonon.cpp b/tools/phonon/phonon.cpp index 43bea111b..065885cf3 100644 --- a/tools/phonon/phonon.cpp +++ b/tools/phonon/phonon.cpp @@ -1,1107 +1,1111 @@ #include #include "string.h" #include "phonon.h" #include "green.h" #include "timer.h" #include "global.h" #ifdef UseSPG extern "C"{ #include "spglib.h" } #endif /* ---------------------------------------------------------------------------- * 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"); + 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]: "); 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"); + 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]); } // 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); memory->create(dos, ndos, "pdos:dos"); for (int i=0; i 0.){ for (int j = 0; j < ndim; ++j){ int idx = int((eigs[iq][j]-offset)*rdf); if (idx>=0 && idx 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; ieml2f*tpi; scale *= scale; 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; 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]);} 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); 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; 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")); 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); 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 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 ){ 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 ){ 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){ 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){ int ipos = j * sysdim; double sum = 0.; fprintf(fp,"%d", j+1); 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]; 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){ dynmat->getDMq(q); dynmat->writeDMq(q, qr, fp); 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; 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){ table[jj] = exp(-double(jj*jj)*fac); fnorm += table[jj]; } fnorm = 1./fnorm; for (int i = 0; i < npt; ++i){ tmp[i] = 0.; 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; 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){ double Utmp = 0., Stmp = 0., Ftmp = 0., Ztmp = 0., Ctmp = 0.; 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; 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 ){ 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 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 (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]: "); 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); #endif memory->destroy(wt); memory->destroy(qpts); #ifdef UseSPG if (method == 1){ #endif nq = nx*ny*nz; double w = 1./double(nq); 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){ 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 } 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); 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]; // 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); // 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); 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){ int iq = map[i]; if (iq == i) iq2idx[iq] = numq++; } for (int iq = 0; iq < nq; ++iq) wt[iq] = 0.; numq = 0; 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; } #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); 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]); 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); 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 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){ 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){ int hit = int((egval[idim] - offset)*rdf); if (hit >= 0 && hit 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"); + 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]; 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; } if (ldos){ 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]; 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; } } 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); memory->create(eigs, nq,ndim,"QMesh_eigs"); 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; memory->create(copy, n, "count_words:copy"); strcpy(copy,line); char *ptr; - if (ptr = strchr(copy,'#')) *ptr = '\0'; + 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/version.h b/tools/phonon/version.h index 8ed0e80aa..decab631b 100644 --- a/tools/phonon/version.h +++ b/tools/phonon/version.h @@ -1 +1 @@ -#define VERSION 7 +#define VERSION 8