jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14366 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2015-12-15 15:59:01 +00:00
parent 75de12f26a
commit 5f802f86b5
10 changed files with 399 additions and 543 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

25
tools/phonon/Makefile
View File

@ -9,23 +9,20 @@ INC = $(LPKINC) $(TCINC) $(SPGINC)
LIB = $(LPKLIB) $(TCLIB) $(SPGLIB) LIB = $(LPKLIB) $(TCLIB) $(SPGLIB)
# #
# cLapack library needed # cLapack library needed
LPKINC = -I/opt/clapack/3.2.1/include LPKINC = -I/opt/libs/clapack/3.2.1/include
LPKLIB = -L/opt/clapack/3.2.1/lib -lclapack -lblas -lf2c #-lm LPKLIB = -L/opt/libs/clapack/3.2.1/lib -lclapack -lblas -lf2c #-lm
# #
# Tricubic library needed # Tricubic library needed
TCINC = -I/opt/tricubic/1.0/include TCINC = -I/opt/libs/tricubic/1.0/include
TCLIB = -L/opt/tricubic/1.0/lib -ltricubic TCLIB = -L/opt/libs/tricubic/1.0/lib -ltricubic
# #
# spglib 0.7.1, used to get the irreducible q-points # spglib 1.8.2, used to get the irreducible q-points
# if UFLAG is not set, spglib won't be used. # if UFLAG is not set, spglib won't be used.
UFLAG = -DUseSPG UFLAG = -DUseSPG
SPGINC = -I/opt/spglib/0.7.1/include SPGINC = -I/opt/libs/spglib/1.8.2/include
SPGLIB = -L/opt/spglib/0.7.1/lib -lsymspg SPGLIB = -L/opt/libs/spglib/1.8.2/lib -lsymspg
# if spglib > 0.7.1 is used, please # if spglib other than version 1.8.2 is used, please
# 1) modify file phonon.cpp, instruction can be found by searching 0.7.1 # modify file phonon.cpp, instruction can be found by searching 1.8.2
# 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 flags
#DEBUG = -g -DDEBUG #DEBUG = -g -DDEBUG
@ -39,7 +36,7 @@ SRC = $(wildcard *.cpp)
OBJ = $(SRC:.cpp=.o) OBJ = $(SRC:.cpp=.o)
#==================================================================== #====================================================================
all: ${EXE} all: ver ${EXE}
${EXE}: $(OBJ) ${EXE}: $(OBJ)
$(LINK) $(OFLAGS) $(OBJ) $(LIB) -o $@ $(LINK) $(OFLAGS) $(OBJ) $(LIB) -o $@
@ -51,7 +48,7 @@ tar:
rm -f ${ROOT}.tar; tar -czvf ${ROOT}.tar.gz *.cpp *.h Makefile README rm -f ${ROOT}.tar; tar -czvf ${ROOT}.tar.gz *.cpp *.h Makefile README
ver: ver:
@echo "#define VERSION `svn info|grep '^Revision'|cut -d: -f2`" > version.h @echo "#define VERSION `git log|grep '^commit'|wc -l`" > version.h
#==================================================================== #====================================================================
.f.o: .f.o:

67
tools/phonon/README
View File

@ -1,31 +1,48 @@
phana = post-processing program for the fix phonon command #-------------------------------------------------------------------------------
phana
#
This program reads the binary file created by fix_phonon and helps to
analyse the phonon related information.
#-------------------------------------------------------------------------------
1. Dependencies
The clapack library is needed to solve the eigen problems,
which could be downloaded from:
http://www.netlib.org/clapack/
This program reads the binary file created by the fix phonon command The tricubic library is also needed to do tricubic interpolations,
and helps to analyse the phonon related info. which could be obtained from:
http://orca.princeton.edu/francois/software/tricubic/
or
http://1drv.ms/1J2WFYk
The clapack library is needed to solve the eigen problems, which could The spglib is optionally needed, enabling one to evaluate the
be downloaded from: http://www.netlib.org/clapack/ 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/
The tricubic library is also needed to to tricubic interpolations, 2. Compilation
which could be obtained from: To compile the code, one needs therefore to install the above
http://orca.princeton.edu/francois/software/tricubic/ or libraries and set the paths correctly in the Makefile.
http://code.google.com/p/fix-phonon/downloads/list Once this is done, by typing
make
will yield the executable "phana".
The spglib (version 0.7.1) is optionally needed, enabling one to 3. Unit system
evaluate the phonon density of states or vibrational thermal The units of the output frequencies by this code is THz for
properties using only the irreducible q-points in the first Brillouin LAMMPS units "real", "si", "metal", and "cgs"; in these cases,
zone. the frequencies are $\nu$ instead of $\omega$.
To compile the code, one needs therefore to install the above 4. Updates
libraries and set the paths correctly in the Makefile. For updates of phana, please check:
http://nes.sjtu.edu.cn/english/research/software.htm
or
https://github.com/lingtikong/phana.git
The units of the output frequencies by this code is THz for LAMMPS 5. Bug report
units "real", "si", "metal", and "cgs"; in these cases, the If any bug found, please drop a line to: konglt(at)sjtu.edu.cn
frequencies are $\nu$ instead of $\omega$. #-------------------------------------------------------------------------------
Author: Ling-Ti Kong, konglt(at)sjtu.edu.cn
One is encouraged to visit http://code.google.com/p/fix-phonon/ to Oct 2015
check out the latest revision on fix-phonon and the post-processing
code.
Author: Ling-Ti Kong (konglt at sjtu.edu.cn)
Feb 2013

197
tools/phonon/dynmat.cpp
View File

@ -1,8 +1,7 @@
#include "dynmat.h" #include "dynmat.h"
#include "math.h" #include "math.h"
#include "version.h" #include "version.h"
#include "global.h"
#define MAXLINE 256
// to intialize the class // to intialize the class
DynMat::DynMat(int narg, char **arg) DynMat::DynMat(int narg, char **arg)
@ -14,6 +13,8 @@ DynMat::DynMat(int narg, char **arg)
DM_q = DM_all = NULL; DM_q = DM_all = NULL;
binfile = funit = dmfile = NULL; binfile = funit = dmfile = NULL;
attyp = NULL;
basis = NULL;
flag_reset_gamma = flag_skip = 0; flag_reset_gamma = flag_skip = 0;
// analyze the command line options // analyze the command line options
@ -74,21 +75,21 @@ DynMat::DynMat(int narg, char **arg)
npt = nx*ny*nz; npt = nx*ny*nz;
// display info related to the read file // display info related to the read file
printf("\n"); for (int i=0; i<80; i++) printf("="); printf("\n"); printf("\n"); for (int i = 0; i < 80; ++i) printf("="); printf("\n");
printf("Dynamical matrix is read from file: %s\n", binfile); printf("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("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("Number of atoms per unit cell : %d\n", nucell);
printf("System dimension : %d\n", sysdim); printf("System dimension : %d\n", sysdim);
printf("Boltzmann constant in used units : %g\n", boltz); 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){ 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); printf("Wrong values read from header of file: %s, please check the binary file!\n", binfile);
fclose(fp); exit(3); fclose(fp); exit(3);
} }
funit = new char[4]; funit = new char[4];
strcpy(funit, "THz"); strcpy(funit, "THz");
if (boltz == 1.){eml2f = 1.; delete funit; funit=new char[22]; strcpy(funit,"sqrt(epsilon/(m.sigma^2))");} 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 == 0.0019872067) eml2f = 3.256576161;
else if (boltz == 8.617343e-5) eml2f = 15.63312493; else if (boltz == 8.617343e-5) eml2f = 15.63312493;
else if (boltz == 1.3806504e-23) eml2f = 1.; else if (boltz == 1.3806504e-23) eml2f = 1.;
@ -100,9 +101,9 @@ DynMat::DynMat(int narg, char **arg)
} }
// now to allocate memory for DM // now to allocate memory for DM
memory = new Memory; memory = new Memory();
DM_all = memory->create(DM_all, npt, fftdim2, "DynMat:DM_all"); memory->create(DM_all, npt, fftdim2, "DynMat:DM_all");
DM_q = memory->create(DM_q, fftdim,fftdim,"DynMat:DM_q"); memory->create(DM_q, fftdim,fftdim,"DynMat:DM_q");
// read all dynamical matrix info into DM_all // read all dynamical matrix info into DM_all
if ( fread(DM_all[0], sizeof(doublecomplex), npt*fftdim2, fp) != size_t(npt*fftdim2)){ if ( fread(DM_all[0], sizeof(doublecomplex), npt*fftdim2, fp) != size_t(npt*fftdim2)){
@ -112,50 +113,38 @@ DynMat::DynMat(int narg, char **arg)
} }
// now try to read unit cell info from the binary file // now try to read unit cell info from the binary file
flag_latinfo = 0; memory->create(basis, nucell, sysdim, "DynMat:basis");
basis = memory->create(basis,nucell,sysdim,"DynMat:basis"); memory->create(attyp, nucell, "DynMat:attyp");
attyp = memory->create(attyp,nucell, "DynMat:attyp"); memory->create(M_inv_sqrt, nucell, "DynMat:M_inv_sqrt");
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(&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) flag_latinfo |= 2; 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) flag_latinfo |= 4; 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) flag_latinfo |= 8; 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) flag_mass_read = 1; if ( fread(&M_inv_sqrt[0], sizeof(double), nucell, fp) != nucell){printf("\nError while reading atomic masses from file: %s\n", binfile); fclose(fp); exit(3);}
fclose(fp); fclose(fp);
if ((flag_latinfo&15) == 15){ car2dir();
car2dir(flag_mass_read);
real2rec(); real2rec();
flag_latinfo = 1;
} else {
Tmeasure = 0.;
flag_latinfo = 0;
}
// initialize interpolation // initialize interpolation
interpolate = new Interpolate(nx,ny,nz,fftdim2,DM_all); interpolate = new Interpolate(nx,ny,nz,fftdim2,DM_all);
if (flag_reset_gamma) interpolate->reset_gamma(); if (flag_reset_gamma) interpolate->reset_gamma();
if ( flag_mass_read ){ // M_inv_sqrt info read, the data stored are force constant matrix instead of dynamical matrix. // Enforcing Austic Sum Rule
EnforceASR(); EnforceASR();
// get the dynamical matrix from force constant matrix: D = 1/M x Phi // get the dynamical matrix from force constant matrix: D = 1/M x Phi
for (int idq=0; idq< npt; idq++){ for (int idq = 0; idq < npt; ++idq){
int ndim =0; int ndim =0;
for (int idim=0; idim<fftdim; idim++){ for (int idim = 0; idim < fftdim; ++idim)
for (int jdim=0; jdim<fftdim; jdim++){ for (int jdim = 0; jdim < fftdim; ++jdim){
double inv_mass = M_inv_sqrt[idim/sysdim]*M_inv_sqrt[jdim/sysdim]; double inv_mass = M_inv_sqrt[idim/sysdim]*M_inv_sqrt[jdim/sysdim];
DM_all[idq][ndim].r *= inv_mass; DM_all[idq][ndim].r *= inv_mass;
DM_all[idq][ndim].i *= inv_mass; DM_all[idq][ndim].i *= inv_mass;
ndim++; ndim++;
} }
} }
}
}
// ask for the interpolation method // ask for the interpolation method
interpolate->set_method(); interpolate->set_method();
@ -191,7 +180,7 @@ void DynMat::writeDMq(double *q)
if (dmfile == NULL){ if (dmfile == NULL){
char str[MAXLINE], *ptr; char str[MAXLINE], *ptr;
printf("\n"); printf("\n");
while (1){ while ( 1 ){
printf("Please input the filename to output the DM at selected q: "); printf("Please input the filename to output the DM at selected q: ");
fgets(str,MAXLINE,stdin); fgets(str,MAXLINE,stdin);
ptr = strtok(str, " \r\t\n\f"); ptr = strtok(str, " \r\t\n\f");
@ -208,8 +197,8 @@ void DynMat::writeDMq(double *q)
} }
fprintf(fp,"# q = [%lg %lg %lg]\n", q[0], q[1], q[2]); fprintf(fp,"# q = [%lg %lg %lg]\n", q[0], q[1], q[2]);
for (int i=0; i<fftdim; i++){ for (int i = 0; i < fftdim; ++i){
for (int j=0; j<fftdim; j++) fprintf(fp,"%lg %lg\t", DM_q[i][j].r, DM_q[i][j].i); 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");
} }
fprintf(fp,"\n"); fprintf(fp,"\n");
@ -225,9 +214,9 @@ void DynMat::writeDMq(double *q, const double qr, FILE *fp)
fprintf(fp, "%lg %lg %lg %lg ", q[0], q[1], q[2], qr); fprintf(fp, "%lg %lg %lg %lg ", q[0], q[1], q[2], qr);
for (int i=0; i<fftdim; i++){ for (int i = 0; i < fftdim; ++i)
for (int j=0; j<fftdim; j++) fprintf(fp,"%lg %lg\t", DM_q[i][j].r, DM_q[i][j].i); 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");
return; return;
} }
@ -254,14 +243,14 @@ int DynMat::geteigen(double *egv, int flag)
liwork = 3 + 5*n; liwork = 3 + 5*n;
lda = n; lda = n;
work = memory->create(work, lwork, "geteigen:work"); memory->create(work, lwork, "geteigen:work");
rwork = memory->create(rwork, lrwork, "geteigen:rwork"); memory->create(rwork, lrwork, "geteigen:rwork");
iwork = memory->create(iwork, liwork, "geteigen:iwork"); memory->create(iwork, liwork, "geteigen:iwork");
zheevd_(&jobz, &uplo, &n, DM_q[0], &lda, w, work, &lwork, rwork, &lrwork, iwork, &liwork, &info); 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) // to get w instead of w^2; and convert w into v (THz hopefully)
for (int i=0; i<n; i++){ for (int i = 0; i < n; ++i){
if (w[i]>= 0.) w[i] = sqrt(w[i]); if (w[i]>= 0.) w[i] = sqrt(w[i]);
else w[i] = -sqrt(-w[i]); else w[i] = -sqrt(-w[i]);
@ -298,24 +287,17 @@ return;
/* ---------------------------------------------------------------------------- /* ----------------------------------------------------------------------------
* private method to convert the cartisan coordinate of basis into fractional * private method to convert the cartisan coordinate of basis into fractional
* ---------------------------------------------------------------------------- */ * ---------------------------------------------------------------------------- */
void DynMat::car2dir(int flag) void DynMat::car2dir()
{ {
if (!flag){ // in newer version, this is done in fix-phonon
for (int i=0; i<3; i++){
basevec[i] /= double(nx);
basevec[i+3] /= double(ny);
basevec[i+6] /= double(nz);
}
}
double mat[9]; double mat[9];
for (int idim=0; idim<9; idim++) mat[idim] = basevec[idim]; for (int idim = 0; idim < 9; ++idim) mat[idim] = basevec[idim];
GaussJordan(3, mat); GaussJordan(3, mat);
for (int i=0; i<nucell; i++){ for (int i = 0; i < nucell; ++i){
double x[3]; double x[3];
x[0] = x[1] = x[2] = 0.; 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++) x[idim] = basis[i][idim];
for (int idim=0; idim<sysdim; 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]; basis[i][idim] = x[0]*mat[idim] + x[1]*mat[3+idim] + x[2]*mat[6+idim];
} }
@ -330,7 +312,7 @@ void DynMat::EnforceASR()
char str[MAXLINE]; char str[MAXLINE];
int nasr = 20; int nasr = 20;
if (nucell <= 1) nasr = 1; if (nucell <= 1) nasr = 1;
printf("\n"); for (int i=0; i<80; i++) printf("="); printf("\n"); for (int i = 0; i < 80; ++i) printf("=");
// compute and display eigenvalues of Phi at gamma before ASR // compute and display eigenvalues of Phi at gamma before ASR
if (nucell > 100){ if (nucell > 100){
@ -339,11 +321,11 @@ void DynMat::EnforceASR()
} }
double egvs[fftdim]; double egvs[fftdim];
for (int i=0; i<fftdim; i++) for (int i = 0; i < fftdim; ++i)
for (int j=0; j<fftdim; j++) DM_q[i][j] = DM_all[0][i*fftdim+j]; for (int j = 0; j < fftdim; ++j) DM_q[i][j] = DM_all[0][i*fftdim+j];
geteigen(egvs, 0); geteigen(egvs, 0);
printf("\nEigenvalues of Phi at gamma before enforcing ASR:\n"); printf("\nEigenvalues of Phi at gamma before enforcing ASR:\n");
for (int i=0; i<fftdim; i++){ for (int i = 0; i < fftdim; ++i){
printf("%lg ", egvs[i]); printf("%lg ", egvs[i]);
if (i%10 == 9) printf("\n"); if (i%10 == 9) printf("\n");
if (i == 99){ printf("...... (%d more skipped)\n", fftdim-100); break;} if (i == 99){ printf("...... (%d more skipped)\n", fftdim-100); break;}
@ -355,20 +337,23 @@ void DynMat::EnforceASR()
fgets(str,MAXLINE,stdin); fgets(str,MAXLINE,stdin);
char *ptr = strtok(str," \t\n\r\f"); char *ptr = strtok(str," \t\n\r\f");
if (ptr) nasr = atoi(ptr); if (ptr) nasr = atoi(ptr);
if (nasr < 1){return; for (int i=0; i<80; i++) printf("="); printf("\n");} if (nasr < 1){
for (int i=0; i<80; i++) printf("="); printf("\n");
return;
}
for (int iit=0; iit<nasr; iit++){ for (int iit = 0; iit < nasr; ++iit){
// simple ASR; the resultant matrix might not be symmetric // simple ASR; the resultant matrix might not be symmetric
for (int a=0; a<sysdim; a++) for (int a = 0; a < sysdim; ++a)
for (int b=0; b<sysdim; b++){ for (int b = 0; b < sysdim; ++b){
for (int k=0; k<nucell; k++){ for (int k = 0; k < nucell; ++k){
double sum = 0.; double sum = 0.;
for (int kp=0; kp<nucell; kp++){ for (int kp = 0; kp < nucell; ++kp){
int idx = (k*sysdim+a)*fftdim+kp*sysdim+b; int idx = (k*sysdim+a)*fftdim+kp*sysdim+b;
sum += DM_all[0][idx].r; sum += DM_all[0][idx].r;
} }
sum /= double(nucell); sum /= double(nucell);
for (int kp=0; kp<nucell; kp++){ for (int kp = 0; kp < nucell; ++kp){
int idx = (k*sysdim+a)*fftdim+kp*sysdim+b; int idx = (k*sysdim+a)*fftdim+kp*sysdim+b;
DM_all[0][idx].r -= sum; DM_all[0][idx].r -= sum;
} }
@ -376,11 +361,11 @@ void DynMat::EnforceASR()
} }
// symmetrize // symmetrize
for (int k=0; k<nucell; k++) for (int k = 0; k < nucell; ++k)
for (int kp=k; kp<nucell; kp++){ for (int kp = k; kp < nucell; ++kp){
double csum = 0.; double csum = 0.;
for (int a=0; a<sysdim; a++) for (int a = 0; a < sysdim; ++a)
for (int b=0; b<sysdim; b++){ for (int b = 0; b < sysdim; ++b){
int idx = (k*sysdim+a)*fftdim+kp*sysdim+b; int idx = (k*sysdim+a)*fftdim+kp*sysdim+b;
int jdx = (kp*sysdim+b)*fftdim+k*sysdim+a; int jdx = (kp*sysdim+b)*fftdim+k*sysdim+a;
csum = (DM_all[0][idx].r + DM_all[0][jdx].r )*0.5; csum = (DM_all[0][idx].r + DM_all[0][jdx].r )*0.5;
@ -390,16 +375,16 @@ void DynMat::EnforceASR()
} }
// symmetric ASR // symmetric ASR
for (int a=0; a<sysdim; a++) for (int a = 0; a < sysdim; ++a)
for (int b=0; b<sysdim; b++){ for (int b = 0; b < sysdim; ++b){
for (int k=0; k<nucell; k++){ for (int k = 0; k < nucell; ++k){
double sum = 0.; double sum = 0.;
for (int kp=0; kp<nucell; kp++){ for (int kp = 0; kp < nucell; ++kp){
int idx = (k*sysdim+a)*fftdim+kp*sysdim+b; int idx = (k*sysdim+a)*fftdim+kp*sysdim+b;
sum += DM_all[0][idx].r; sum += DM_all[0][idx].r;
} }
sum /= double(nucell-k); sum /= double(nucell-k);
for (int kp=k; kp<nucell; kp++){ for (int kp = k; kp < nucell; ++kp){
int idx = (k*sysdim+a)*fftdim+kp*sysdim+b; int idx = (k*sysdim+a)*fftdim+kp*sysdim+b;
int jdx = (kp*sysdim+b)*fftdim+k*sysdim+a; int jdx = (kp*sysdim+b)*fftdim+k*sysdim+a;
DM_all[0][idx].r -= sum; DM_all[0][idx].r -= sum;
@ -409,16 +394,17 @@ void DynMat::EnforceASR()
} }
// compute and display eigenvalues of Phi at gamma after ASR // compute and display eigenvalues of Phi at gamma after ASR
for (int i=0; i<fftdim; i++) for (int i = 0; i < fftdim; ++i)
for (int j=0; j<fftdim; j++) DM_q[i][j] = DM_all[0][i*fftdim+j]; for (int j = 0; j < fftdim; ++j) DM_q[i][j] = DM_all[0][i*fftdim+j];
geteigen(egvs, 0); geteigen(egvs, 0);
printf("Eigenvalues of Phi at gamma after enforcing ASR:\n"); printf("Eigenvalues of Phi at gamma after enforcing ASR:\n");
for (int i=0; i<fftdim; i++){ for (int i = 0; i < fftdim; ++i){
printf("%lg ", egvs[i]); printf("%lg ", egvs[i]);
if (i%10 == 9) printf("\n"); if (i%10 == 9) printf("\n");
if (i == 99){ printf("...... (%d more skiped)", fftdim-100); break;} if (i == 99){ printf("...... (%d more skiped)", fftdim-100); break;}
} }
printf("\n"); for (int i=0; i<80; i++) printf("="); printf("\n\n"); printf("\n");
for (int i = 0; i < 80; ++i) printf("="); printf("\n\n");
return; return;
} }
@ -441,23 +427,23 @@ void DynMat::real2rec()
ibasevec[8] = basevec[0]*basevec[4] - basevec[1]*basevec[3]; ibasevec[8] = basevec[0]*basevec[4] - basevec[1]*basevec[3];
double vol = 0.; double vol = 0.;
for (int i=0; i<sysdim; i++) vol += ibasevec[i] * basevec[i]; for (int i = 0; i < sysdim; ++i) vol += ibasevec[i] * basevec[i];
vol = 8.*atan(1.)/vol; vol = 8.*atan(1.)/vol;
for (int i=0; i<9; i++) ibasevec[i] *= vol; for (int i = 0; i < 9; ++i) ibasevec[i] *= vol;
printf("\n"); for (int i=0; i<80; i++) printf("="); printf("\n"); for (int i = 0; i < 80; ++i) printf("=");
printf("\nBasis vectors of the unit cell in real space:"); printf("\nBasis vectors of the unit cell in real space:");
for (int i=0; i<sysdim; i++){ for (int i = 0; i < sysdim; ++i){
printf("\n A%d: ", i+1); printf("\n A%d: ", i+1);
for (int j=0; j<sysdim; j++) printf("%8.4f ", basevec[i*3+j]); for (int j = 0; j < sysdim; ++j) printf("%8.4f ", basevec[i*3+j]);
} }
printf("\nBasis vectors of the corresponding reciprocal cell:"); printf("\nBasis vectors of the corresponding reciprocal cell:");
for (int i=0; i<sysdim; i++){ for (int i = 0; i < sysdim; ++i){
printf("\n B%d: ", i+1); printf("\n B%d: ", i+1);
for (int j=0; j<sysdim; j++) printf("%8.4f ", ibasevec[i*3+j]); for (int j = 0; j < sysdim; ++j) printf("%8.4f ", ibasevec[i*3+j]);
} }
printf("\n"); for (int i=0; i<80; i++) printf("="); printf("\n"); printf("\n"); for (int i = 0; i < 80; ++i) printf("="); printf("\n");
return; return;
} }
@ -479,21 +465,21 @@ void DynMat::GaussJordan(int n, double *Mat)
indxr = new int[n]; indxr = new int[n];
ipiv = new int[n]; ipiv = new int[n];
for (i=0; i<n; i++) ipiv[i] = 0; for (i = 0; i < n; ++i) ipiv[i] = 0;
for (i=0; i<n; i++){ for (i = 0; i < n; ++i){
big = 0.; big = 0.;
for (j=0; j<n; j++){ for (j = 0; j < n; ++j){
if (ipiv[j] != 1){ if (ipiv[j] != 1){
for (k=0; k<n; k++){ for (k = 0; k < n; ++k){
if (ipiv[k] == 0){ if (ipiv[k] == 0){
idr = j*n+k; idr = j * n + k;
nmjk = abs(Mat[idr]); nmjk = abs(Mat[idr]);
if (nmjk >= big){ if (nmjk >= big){
big = nmjk; big = nmjk;
irow = j; irow = j;
icol = k; icol = k;
} }
}else if (ipiv[k]>1){ } else if (ipiv[k] > 1){
printf("DynMat: Singular matrix in double GaussJordan!\n"); exit(1); printf("DynMat: Singular matrix in double GaussJordan!\n"); exit(1);
} }
} }
@ -501,7 +487,7 @@ void DynMat::GaussJordan(int n, double *Mat)
} }
ipiv[icol] += 1; ipiv[icol] += 1;
if (irow != icol){ if (irow != icol){
for (l=0; l<n; l++){ for (l = 0; l < n; ++l){
idr = irow*n+l; idr = irow*n+l;
idc = icol*n+l; idc = icol*n+l;
dum = Mat[idr]; dum = Mat[idr];
@ -511,7 +497,7 @@ void DynMat::GaussJordan(int n, double *Mat)
} }
indxr[i] = irow; indxr[i] = irow;
indxc[i] = icol; indxc[i] = icol;
idr = icol*n+icol; idr = icol * n + icol;
if (Mat[idr] == 0.){ if (Mat[idr] == 0.){
printf("DynMat: Singular matrix in double GaussJordan!"); printf("DynMat: Singular matrix in double GaussJordan!");
exit(1); exit(1);
@ -520,24 +506,24 @@ void DynMat::GaussJordan(int n, double *Mat)
pivinv = 1./ Mat[idr]; pivinv = 1./ Mat[idr];
Mat[idr] = 1.; Mat[idr] = 1.;
idr = icol*n; idr = icol*n;
for (l=0; l<n; l++) Mat[idr+l] *= pivinv; for (l = 0; l < n; ++l) Mat[idr+l] *= pivinv;
for (ll=0; ll<n; ll++){ for (ll = 0; ll < n; ++ll){
if (ll != icol){ if (ll != icol){
idc = ll*n+icol; idc = ll * n + icol;
dum = Mat[idc]; dum = Mat[idc];
Mat[idc] = 0.; Mat[idc] = 0.;
idc -= icol; idc -= icol;
for (l=0; l<n; l++) Mat[idc+l] -= Mat[idr+l]*dum; for (l = 0; l < n; ++l) Mat[idc+l] -= Mat[idr+l]*dum;
} }
} }
} }
for (l=n-1; l>=0; l--){ for (l = n-1; l >= 0; --l){
int rl = indxr[l]; int rl = indxr[l];
int cl = indxc[l]; int cl = indxc[l];
if (rl != cl){ if (rl != cl){
for (k=0; k<n; k++){ for (k = 0; k < n; ++k){
idr = k*n+rl; idr = k * n + rl;
idc = k*n+cl; idc = k * n + cl;
dum = Mat[idr]; dum = Mat[idr];
Mat[idr] = Mat[idc]; Mat[idr] = Mat[idc];
Mat[idc] = dum; Mat[idc] = dum;
@ -547,6 +533,7 @@ void DynMat::GaussJordan(int n, double *Mat)
delete []indxr; delete []indxr;
delete []indxc; delete []indxc;
delete []ipiv; delete []ipiv;
return; return;
} }
@ -595,7 +582,7 @@ void DynMat::ShowVersion()
printf(" ( _ \\( )_( ) /__\\ ( \\( ) /__\\ \n"); printf(" ( _ \\( )_( ) /__\\ ( \\( ) /__\\ \n");
printf(" )___/ ) _ ( /(__)\\ ) ( /(__)\\ \n"); printf(" )___/ ) _ ( /(__)\\ ) ( /(__)\\ \n");
printf(" (__) (_) (_)(__)(__)(_)\\_)(__)(__)\n"); printf(" (__) (_) (_)(__)(__)(_)\\_)(__)(__)\n");
printf("\nPHonon ANAlyzer for Fix-Phonon, version 2.%d, compiled on %s.\n", VERSION, __DATE__); printf("\nPHonon ANAlyzer for Fix-Phonon, version 2.%02d, compiled on %s.\n", VERSION, __DATE__);
return; return;
} }

2
tools/phonon/dynmat.h
View File

@ -56,7 +56,7 @@ private:
doublecomplex **DM_all; doublecomplex **DM_all;
void car2dir(int); // to convert basis from cartisian coordinate into factional. void car2dir(); // to convert basis from cartisian coordinate into factional.
void real2rec(); void real2rec();
void GaussJordan(int, double *); void GaussJordan(int, double *);

47
tools/phonon/green.cpp
View File

@ -4,8 +4,7 @@
#include <math.h> #include <math.h>
#include "green.h" #include "green.h"
#include <complex> #include <complex>
#include "global.h"
#define MAXLINE 256
/******************************************************************************* /*******************************************************************************
* The class of Green is designed to evaluate the LDOS via the Green's Function * The class of Green is designed to evaluate the LDOS via the Green's Function
@ -45,7 +44,7 @@ Green::Green(const int ntm, const int sdim, const int niter, const double min, c
H = Hessian; iatom = itm; H = Hessian; iatom = itm;
ldos = lpdos; ldos = lpdos;
memory = new Memory; memory = new Memory();
if (natom < 1 || iatom < 0 || iatom >= natom){ if (natom < 1 || iatom < 0 || iatom >= natom){
printf("\nError: Wrong number of total atoms or wrong index of interested atom!\n"); printf("\nError: Wrong number of total atoms or wrong index of interested atom!\n");
return; return;
@ -58,9 +57,9 @@ Green::Green(const int ntm, const int sdim, const int niter, const double min, c
// initialize variables and allocate local memories // initialize variables and allocate local memories
dw = (wmax - wmin)/double(nw-1); dw = (wmax - wmin)/double(nw-1);
alpha = memory->create(alpha, sysdim,nit, "Green_Green:alpha"); memory->create(alpha, sysdim,nit, "Green_Green:alpha");
beta = memory->create(beta, sysdim,nit+1,"Green_Green:beta"); memory->create(beta, sysdim,nit+1,"Green_Green:beta");
//ldos = memory->create(ldos, nw,sysdim, "Green_Green:ldos"); //memory->create(ldos, nw,sysdim, "Green_Green:ldos");
// use Lanczos algorithm to diagonalize the Hessian // use Lanczos algorithm to diagonalize the Hessian
Lanczos(); Lanczos();
@ -100,35 +99,35 @@ void Green::Lanczos()
int ipos = iatom*sysdim; int ipos = iatom*sysdim;
// Loop over dimension // Loop over dimension
for (int idim=0; idim<sysdim; idim++){ for (int idim = 0; idim < sysdim; ++idim){
beta[idim][0] = 0.; beta[idim][0] = 0.;
for (int i=0; i<ndim; i++){vp[i] = v[i] = 0.;} for (int i = 0; i < ndim; ++i) vp[i] = v[i] = 0.;
v[ipos+idim] = 1.; v[ipos+idim] = 1.;
// Loop on fraction levels // Loop on fraction levels
for (int i=0; i<nit; i++){ for (int i = 0; i < nit; ++i){
double sum_a = 0.; double sum_a = 0.;
for (int j=0; j<ndim; j++){ for (int j = 0; j < ndim; ++j){
double sumHv = 0.; double sumHv = 0.;
for (int k=0; k<ndim; k++) sumHv += H[j][k]*v[k]; for (int k = 0; k < ndim; ++k) sumHv += H[j][k]*v[k];
w[j] = sumHv - beta[idim][i]*vp[j]; w[j] = sumHv - beta[idim][i]*vp[j];
sum_a += w[j]*v[j]; sum_a += w[j]*v[j];
} }
alpha[idim][i] = sum_a; alpha[idim][i] = sum_a;
for (int k=0; k<ndim; k++) w[k] -= alpha[idim][i]*v[k]; for (int k = 0; k < ndim; ++k) w[k] -= alpha[idim][i]*v[k];
double gamma = 0.; double gamma = 0.;
for (int k=0; k<ndim; k++) gamma += w[k]*v[k]; for (int k = 0; k < ndim; ++k) gamma += w[k]*v[k];
for (int k=0; k<ndim; k++) w[k] -= gamma*v[k]; for (int k = 0; k < ndim; ++k) w[k] -= gamma*v[k];
double sum_b = 0.; double sum_b = 0.;
for (int k=0; k<ndim; k++) sum_b += w[k]*w[k]; for (int k = 0; k < ndim; ++k) sum_b += w[k]*w[k];
beta[idim][i+1] = sqrt(sum_b); beta[idim][i+1] = sqrt(sum_b);
ptr = vp; vp = v; v = ptr; ptr = vp; vp = v; v = ptr;
double tmp = 1./beta[idim][i+1]; double tmp = 1./beta[idim][i+1];
for (int k=0; k<ndim; k++) v[k] = w[k]*tmp; for (int k = 0; k < ndim; ++k) v[k] = w[k] * tmp;
} }
} }
@ -154,10 +153,10 @@ void Green::Recursion()
xmax = new double [sysdim]; xmax = new double [sysdim];
int nave = nit/4; int nave = nit/4;
for (int idim=0; idim<sysdim; idim++){ for (int idim = 0; idim < sysdim; ++idim){
alpha_inf[idim] = beta_inf[idim] = 0.; alpha_inf[idim] = beta_inf[idim] = 0.;
for (int i= nit-nave; i<nit; i++){ for (int i = nit-nave; i < nit; ++i){
alpha_inf[idim] += alpha[idim][i]; alpha_inf[idim] += alpha[idim][i];
beta_inf[idim] += beta[idim][i+1]; beta_inf[idim] += beta[idim][i+1];
} }
@ -173,11 +172,11 @@ void Green::Recursion()
double sr, si; double sr, si;
double w = wmin; double w = wmin;
for (int i=0; i<nw; i++){ for (int i = 0; i < nw; ++i){
double a = w*w, ax, bx; double a = w*w, ax, bx;
Z = std::complex<double>(w*w, epson); Z = std::complex<double>(w*w, epson);
for (int idim=0; idim<sysdim; idim++){ for (int idim = 0; idim < sysdim; ++idim){
double two_b = 2.*beta_inf[idim]*beta_inf[idim]; double two_b = 2.*beta_inf[idim]*beta_inf[idim];
double rtwob = 1./two_b; double rtwob = 1./two_b;
@ -201,7 +200,7 @@ void Green::Recursion()
si = epson * rtwob + bx; si = epson * rtwob + bx;
rec_x = std::complex<double> (sr, si); rec_x = std::complex<double> (sr, si);
for (int j=0; j<nit; j++){ for (int j = 0; j < nit; ++j){
rec_x_inv = Z - alpha[idim][nit-j-1] - beta[idim][nit-j]*beta[idim][nit-j]*rec_x; rec_x_inv = Z - alpha[idim][nit-j-1] - beta[idim][nit-j]*beta[idim][nit-j]*rec_x;
rec_x = 1./rec_x_inv; rec_x = 1./rec_x_inv;
} }
@ -229,13 +228,13 @@ void Green::recursion()
double w = wmin; double w = wmin;
for (int i=0; i<nw; i++){ for (int i = 0; i < nw; ++i){
Z = std::complex<double>(w*w, epson); Z = std::complex<double>(w*w, epson);
for (int idim=0; idim<sysdim; idim++){ for (int idim = 0; idim < sysdim; ++idim){
rec_x = std::complex<double>(0.,0.); rec_x = std::complex<double>(0.,0.);
for (int j=0; j<nit; j++){ for (int j = 0; j < nit; ++j){
rec_x_inv = Z - alpha[idim][nit-j-1] - beta[idim][nit-j]*beta[idim][nit-j]*rec_x; rec_x_inv = Z - alpha[idim][nit-j-1] - beta[idim][nit-j]*beta[idim][nit-j]*rec_x;
rec_x = 1./rec_x_inv; rec_x = 1./rec_x_inv;
} }

51
tools/phonon/interpolate.cpp
View File

@ -1,9 +1,6 @@
#include "interpolate.h" #include "interpolate.h"
#include "math.h" #include "math.h"
#include "global.h"
#define MAXLINE 256
#define MIN(a,b) ((a)>(b)?(b):(a))
#define MAX(a,b) ((a)>(b)?(a):(b))
/* ---------------------------------------------------------------------------- /* ----------------------------------------------------------------------------
* Constructor used to get info from caller, and prepare other necessary data * Constructor used to get info from caller, and prepare other necessary data
@ -15,7 +12,7 @@ Interpolate::Interpolate(int nx, int ny, int nz, int ndm, doublecomplex **DM)
Nz = nz; Nz = nz;
Npt = Nx*Ny*Nz; Npt = Nx*Ny*Nz;
ndim = ndm; ndim = ndm;
memory = new Memory; memory = new Memory();
which = UseGamma = 0; which = UseGamma = 0;
@ -33,13 +30,13 @@ void Interpolate::tricubic_init()
{ {
// prepare necessary data for tricubic // prepare necessary data for tricubic
if (flag_allocated_dfs == 0){ if (flag_allocated_dfs == 0){
Dfdx = memory->create(Dfdx, Npt, ndim, "Interpolate_Interpolate:Dfdx"); memory->create(Dfdx, Npt, ndim, "Interpolate_Interpolate:Dfdx");
Dfdy = memory->create(Dfdy, Npt, ndim, "Interpolate_Interpolate:Dfdy"); memory->create(Dfdy, Npt, ndim, "Interpolate_Interpolate:Dfdy");
Dfdz = memory->create(Dfdz, Npt, ndim, "Interpolate_Interpolate:Dfdz"); memory->create(Dfdz, Npt, ndim, "Interpolate_Interpolate:Dfdz");
D2fdxdy = memory->create(D2fdxdy, Npt, ndim, "Interpolate_Interpolate:D2fdxdy"); memory->create(D2fdxdy, Npt, ndim, "Interpolate_Interpolate:D2fdxdy");
D2fdxdz = memory->create(D2fdxdz, Npt, ndim, "Interpolate_Interpolate:D2fdxdz"); memory->create(D2fdxdz, Npt, ndim, "Interpolate_Interpolate:D2fdxdz");
D2fdydz = memory->create(D2fdydz, Npt, ndim, "Interpolate_Interpolate:D2fdydz"); memory->create(D2fdydz, Npt, ndim, "Interpolate_Interpolate:D2fdydz");
D3fdxdydz = memory->create(D3fdxdydz, Npt, ndim, "Interpolate_Interpolate:D2fdxdydz"); memory->create(D3fdxdydz, Npt, ndim, "Interpolate_Interpolate:D2fdxdydz");
flag_allocated_dfs = 1; flag_allocated_dfs = 1;
} }
@ -47,9 +44,9 @@ void Interpolate::tricubic_init()
// get the derivatives // get the derivatives
int n=0; int n=0;
const double half = 0.5, one4 = 0.25, one8 = 0.125; const double half = 0.5, one4 = 0.25, one8 = 0.125;
for (int ii=0; ii<Nx; ii++) for (int ii = 0; ii < Nx; ++ii)
for (int jj=0; jj<Ny; jj++) for (int jj = 0; jj < Ny; ++jj)
for (int kk=0; kk<Nz; kk++){ for (int kk = 0; kk < Nz; ++kk){
int ip = (ii+1)%Nx, jp = (jj+1)%Ny, kp = (kk+1)%Nz; 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 im = (ii-1+Nx)%Nx, jm = (jj-1+Ny)%Ny, km = (kk-1+Nz)%Nz;
@ -127,8 +124,8 @@ void Interpolate::tricubic(double *qin, doublecomplex *DMq)
{ {
// qin should be in unit of 2*pi/L // qin should be in unit of 2*pi/L
double q[3]; double q[3];
for (int i=0; i<3; i++) q[i] = qin[i]; for (int i = 0; i < 3; ++i) q[i] = qin[i];
for (int i=0; i<3; i++){ for (int i = 0; i < 3; ++i){
while (q[i] < 0.) q[i] += 1.; while (q[i] < 0.) q[i] += 1.;
while (q[i] >= 1.) q[i] -= 1.; while (q[i] >= 1.) q[i] -= 1.;
} }
@ -150,8 +147,8 @@ void Interpolate::tricubic(double *qin, doublecomplex *DMq)
vidx[7] = (ixp*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 i=0; i<8; i++) if (vidx[i] == 0) UseGamma = 1;
for (int idim=0; idim<ndim; idim++){ for (int idim = 0; idim < ndim; ++idim){
for (int i=0; i<8; i++){ for (int i = 0; i < 8; ++i){
f[i] = data[vidx[i]][idim].r; f[i] = data[vidx[i]][idim].r;
dfdx[i] = Dfdx[vidx[i]][idim].r; dfdx[i] = Dfdx[vidx[i]][idim].r;
dfdy[i] = Dfdy[vidx[i]][idim].r; dfdy[i] = Dfdy[vidx[i]][idim].r;
@ -164,7 +161,7 @@ void Interpolate::tricubic(double *qin, doublecomplex *DMq)
tricubic_get_coeff(&a[0],&f[0],&dfdx[0],&dfdy[0],&dfdz[0],&d2fdxdy[0],&d2fdxdz[0],&d2fdydz[0],&d3fdxdydz[0]); 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); DMq[idim].r = tricubic_eval(&a[0],x,y,z);
for (int i=0; i<8; i++){ for (int i = 0; i < 8; ++i){
f[i] = data[vidx[i]][idim].i; f[i] = data[vidx[i]][idim].i;
dfdx[i] = Dfdx[vidx[i]][idim].i; dfdx[i] = Dfdx[vidx[i]][idim].i;
dfdy[i] = Dfdy[vidx[i]][idim].i; dfdy[i] = Dfdy[vidx[i]][idim].i;
@ -177,6 +174,8 @@ void Interpolate::tricubic(double *qin, doublecomplex *DMq)
tricubic_get_coeff(&a[0],&f[0],&dfdx[0],&dfdy[0],&dfdz[0],&d2fdxdy[0],&d2fdxdz[0],&d2fdydz[0],&d3fdxdydz[0]); 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); DMq[idim].i = tricubic_eval(&a[0],x,y,z);
} }
return;
} }
/* ---------------------------------------------------------------------------- /* ----------------------------------------------------------------------------
@ -188,8 +187,8 @@ void Interpolate::trilinear(double *qin, doublecomplex *DMq)
{ {
// rescale q[i] into [0 1) // rescale q[i] into [0 1)
double q[3]; double q[3];
for (int i=0; i<3; i++) q[i] = qin[i]; for (int i = 0; i < 3; ++i) q[i] = qin[i];
for (int i=0; i<3; i++){ for (int i = 0; i < 3; ++i){
while (q[i] < 0.) q[i] += 1.; while (q[i] < 0.) q[i] += 1.;
while (q[i] >= 1.) q[i] -= 1.; while (q[i] >= 1.) q[i] -= 1.;
} }
@ -220,7 +219,7 @@ void Interpolate::trilinear(double *qin, doublecomplex *DMq)
vidx[5] = ((ix*Ny)+iyp)*Nz + izp; vidx[5] = ((ix*Ny)+iyp)*Nz + izp;
vidx[6] = ((ixp*Ny)+iyp)*Nz + iz; vidx[6] = ((ixp*Ny)+iyp)*Nz + iz;
vidx[7] = ((ixp*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 i = 0; i < 8; ++i) if (vidx[i] == 0) UseGamma = 1;
double fac[8]; double fac[8];
fac[0] = (1.-x)*(1.-y)*(1.-z); fac[0] = (1.-x)*(1.-y)*(1.-z);
@ -233,10 +232,10 @@ void Interpolate::trilinear(double *qin, doublecomplex *DMq)
fac[7] = x*y*z; fac[7] = x*y*z;
// now to do the interpolation // now to do the interpolation
for (int idim=0; idim<ndim; idim++){ for (int idim = 0; idim < ndim; ++idim){
DMq[idim].r = 0.; DMq[idim].r = 0.;
DMq[idim].i = 0.; DMq[idim].i = 0.;
for (int i=0; i<8; i++){ for (int i = 0; i < 8; ++i){
DMq[idim].r += data[vidx[i]][idim].r*fac[i]; DMq[idim].r += data[vidx[i]][idim].r*fac[i];
DMq[idim].i += data[vidx[i]][idim].i*fac[i]; DMq[idim].i += data[vidx[i]][idim].i*fac[i];
} }
@ -268,7 +267,7 @@ void Interpolate::set_method()
printf("\n");for(int i=0; i<80; i++) printf("="); printf("\n");for(int i=0; i<80; i++) printf("=");
printf("\nWhich interpolation method would you like to use?\n"); printf("\nWhich interpolation method would you like to use?\n");
printf(" 1. Tricubic;\n 2. Trilinear;\n"); printf(" 1. Tricubic;\n 2. Trilinear;\n");
printf("Your choice[1]: "); printf("Your choice [1]: ");
fgets(str,MAXLINE,stdin); fgets(str,MAXLINE,stdin);
char *ptr = strtok(str," \t\n\r\f"); char *ptr = strtok(str," \t\n\r\f");
if (ptr) im = atoi(ptr); if (ptr) im = atoi(ptr);

481
tools/phonon/phonon.cpp
View File

@ -3,6 +3,7 @@
#include "phonon.h" #include "phonon.h"
#include "green.h" #include "green.h"
#include "timer.h" #include "timer.h"
#include "global.h"
#ifdef UseSPG #ifdef UseSPG
extern "C"{ extern "C"{
@ -10,10 +11,6 @@ extern "C"{
} }
#endif #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 * Class Phonon is the main driver to calculate phonon DOS, phonon
* dispersion curve and some other things. * dispersion curve and some other things.
@ -42,7 +39,10 @@ Phonon::Phonon(DynMat *dm)
// display the menu // display the menu
char str[MAXLINE]; char str[MAXLINE];
while ( 1 ){ while ( 1 ){
printf("\n"); for (int i=0; i<37;i++) printf("="); printf(" Menu "); for (int i=0; i<37;i++) printf("="); printf("\n"); printf("\n");
for (int i = 0; i < 37; ++i) printf("=");
printf(" Menu ");
for (int i = 0; i < 37; ++i) printf("="); printf("\n");
printf(" 1. Phonon DOS evaluation;\n"); printf(" 1. Phonon DOS evaluation;\n");
printf(" 2. Phonon dispersion curves;\n"); printf(" 2. Phonon dispersion curves;\n");
printf(" 3. Dynamical matrix at arbitrary q;\n"); printf(" 3. Dynamical matrix at arbitrary q;\n");
@ -52,14 +52,15 @@ Phonon::Phonon(DynMat *dm)
printf(" 7. Local phonon DOS from eigenvectors;\n"); printf(" 7. Local phonon DOS from eigenvectors;\n");
printf(" 8. Local phonon DOS by RSGF method;\n"); printf(" 8. Local phonon DOS by RSGF method;\n");
printf(" 9. Freqs and eigenvectors at arbitrary q;\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(" -1. Reset the interpolation method;\n");
printf(" 0. Exit.\n"); printf(" 0. Exit.\n");
// read user choice // read user choice
int job = 0; int job = 0;
printf("Your choice[0]: "); printf("Your choice [0]: ");
if (count_words(fgets(str,MAXLINE,stdin)) > 0) job = atoi(strtok(str," \t\n\r\f")); if (count_words(fgets(str,MAXLINE,stdin)) > 0) job = atoi(strtok(str," \t\n\r\f"));
printf("\nYour selection: %d\n", job); 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 // now to do the job according to user's choice
if (job == 1) pdos(); if (job == 1) pdos();
@ -71,6 +72,7 @@ Phonon::Phonon(DynMat *dm)
else if (job == 7) ldos_egv(); else if (job == 7) ldos_egv();
else if (job == 8) ldos_rsgf(); else if (job == 8) ldos_rsgf();
else if (job == 9) vecanyq(); else if (job == 9) vecanyq();
else if (job ==10) ShowCell();
else if (job ==-1) dynmat->reset_interp_method(); else if (job ==-1) dynmat->reset_interp_method();
else break; else break;
} }
@ -113,12 +115,11 @@ void Phonon::pdos()
char str[MAXLINE]; char str[MAXLINE];
fmin = fmax = eigs[0][0]; fmin = fmax = eigs[0][0];
for (int iq=0; iq<nq; iq++){ for (int iq = 0; iq < nq; ++iq)
for (int j=0; j<ndim; j++){ for (int j = 0; j < ndim; ++j){
fmin = MIN(fmin, eigs[iq][j]); fmin = MIN(fmin, eigs[iq][j]);
fmax = MAX(fmax, eigs[iq][j]); fmax = MAX(fmax, eigs[iq][j]);
} }
}
// Now to ask for the output frequency range // Now to ask for the output frequency range
printf("\nThe frequency range of all q-points are: [%g %g]\n", fmin, fmax); printf("\nThe frequency range of all q-points are: [%g %g]\n", fmin, fmax);
@ -139,14 +140,14 @@ void Phonon::pdos()
df = (fmax-fmin)/double(ndos-1); df = (fmax-fmin)/double(ndos-1);
rdf = 1./df; rdf = 1./df;
memory->destroy(dos); memory->destroy(dos);
dos = memory->create(dos, ndos, "pdos:dos"); memory->create(dos, ndos, "pdos:dos");
for (int i=0; i<ndos; i++) dos[i] = 0.; for (int i=0; i<ndos; i++) dos[i] = 0.;
// now to calculate the DOS // now to calculate the DOS
double offset = fmin-0.5*df; double offset = fmin-0.5*df;
for (int iq=0; iq<nq; iq++){ for (int iq = 0; iq < nq; ++iq){
if (wt[iq] > 0.){ if (wt[iq] > 0.){
for (int j=0; j<ndim; j++){ for (int j = 0; j < ndim; ++j){
int idx = int((eigs[iq][j]-offset)*rdf); int idx = int((eigs[iq][j]-offset)*rdf);
if (idx>=0 && idx<ndos) dos[idx] += wt[iq]; if (idx>=0 && idx<ndos) dos[idx] += wt[iq];
} }
@ -220,17 +221,20 @@ void Phonon::writeLDOS()
printf("The phonon LDOSs will be written to file(s) : pldos_?.dat\n\n"); printf("The phonon LDOSs will be written to file(s) : pldos_?.dat\n\n");
const double one3 = 1./double(sysdim); const double one3 = 1./double(sysdim);
char str[MAXLINE]; char str[MAXLINE];
for (int ilocal=0; ilocal<nlocal; ilocal++){ for (int ilocal = 0; ilocal < nlocal; ++ilocal){
sprintf(str,"pldos_%d.dat", locals[ilocal]); sprintf(str,"pldos_%d.dat", locals[ilocal]);
char *fname = strtok(str," \t\n\r\f"); char *fname = strtok(str," \t\n\r\f");
FILE *fp = fopen(fname, "w"); fname = NULL; FILE *fp = fopen(fname, "w"); fname = NULL;
fprintf(fp,"#freq xDOS yDOS zDOS total\n"); fprintf(fp,"#freq xDOS yDOS zDOS total\n");
double freq = fmin; double freq = fmin;
for (int i=0; i<ndos; i++){ for (int i = 0; i < ndos; ++i){
fprintf(fp,"%lg", freq); fprintf(fp,"%lg", freq);
double total = 0.; double total = 0.;
for (int idim=0; idim<sysdim; idim++) {fprintf(fp," %lg",ldos[ilocal][i][idim]); total += ldos[ilocal][i][idim];} for (int idim = 0; idim < sysdim; ++idim){
fprintf(fp," %lg",ldos[ilocal][i][idim]);
total += ldos[ilocal][i][idim];
}
fprintf(fp," %lg\n", total*one3); fprintf(fp," %lg\n", total*one3);
freq += df; freq += df;
} }
@ -250,21 +254,21 @@ void Phonon::ldos_rsgf()
const double tpi = 8.*atan(1.); const double tpi = 8.*atan(1.);
double **Hessian, scale; double **Hessian, scale;
scale = dynmat->eml2f*tpi; scale *= scale; scale = dynmat->eml2f*tpi; scale *= scale;
Hessian = memory->create(Hessian, ndim, ndim, "phonon_ldos:Hessian"); memory->create(Hessian, ndim, ndim, "phonon_ldos:Hessian");
double q0[3]; double q0[3];
q0[0] = q0[1] = q0[2] = 0.; q0[0] = q0[1] = q0[2] = 0.;
dynmat->getDMq(q0); dynmat->getDMq(q0);
for (int i=0; i<ndim; i++) for (int i = 0; i < ndim; ++i)
for (int j=0; j<ndim; j++) Hessian[i][j] = dynmat->DM_q[i][j].r*scale; for (int j = 0; j < ndim; ++j) Hessian[i][j] = dynmat->DM_q[i][j].r*scale;
if (ndim < 300){ if (ndim < 300){
double *egvs = new double [ndim]; double *egvs = new double [ndim];
dynmat->geteigen(egvs, 0); dynmat->geteigen(egvs, 0);
fmin = fmax = egvs[0]; fmin = fmax = egvs[0];
for (int i=1; i<ndim; i++){fmin = MIN(fmin, egvs[i]); fmax = MAX(fmax, egvs[i]);} for (int i = 1; i < ndim; ++i){fmin = MIN(fmin, egvs[i]); fmax = MAX(fmax, egvs[i]);}
delete []egvs; delete []egvs;
} else { } else {
fmin = 0.; fmax = 20.; fmin = 0.; fmax = 20.;
@ -284,24 +288,29 @@ void Phonon::ldos_rsgf()
if (nr < 1){ if (nr < 1){
istr = iend = ik; istr = iend = ik;
iinc = 1; iinc = 1;
} else if (nr == 1) { } else if (nr == 1) {
char *ptr = strtok(str," \t\n\r\f"); char *ptr = strtok(str," \t\n\r\f");
if (strcmp(ptr,"q") == 0) break; if (strcmp(ptr,"q") == 0) break;
ik = atoi(ptr); ik = atoi(ptr);
if (ik < 0 || ik >= dynmat->nucell) break; ik = MAX(0, MIN(ik, dynmat->nucell-1));
istr = iend = ik; istr = iend = ik;
iinc = 1; iinc = 1;
} else if (nr == 2) { } else if (nr == 2) {
istr = atoi(strtok(str," \t\n\r\f")); istr = atoi(strtok(str," \t\n\r\f"));
iend = atoi(strtok(NULL," \t\n\r\f")); iend = atoi(strtok(NULL," \t\n\r\f"));
iinc = 1; iinc = 1;
if (istr < 0||iend >= dynmat->nucell||istr > iend) break; istr = MAX(0, MIN(istr, dynmat->nucell-1));
iend = MAX(0, MIN(iend, dynmat->nucell-1));
} else if (nr >= 3) { } else if (nr >= 3) {
istr = atoi(strtok(str," \t\n\r\f")); istr = atoi(strtok(str," \t\n\r\f"));
iend = atoi(strtok(NULL," \t\n\r\f")); iend = atoi(strtok(NULL," \t\n\r\f"));
iinc = atoi(strtok(NULL," \t\n\r\f")); iinc = atoi(strtok(NULL," \t\n\r\f"));
if (istr<0 || iend >= dynmat->nucell || istr > iend || iinc<1) break; istr = MAX(0, MIN(istr, dynmat->nucell-1));
iend = MAX(0, MIN(iend, dynmat->nucell-1));
} }
printf("Please input the frequency range to evaluate LDOS [%g %g]: ", fmin, fmax); printf("Please input the frequency range to evaluate LDOS [%g %g]: ", fmin, fmax);
@ -332,7 +341,7 @@ void Phonon::ldos_rsgf()
ldos = memory->create(ldos,nlocal,ndos,dynmat->sysdim,"ldos_rsgf:ldos"); ldos = memory->create(ldos,nlocal,ndos,dynmat->sysdim,"ldos_rsgf:ldos");
memory->destroy(locals); memory->destroy(locals);
locals = memory->create(locals, nlocal, "ldos_rsgf:locals"); memory->create(locals, nlocal, "ldos_rsgf:locals");
df = (fmax-fmin)/double(ndos-1); df = (fmax-fmin)/double(ndos-1);
rdf = 1./df; rdf = 1./df;
@ -366,104 +375,6 @@ void Phonon::ldos_rsgf()
return; return;
} }
/*------------------------------------------------------------------------------
* Private method to evaluate the phonon dispersion curves
*----------------------------------------------------------------------------*/
void Phonon::pdisp()
{
// ask the output file name and write the header.
char str[MAXLINE];
printf("Please input the filename to output the dispersion data [pdisp.dat]:");
if (count_words(fgets(str,MAXLINE,stdin)) < 1) strcpy(str, "pdisp.dat");
char *ptr = strtok(str," \t\n\r\f");
char *fname = new char[strlen(ptr)+1];
strcpy(fname,ptr);
FILE *fp = fopen(fname, "w");
fprintf(fp,"# q qr freq\n");
fprintf(fp,"# 2pi/L 2pi/L %s\n", dynmat->funit);
// to store the nodes of the dispersion curve
std::vector<double> nodes; nodes.clear();
// now the calculate the dispersion curve
double qstr[3], qend[3], q[3], qinc[3], qr=0., dq;
int nq = MAX(MAX(dynmat->nx,dynmat->ny),dynmat->nz)/2+1;
qend[0] = qend[1] = qend[2] = 0.;
double *egvs = new double [ndim];
while (1){
for (int i=0; i<3; i++) qstr[i] = qend[i];
int quit = 0;
printf("\nPlease input the start q-point in unit of B1->B3, q to exit [%g %g %g]: ", qstr[0], qstr[1], qstr[2]);
int n = count_words(fgets(str,MAXLINE,stdin));
ptr = strtok(str," \t\n\r\f");
if ((n == 1) && (strcmp(ptr,"q") == 0)) break;
else if (n >= 3){
qstr[0] = atof(ptr);
qstr[1] = atof(strtok(NULL," \t\n\r\f"));
qstr[2] = atof(strtok(NULL," \t\n\r\f"));
}
do printf("Please input the end q-point in unit of B1->B3: ");
while (count_words(fgets(str,MAXLINE,stdin)) < 3);
qend[0] = atof(strtok(str," \t\n\r\f"));
qend[1] = atof(strtok(NULL," \t\n\r\f"));
qend[2] = atof(strtok(NULL," \t\n\r\f"));
printf("Please input the # of points along the line [%d]: ", nq);
if (count_words(fgets(str,MAXLINE,stdin)) > 0) nq = atoi(strtok(str," \t\n\r\f"));
nq = MAX(nq,2);
for (int i=0; i<3; i++) qinc[i] = (qend[i]-qstr[i])/double(nq-1);
dq = sqrt(qinc[0]*qinc[0]+qinc[1]*qinc[1]+qinc[2]*qinc[2]);
nodes.push_back(qr);
for (int i=0; i<3; i++) q[i] = qstr[i];
for (int ii=0; ii<nq; ii++){
double wii = 1.;
dynmat->getDMq(q, &wii);
if (wii > 0.){
dynmat->geteigen(egvs, 0);
fprintf(fp,"%lg %lg %lg %lg ", q[0], q[1], q[2], qr);
for (int i=0; i<ndim; i++) fprintf(fp," %lg", egvs[i]);
}
fprintf(fp,"\n");
for (int i=0; i<3; i++) q[i] += qinc[i];
qr += dq;
}
qr -= dq;
}
if (qr > 0.) nodes.push_back(qr);
fclose(fp);
delete []egvs;
// write the gnuplot script which helps to visualize the result
int nnd = nodes.size();
if (nnd > 1){
fp = fopen("pdisp.gnuplot", "w");
fprintf(fp,"set term post enha colo 20\nset out %cpdisp.eps%c\n\n",char(34),char(34));
fprintf(fp,"set xlabel %cq%c\n",char(34),char(34));
fprintf(fp,"set ylabel %cfrequency (THz)%c\n\n",char(34),char(34));
fprintf(fp,"set xrange [0:%lg]\nset yrange [0:*]\n\n", nodes[nnd-1]);
fprintf(fp,"set grid xtics\n");
fprintf(fp,"# {/Symbol G} will give you letter gamma in the label\nset xtics (");
for (int i=0; i<nnd-1; i++) fprintf(fp,"%c%c %lg, ", char(34),char(34),nodes[i]);
fprintf(fp,"%c%c %lg)\n\n",char(34),char(34),nodes[nnd-1]);
fprintf(fp,"unset key\n\n");
fprintf(fp,"plot %c%s%c u 4:5 w l lt 1",char(34),fname,char(34));
for (int i=1; i<ndim; i++) fprintf(fp,",\\\n%c%c u 4:%d w l lt 1",char(34),char(34),i+5);
fclose(fp);
}
delete []fname;
nodes.clear();
return;
}
/* ---------------------------------------------------------------------------- /* ----------------------------------------------------------------------------
* Private method to write out the dynamical matrix at selected q * Private method to write out the dynamical matrix at selected q
* ---------------------------------------------------------------------------- */ * ---------------------------------------------------------------------------- */
@ -491,7 +402,7 @@ void Phonon::vfanyq()
char str[MAXLINE]; char str[MAXLINE];
double q[3], egvs[ndim]; double q[3], egvs[ndim];
while (1){ while ( 1 ){
printf("Please input the q-point to compute the frequencies, q to exit: "); printf("Please input the q-point to compute the frequencies, q to exit: ");
if (count_words(fgets(str,MAXLINE,stdin)) < 3) break; if (count_words(fgets(str,MAXLINE,stdin)) < 3) break;
@ -503,7 +414,7 @@ void Phonon::vfanyq()
dynmat->geteigen(egvs, 0); dynmat->geteigen(egvs, 0);
printf("q-point: [%lg %lg %lg], ", q[0], q[1], q[2]); printf("q-point: [%lg %lg %lg], ", q[0], q[1], q[2]);
printf("vibrational frequencies at this q-point:\n"); 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; return;
@ -521,7 +432,7 @@ void Phonon::vecanyq()
if (count_words(fgets(str,MAXLINE,stdin)) < 1) strcpy(str,"eigvec.dat"); if (count_words(fgets(str,MAXLINE,stdin)) < 1) strcpy(str,"eigvec.dat");
FILE *fp = fopen(strtok(str," \t\n\r\f"), "w"); FILE *fp = fopen(strtok(str," \t\n\r\f"), "w");
while (1){ while ( 1 ){
printf("Please input the q-point to compute the frequencies, q to exit: "); printf("Please input the q-point to compute the frequencies, q to exit: ");
if (count_words(fgets(str,MAXLINE,stdin)) < 3) break; if (count_words(fgets(str,MAXLINE,stdin)) < 3) break;
@ -533,14 +444,14 @@ void Phonon::vecanyq()
dynmat->geteigen(egvs, 1); dynmat->geteigen(egvs, 1);
fprintf(fp,"# q-point: [%lg %lg %lg], sysdim: %d, # of atoms per cell: %d\n", fprintf(fp,"# q-point: [%lg %lg %lg], sysdim: %d, # of atoms per cell: %d\n",
q[0],q[1],q[2], sysdim, dynmat->nucell); q[0],q[1],q[2], sysdim, dynmat->nucell);
for (int i=0; i<ndim; i++){ for (int i = 0; i < ndim; ++i){
fprintf(fp,"# frequency %d at [%lg %lg %lg]: %lg\n",i+1,q[0],q[1],q[2],egvs[i]); fprintf(fp,"# frequency %d at [%lg %lg %lg]: %lg\n",i+1,q[0],q[1],q[2],egvs[i]);
fprintf(fp,"# atom eigenvector : |e|\n"); fprintf(fp,"# atom eigenvector : |e|\n");
for (int j=0; j<dynmat->nucell; j++){ for (int j = 0; j < dynmat->nucell; ++j){
int ipos = j * sysdim; int ipos = j * sysdim;
double sum = 0.; double sum = 0.;
fprintf(fp,"%d", j+1); fprintf(fp,"%d", j+1);
for (int idim=0; idim<sysdim; idim++){ for (int idim = 0; idim < sysdim; ++idim){
fprintf(fp," %lg %lg", eigvec[i][ipos+idim].r, eigvec[i][ipos+idim].i); 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; sum += eigvec[i][ipos+idim].r * eigvec[i][ipos+idim].r + eigvec[i][ipos+idim].i * eigvec[i][ipos+idim].i;
} }
@ -575,9 +486,8 @@ void Phonon::DMdisp()
int nq = MAX(MAX(dynmat->nx,dynmat->ny),dynmat->nz)/2; int nq = MAX(MAX(dynmat->nx,dynmat->ny),dynmat->nz)/2;
qend[0] = qend[1] = qend[2] = 0.; qend[0] = qend[1] = qend[2] = 0.;
while (1){ while ( 1 ){
for (int i = 0; i < 3; ++i) qstr[i] = qend[i];
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]); 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)); int n = count_words(fgets(str,MAXLINE,stdin));
@ -602,16 +512,17 @@ void Phonon::DMdisp()
for (int i=0; i<3; i++) qinc[i] = (qend[i]-qstr[i])/double(nq-1); 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]); 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 i = 0; i < 3; ++i) q[i] = qstr[i];
for (int ii=0; ii<nq; ii++){ for (int ii = 0; ii < nq; ++ii){
dynmat->getDMq(q); dynmat->getDMq(q);
dynmat->writeDMq(q, qr, fp); dynmat->writeDMq(q, qr, fp);
for (int i=0; i<3; i++) q[i] += qinc[i]; for (int i = 0; i < 3; ++i) q[i] += qinc[i];
qr += dq; qr += dq;
} }
qr -= dq; qr -= dq;
} }
fclose(fp); fclose(fp);
return; return;
} }
@ -625,26 +536,26 @@ void Phonon::smooth(double *array, const int npt)
int nlag = npt/4; int nlag = npt/4;
double *tmp, *table; double *tmp, *table;
tmp = memory->create(tmp, npt, "smooth:tmp"); memory->create(tmp, npt, "smooth:tmp");
table = memory->create(table, nlag+1, "smooth:table"); memory->create(table, nlag+1, "smooth:table");
double fnorm = -1.; double fnorm = -1.;
double sigma = 4., fac = 1./(sigma*sigma); double sigma = 4., fac = 1./(sigma*sigma);
for (int jj=0; jj<= nlag; jj++){ for (int jj = 0; jj <= nlag; ++jj){
table[jj] = exp(-double(jj*jj)*fac); table[jj] = exp(-double(jj*jj)*fac);
fnorm += table[jj]; fnorm += table[jj];
} }
fnorm = 1./fnorm; fnorm = 1./fnorm;
for (int i=0; i<npt; i++){ for (int i = 0; i < npt; ++i){
tmp[i] = 0.; tmp[i] = 0.;
for (int jj=-nlag; jj<= nlag; jj++){ for (int jj = -nlag; jj <= nlag; ++jj){
int j = (i+jj+npt)%npt; // assume periodical data int j = (i+jj+npt)%npt; // assume periodical data
tmp [i] += array[j]*table[abs(jj)]; tmp [i] += array[j]*table[abs(jj)];
} }
} }
for (int i=0; i<npt; i++) array[i] = tmp[i]*fnorm; for (int i = 0; i < npt; ++i) array[i] = tmp[i]*fnorm;
memory->destroy(tmp); memory->destroy(tmp);
memory->destroy(table); memory->destroy(table);
@ -682,9 +593,9 @@ void Phonon::therm()
double h_o_KbT = h/(Kb*T)*1.e12, KbT_in_eV = Kb*T/eV; double 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.; double Uvib = 0., Svib = 0., Fvib = 0., Cvib = 0., ZPE = 0.;
for (int iq=0; iq<nq; iq++){ for (int iq = 0; iq < nq; ++iq){
double Utmp = 0., Stmp = 0., Ftmp = 0., Ztmp = 0., Ctmp = 0.; double Utmp = 0., Stmp = 0., Ftmp = 0., Ztmp = 0., Ctmp = 0.;
for (int i=0; i<ndim; i++){ for (int i = 0; i < ndim; ++i){
if (eigs[iq][i] <= 0.) continue; if (eigs[iq][i] <= 0.) continue;
double x = eigs[iq][i] * h_o_KbT; double x = eigs[iq][i] * h_o_KbT;
double expterm = 1./(exp(x)-1.); double expterm = 1./(exp(x)-1.);
@ -709,6 +620,7 @@ void Phonon::therm()
printf("Please input the desired temperature (K), enter to exit: "); printf("Please input the desired temperature (K), enter to exit: ");
if (count_words(fgets(str,MAXLINE,stdin)) < 1) break; if (count_words(fgets(str,MAXLINE,stdin)) < 1) break;
T = atof(strtok(str," \t\n\r\f")); T = atof(strtok(str," \t\n\r\f"));
} while (T > 0.); } while (T > 0.);
fclose(fp); fclose(fp);
@ -737,15 +649,15 @@ void Phonon::local_therm()
fprintf(fp,"#-------------------------------------------------------------------------------\n"); fprintf(fp,"#-------------------------------------------------------------------------------\n");
double **Uvib, **Svib, **Fvib, **Cvib, **ZPE; double **Uvib, **Svib, **Fvib, **Cvib, **ZPE;
Uvib = memory->create(Uvib,nlocal,sysdim,"local_therm:Uvib"); memory->create(Uvib,nlocal,sysdim,"local_therm:Uvib");
Svib = memory->create(Svib,nlocal,sysdim,"local_therm:Svib"); memory->create(Svib,nlocal,sysdim,"local_therm:Svib");
Fvib = memory->create(Fvib,nlocal,sysdim,"local_therm:Fvib"); memory->create(Fvib,nlocal,sysdim,"local_therm:Fvib");
Cvib = memory->create(Cvib,nlocal,sysdim,"local_therm:Cvib"); memory->create(Cvib,nlocal,sysdim,"local_therm:Cvib");
ZPE = memory->create(ZPE ,nlocal,sysdim,"local_therm:ZPE"); memory->create(ZPE ,nlocal,sysdim,"local_therm:ZPE");
// constants J.s J/K J // constants J.s J/K J
const double h = 6.62606896e-34, Kb = 1.380658e-23, eV = 1.60217733e-19; const double h = 6.62606896e-34, Kb = 1.380658e-23, eV = 1.60217733e-19;
double T = dynmat->Tmeasure; double T = dynmat->Tmeasure;
while (1){ while ( 1 ){
printf("\nPlease input the temperature at which to evaluate the local vibrational\n"); printf("\nPlease input the temperature at which to evaluate the local vibrational\n");
printf("thermal properties, non-positive number to exit [%g]: ", T); printf("thermal properties, non-positive number to exit [%g]: ", T);
if (count_words(fgets(str,MAXLINE,stdin)) > 0){ if (count_words(fgets(str,MAXLINE,stdin)) > 0){
@ -755,11 +667,11 @@ void Phonon::local_therm()
// constants under the same temperature; assuming angular frequency in THz // 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 h_o_KbT = h/(Kb*T)*1.e12, KbT_in_eV = Kb*T/eV;
for (int i=0; i<nlocal; i++) for (int i = 0; i <nlocal; ++i)
for (int j=0; j<sysdim; j++) Uvib[i][j] = Svib[i][j] = Fvib[i][j] = Cvib[i][j] = ZPE[i][j] = 0.; for (int j = 0; j <sysdim; ++j) Uvib[i][j] = Svib[i][j] = Fvib[i][j] = Cvib[i][j] = ZPE[i][j] = 0.;
double freq = fmin-df; double freq = fmin-df;
for (int i=0; i<ndos; i++){ for (int i = 0; i < ndos; ++i){
freq += df; freq += df;
if (freq <= 0.) continue; if (freq <= 0.) continue;
@ -772,8 +684,8 @@ void Phonon::local_therm()
double Ctmp = x*x*exp(x)*expterm*expterm; double Ctmp = x*x*exp(x)*expterm*expterm;
double Ztmp = 0.5*h*freq; double Ztmp = 0.5*h*freq;
for (int il=0; il<nlocal; il++) for (int il = 0; il < nlocal; ++il)
for (int idim=0; idim<sysdim; idim++){ for (int idim = 0; idim < sysdim; ++idim){
Uvib[il][idim] += ldos[il][i][idim]*Utmp; Uvib[il][idim] += ldos[il][i][idim]*Utmp;
Svib[il][idim] += ldos[il][i][idim]*Stmp; Svib[il][idim] += ldos[il][i][idim]*Stmp;
Fvib[il][idim] += ldos[il][i][idim]*Ftmp; Fvib[il][idim] += ldos[il][i][idim]*Ftmp;
@ -781,8 +693,8 @@ void Phonon::local_therm()
ZPE [il][idim] += ldos[il][i][idim]*Ztmp; ZPE [il][idim] += ldos[il][i][idim]*Ztmp;
} }
} }
for (int il=0; il<nlocal; il++) for (int il = 0; il < nlocal; ++il)
for (int idim=0; idim<sysdim; idim++){ for (int idim = 0; idim < sysdim; ++idim){
Uvib[il][idim] *= KbT_in_eV*df; Uvib[il][idim] *= KbT_in_eV*df;
Svib[il][idim] *= df; Svib[il][idim] *= df;
Fvib[il][idim] *= KbT_in_eV*df; Fvib[il][idim] *= KbT_in_eV*df;
@ -791,34 +703,34 @@ void Phonon::local_therm()
} }
// output result under current temperature // output result under current temperature
for (int il=0; il<nlocal; il++){ for (int il = 0; il < nlocal; ++il){
fprintf(fp,"%d %g ", locals[il], T); fprintf(fp,"%d %g ", locals[il], T);
double total = 0.; double total = 0.;
for (int idim=0; idim<sysdim; idim++){ for (int idim = 0; idim < sysdim; ++idim){
fprintf(fp,"%g ", Uvib[il][idim]); fprintf(fp,"%g ", Uvib[il][idim]);
total += Uvib[il][idim]; total += Uvib[il][idim];
} }
fprintf(fp,"%g ", total); total = 0.; fprintf(fp,"%g ", total); total = 0.;
for (int idim=0; idim<sysdim; idim++){ for (int idim = 0; idim < sysdim; ++idim){
fprintf(fp,"%g ", Svib[il][idim]); fprintf(fp,"%g ", Svib[il][idim]);
total += Svib[il][idim]; total += Svib[il][idim];
} }
fprintf(fp,"%g ", total); total = 0.; fprintf(fp,"%g ", total); total = 0.;
for (int idim=0; idim<sysdim; idim++){ for (int idim = 0; idim < sysdim; ++idim){
fprintf(fp,"%g ", Fvib[il][idim]); fprintf(fp,"%g ", Fvib[il][idim]);
total += Fvib[il][idim]; total += Fvib[il][idim];
} }
fprintf(fp,"%g ", total); total = 0.; fprintf(fp,"%g ", total); total = 0.;
for (int idim=0; idim<sysdim; idim++){ for (int idim = 0; idim < sysdim; ++idim){
fprintf(fp,"%g ", Cvib[il][idim]); fprintf(fp,"%g ", Cvib[il][idim]);
total += Cvib[il][idim]; total += Cvib[il][idim];
} }
fprintf(fp,"%g ", total); total = 0.; fprintf(fp,"%g ", total); total = 0.;
for (int idim=0; idim<sysdim; idim++){ for (int idim = 0; idim < sysdim; ++idim){
fprintf(fp,"%g ", ZPE[il][idim]); fprintf(fp,"%g ", ZPE[il][idim]);
total += ZPE[il][idim]; total += ZPE[il][idim];
} }
@ -846,7 +758,7 @@ void Phonon::QMesh()
ny = atoi(strtok(NULL," \t\n\r\f")); ny = atoi(strtok(NULL," \t\n\r\f"));
nz = atoi(strtok(NULL," \t\n\r\f")); nz = atoi(strtok(NULL," \t\n\r\f"));
} }
if (nx<1||ny<1||nz<1) return; if (nx < 1||ny < 1||nz < 1) return;
if (dynmat->nx == 1) nx = 1; if (dynmat->nx == 1) nx = 1;
if (dynmat->ny == 1) ny = 1; if (dynmat->ny == 1) ny = 1;
if (dynmat->nz == 1) nz = 1; if (dynmat->nz == 1) nz = 1;
@ -855,9 +767,9 @@ void Phonon::QMesh()
// ask method to generate q-points // ask method to generate q-points
int method = 2; int method = 2;
printf("Please select your method to generate the q-points:\n"); printf("Please select your method to generate the q-points:\n");
printf(" 1. uniform;\n 2. Monkhost-Pack mesh;\nYour choice[2]: "); printf(" 1. uniform;\n 2. Monkhost-Pack mesh;\nYour choice [2]: ");
if (count_words(fgets(str,MAXLINE,stdin)) > 0) method = atoi(strtok(str," \t\n\r\f")); if (count_words(fgets(str,MAXLINE,stdin)) > 0) method = atoi(strtok(str," \t\n\r\f"));
method = 2-method%2; method = 2 - method%2;
printf("Your selection: %d\n", method); printf("Your selection: %d\n", method);
#endif #endif
@ -869,149 +781,52 @@ void Phonon::QMesh()
#endif #endif
nq = nx*ny*nz; nq = nx*ny*nz;
double w = 1./double(nq); double w = 1./double(nq);
wt = memory->create(wt, nq, "QMesh:wt"); memory->create(wt, nq, "QMesh:wt");
qpts = memory->create(qpts, nq, 3, "QMesh:qpts"); memory->create(qpts, nq, 3, "QMesh:qpts");
int iq = 0; int iq = 0;
for (int i=0; i<nx; i++) for (int i = 0; i < nx; ++i)
for (int j=0; j<ny; j++) for (int j = 0; j < ny; ++j)
for (int k=0; k<nz; k++){ for (int k = 0; k < nz; ++k){
qpts[iq][0] = double(i)/double(nx); qpts[iq][0] = double(i)/double(nx);
qpts[iq][1] = double(j)/double(ny); qpts[iq][1] = double(j)/double(ny);
qpts[iq][2] = double(k)/double(nz); qpts[iq][2] = double(k)/double(nz);
wt[iq++] = w; wt[iq++] = w;
} }
#ifdef UseSPG #ifdef UseSPG
} } else {
if ((method == 2) && (atpos == NULL)){ if (atpos == NULL) memory->create(atpos, dynmat->nucell,3,"QMesh:atpos");
atpos = memory->create(atpos, dynmat->nucell,3,"QMesh:atpos"); if (attyp == NULL) memory->create(attyp, dynmat->nucell, "QMesh:attyp");
attyp = memory->create(attyp, dynmat->nucell, "QMesh:attyp");
for (int i=0; i<dynmat->nucell; i++)
for (int idim=0; idim<3; idim++) atpos[i][idim] = 0.;
for (int i=0; i<3; i++) latvec[i][i] = 1.;
int flag_lat_info_read = dynmat->flag_latinfo;
if ( flag_lat_info_read ){ // get unit cell info from binary file; done by dynmat
num_atom = dynmat->nucell; num_atom = dynmat->nucell;
// set default, in case system dimension under study is not 3. // set default, in case system dimension under study is not 3.
for (int i=0; i<dynmat->nucell; i++) for (int i = 0; i < dynmat->nucell; ++i)
for (int idim=0; idim<3; idim++) atpos[i][idim] = 0.; for (int idim = 0; idim < 3; ++idim) atpos[i][idim] = 0.;
for (int i=0; i<3; i++) latvec[i][i] = 1.; 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 // get atomic type info
for (int i=0; i<num_atom; i++) attyp[i] = dynmat->attyp[i]; for (int i = 0; i < num_atom; ++i) attyp[i] = dynmat->attyp[i];
// get unit cell vector info // get unit cell vector info
int ndim = 0; int ndim = 0;
for (int idim=0; idim<3; idim++) for (int idim = 0; idim < 3; ++idim)
for (int jdim=0; jdim<3; jdim++) latvec[jdim][idim] = dynmat->basevec[ndim++]; for (int jdim = 0; jdim < 3; ++jdim) latvec[jdim][idim] = dynmat->basevec[ndim++];
// get atom position in unit cell; fractional // get atom position in unit cell; fractional
for (int i=0; i<num_atom; i++) for (int i = 0; i < num_atom; ++i)
for (int idim=0; idim<sysdim; idim++) atpos[i][idim] = dynmat->basis[i][idim]; for (int idim = 0; idim < sysdim; ++idim) atpos[i][idim] = dynmat->basis[i][idim];
// display the unit cell info read // display the unit cell info read
printf("\n");for (int ii=0; ii<80; ii++) printf("="); printf("\n"); printf("\n");for (int ii = 0; ii < 80; ++ii) printf("="); printf("\n");
printf("The basis vectors of the unit cell:\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]); 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("Atom(s) in the unit cell:\n");
printf(" No. type sx sy sz\n"); printf(" No. type sx sy sz\n");
for (int i=0; i<num_atom; i++) printf(" %d %d %lg %lg %lg\n", i+1, attyp[i], atpos[i][0], atpos[i][1], atpos[i][2]); 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]);
printf("\nIs the above info correct? (y/n)[y]: "); if (num_atom > NUMATOM) printf(" ... (%d atoms omitted.)\n", num_atom - NUMATOM);
fgets(str,MAXLINE,stdin);
char *ptr = strtok(str," \t\n\r\f");
if ( (ptr) && ( (strcmp(ptr,"y") != 0) && (strcmp(ptr,"Y") != 0)) ) flag_lat_info_read = 0;
}
if (flag_lat_info_read == 0) { // get unit cell info from file or user input
int latsrc = 1;
printf("\nPlease select the way to provide the unit cell info:\n");
printf(" 1. By file;\n 2. Read in.\nYour choice [1]: ");
if (count_words(fgets(str,MAXLINE,stdin)) > 0) latsrc = atoi(strtok(str," \t\n\r\f"));
latsrc = 2-latsrc%2;
/*----------------------------------------------------------------
* Ask for lattice info from the user; the format of the file is:
* A1_x A1_y A1_z
* A2_x A2_y A2_z
* A3_x A3_y A3_z
* natom
* Type_1 sx_1 sy_1 sz_1
* ...
* Type_n sx_n sy_n sz_n
*----------------------------------------------------------------*/
if (latsrc == 1){ // to read unit cell info from file; get file name first
do printf("Please input the file name containing the unit cell info: ");
while (count_words(fgets(str,MAXLINE,stdin)) < 1);
char *fname = strtok(str," \t\n\r\f");
FILE *fp = fopen(fname,"r"); fname = NULL;
if (fp == NULL) latsrc = 2;
else {
for (int i=0; i<3; i++){ // read unit cell vector info; # of atoms per unit cell
if (count_words(fgets(str,MAXLINE,fp)) < 3){
latsrc = 2;
break;
}
latvec[0][i] = atof(strtok(str, " \t\n\r\f"));
latvec[1][i] = atof(strtok(NULL," \t\n\r\f"));
latvec[2][i] = atof(strtok(NULL," \t\n\r\f"));
}
if (count_words(fgets(str,MAXLINE,fp)) < 1) latsrc = 2;
else {
num_atom = atoi(strtok(str," \t\n\r\f"));
if (num_atom > dynmat->nucell){
printf("\nError: # of atoms read from file (%d) is bigger than that given by the dynamical matrix (%d)!\n", num_atom, dynmat->nucell);
return;
}
for (int i=0; i<num_atom; i++){ // read atomic type and fractional positions
if (count_words(fgets(str,MAXLINE,fp)) < 4){
latsrc = 2;
break;
} else {
attyp[i] = atoi(strtok(str," \t\n\r\f"));
atpos[i][0] = atof(strtok(NULL," \t\n\r\f"));
atpos[i][1] = atof(strtok(NULL," \t\n\r\f"));
atpos[i][2] = atof(strtok(NULL," \t\n\r\f"));
}
}
}
}
fclose(fp);
}
if (latsrc == 2){
for (int i=0; i<3; i++){
do printf("Please input the vector A%d: ", i+1);
while (count_words(fgets(str,MAXLINE,stdin)) < 3);
latvec[0][i] = atof(strtok(str," \t\n\r\f"));
latvec[1][i] = atof(strtok(NULL," \t\n\r\f"));
latvec[2][i] = atof(strtok(NULL," \t\n\r\f"));
}
do printf("please input the number of atoms per unit cell: ");
while (count_words(fgets(str,MAXLINE,stdin)) < 1);
num_atom = atoi(strtok(str," \t\n\r\f"));
if (num_atom > dynmat->nucell){
printf("\nError: # of atoms input (%d) is bigger than that given by the dynamical matrix (%d)!\n", num_atom, dynmat->nucell);
return;
}
for (int i=0; i<num_atom; i++){
do printf("Please input the type, and fractional coordinate of atom No.%d: ", i+1);
while (count_words(fgets(str,MAXLINE,stdin)) < 4);
attyp[i] = atoi(strtok(str," \t\n\r\f"));
atpos[i][0] = atof(strtok(NULL," \t\n\r\f"));
atpos[i][1] = atof(strtok(NULL," \t\n\r\f"));
atpos[i][2] = atof(strtok(NULL," \t\n\r\f"));
}
}
} // end of read from file or input
} // end of if (method == 2 && ...
if (method == 2){
int mesh[3], shift[3], is_time_reversal = 0; int mesh[3], shift[3], is_time_reversal = 0;
mesh[0] = nx; mesh[1] = ny; mesh[2] = nz; mesh[0] = nx; mesh[1] = ny; mesh[2] = nz;
shift[0] = shift[1] = shift[2] = 0; shift[0] = shift[1] = shift[2] = 0;
@ -1019,29 +834,28 @@ void Phonon::QMesh()
int grid_point[num_grid][3], map[num_grid]; int grid_point[num_grid][3], map[num_grid];
double symprec = 1.e-4, pos[num_atom][3]; double symprec = 1.e-4, pos[num_atom][3];
for (int i=0; i<num_atom; i++) for (int i = 0; i < num_atom; ++i)
for (int j=0; j<3; j++) pos[i][j] = atpos[i][j]; for (int j = 0; j < 3; ++j) pos[i][j] = atpos[i][j];
//spg_show_symmetry(latvec, pos, attyp, num_atom, symprec);
// if spglib 0.7.1 is used // if spglib 0.7.1 is used
nq = spg_get_ir_reciprocal_mesh(grid_point, map, num_grid, mesh, shift, is_time_reversal, latvec, pos, attyp, num_atom, symprec); // nq = spg_get_ir_reciprocal_mesh(grid_point, map, num_grid, mesh, shift, is_time_reversal, latvec, pos, attyp, num_atom, symprec);
// if spglib >= 1.0.3 is used // if spglib >= 1.0.3 is used
//nq = spg_get_ir_reciprocal_mesh(grid_point, map, mesh, shift, is_time_reversal, latvec, pos, attyp, num_atom, symprec); nq = spg_get_ir_reciprocal_mesh(grid_point, map, mesh, shift, is_time_reversal, latvec, pos, attyp, num_atom, symprec);
wt = memory->create(wt, nq, "QMesh:wt"); memory->create(wt, nq, "QMesh:wt");
qpts = memory->create(qpts, nq,3,"QMesh:qpts"); memory->create(qpts, nq,3,"QMesh:qpts");
int *iq2idx = new int[num_grid]; int *iq2idx = new int[num_grid];
int numq = 0; int numq = 0;
for (int i=0; i<num_grid; i++){ for (int i = 0; i < num_grid; ++i){
int iq = map[i]; int iq = map[i];
if (iq == i) iq2idx[iq] = numq++; if (iq == i) iq2idx[iq] = numq++;
} }
for (int iq=0; iq<nq; iq++) wt[iq] = 0.; for (int iq = 0; iq < nq; ++iq) wt[iq] = 0.;
numq = 0; numq = 0;
for (int i=0; i<num_grid; i++){ for (int i = 0; i < num_grid; ++i){
int iq = map[i]; int iq = map[i];
if (iq == i){ if (iq == i){
qpts[numq][0] = double(grid_point[i][0])/double(mesh[0]); qpts[numq][0] = double(grid_point[i][0])/double(mesh[0]);
@ -1054,8 +868,8 @@ void Phonon::QMesh()
delete []iq2idx; delete []iq2idx;
double wsum = 0.; double wsum = 0.;
for (int iq=0; iq<nq; iq++) wsum += wt[iq]; for (int iq = 0; iq < nq; ++iq) wsum += wt[iq];
for (int iq=0; iq<nq; iq++) wt[iq] /= wsum; for (int iq = 0; iq < nq; ++iq) wt[iq] /= wsum;
} }
#endif #endif
@ -1078,7 +892,7 @@ void Phonon::ldos_egv()
if (nmax < 1) return; if (nmax < 1) return;
memory->destroy(locals); memory->destroy(locals);
locals = memory->create(locals, nmax, "ldos_egv:locals"); memory->create(locals, nmax, "ldos_egv:locals");
nlocal = 0; nlocal = 0;
ptr = strtok(str," \t\n\r\f"); ptr = strtok(str," \t\n\r\f");
@ -1091,7 +905,7 @@ void Phonon::ldos_egv()
if (nlocal < 1) return; if (nlocal < 1) return;
printf("Local PDOS for atom(s):"); printf("Local PDOS for atom(s):");
for (int i=0; i<nlocal; i++) printf(" %d", locals[i]); for (int i = 0; i < nlocal; ++i) printf(" %d", locals[i]);
printf(" will be computed.\n"); printf(" will be computed.\n");
fmin = 0.; fmax = 10.; fmin = 0.; fmax = 10.;
@ -1118,14 +932,14 @@ void Phonon::ldos_egv()
memory->destroy(dos); memory->destroy(dos);
memory->destroy(ldos); memory->destroy(ldos);
dos = memory->create(dos, ndos,"ldos_egv:dos"); memory->create(dos, ndos,"ldos_egv:dos");
ldos = memory->create(ldos,nlocal,ndos,sysdim,"ldos_egv:ldos"); memory->create(ldos,nlocal,ndos,sysdim,"ldos_egv:ldos");
for (int i=0; i<ndos; i++) dos[i] = 0.; for (int i = 0; i < ndos; ++i) dos[i] = 0.;
for (int ilocal=0; ilocal<nlocal; ilocal++) for (int ilocal = 0; ilocal < nlocal; ++ilocal)
for (int i=0; i<ndos; i++) for (int i = 0; i < ndos; ++i)
for (int idim=0; idim<sysdim; idim++) ldos[ilocal][i][idim] = 0.; for (int idim = 0; idim < sysdim; ++idim) ldos[ilocal][i][idim] = 0.;
int nprint; int nprint;
if (nq > 10) nprint = nq/10; if (nq > 10) nprint = nq/10;
@ -1137,7 +951,7 @@ void Phonon::ldos_egv()
doublecomplex **egvec = dynmat->DM_q; doublecomplex **egvec = dynmat->DM_q;
printf("\nNow to compute the phonons and DOSs "); fflush(stdout); printf("\nNow to compute the phonons and DOSs "); fflush(stdout);
for (int iq=0; iq<nq; iq++){ for (int iq = 0; iq < nq; ++iq){
if ((iq+1)%nprint == 0) {printf("."); fflush(stdout);} if ((iq+1)%nprint == 0) {printf("."); fflush(stdout);}
dynmat->getDMq(qpts[iq], &wt[iq]); dynmat->getDMq(qpts[iq], &wt[iq]);
@ -1145,14 +959,14 @@ void Phonon::ldos_egv()
dynmat->geteigen(&egval[0], 1); dynmat->geteigen(&egval[0], 1);
for (int idim=0; idim<ndim; idim++){ for (int idim = 0; idim < ndim; ++idim){
int hit = int((egval[idim] - offset)*rdf); int hit = int((egval[idim] - offset)*rdf);
if (hit >= 0 && hit <ndos){ if (hit >= 0 && hit <ndos){
dos[hit] += wt[iq]; dos[hit] += wt[iq];
for (int ilocal=0; ilocal<nlocal; ilocal++){ for (int ilocal = 0; ilocal < nlocal; ++ilocal){
int ipos = locals[ilocal]*sysdim; int ipos = locals[ilocal]*sysdim;
for (int jdim=0; jdim<sysdim; jdim++){ for (int jdim = 0; jdim < sysdim; ++jdim){
double dr = egvec[idim][ipos+jdim].r, di = egvec[idim][ipos+jdim].i; double dr = egvec[idim][ipos+jdim].r, di = egvec[idim][ipos+jdim].i;
double norm = dr * dr + di * di; double norm = dr * dr + di * di;
ldos[ilocal][hit][jdim] += wt[iq] * norm; ldos[ilocal][hit][jdim] += wt[iq] * norm;
@ -1179,6 +993,33 @@ void Phonon::ldos_egv()
return; return;
} }
/* ----------------------------------------------------------------------------
* Private method to show the unit cell info
* ---------------------------------------------------------------------------- */
void Phonon::ShowCell()
{
printf("\n");
for (int i = 0; i < 30; ++i) printf("=");
printf(" Unit Cell Info ");
for (int i = 0; i < 30; ++i) printf("="); printf("\n");
printf("Number of atoms in the unit cell: %d\n", dynmat->nucell);
printf("Basis vectors of the unit cell:\n");
printf(" %15.8f %15.8f %15.8f\n", dynmat->basevec[0], dynmat->basevec[1], dynmat->basevec[2]);
printf(" %15.8f %15.8f %15.8f\n", dynmat->basevec[3], dynmat->basevec[4], dynmat->basevec[5]);
printf(" %15.8f %15.8f %15.8f\n", dynmat->basevec[6], dynmat->basevec[7], dynmat->basevec[8]);
printf("Basis vectors of the reciprocal:\n");
printf(" %15.8f %15.8f %15.8f\n", dynmat->ibasevec[0], dynmat->ibasevec[1], dynmat->ibasevec[2]);
printf(" %15.8f %15.8f %15.8f\n", dynmat->ibasevec[3], dynmat->ibasevec[4], dynmat->ibasevec[5]);
printf(" %15.8f %15.8f %15.8f\n", dynmat->ibasevec[6], dynmat->ibasevec[7], dynmat->ibasevec[8]);
printf("Atomic type and fractional coordinates:\n");
for (int i = 0; i < dynmat->nucell; ++i)
printf("%4d %12.8f %12.8f %12.8f\n", dynmat->attyp[i], dynmat->basis[i][0], dynmat->basis[i][1], dynmat->basis[i][2]);
for (int i = 0; i < 80; ++i) printf("=");
printf("\n");
return;
}
/* ---------------------------------------------------------------------------- /* ----------------------------------------------------------------------------
* Private method to normalize the DOS and/or Local DOS. * Private method to normalize the DOS and/or Local DOS.
* Simpson's rule is used for the integration. * Simpson's rule is used for the integration.
@ -1188,24 +1029,24 @@ void Phonon::Normalize()
double odd, even, sum; double odd, even, sum;
if (dos){ if (dos){
odd = even = 0.; odd = even = 0.;
for (int i=1; i<ndos-1; i +=2) odd += dos[i]; for (int i = 1; i < ndos-1; i +=2) odd += dos[i];
for (int i=2; i<ndos-1; i +=2) even += dos[i]; for (int i = 2; i < ndos-1; i +=2) even += dos[i];
sum = dos[0] + dos[ndos-1]; sum = dos[0] + dos[ndos-1];
sum += 4.*odd + 2.*even; sum += 4.*odd + 2.*even;
sum = 3.*rdf/sum; sum = 3.*rdf/sum;
for (int i=0; i<ndos; i++) dos[i] *= sum; for (int i = 0; i < ndos; ++i) dos[i] *= sum;
} }
if (ldos){ if (ldos){
for (int ilocal=0; ilocal<nlocal; ilocal++) for (int ilocal = 0; ilocal < nlocal; ++ilocal)
for (int idim=0; idim<sysdim; idim++){ for (int idim = 0; idim < sysdim; ++idim){
odd = even = 0.; odd = even = 0.;
for (int i=1; i<ndos-1; i +=2) odd += ldos[ilocal][i][idim]; for (int i = 1; i < ndos-1; i += 2) odd += ldos[ilocal][i][idim];
for (int i=2; i<ndos-1; i +=2) even += ldos[ilocal][i][idim]; 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 = ldos[ilocal][0][idim] + ldos[ilocal][ndos-1][idim];
sum += 4.*odd + 2.*even; sum += 4.*odd + 2.*even;
sum = 3.*rdf/sum; sum = 3.*rdf/sum;
for (int i=0; i<ndos; i++) ldos[ilocal][i][idim] *= sum; for (int i = 0; i < ndos; ++i) ldos[ilocal][i][idim] *= sum;
} }
} }
@ -1225,9 +1066,9 @@ void Phonon::ComputeAll()
printf("\nNow to compute the phonons "); fflush(stdout); printf("\nNow to compute the phonons "); fflush(stdout);
// now to calculate the frequencies at all q-points // now to calculate the frequencies at all q-points
memory->destroy(eigs); memory->destroy(eigs);
eigs = memory->create(eigs, nq,ndim,"QMesh_eigs"); memory->create(eigs, nq,ndim,"QMesh_eigs");
for (int iq=0; iq<nq; iq++){ for (int iq = 0; iq < nq; ++iq){
if ((iq+1)%nprint == 0) {printf("."); fflush(stdout);} if ((iq+1)%nprint == 0) {printf("."); fflush(stdout);}
dynmat->getDMq(qpts[iq], &wt[iq]); dynmat->getDMq(qpts[iq], &wt[iq]);
@ -1246,7 +1087,7 @@ int Phonon::count_words(const char *line)
{ {
int n = strlen(line) + 1; int n = strlen(line) + 1;
char *copy; char *copy;
copy = memory->create(copy, n, "count_words:copy"); memory->create(copy, n, "count_words:copy");
strcpy(copy,line); strcpy(copy,line);
char *ptr; char *ptr;

2
tools/phonon/phonon.h
View File

@ -43,6 +43,8 @@ private:
void DMdisp(); void DMdisp();
void vecanyq(); void vecanyq();
void ShowCell();
void smooth(double *, const int); void smooth(double *, const int);
void writeDOS(); void writeDOS();
void writeLDOS(); void writeLDOS();

16
tools/phonon/timer.cpp
View File

@ -1,5 +1,7 @@
#include "timer.h" #include "timer.h"
/* -----------------------------------------------------------------------------
* Initialization of time
* -------------------------------------------------------------------------- */
Timer::Timer() Timer::Timer()
{ {
flag = 0; flag = 0;
@ -7,6 +9,9 @@ Timer::Timer()
return; return;
} }
/* -----------------------------------------------------------------------------
* public function, start the timer
* -------------------------------------------------------------------------- */
void Timer::start() void Timer::start()
{ {
t1 = clock(); t1 = clock();
@ -15,6 +20,9 @@ void Timer::start()
return; return;
} }
/* -----------------------------------------------------------------------------
* public function, stop the timer
* -------------------------------------------------------------------------- */
void Timer::stop() void Timer::stop()
{ {
if ( flag&1 ) { if ( flag&1 ) {
@ -24,6 +32,9 @@ void Timer::stop()
return; return;
} }
/* -----------------------------------------------------------------------------
* public function, print the total time used after timer stops
* -------------------------------------------------------------------------- */
void Timer::print() void Timer::print()
{ {
if ( (flag&3) != 3) return; if ( (flag&3) != 3) return;
@ -34,6 +45,9 @@ void Timer::print()
return; return;
} }
/* -----------------------------------------------------------------------------
* public function, return the total time used up to now, in seconds
* -------------------------------------------------------------------------- */
double Timer::elapse() double Timer::elapse()
{ {
if ( (flag&3) != 3) return 0.; if ( (flag&3) != 3) return 0.;

2
tools/phonon/version.h
View File

@ -1 +1 @@
#define VERSION 75 #define VERSION 7