diff --git a/tools/phonon/Makefile b/tools/phonon/Makefile new file mode 100644 index 0000000000..6c39f55cea --- /dev/null +++ b/tools/phonon/Makefile @@ -0,0 +1,64 @@ +.SUFFIXES : .o .cpp +# compiler and flags +CC = g++ -Wno-unused-result +LINK = $(CC) -static +CFLAGS = -O3 $(DEBUG) $(UFLAG) +# +OFLAGS = -O3 $(DEBUG) +INC = $(LPKINC) $(TCINC) $(SPGINC) +LIB = $(LPKLIB) $(TCLIB) $(SPGLIB) +# +# cLapack library needed +LPKINC = -I/opt/clapack/3.2.1/include +LPKLIB = -L/opt/clapack/3.2.1/lib -lclapack -lblas -lf2c #-lm +# +# Tricubic library needed +TCINC = -I/opt/tricubic/1.0/include +TCLIB = -L/opt/tricubic/1.0/lib -ltricubic +# +# spglib 0.7.1, used to get the irreducible q-points +# if UFLAG is not set, spglib won't be used. +UFLAG = -DUseSPG +SPGINC = -I/opt/spglib/0.7.1/include +SPGLIB = -L/opt/spglib/0.7.1/lib -lsymspg +# if spglib > 0.7.1 is used, please +# 1) modify file phonon.cpp, instruction can be found by searching 0.7.1 +# 2) uncomment the following two lines +#SPGINC = -I/opt/spglib/1.1.2/include +#SPGLIB = -L/opt/spglib/1.1.2/lib -lsymspg + +# Debug flags +#DEBUG = -g -DDEBUG +#==================================================================== +ROOT = phana +# executable name +EXE = $(ROOT) +#==================================================================== +# source and rules +SRC = $(wildcard *.cpp) +OBJ = $(SRC:.cpp=.o) + +#==================================================================== +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 `svn info|grep '^Revision'|cut -d: -f2`" > version.h + +#==================================================================== +.f.o: + $(FC) $(FFLAGS) $(FREE) $(MPI) ${INC} -c $< +.f90.o: + $(FC) $(FFLAGS) $(FREE) $(MPI) ${INC} -c $< +.c.o: + $(CC) $(CFLAGS) -c $< +.cpp.o: + $(CC) $(CFLAGS) $(INC) -c $< diff --git a/tools/phonon/README b/tools/phonon/README new file mode 100644 index 0000000000..042cc13c7d --- /dev/null +++ b/tools/phonon/README @@ -0,0 +1,31 @@ +Phonon post-processing tool + +This program reads the binary file created by fix_phonon +and helps to analyse the phonon related info. + +The clapack library is needed to solve the eigen problems, +which could be downloaded from: +http://www.netlib.org/clapack/ + +The tricubic library is also needed to to tricubic interpolations, +which could be obtained from: +http://orca.princeton.edu/francois/software/tricubic/ + +The spglib (version 0.7.1) is optionally needed, enabling one to +evaluate the phonon density of states or vibrational thermal +properties using only the irreducible q-points in the first +Brillouin zone. + +To compile the code, one needs therefore to install the above +libraries and set the paths correctly in the Makefile. + +The units of the output frequencies by this code is THz for +LAMMPS units "real", "si", "metal", and "cgs"; in these cases, +the frequencies are $\nu$ instead of $\omega$. + +One is encouraged to visit http://code.google.com/p/fix-phonon/ +to check out the latest revision on fix-phonon and the post-processing +code. + +Author: Ling-Ti Kong +Feb 2013 diff --git a/tools/phonon/dynmat.cpp b/tools/phonon/dynmat.cpp new file mode 100644 index 0000000000..1886901647 --- /dev/null +++ b/tools/phonon/dynmat.cpp @@ -0,0 +1,602 @@ +#include "dynmat.h" +#include "math.h" +#include "version.h" + +#define MAXLINE 256 + +// to intialize the class +DynMat::DynMat(int narg, char **arg) +{ + attyp = NULL; + memory = NULL; + M_inv_sqrt = NULL; + interpolate = NULL; + DM_q = DM_all = NULL; + binfile = funit = dmfile = NULL; + + 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"); + if (sysdim<1||sysdim>3||nx<1||ny<1||nz<1||nucell<1){ + printf("Wrong values read from header of file: %s, please check the binary file!\n", binfile); + fclose(fp); exit(3); + } + + funit = new char[4]; + strcpy(funit, "THz"); + if (boltz == 1.){eml2f = 1.; delete funit; funit=new char[22]; strcpy(funit,"sqrt(epsilon/(m.sigma^2))");} + else if (boltz == 0.0019872067) eml2f = 3.256576161; + else if (boltz == 8.617343e-5) eml2f = 15.63312493; + else if (boltz == 1.3806504e-23) eml2f = 1.; + else if (boltz == 1.3806504e-16) eml2f = 1.591549431e-14; + else { + printf("WARNING: Because of float precision, I cannot get the factor to convert sqrt(E/ML^2)\n"); + printf("into THz, instead, I set it to be 1; you should check the unit used by LAMMPS.\n"); + eml2f = 1.; + } + + // now to allocate memory for DM + memory = new Memory; + DM_all = memory->create(DM_all, npt, fftdim2, "DynMat:DM_all"); + DM_q = memory->create(DM_q, fftdim,fftdim,"DynMat:DM_q"); + + // read all dynamical matrix info into DM_all + if ( fread(DM_all[0], sizeof(doublecomplex), npt*fftdim2, fp) != size_t(npt*fftdim2)){ + printf("\nError while reading the DM from file: %s\n", binfile); + fclose(fp); + exit(1); + } + + // now try to read unit cell info from the binary file + flag_latinfo = 0; + basis = memory->create(basis,nucell,sysdim,"DynMat:basis"); + attyp = memory->create(attyp,nucell, "DynMat:attyp"); + M_inv_sqrt = memory->create(M_inv_sqrt, nucell, "DynMat:M_inv_sqrt"); + int flag_mass_read = 0; + + if ( fread(&Tmeasure, sizeof(double), 1, fp) == 1) flag_latinfo |= 1; + if ( fread(&basevec[0], sizeof(double), 9, fp) == 9) flag_latinfo |= 2; + if ( fread(basis[0], sizeof(double), fftdim, fp) == fftdim) flag_latinfo |= 4; + if ( fread(&attyp[0], sizeof(int), nucell, fp) == nucell) flag_latinfo |= 8; + if ( fread(&M_inv_sqrt[0], sizeof(double), nucell, fp) == nucell) flag_mass_read = 1; + fclose(fp); + + if ((flag_latinfo&15) == 15){ + car2dir(flag_mass_read); + real2rec(); + + flag_latinfo = 1; + } else { + Tmeasure = 0.; + flag_latinfo = 0; + } + + // initialize interpolation + interpolate = new Interpolate(nx,ny,nz,fftdim2,DM_all); + if (flag_reset_gamma) interpolate->reset_gamma(); + + if ( flag_mass_read ){ // M_inv_sqrt info read, the data stored are force constant matrix instead of dynamical matrix. + + EnforceASR(); + + // get the dynamical matrix from force constant matrix: D = 1/M x Phi + for (int idq=0; idq< npt; idq++){ + int ndim =0; + for (int idim=0; idimset_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; icreate(work, lwork, "geteigen:work"); + rwork = memory->create(rwork, lrwork, "geteigen:rwork"); + iwork = memory->create(iwork, liwork, "geteigen:iwork"); + + zheevd_(&jobz, &uplo, &n, DM_q[0], &lda, w, work, &lwork, rwork, &lrwork, iwork, &liwork, &info); + + // to get w instead of w^2; and convert w into v (THz hopefully) + for (int i=0; i= 0.) w[i] = sqrt(w[i]); + else w[i] = -sqrt(-w[i]); + + w[i] *= eml2f; + } + + memory->destroy(work); + memory->destroy(rwork); + memory->destroy(iwork); + +return info; +} + +/* ---------------------------------------------------------------------------- + * method to get the Dynamical Matrix at q + * ---------------------------------------------------------------------------- */ +void DynMat::getDMq(double *q) +{ + interpolate->execute(q, DM_q[0]); +return; +} + +/* ---------------------------------------------------------------------------- + * method to get the Dynamical Matrix at q + * ---------------------------------------------------------------------------- */ +void DynMat::getDMq(double *q, double *wt) +{ + interpolate->execute(q, DM_q[0]); + + if (flag_skip && interpolate->UseGamma ) wt[0] = 0.; +return; +} + +/* ---------------------------------------------------------------------------- + * private method to convert the cartisan coordinate of basis into fractional + * ---------------------------------------------------------------------------- */ +void DynMat::car2dir(int flag) +{ + if (!flag){ // in newer version, this is done in fix-phonon + for (int i=0; i<3; i++){ + basevec[i] /= double(nx); + basevec[i+3] /= double(ny); + basevec[i+6] /= double(nz); + } + } + double mat[9]; + for (int idim=0; idim<9; idim++) mat[idim] = basevec[idim]; + GaussJordan(3, mat); + + for (int i=0; i= 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=0; l--){ + int rl = indxr[l]; + int cl = indxc[l]; + if (rl != cl){ + for (k=0; kset_method(); + +return; +} + +/* ---------------------------------------------------------------------------- + * Private method to display help info + * ---------------------------------------------------------------------------- */ +void DynMat::help() +{ + ShowVersion(); + printf("\nUsage:\n phana [options] [file]\n\n"); + printf("Available options:\n"); + printf(" -r To reset the dynamical matrix at the gamma point by a 4th order\n"); + printf(" polynomial interpolation along the [100] direction; this might be\n"); + printf(" useful for systems with charges. As for charged system, the dynamical\n"); + printf(" matrix at Gamma is far from accurate because of the long range nature\n"); + printf(" of Coulombic interaction. By reset it by interpolation, will partially\n"); + printf(" elliminate the unwanted behavior, but the result is still inaccurate.\n"); + printf(" By default, this is not set; and not expected for uncharged systems.\n\n"); + printf(" -s This will reset the dynamical matrix at the gamma point, too, but it\n"); + printf(" will also inform the code to skip all q-points that is in the vicinity\n"); + printf(" of the gamma point when evaluating phonon DOS and/or phonon dispersion.\n\n"); + printf(" By default, this is not set; and not expected for uncharged systems.\n\n"); + printf(" -h To print out this help info.\n\n"); + printf(" file To define the filename that carries the binary dynamical matrice generated\n"); + printf(" by fix-phonon. If not provided, the code will ask for it.\n"); + printf("\n\n"); + exit(0); +} + +/* ---------------------------------------------------------------------------- + * Private method to display the version info + * ---------------------------------------------------------------------------- */ +void DynMat::ShowVersion() +{ + printf(" ____ _ _ __ _ _ __ \n"); + printf(" ( _ \\( )_( ) /__\\ ( \\( ) /__\\ \n"); + printf(" )___/ ) _ ( /(__)\\ ) ( /(__)\\ \n"); + printf(" (__) (_) (_)(__)(__)(_)\\_)(__)(__)\n"); + printf("\nPHonon ANAlyzer for Fix-Phonon, version 2.%d, compiled on %s.\n", VERSION, __DATE__); + +return; +} +/* --------------------------------------------------------------------*/ diff --git a/tools/phonon/dynmat.h b/tools/phonon/dynmat.h new file mode 100644 index 0000000000..95cb2d7a32 --- /dev/null +++ b/tools/phonon/dynmat.h @@ -0,0 +1,66 @@ +#ifndef DYNMAT_H +#define DYNMAT_H + +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "memory.h" +#include "interpolate.h" + +extern "C"{ +#include "f2c.h" +#include "clapack.h" +} + +using namespace std; + +class DynMat { +public: + + DynMat(int, char**); + ~DynMat(); + + int nx, ny, nz, nucell; + int sysdim, fftdim; + double eml2f; + char *funit; + + void getDMq(double *); + void getDMq(double *, double *); + void writeDMq(double *); + void writeDMq(double *, const double, FILE *fp); + int geteigen(double *, int); + void reset_interp_method(); + + doublecomplex **DM_q; + + int flag_latinfo; + double Tmeasure, basevec[9], ibasevec[9]; + double **basis; + int *attyp; + +private: + + int flag_skip, flag_reset_gamma; + Interpolate *interpolate; + + Memory *memory; + int npt, fftdim2; + + int nasr; + void EnforceASR(); + + char *binfile, *dmfile; + double boltz, q[3]; + double *M_inv_sqrt; + + doublecomplex **DM_all; + + void car2dir(int); // to convert basis from cartisian coordinate into factional. + void real2rec(); + void GaussJordan(int, double *); + + void help(); + void ShowVersion(); +}; +#endif diff --git a/tools/phonon/green.cpp b/tools/phonon/green.cpp new file mode 100644 index 0000000000..5136ecc65c --- /dev/null +++ b/tools/phonon/green.cpp @@ -0,0 +1,250 @@ +#include +#include +#include +#include +#include "green.h" +#include + +#define MAXLINE 256 + +/******************************************************************************* + * 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); + alpha = memory->create(alpha, sysdim,nit, "Green_Green:alpha"); + beta = memory->create(beta, sysdim,nit+1,"Green_Green:beta"); + //ldos = memory->create(ldos, nw,sysdim, "Green_Green:ldos"); + + // 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 Z, z_m_a, r_x, rec_x, rec_x_inv; + double sr, si; + + double w = wmin; + for (int i=0; i(w*w, epson); + + for (int idim=0; idim 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 Z, rec_x, rec_x_inv; + std::complex cunit = std::complex(0.,1.); + + double w = wmin; + + for (int i=0; i(w*w, epson); + + for (int idim=0; idim(0.,0.); + + for (int j=0; j(b)?(b):(a)) +#define MAX(a,b) ((a)>(b)?(a):(b)) + +/* ---------------------------------------------------------------------------- + * 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){ + Dfdx = memory->create(Dfdx, Npt, ndim, "Interpolate_Interpolate:Dfdx"); + Dfdy = memory->create(Dfdy, Npt, ndim, "Interpolate_Interpolate:Dfdy"); + Dfdz = memory->create(Dfdz, Npt, ndim, "Interpolate_Interpolate:Dfdz"); + D2fdxdy = memory->create(D2fdxdy, Npt, ndim, "Interpolate_Interpolate:D2fdxdy"); + D2fdxdz = memory->create(D2fdxdz, Npt, ndim, "Interpolate_Interpolate:D2fdxdz"); + D2fdydz = memory->create(D2fdydz, Npt, ndim, "Interpolate_Interpolate:D2fdydz"); + D3fdxdydz = memory->create(D3fdxdydz, Npt, ndim, "Interpolate_Interpolate:D2fdxdydz"); + + 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; iidestroy(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= 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 +extern "C"{ +#include "f2c.h" +#include "clapack.h" +} + +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/main.cpp b/tools/phonon/main.cpp new file mode 100644 index 0000000000..d7d5baa1cc --- /dev/null +++ b/tools/phonon/main.cpp @@ -0,0 +1,18 @@ +#include "stdio.h" +#include "stdlib.h" +#include "dynmat.h" +#include "phonon.h" + +using namespace std; + +int main(int argc, char** argv) +{ + + DynMat *dynmat = new DynMat(argc, argv); + Phonon *phonon = new Phonon(dynmat); + + delete phonon; + delete dynmat; + +return 0; +} diff --git a/tools/phonon/memory.cpp b/tools/phonon/memory.cpp new file mode 100644 index 0000000000..4d65e83ef8 --- /dev/null +++ b/tools/phonon/memory.cpp @@ -0,0 +1,51 @@ +#include "stdio.h" +#include "stdlib.h" +#include "string.h" +#include "memory.h" + +/* ---------------------------------------------------------------------- + safe malloc +------------------------------------------------------------------------- */ + +void *Memory::smalloc(bigint nbytes, const char *name) +{ + if (nbytes == 0) return NULL; + + void *ptr = malloc(nbytes); + if (ptr == NULL) printf("Failed to allocate " BIGINT_FORMAT "bytes for array %s", nbytes,name); + return ptr; +} + +/* ---------------------------------------------------------------------- + safe realloc +------------------------------------------------------------------------- */ +void *Memory::srealloc(void *ptr, bigint nbytes, const char *name) +{ + if (nbytes == 0) { + destroy(ptr); + return NULL; + } + + ptr = realloc(ptr,nbytes); + if (ptr == NULL) printf("Failed to reallocate " BIGINT_FORMAT "bytes for array %s", nbytes,name); + return ptr; +} + +/* ---------------------------------------------------------------------- + safe free +------------------------------------------------------------------------- */ + +void Memory::sfree(void *ptr) +{ + if (ptr == NULL) return; + free(ptr); +} + +/* ---------------------------------------------------------------------- + erroneous usage of templated create/grow functions +------------------------------------------------------------------------- */ + +void Memory::fail(const char *name) +{ + printf("Cannot create/grow a vector/array of pointers for %s",name); +} diff --git a/tools/phonon/memory.h b/tools/phonon/memory.h new file mode 100644 index 0000000000..ae2feceba3 --- /dev/null +++ b/tools/phonon/memory.h @@ -0,0 +1,430 @@ +#ifndef LMP_MEMORY_H +#define LMP_MEMORY_H + +#define __STDC_LIMIT_MACROS +#define __STDC_FORMAT_MACROS + +#include "stdio.h" +#include "stdlib.h" +#include "limits.h" +#include "stdint.h" +#include "inttypes.h" + +typedef int64_t bigint; +#define BIGINT_FORMAT "%" PRId64 +#define ATOBIGINT atoll + +class Memory { + public: + Memory(){}; + + void *smalloc(bigint n, const char *); + void *srealloc(void *, bigint n, const char *); + void sfree(void *); + void fail(const char *); + +/* ---------------------------------------------------------------------- + create a 1d array +------------------------------------------------------------------------- */ + + template + TYPE *create(TYPE *&array, int n, const char *name) + { + bigint nbytes = sizeof(TYPE) * n; + array = (TYPE *) smalloc(nbytes,name); + return array; + }; + + template + TYPE **create(TYPE **&array, int n, const char *name) {fail(name);} + +/* ---------------------------------------------------------------------- + grow or shrink 1d array +------------------------------------------------------------------------- */ + + template + TYPE *grow(TYPE *&array, int n, const char *name) + { + if (array == NULL) return create(array,n,name); + + bigint nbytes = sizeof(TYPE) * n; + array = (TYPE *) srealloc(array,nbytes,name); + return array; + }; + + template + TYPE **grow(TYPE **&array, int n, const char *name) {fail(name);} + +/* ---------------------------------------------------------------------- + destroy a 1d array +------------------------------------------------------------------------- */ + + template + void destroy(TYPE *array) + { + sfree(array); + }; + +/* ---------------------------------------------------------------------- + create a 1d array with index from nlo to nhi inclusive + cannot grow it +------------------------------------------------------------------------- */ + + template + TYPE *create1d_offset(TYPE *&array, int nlo, int nhi, const char *name) + { + bigint nbytes = sizeof(TYPE) * (nhi-nlo+1); + array = (TYPE *) smalloc(nbytes,name); + array -= nlo; + return array; + } + + template + TYPE **create1d_offset(TYPE **&array, int nlo, int nhi, const char *name) + {fail(name);} + +/* ---------------------------------------------------------------------- + destroy a 1d array with index offset +------------------------------------------------------------------------- */ + + template + void destroy1d_offset(TYPE *array, int offset) + { + if (array) sfree(&array[offset]); + } + +/* ---------------------------------------------------------------------- + create a 2d array +------------------------------------------------------------------------- */ + + template + TYPE **create(TYPE **&array, int n1, int n2, const char *name) + { + bigint nbytes = sizeof(TYPE) * n1*n2; + TYPE *data = (TYPE *) smalloc(nbytes,name); + nbytes = sizeof(TYPE *) * n1; + array = (TYPE **) smalloc(nbytes,name); + + int n = 0; + for (int i = 0; i < n1; i++) { + array[i] = &data[n]; + n += n2; + } + return array; + } + + template + TYPE ***create(TYPE ***&array, int n1, int n2, const char *name) + {fail(name);} + +/* ---------------------------------------------------------------------- + grow or shrink 1st dim of a 2d array + last dim must stay the same +------------------------------------------------------------------------- */ + + template + TYPE **grow(TYPE **&array, int n1, int n2, const char *name) + { + if (array == NULL) return create(array,n1,n2,name); + + bigint nbytes = sizeof(TYPE) * n1*n2; + TYPE *data = (TYPE *) srealloc(array[0],nbytes,name); + nbytes = sizeof(TYPE *) * n1; + array = (TYPE **) srealloc(array,nbytes,name); + + int n = 0; + for (int i = 0; i < n1; i++) { + array[i] = &data[n]; + n += n2; + } + return array; + } + + template + TYPE ***grow(TYPE ***&array, int n1, int n2, const char *name) + {fail(name);} + +/* ---------------------------------------------------------------------- + destroy a 2d array +------------------------------------------------------------------------- */ + + template + void destroy(TYPE **array) + { + if (array == NULL) return; + sfree(array[0]); + sfree(array); + } + +/* ---------------------------------------------------------------------- + create a 2d array with 2nd index from n2lo to n2hi inclusive + cannot grow it +------------------------------------------------------------------------- */ + + template + TYPE **create2d_offset(TYPE **&array, int n1, int n2lo, int n2hi, + const char *name) + { + int n2 = n2hi - n2lo + 1; + create(array,n1,n2,name); + for (int i = 0; i < n1; i++) array[i] -= n2lo; + return array; + } + + template + TYPE ***create2d_offset(TYPE ***&array, int n1, int n2lo, int n2hi, + const char *name) {fail(name);} + +/* ---------------------------------------------------------------------- + destroy a 2d array with 2nd index offset +------------------------------------------------------------------------- */ + + template + void destroy2d_offset(TYPE **array, int offset) + { + if (array == NULL) return; + sfree(&array[0][offset]); + sfree(array); + } + +/* ---------------------------------------------------------------------- + create a 3d array +------------------------------------------------------------------------- */ + + template + TYPE ***create(TYPE ***&array, int n1, int n2, int n3, const char *name) + { + bigint nbytes = sizeof(TYPE) * n1*n2*n3; + TYPE *data = (TYPE *) smalloc(nbytes,name); + nbytes = sizeof(TYPE *) * n1*n2; + TYPE **plane = (TYPE **) smalloc(nbytes,name); + nbytes = sizeof(TYPE **) * n1; + array = (TYPE ***) smalloc(nbytes,name); + + int i,j; + int n = 0; + for (i = 0; i < n1; i++) { + array[i] = &plane[i*n2]; + for (j = 0; j < n2; j++) { + plane[i*n2+j] = &data[n]; + n += n3; + } + } + return array; + } + + template + TYPE ****create(TYPE ****&array, int n1, int n2, int n3, const char *name) + {fail(name);} + +/* ---------------------------------------------------------------------- + grow or shrink 1st dim of a 3d array + last 2 dims must stay the same +------------------------------------------------------------------------- */ + + template + TYPE ***grow(TYPE ***&array, int n1, int n2, int n3, const char *name) + { + if (array == NULL) return create(array,n1,n2,n3,name); + + bigint nbytes = sizeof(TYPE) * n1*n2*n3; + TYPE *data = (TYPE *) srealloc(array[0][0],nbytes,name); + nbytes = sizeof(TYPE *) * n1*n2; + TYPE **plane = (TYPE **) srealloc(array[0],nbytes,name); + nbytes = sizeof(TYPE **) * n1; + array = (TYPE ***) srealloc(array,nbytes,name); + + int i,j; + int n = 0; + for (i = 0; i < n1; i++) { + array[i] = &plane[i*n2]; + for (j = 0; j < n2; j++) { + plane[i*n2+j] = &data[n]; + n += n3; + } + } + return array; + } + + template + TYPE ****grow(TYPE ****&array, int n1, int n2, int n3, const char *name) + {fail(name);} + +/* ---------------------------------------------------------------------- + destroy a 3d array +------------------------------------------------------------------------- */ + + template + void destroy(TYPE ***array) + { + if (array == NULL) return; + sfree(array[0][0]); + sfree(array[0]); + sfree(array); + } + +/* ---------------------------------------------------------------------- + create a 3d array with 1st index from n1lo to n1hi inclusive + cannot grow it +------------------------------------------------------------------------- */ + + template + TYPE ***create3d_offset(TYPE ***&array, int n1lo, int n1hi, + int n2, int n3, const char *name) + { + int n1 = n1hi - n1lo + 1; + create(array,n1,n2,n3,name); + array -= n1lo; + return array; + } + + template + TYPE ****create3d_offset(TYPE ****&array, int n1lo, int n1hi, + int n2, int n3, const char *name) + {fail(name);} + +/* ---------------------------------------------------------------------- + free a 3d array with 1st index offset +------------------------------------------------------------------------- */ + + template + void destroy3d_offset(TYPE ***array, int offset) + { + if (array) destroy(&array[offset]); + } + +/* ---------------------------------------------------------------------- + create a 3d array with + 1st index from n1lo to n1hi inclusive, + 2nd index from n2lo to n2hi inclusive, + 3rd index from n3lo to n3hi inclusive + cannot grow it +------------------------------------------------------------------------- */ + + template + TYPE ***create3d_offset(TYPE ***&array, int n1lo, int n1hi, + int n2lo, int n2hi, int n3lo, int n3hi, + const char *name) + { + int n1 = n1hi - n1lo + 1; + int n2 = n2hi - n2lo + 1; + int n3 = n3hi - n3lo + 1; + create(array,n1,n2,n3,name); + + for (int i = 0; i < n1*n2; i++) array[0][i] -= n3lo; + for (int i = 0; i < n1; i++) array[i] -= n2lo; + array -= n1lo; + return array; + } + + template + TYPE ****create3d_offset(TYPE ****&array, int n1lo, int n1hi, + int n2lo, int n2hi, int n3lo, int n3hi, + const char *name) + {fail(name);} + +/* ---------------------------------------------------------------------- + free a 3d array with all 3 indices offset +------------------------------------------------------------------------- */ + + template + void destroy3d_offset(TYPE ***array, + int n1_offset, int n2_offset, int n3_offset) + { + if (array == NULL) return; + sfree(&array[n1_offset][n2_offset][n3_offset]); + sfree(&array[n1_offset][n2_offset]); + sfree(&array[n1_offset]); + } + +/* ---------------------------------------------------------------------- + create a 4d array +------------------------------------------------------------------------- */ + + template + TYPE ****create(TYPE ****&array, int n1, int n2, int n3, int n4, + const char *name) + { + bigint nbytes = sizeof(TYPE) * n1*n2*n3*n4; + TYPE *data = (double *) smalloc(nbytes,name); + nbytes = sizeof(TYPE *) * n1*n2*n3; + TYPE **cube = (double **) smalloc(nbytes,name); + nbytes = sizeof(TYPE **) * n1*n2; + TYPE ***plane = (double ***) smalloc(nbytes,name); + nbytes = sizeof(TYPE ***) * n1; + array = (double ****) smalloc(nbytes,name); + + int i,j,k; + int n = 0; + for (i = 0; i < n1; i++) { + array[i] = &plane[i*n2]; + for (j = 0; j < n2; j++) { + plane[i*n2+j] = &cube[i*n2*n3+j*n3]; + for (k = 0; k < n3; k++) { + cube[i*n2*n3+j*n3+k] = &data[n]; + n += n4; + } + } + } + return array; + } + + template + TYPE *****create(TYPE *****&array, int n1, int n2, int n3, int n4, + const char *name) + {fail(name);} + +/* ---------------------------------------------------------------------- + destroy a 4d array +------------------------------------------------------------------------- */ + + template + void destroy(TYPE ****array) + { + if (array == NULL) return; + sfree(array[0][0][0]); + sfree(array[0][0]); + sfree(array[0]); + sfree(array); + } + +/* ---------------------------------------------------------------------- + memory usage of arrays, including pointers +------------------------------------------------------------------------- */ + + template + bigint usage(TYPE *array, int n) + { + bigint bytes = sizeof(TYPE) * n; + return bytes; + } + + template + bigint usage(TYPE **array, int n1, int n2) + { + bigint bytes = sizeof(TYPE) * n1*n2; + bytes += sizeof(TYPE *) * n1; + return bytes; + } + + template + bigint usage(TYPE ***array, int n1, int n2, int n3) + { + bigint bytes = sizeof(TYPE) * n1*n2*n3; + bytes += sizeof(TYPE *) * n1*n2; + bytes += sizeof(TYPE **) * n1; + return bytes; + } + + template + bigint usage(TYPE ****array, int n1, int n2, int n3, int n4) + { + bigint bytes = sizeof(TYPE) * n1*n2*n3*n4; + bytes += sizeof(TYPE *) * n1*n2*n3; + bytes += sizeof(TYPE **) * n1*n2; + bytes += sizeof(TYPE ***) * n1; + return bytes; + } +}; + +#endif diff --git a/tools/phonon/phonon.cpp b/tools/phonon/phonon.cpp new file mode 100644 index 0000000000..5ee34a6ce8 --- /dev/null +++ b/tools/phonon/phonon.cpp @@ -0,0 +1,1266 @@ +#include +#include "string.h" +#include "phonon.h" +#include "green.h" +#include "timer.h" + +#ifdef UseSPG +extern "C"{ +#include "spglib.h" +} +#endif + +#define MAXLINE 256 +#define MIN(a,b) ((a)>(b)?(b):(a)) +#define MAX(a,b) ((a)>(b)?(a):(b)) + +/* ---------------------------------------------------------------------------- + * Class Phonon is the main driver to calculate phonon DOS, phonon + * dispersion curve and some other things. + * ---------------------------------------------------------------------------- */ +Phonon::Phonon(DynMat *dm) +{ + // create memory + memory = new Memory(); + + // pass the class from main + dynmat = dm; + sysdim = dynmat->sysdim; + ndim = dynmat->fftdim; + dos = NULL; + ldos = NULL; + qpts = NULL; + wt = NULL; + eigs = NULL; + locals = NULL; + +#ifdef UseSPG + attyp = NULL; + atpos = NULL; +#endif + + // display the menu + char str[MAXLINE]; + while ( 1 ){ + printf("\n"); for (int i=0; i<37;i++) printf("="); printf(" Menu "); for (int i=0; i<37;i++) printf("="); printf("\n"); + printf(" 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(" -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"); + + // 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 ==-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= 2){ + fmin = atof(strtok(str," \t\n\r\f")); + fmax = atof(strtok(NULL," \t\n\r\f")); + } + if (fmin > fmax){double swap = fmin; fmin = fmax; fmax = swap;} + printf("The fequency range for your phonon DOS is [%g %g].\n", fmin, fmax); + + ndos = 201; + printf("Please input the number of intervals [%d]: ", ndos); + if (count_words(fgets(str,MAXLINE,stdin)) > 0) ndos = atoi(strtok(str," \t\n\r\f")); + ndos += (ndos+1)%2; ndos = MAX(2,ndos); + + + df = (fmax-fmin)/double(ndos-1); + rdf = 1./df; + memory->destroy(dos); + dos = memory->create(dos, ndos, "pdos:dos"); + for (int i=0; i 0.){ + for (int j=0; j=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; + Hessian = memory->create(Hessian, ndim, ndim, "phonon_ldos:Hessian"); + + double q0[3]; + q0[0] = q0[1] = q0[2] = 0.; + + dynmat->getDMq(q0); + for (int i=0; iDM_q[i][j].r*scale; + + if (ndim < 300){ + double *egvs = new double [ndim]; + dynmat->geteigen(egvs, 0); + + fmin = fmax = egvs[0]; + for (int i=1; inucell); + printf("Please input the index/index range/index range and increment of atom(s)\n"); + printf("in the unit cell to evaluate LDOS, q to exit [%d]: ", ik); + int nr = count_words( fgets(str,MAXLINE,stdin) ); + if (nr < 1){ + istr = iend = ik; + iinc = 1; + } else if (nr == 1) { + char *ptr = strtok(str," \t\n\r\f"); + if (strcmp(ptr,"q") == 0) break; + + ik = atoi(ptr); + if (ik < 0 || ik >= dynmat->nucell) break; + istr = iend = ik; + iinc = 1; + } else if (nr == 2) { + istr = atoi(strtok(str," \t\n\r\f")); + iend = atoi(strtok(NULL," \t\n\r\f")); + iinc = 1; + if (istr < 0||iend >= dynmat->nucell||istr > iend) break; + } else if (nr >= 3) { + istr = atoi(strtok(str," \t\n\r\f")); + iend = atoi(strtok(NULL," \t\n\r\f")); + iinc = atoi(strtok(NULL," \t\n\r\f")); + if (istr<0 || iend >= dynmat->nucell || istr > iend || iinc<1) break; + } + + printf("Please input the frequency range to evaluate LDOS [%g %g]: ", fmin, fmax); + if (count_words(fgets(str,MAXLINE,stdin)) >= 2){ + fmin = atof(strtok(str," \t\n\r\f")); + fmax = atof(strtok(NULL," \t\n\r\f")); + } + if (fmax < fmin) break; + printf("The frequency range for your LDOS is [%g %g].\n", fmin, fmax); + + printf("Please input the desired number of points in LDOS [%d]: ", ndos); + if (count_words(fgets(str,MAXLINE,stdin)) > 0) ndos = atoi(strtok(str," \t\n\r\f")); + if (ndos < 2) break; + ndos += (ndos+1)%2; + + printf("Please input the maximum # of Lanczos iterations [%d]: ", nit); + if (count_words(fgets(str,MAXLINE,stdin)) > 0) nit = atoi(strtok(str," \t\n\r\f")); + if (nit < 1) break; + + printf("Please input the value of epsilon for delta-function [%g]: ", eps); + if (count_words(fgets(str,MAXLINE,stdin)) > 0) eps = atof(strtok(str," \t\n\r\f")); + if (eps <= 0.) break; + + // prepare array for local pdos + nlocal = 0; + for (ik = istr; ik <= iend; ik += iinc) nlocal++; + memory->destroy(ldos); + ldos = memory->create(ldos,nlocal,ndos,dynmat->sysdim,"ldos_rsgf:ldos"); + + memory->destroy(locals); + locals = memory->create(locals, nlocal, "ldos_rsgf:locals"); + + df = (fmax-fmin)/double(ndos-1); + rdf = 1./df; + + // to measure the LDOS via real space Green's function method + int ilocal = 0; + for (ik = istr; ik <= iend; ik += iinc){ + locals[ilocal] = ik; + + // time info + Timer *time = new Timer(); + printf("\nNow to compute the LDOS for atom %d by Real Space Greens function method ...\n", ik); + fflush(stdout); + + // run real space green's function calculation + Green *green = new Green(dynmat->nucell, dynmat->sysdim, nit, fmin, fmax, ndos, eps, Hessian, ik, ldos[ilocal++]); + delete green; + + time->stop(); time->print(); delete time; + } + + Normalize(); + writeLDOS(); + + // evaluate the local vibrational thermal properties optionally + local_therm(); + + } + memory->destroy(Hessian); + +return; +} + +/*------------------------------------------------------------------------------ + * Private method to evaluate the phonon dispersion curves + *----------------------------------------------------------------------------*/ +void Phonon::pdisp() +{ + // ask the output file name and write the header. + char str[MAXLINE]; + printf("Please input the filename to output the dispersion data [pdisp.dat]:"); + if (count_words(fgets(str,MAXLINE,stdin)) < 1) strcpy(str, "pdisp.dat"); + char *ptr = strtok(str," \t\n\r\f"); + char *fname = new char[strlen(ptr)+1]; + strcpy(fname,ptr); + + FILE *fp = fopen(fname, "w"); + fprintf(fp,"# q qr freq\n"); + fprintf(fp,"# 2pi/L 2pi/L %s\n", dynmat->funit); + + // to store the nodes of the dispersion curve + std::vector nodes; nodes.clear(); + + // now the calculate the dispersion curve + double qstr[3], qend[3], q[3], qinc[3], qr=0., dq; + int nq = MAX(MAX(dynmat->nx,dynmat->ny),dynmat->nz)/2+1; + qend[0] = qend[1] = qend[2] = 0.; + + double *egvs = new double [ndim]; + while (1){ + for (int i=0; i<3; i++) qstr[i] = qend[i]; + + int quit = 0; + printf("\nPlease input the start q-point in unit of B1->B3, q to exit [%g %g %g]: ", qstr[0], qstr[1], qstr[2]); + int n = count_words(fgets(str,MAXLINE,stdin)); + ptr = strtok(str," \t\n\r\f"); + if ((n == 1) && (strcmp(ptr,"q") == 0)) break; + else if (n >= 3){ + qstr[0] = atof(ptr); + qstr[1] = atof(strtok(NULL," \t\n\r\f")); + qstr[2] = atof(strtok(NULL," \t\n\r\f")); + } + + do printf("Please input the end q-point in unit of B1->B3: "); + while (count_words(fgets(str,MAXLINE,stdin)) < 3); + qend[0] = atof(strtok(str," \t\n\r\f")); + qend[1] = atof(strtok(NULL," \t\n\r\f")); + qend[2] = atof(strtok(NULL," \t\n\r\f")); + + printf("Please input the # of points along the line [%d]: ", nq); + if (count_words(fgets(str,MAXLINE,stdin)) > 0) nq = atoi(strtok(str," \t\n\r\f")); + nq = MAX(nq,2); + + for (int i=0; i<3; i++) qinc[i] = (qend[i]-qstr[i])/double(nq-1); + dq = sqrt(qinc[0]*qinc[0]+qinc[1]*qinc[1]+qinc[2]*qinc[2]); + + nodes.push_back(qr); + for (int i=0; i<3; i++) q[i] = qstr[i]; + for (int ii=0; iigetDMq(q, &wii); + if (wii > 0.){ + dynmat->geteigen(egvs, 0); + fprintf(fp,"%lg %lg %lg %lg ", q[0], q[1], q[2], qr); + for (int i=0; i 0.) nodes.push_back(qr); + fclose(fp); + delete []egvs; + + // write the gnuplot script which helps to visualize the result + int nnd = nodes.size(); + if (nnd > 1){ + fp = fopen("pdisp.gnuplot", "w"); + fprintf(fp,"set term post enha colo 20\nset out %cpdisp.eps%c\n\n",char(34),char(34)); + fprintf(fp,"set xlabel %cq%c\n",char(34),char(34)); + fprintf(fp,"set ylabel %cfrequency (THz)%c\n\n",char(34),char(34)); + fprintf(fp,"set xrange [0:%lg]\nset yrange [0:*]\n\n", nodes[nnd-1]); + fprintf(fp,"set grid xtics\n"); + fprintf(fp,"# {/Symbol G} will give you letter gamma in the label\nset xtics ("); + for (int i=0; igetDMq(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; iDM_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; inucell; j++){ + int ipos = j * sysdim; + double sum = 0.; + fprintf(fp,"%d", j+1); + for (int idim=0; idimnx,dynmat->ny),dynmat->nz)/2; + qend[0] = qend[1] = qend[2] = 0.; + + while (1){ + + for (int i=0; i<3; i++) qstr[i] = qend[i]; + + 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; iigetDMq(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; + tmp = memory->create(tmp, npt, "smooth:tmp"); + table = 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; idestroy(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 0.); + fclose(fp); + +return; +} + +/* ---------------------------------------------------------------------------- + * Private method to calculate the local thermal properties + * ---------------------------------------------------------------------------- */ +void Phonon::local_therm() +{ + char str[MAXLINE]; + printf("\nWould you like to compute the local thermal properties (y/n)[n]: "); + if (count_words(fgets(str,MAXLINE,stdin)) < 1) return; + char *ptr = strtok(str," \t\n\r\f"); + if (strcmp(ptr,"y") != 0 && strcmp(ptr, "Y") != 0 && strcmp(ptr, "yes") != 0) return; + + printf("Please input the filename to output vibrational thermal info [localtherm.dat]: "); + if (count_words(fgets(str,MAXLINE,stdin)) < 1) strcpy(str, "localtherm.dat"); + + FILE *fp = fopen(strtok(str," \t\n\r\f"), "w"); + fprintf(fp,"# atom Temp U_vib (eV) S_vib (kB) F_vib (eV) C_vib (kB) ZPE (eV)\n"); + fprintf(fp,"# ------------ ------------ ----------- ----------- ------------\n"); + fprintf(fp,"# Ux Uy Uz Ut Sx Sy Sz St Fx Fy Fz Ft Cx Cy Cz Ct Zx Zy Zz Zt\n"); + fprintf(fp,"# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22\n"); + fprintf(fp,"#-------------------------------------------------------------------------------\n"); + + double **Uvib, **Svib, **Fvib, **Cvib, **ZPE; + Uvib = memory->create(Uvib,nlocal,sysdim,"local_therm:Uvib"); + Svib = memory->create(Svib,nlocal,sysdim,"local_therm:Svib"); + Fvib = memory->create(Fvib,nlocal,sysdim,"local_therm:Fvib"); + Cvib = memory->create(Cvib,nlocal,sysdim,"local_therm:Cvib"); + ZPE = memory->create(ZPE ,nlocal,sysdim,"local_therm:ZPE"); + // 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; inx, 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); + wt = memory->create(wt, nq, "QMesh:wt"); + qpts = memory->create(qpts, nq, 3, "QMesh:qpts"); + + int iq = 0; + for (int i=0; icreate(atpos, dynmat->nucell,3,"QMesh:atpos"); + attyp = memory->create(attyp, dynmat->nucell, "QMesh:attyp"); + + for (int i=0; inucell; i++) + for (int idim=0; idim<3; idim++) atpos[i][idim] = 0.; + for (int i=0; i<3; i++) latvec[i][i] = 1.; + + int flag_lat_info_read = dynmat->flag_latinfo; + + if ( flag_lat_info_read ){ // get unit cell info from binary file; done by dynmat + num_atom = dynmat->nucell; + // set default, in case system dimension under study is not 3. + for (int i=0; inucell; i++) + for (int idim=0; idim<3; idim++) atpos[i][idim] = 0.; + for (int i=0; i<3; i++) latvec[i][i] = 1.; + + // get atomic type info + for (int i=0; iattyp[i]; + // get unit cell vector info + int ndim = 0; + for (int idim=0; idim<3; idim++) + for (int jdim=0; jdim<3; jdim++) latvec[jdim][idim] = dynmat->basevec[ndim++]; + // get atom position in unit cell; fractional + for (int i=0; ibasis[i][idim]; + + // 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 0) latsrc = atoi(strtok(str," \t\n\r\f")); + latsrc = 2-latsrc%2; + /*---------------------------------------------------------------- + * Ask for lattice info from the user; the format of the file is: + * A1_x A1_y A1_z + * A2_x A2_y A2_z + * A3_x A3_y A3_z + * natom + * Type_1 sx_1 sy_1 sz_1 + * ... + * Type_n sx_n sy_n sz_n + *----------------------------------------------------------------*/ + if (latsrc == 1){ // to read unit cell info from file; get file name first + do printf("Please input the file name containing the unit cell info: "); + while (count_words(fgets(str,MAXLINE,stdin)) < 1); + char *fname = strtok(str," \t\n\r\f"); + FILE *fp = fopen(fname,"r"); fname = NULL; + + if (fp == NULL) latsrc = 2; + else { + for (int i=0; i<3; i++){ // read unit cell vector info; # of atoms per unit cell + if (count_words(fgets(str,MAXLINE,fp)) < 3){ + latsrc = 2; + break; + } + latvec[0][i] = atof(strtok(str, " \t\n\r\f")); + latvec[1][i] = atof(strtok(NULL," \t\n\r\f")); + latvec[2][i] = atof(strtok(NULL," \t\n\r\f")); + } + if (count_words(fgets(str,MAXLINE,fp)) < 1) latsrc = 2; + else { + num_atom = atoi(strtok(str," \t\n\r\f")); + if (num_atom > dynmat->nucell){ + printf("\nError: # of atoms read from file (%d) is bigger than that given by the dynamical matrix (%d)!\n", num_atom, dynmat->nucell); + return; + } + + for (int i=0; i dynmat->nucell){ + printf("\nError: # of atoms input (%d) is bigger than that given by the dynamical matrix (%d)!\n", num_atom, dynmat->nucell); + return; + } + + for (int i=0; i= 1.0.3 is used + //nq = spg_get_ir_reciprocal_mesh(grid_point, map, mesh, shift, is_time_reversal, latvec, pos, attyp, num_atom, symprec); + + wt = memory->create(wt, nq, "QMesh:wt"); + qpts = memory->create(qpts, nq,3,"QMesh:qpts"); + + int *iq2idx = new int[num_grid]; + int numq = 0; + for (int i=0; i %d points\n", nx,ny,nz,nq); + +return; +} + +/* ---------------------------------------------------------------------------- + * Private method to calculate the local phonon DOS and total phonon DOS based + * on the eigenvectors + * ---------------------------------------------------------------------------- */ +void Phonon::ldos_egv() +{ + // get local position info + char str[MAXLINE], *ptr; + printf("\nThe # of atoms per cell is: %d, please input the atom IDs to compute\n", dynmat->nucell); + printf("local PDOS, IDs begin with 0: "); + int nmax = count_words(fgets(str,MAXLINE,stdin)); + if (nmax < 1) return; + + memory->destroy(locals); + locals = memory->create(locals, nmax, "ldos_egv:locals"); + + 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= 2) { + fmin = atof(strtok(str," \t\n\r\f")); + fmax = atof(strtok(NULL," \t\n\r\f")); + } + if (fmax < 0. || fmax < fmin) return; + + ndos = 201; + printf("Please input your desired # of points in PDOS [%d]: ", ndos); + if (count_words(fgets(str,MAXLINE,stdin)) > 0) ndos = atoi(strtok(str," \t\n\r\f")); + if (ndos < 2) return; + ndos += (ndos+1)%2; + + df = (fmax-fmin)/double(ndos-1); + rdf = 1./df; + + // get the q-points + QMesh(); + + // allocate memory for DOS and LDOS + memory->destroy(dos); + memory->destroy(ldos); + + dos = memory->create(dos, ndos,"ldos_egv:dos"); + ldos = memory->create(ldos,nlocal,ndos,sysdim,"ldos_egv:ldos"); + + for (int i=0; i 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; iqgetDMq(qpts[iq], &wt[iq]); + if (wt[iq] <= 0.) continue; + + dynmat->geteigen(&egval[0], 1); + + for (int idim=0; idim= 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 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 10) nprint = nq/10; + else nprint = 1; + Timer *time = new Timer(); + + printf("\nNow to compute the phonons "); fflush(stdout); + // now to calculate the frequencies at all q-points + memory->destroy(eigs); + eigs = memory->create(eigs, nq,ndim,"QMesh_eigs"); + + for (int iq=0; iqgetDMq(qpts[iq], &wt[iq]); + if (wt[iq] > 0.) dynmat->geteigen(eigs[iq], 0); + } + printf("Done!\n"); + time->stop(); time->print(); delete time; + +return; +} + +/*------------------------------------------------------------------------------ + * Method to count # of words in a string, without destroying the string + *----------------------------------------------------------------------------*/ +int Phonon::count_words(const char *line) +{ + int n = strlen(line) + 1; + char *copy; + copy = memory->create(copy, n, "count_words:copy"); + strcpy(copy,line); + + char *ptr; + if (ptr = strchr(copy,'#')) *ptr = '\0'; + + if (strtok(copy," \t\n\r\f") == NULL) { + memory->destroy(copy); + return 0; + } + n = 1; + while (strtok(NULL," \t\n\r\f")) n++; + + memory->destroy(copy); + return n; +} + +/*----------------------------------------------------------------------------*/ diff --git a/tools/phonon/phonon.h b/tools/phonon/phonon.h new file mode 100644 index 0000000000..02552fc409 --- /dev/null +++ b/tools/phonon/phonon.h @@ -0,0 +1,59 @@ +#ifndef PHONON_H +#define PHONON_H + +#include "stdio.h" +#include "stdlib.h" +#include +#include "dynmat.h" +#include "memory.h" + +using namespace std; + +class Phonon{ +public: + Phonon(DynMat *); + ~Phonon(); + + DynMat *dynmat; + +private: + int nq, ndim, sysdim; + double **qpts, *wt; + double **eigs; + + int ndos, nlocal, *locals; + double *dos, fmin, fmax, df, rdf; + double ***ldos; + + Memory *memory; + + void QMesh(); + void ComputeAll(); + + void pdos(); + void pdisp(); + void therm(); + + void ldos_egv(); + void ldos_rsgf(); + void local_therm(); + + void dmanyq(); + void vfanyq(); + void DMdisp(); + void vecanyq(); + + void smooth(double *, const int); + void writeDOS(); + void writeLDOS(); + void Normalize(); + + int count_words(const char *); + +#ifdef UseSPG + int num_atom, *attyp; + double latvec[3][3], **atpos; +#endif +}; + +#endif diff --git a/tools/phonon/timer.cpp b/tools/phonon/timer.cpp new file mode 100644 index 0000000000..0b66a9572e --- /dev/null +++ b/tools/phonon/timer.cpp @@ -0,0 +1,41 @@ +#include "timer.h" + +Timer::Timer() +{ + flag = 0; + start(); +return; +} + +void Timer::start() +{ + t1 = clock(); + flag |= 1; + +return; +} + +void Timer::stop() +{ + if ( flag&1 ) { + t2 = clock(); + flag |= 2; + } +return; +} + +void Timer::print() +{ + if ( (flag&3) != 3) return; + + cpu_time_used = ((double) (t2 - t1)) / CLOCKS_PER_SEC; + printf("Total CPU time used: %g seconds.\n", cpu_time_used); + +return; +} + +double Timer::elapse() +{ + if ( (flag&3) != 3) return 0.; + else return ((double) (t2 - t1)) / CLOCKS_PER_SEC; +} diff --git a/tools/phonon/timer.h b/tools/phonon/timer.h new file mode 100644 index 0000000000..cd04dea56d --- /dev/null +++ b/tools/phonon/timer.h @@ -0,0 +1,24 @@ +#ifndef TIMER_H +#define TIMER_H + +#include "stdio.h" +#include "stdlib.h" +#include "time.h" + +class Timer { +public: + Timer(); + + void start(); + void stop(); + void print(); + double elapse(); + +private: + clock_t t1, t2; + double cpu_time_used; + + int flag; +}; + +#endif