git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@9633 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
sjplimp
64
tools/phonon/Makefile
Normal file
64
tools/phonon/Makefile
Normal file
@ -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 $<
|
||||||
31
tools/phonon/README
Normal file
31
tools/phonon/README
Normal file
@ -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
|
||||||
602
tools/phonon/dynmat.cpp
Normal file
602
tools/phonon/dynmat.cpp
Normal file
@ -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; idim<fftdim; idim++){
|
||||||
|
for (int jdim=0; jdim<fftdim; jdim++){
|
||||||
|
double inv_mass = M_inv_sqrt[idim/sysdim]*M_inv_sqrt[jdim/sysdim];
|
||||||
|
DM_all[idq][ndim].r *= inv_mass;
|
||||||
|
DM_all[idq][ndim].i *= inv_mass;
|
||||||
|
ndim++;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// ask for the interpolation method
|
||||||
|
interpolate->set_method();
|
||||||
|
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
|
||||||
|
// to destroy the class
|
||||||
|
DynMat::~DynMat()
|
||||||
|
{
|
||||||
|
// destroy all memory allocated
|
||||||
|
if (funit) delete []funit;
|
||||||
|
if (dmfile) delete []dmfile;
|
||||||
|
if (binfile) delete []binfile;
|
||||||
|
if (interpolate) delete interpolate;
|
||||||
|
|
||||||
|
memory->destroy(DM_q);
|
||||||
|
memory->destroy(attyp);
|
||||||
|
memory->destroy(basis);
|
||||||
|
memory->destroy(DM_all);
|
||||||
|
memory->destroy(M_inv_sqrt);
|
||||||
|
if (memory) delete memory;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------------
|
||||||
|
* method to write DM_q to file, single point
|
||||||
|
* ---------------------------------------------------------------------------- */
|
||||||
|
void DynMat::writeDMq(double *q)
|
||||||
|
{
|
||||||
|
FILE *fp;
|
||||||
|
// only ask for file name for the first time
|
||||||
|
// other calls will append the result to the file.
|
||||||
|
if (dmfile == NULL){
|
||||||
|
char str[MAXLINE], *ptr;
|
||||||
|
printf("\n");
|
||||||
|
while (1){
|
||||||
|
printf("Please input the filename to output the DM at selected q: ");
|
||||||
|
fgets(str,MAXLINE,stdin);
|
||||||
|
ptr = strtok(str, " \r\t\n\f");
|
||||||
|
if (ptr) break;
|
||||||
|
}
|
||||||
|
|
||||||
|
int n = strlen(ptr) + 1;
|
||||||
|
dmfile = new char[n];
|
||||||
|
strcpy(dmfile, ptr);
|
||||||
|
fp = fopen(dmfile,"w");
|
||||||
|
|
||||||
|
} else {
|
||||||
|
fp = fopen(dmfile,"a");
|
||||||
|
}
|
||||||
|
fprintf(fp,"# q = [%lg %lg %lg]\n", q[0], q[1], q[2]);
|
||||||
|
|
||||||
|
for (int i=0; i<fftdim; i++){
|
||||||
|
for (int j=0; j<fftdim; j++) fprintf(fp,"%lg %lg\t", DM_q[i][j].r, DM_q[i][j].i);
|
||||||
|
fprintf(fp,"\n");
|
||||||
|
}
|
||||||
|
fprintf(fp,"\n");
|
||||||
|
fclose(fp);
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------------
|
||||||
|
* method to write DM_q to file, dispersion-like
|
||||||
|
* ---------------------------------------------------------------------------- */
|
||||||
|
void DynMat::writeDMq(double *q, const double qr, FILE *fp)
|
||||||
|
{
|
||||||
|
|
||||||
|
fprintf(fp, "%lg %lg %lg %lg ", q[0], q[1], q[2], qr);
|
||||||
|
|
||||||
|
for (int i=0; i<fftdim; i++){
|
||||||
|
for (int j=0; j<fftdim; j++) fprintf(fp,"%lg %lg\t", DM_q[i][j].r, DM_q[i][j].i);
|
||||||
|
}
|
||||||
|
fprintf(fp,"\n");
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------------
|
||||||
|
* method to evaluate the eigenvalues of current q-point;
|
||||||
|
* return the eigenvalues in egv.
|
||||||
|
* cLapack subroutine zheevd is employed.
|
||||||
|
* ---------------------------------------------------------------------------- */
|
||||||
|
int DynMat::geteigen(double *egv, int flag)
|
||||||
|
{
|
||||||
|
char jobz, uplo;
|
||||||
|
integer n, lda, lwork, lrwork, *iwork, liwork, info;
|
||||||
|
doublecomplex *work;
|
||||||
|
doublereal *w = &egv[0], *rwork;
|
||||||
|
|
||||||
|
n = fftdim;
|
||||||
|
if (flag) jobz = 'V';
|
||||||
|
else jobz = 'N';
|
||||||
|
|
||||||
|
uplo = 'U';
|
||||||
|
lwork = (n+2)*n;
|
||||||
|
lrwork = 1 + (5+n+n)*n;
|
||||||
|
liwork = 3 + 5*n;
|
||||||
|
lda = n;
|
||||||
|
|
||||||
|
work = memory->create(work, lwork, "geteigen:work");
|
||||||
|
rwork = memory->create(rwork, lrwork, "geteigen:rwork");
|
||||||
|
iwork = memory->create(iwork, liwork, "geteigen:iwork");
|
||||||
|
|
||||||
|
zheevd_(&jobz, &uplo, &n, DM_q[0], &lda, w, work, &lwork, rwork, &lrwork, iwork, &liwork, &info);
|
||||||
|
|
||||||
|
// to get w instead of w^2; and convert w into v (THz hopefully)
|
||||||
|
for (int i=0; i<n; i++){
|
||||||
|
if (w[i]>= 0.) w[i] = sqrt(w[i]);
|
||||||
|
else w[i] = -sqrt(-w[i]);
|
||||||
|
|
||||||
|
w[i] *= eml2f;
|
||||||
|
}
|
||||||
|
|
||||||
|
memory->destroy(work);
|
||||||
|
memory->destroy(rwork);
|
||||||
|
memory->destroy(iwork);
|
||||||
|
|
||||||
|
return info;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------------
|
||||||
|
* method to get the Dynamical Matrix at q
|
||||||
|
* ---------------------------------------------------------------------------- */
|
||||||
|
void DynMat::getDMq(double *q)
|
||||||
|
{
|
||||||
|
interpolate->execute(q, DM_q[0]);
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------------
|
||||||
|
* method to get the Dynamical Matrix at q
|
||||||
|
* ---------------------------------------------------------------------------- */
|
||||||
|
void DynMat::getDMq(double *q, double *wt)
|
||||||
|
{
|
||||||
|
interpolate->execute(q, DM_q[0]);
|
||||||
|
|
||||||
|
if (flag_skip && interpolate->UseGamma ) wt[0] = 0.;
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------------
|
||||||
|
* private method to convert the cartisan coordinate of basis into fractional
|
||||||
|
* ---------------------------------------------------------------------------- */
|
||||||
|
void DynMat::car2dir(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<nucell; i++){
|
||||||
|
double x[3];
|
||||||
|
x[0] = x[1] = x[2] = 0.;
|
||||||
|
for (int idim=0; idim<sysdim; idim++) x[idim] = basis[i][idim];
|
||||||
|
for (int idim=0; idim<sysdim; idim++)
|
||||||
|
basis[i][idim] = x[0]*mat[idim] + x[1]*mat[3+idim] + x[2]*mat[6+idim];
|
||||||
|
}
|
||||||
|
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------------
|
||||||
|
* private method to enforce the acoustic sum rule on force constant matrix at G
|
||||||
|
* ---------------------------------------------------------------------------- */
|
||||||
|
void DynMat::EnforceASR()
|
||||||
|
{
|
||||||
|
char str[MAXLINE];
|
||||||
|
int nasr = 20;
|
||||||
|
if (nucell <= 1) nasr = 1;
|
||||||
|
printf("\n"); for (int i=0; i<80; i++) printf("=");
|
||||||
|
|
||||||
|
// compute and display eigenvalues of Phi at gamma before ASR
|
||||||
|
if (nucell > 100){
|
||||||
|
printf("\nYour unit cell is rather large, eigenvalue evaluation takes some time...");
|
||||||
|
fflush(stdout);
|
||||||
|
}
|
||||||
|
|
||||||
|
double egvs[fftdim];
|
||||||
|
for (int i=0; i<fftdim; i++)
|
||||||
|
for (int j=0; j<fftdim; j++) DM_q[i][j] = DM_all[0][i*fftdim+j];
|
||||||
|
geteigen(egvs, 0);
|
||||||
|
printf("\nEigenvalues of Phi at gamma before enforcing ASR:\n");
|
||||||
|
for (int i=0; i<fftdim; i++){
|
||||||
|
printf("%lg ", egvs[i]);
|
||||||
|
if (i%10 == 9) printf("\n");
|
||||||
|
if (i == 99){ printf("...... (%d more skipped)\n", fftdim-100); break;}
|
||||||
|
}
|
||||||
|
printf("\n\n");
|
||||||
|
|
||||||
|
// ask for iterations to enforce ASR
|
||||||
|
printf("Please input the # of iterations to enforce ASR [%d]: ", nasr);
|
||||||
|
fgets(str,MAXLINE,stdin);
|
||||||
|
char *ptr = strtok(str," \t\n\r\f");
|
||||||
|
if (ptr) nasr = atoi(ptr);
|
||||||
|
if (nasr < 1){return; for (int i=0; i<80; i++) printf("="); printf("\n");}
|
||||||
|
|
||||||
|
for (int iit=0; iit<nasr; iit++){
|
||||||
|
// simple ASR; the resultant matrix might not be symmetric
|
||||||
|
for (int a=0; a<sysdim; a++)
|
||||||
|
for (int b=0; b<sysdim; b++){
|
||||||
|
for (int k=0; k<nucell; k++){
|
||||||
|
double sum = 0.;
|
||||||
|
for (int kp=0; kp<nucell; kp++){
|
||||||
|
int idx = (k*sysdim+a)*fftdim+kp*sysdim+b;
|
||||||
|
sum += DM_all[0][idx].r;
|
||||||
|
}
|
||||||
|
sum /= double(nucell);
|
||||||
|
for (int kp=0; kp<nucell; kp++){
|
||||||
|
int idx = (k*sysdim+a)*fftdim+kp*sysdim+b;
|
||||||
|
DM_all[0][idx].r -= sum;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// symmetrize
|
||||||
|
for (int k=0; k<nucell; k++)
|
||||||
|
for (int kp=k; kp<nucell; kp++){
|
||||||
|
double csum = 0.;
|
||||||
|
for (int a=0; a<sysdim; a++)
|
||||||
|
for (int b=0; b<sysdim; b++){
|
||||||
|
int idx = (k*sysdim+a)*fftdim+kp*sysdim+b;
|
||||||
|
int jdx = (kp*sysdim+b)*fftdim+k*sysdim+a;
|
||||||
|
csum = (DM_all[0][idx].r + DM_all[0][jdx].r )*0.5;
|
||||||
|
DM_all[0][idx].r = DM_all[0][jdx].r = csum;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// symmetric ASR
|
||||||
|
for (int a=0; a<sysdim; a++)
|
||||||
|
for (int b=0; b<sysdim; b++){
|
||||||
|
for (int k=0; k<nucell; k++){
|
||||||
|
double sum = 0.;
|
||||||
|
for (int kp=0; kp<nucell; kp++){
|
||||||
|
int idx = (k*sysdim+a)*fftdim+kp*sysdim+b;
|
||||||
|
sum += DM_all[0][idx].r;
|
||||||
|
}
|
||||||
|
sum /= double(nucell-k);
|
||||||
|
for (int kp=k; kp<nucell; kp++){
|
||||||
|
int idx = (k*sysdim+a)*fftdim+kp*sysdim+b;
|
||||||
|
int jdx = (kp*sysdim+b)*fftdim+k*sysdim+a;
|
||||||
|
DM_all[0][idx].r -= sum;
|
||||||
|
DM_all[0][jdx].r = DM_all[0][idx].r;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// compute and display eigenvalues of Phi at gamma after ASR
|
||||||
|
for (int i=0; i<fftdim; i++)
|
||||||
|
for (int j=0; j<fftdim; j++) DM_q[i][j] = DM_all[0][i*fftdim+j];
|
||||||
|
geteigen(egvs, 0);
|
||||||
|
printf("Eigenvalues of Phi at gamma after enforcing ASR:\n");
|
||||||
|
for (int i=0; i<fftdim; i++){
|
||||||
|
printf("%lg ", egvs[i]);
|
||||||
|
if (i%10 == 9) printf("\n");
|
||||||
|
if (i == 99){ printf("...... (%d more skiped)", fftdim-100); break;}
|
||||||
|
}
|
||||||
|
printf("\n"); for (int i=0; i<80; i++) printf("="); printf("\n\n");
|
||||||
|
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------------
|
||||||
|
* private method to get the reciprocal lattice vectors from the real space ones
|
||||||
|
* ---------------------------------------------------------------------------- */
|
||||||
|
void DynMat::real2rec()
|
||||||
|
{
|
||||||
|
ibasevec[0] = basevec[4]*basevec[8] - basevec[5]*basevec[7];
|
||||||
|
ibasevec[1] = basevec[5]*basevec[6] - basevec[3]*basevec[8];
|
||||||
|
ibasevec[2] = basevec[3]*basevec[7] - basevec[4]*basevec[6];
|
||||||
|
|
||||||
|
ibasevec[3] = basevec[7]*basevec[2] - basevec[8]*basevec[1];
|
||||||
|
ibasevec[4] = basevec[8]*basevec[0] - basevec[6]*basevec[2];
|
||||||
|
ibasevec[5] = basevec[6]*basevec[1] - basevec[7]*basevec[0];
|
||||||
|
|
||||||
|
ibasevec[6] = basevec[1]*basevec[5] - basevec[2]*basevec[4];
|
||||||
|
ibasevec[7] = basevec[2]*basevec[3] - basevec[0]*basevec[5];
|
||||||
|
ibasevec[8] = basevec[0]*basevec[4] - basevec[1]*basevec[3];
|
||||||
|
|
||||||
|
double vol = 0.;
|
||||||
|
for (int i=0; i<sysdim; i++) vol += ibasevec[i] * basevec[i];
|
||||||
|
vol = 8.*atan(1.)/vol;
|
||||||
|
|
||||||
|
for (int i=0; i<9; i++) ibasevec[i] *= vol;
|
||||||
|
|
||||||
|
printf("\n"); for (int i=0; i<80; i++) printf("=");
|
||||||
|
printf("\nBasis vectors of the unit cell in real space:");
|
||||||
|
for (int i=0; i<sysdim; i++){
|
||||||
|
printf("\n A%d: ", i+1);
|
||||||
|
for (int j=0; j<sysdim; j++) printf("%8.4f ", basevec[i*3+j]);
|
||||||
|
}
|
||||||
|
printf("\nBasis vectors of the corresponding reciprocal cell:");
|
||||||
|
for (int i=0; i<sysdim; i++){
|
||||||
|
printf("\n B%d: ", i+1);
|
||||||
|
for (int j=0; j<sysdim; j++) printf("%8.4f ", ibasevec[i*3+j]);
|
||||||
|
}
|
||||||
|
printf("\n"); for (int i=0; i<80; i++) printf("="); printf("\n");
|
||||||
|
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
* private method, to get the inverse of a double matrix by means of
|
||||||
|
* Gaussian-Jordan Elimination with full pivoting; square matrix required.
|
||||||
|
*
|
||||||
|
* Adapted from the Numerical Recipes in Fortran.
|
||||||
|
* --------------------------------------------------------------------*/
|
||||||
|
void DynMat::GaussJordan(int n, double *Mat)
|
||||||
|
{
|
||||||
|
int i,icol,irow,j,k,l,ll,idr,idc;
|
||||||
|
int *indxc,*indxr,*ipiv;
|
||||||
|
double big, nmjk;
|
||||||
|
double dum, pivinv;
|
||||||
|
|
||||||
|
indxc = new int[n];
|
||||||
|
indxr = new int[n];
|
||||||
|
ipiv = new int[n];
|
||||||
|
|
||||||
|
for (i=0; i<n; i++) ipiv[i] = 0;
|
||||||
|
for (i=0; i<n; i++){
|
||||||
|
big = 0.;
|
||||||
|
for (j=0; j<n; j++){
|
||||||
|
if (ipiv[j] != 1){
|
||||||
|
for (k=0; k<n; k++){
|
||||||
|
if (ipiv[k] == 0){
|
||||||
|
idr = j*n+k;
|
||||||
|
nmjk = abs(Mat[idr]);
|
||||||
|
if (nmjk >= big){
|
||||||
|
big = nmjk;
|
||||||
|
irow = j;
|
||||||
|
icol = k;
|
||||||
|
}
|
||||||
|
}else if (ipiv[k]>1){
|
||||||
|
printf("DynMat: Singular matrix in double GaussJordan!\n"); exit(1);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
ipiv[icol] += 1;
|
||||||
|
if (irow != icol){
|
||||||
|
for (l=0; l<n; l++){
|
||||||
|
idr = irow*n+l;
|
||||||
|
idc = icol*n+l;
|
||||||
|
dum = Mat[idr];
|
||||||
|
Mat[idr] = Mat[idc];
|
||||||
|
Mat[idc] = dum;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
indxr[i] = irow;
|
||||||
|
indxc[i] = icol;
|
||||||
|
idr = icol*n+icol;
|
||||||
|
if (Mat[idr] == 0.){
|
||||||
|
printf("DynMat: Singular matrix in double GaussJordan!");
|
||||||
|
exit(1);
|
||||||
|
}
|
||||||
|
|
||||||
|
pivinv = 1./ Mat[idr];
|
||||||
|
Mat[idr] = 1.;
|
||||||
|
idr = icol*n;
|
||||||
|
for (l=0; l<n; l++) Mat[idr+l] *= pivinv;
|
||||||
|
for (ll=0; ll<n; ll++){
|
||||||
|
if (ll != icol){
|
||||||
|
idc = ll*n+icol;
|
||||||
|
dum = Mat[idc];
|
||||||
|
Mat[idc] = 0.;
|
||||||
|
idc -= icol;
|
||||||
|
for (l=0; l<n; l++) Mat[idc+l] -= Mat[idr+l]*dum;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
for (l=n-1; l>=0; l--){
|
||||||
|
int rl = indxr[l];
|
||||||
|
int cl = indxc[l];
|
||||||
|
if (rl != cl){
|
||||||
|
for (k=0; k<n; k++){
|
||||||
|
idr = k*n+rl;
|
||||||
|
idc = k*n+cl;
|
||||||
|
dum = Mat[idr];
|
||||||
|
Mat[idr] = Mat[idc];
|
||||||
|
Mat[idc] = dum;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
delete []indxr;
|
||||||
|
delete []indxc;
|
||||||
|
delete []ipiv;
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------------
|
||||||
|
* Public method to reset the interpolation method
|
||||||
|
* ---------------------------------------------------------------------------- */
|
||||||
|
void DynMat::reset_interp_method()
|
||||||
|
{
|
||||||
|
interpolate->set_method();
|
||||||
|
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------------
|
||||||
|
* Private method to display help info
|
||||||
|
* ---------------------------------------------------------------------------- */
|
||||||
|
void DynMat::help()
|
||||||
|
{
|
||||||
|
ShowVersion();
|
||||||
|
printf("\nUsage:\n phana [options] [file]\n\n");
|
||||||
|
printf("Available options:\n");
|
||||||
|
printf(" -r To reset the dynamical matrix at the gamma point by a 4th order\n");
|
||||||
|
printf(" polynomial interpolation along the [100] direction; this might be\n");
|
||||||
|
printf(" useful for systems with charges. As for charged system, the dynamical\n");
|
||||||
|
printf(" matrix at Gamma is far from accurate because of the long range nature\n");
|
||||||
|
printf(" of Coulombic interaction. By reset it by interpolation, will partially\n");
|
||||||
|
printf(" elliminate the unwanted behavior, but the result is still inaccurate.\n");
|
||||||
|
printf(" By default, this is not set; and not expected for uncharged systems.\n\n");
|
||||||
|
printf(" -s This will reset the dynamical matrix at the gamma point, too, but it\n");
|
||||||
|
printf(" will also inform the code to skip all q-points that is in the vicinity\n");
|
||||||
|
printf(" of the gamma point when evaluating phonon DOS and/or phonon dispersion.\n\n");
|
||||||
|
printf(" By default, this is not set; and not expected for uncharged systems.\n\n");
|
||||||
|
printf(" -h To print out this help info.\n\n");
|
||||||
|
printf(" file To define the filename that carries the binary dynamical matrice generated\n");
|
||||||
|
printf(" by fix-phonon. If not provided, the code will ask for it.\n");
|
||||||
|
printf("\n\n");
|
||||||
|
exit(0);
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------------
|
||||||
|
* Private method to display the version info
|
||||||
|
* ---------------------------------------------------------------------------- */
|
||||||
|
void DynMat::ShowVersion()
|
||||||
|
{
|
||||||
|
printf(" ____ _ _ __ _ _ __ \n");
|
||||||
|
printf(" ( _ \\( )_( ) /__\\ ( \\( ) /__\\ \n");
|
||||||
|
printf(" )___/ ) _ ( /(__)\\ ) ( /(__)\\ \n");
|
||||||
|
printf(" (__) (_) (_)(__)(__)(_)\\_)(__)(__)\n");
|
||||||
|
printf("\nPHonon ANAlyzer for Fix-Phonon, version 2.%d, compiled on %s.\n", VERSION, __DATE__);
|
||||||
|
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
/* --------------------------------------------------------------------*/
|
||||||
66
tools/phonon/dynmat.h
Normal file
66
tools/phonon/dynmat.h
Normal file
@ -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
|
||||||
250
tools/phonon/green.cpp
Normal file
250
tools/phonon/green.cpp
Normal file
@ -0,0 +1,250 @@
|
|||||||
|
#include <stdio.h>
|
||||||
|
#include <stdlib.h>
|
||||||
|
#include <string.h>
|
||||||
|
#include <math.h>
|
||||||
|
#include "green.h"
|
||||||
|
#include <complex>
|
||||||
|
|
||||||
|
#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<sysdim; idim++){
|
||||||
|
beta[idim][0] = 0.;
|
||||||
|
for (int i=0; i<ndim; i++){vp[i] = v[i] = 0.;}
|
||||||
|
v[ipos+idim] = 1.;
|
||||||
|
|
||||||
|
// Loop on fraction levels
|
||||||
|
for (int i=0; i<nit; i++){
|
||||||
|
double sum_a = 0.;
|
||||||
|
for (int j=0; j<ndim; j++){
|
||||||
|
double sumHv = 0.;
|
||||||
|
for (int k=0; k<ndim; k++) sumHv += H[j][k]*v[k];
|
||||||
|
w[j] = sumHv - beta[idim][i]*vp[j];
|
||||||
|
sum_a += w[j]*v[j];
|
||||||
|
}
|
||||||
|
alpha[idim][i] = sum_a;
|
||||||
|
|
||||||
|
for (int k=0; k<ndim; k++) w[k] -= alpha[idim][i]*v[k];
|
||||||
|
|
||||||
|
double gamma = 0.;
|
||||||
|
for (int k=0; k<ndim; k++) gamma += w[k]*v[k];
|
||||||
|
for (int k=0; k<ndim; k++) w[k] -= gamma*v[k];
|
||||||
|
|
||||||
|
double sum_b = 0.;
|
||||||
|
for (int k=0; k<ndim; k++) sum_b += w[k]*w[k];
|
||||||
|
beta[idim][i+1] = sqrt(sum_b);
|
||||||
|
|
||||||
|
ptr = vp; vp = v; v = ptr;
|
||||||
|
double tmp = 1./beta[idim][i+1];
|
||||||
|
for (int k=0; k<ndim; k++) v[k] = w[k]*tmp;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
ptr = NULL;
|
||||||
|
delete []vp;
|
||||||
|
delete []v;
|
||||||
|
delete []w;
|
||||||
|
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
|
||||||
|
/*------------------------------------------------------------------------------
|
||||||
|
* Private method to compute the LDOS via the recusive method for system with
|
||||||
|
* many atoms
|
||||||
|
*----------------------------------------------------------------------------*/
|
||||||
|
void Green::Recursion()
|
||||||
|
{
|
||||||
|
// local variables
|
||||||
|
double *alpha_inf, *beta_inf, *xmin, *xmax;
|
||||||
|
alpha_inf = new double [sysdim];
|
||||||
|
beta_inf = new double [sysdim];
|
||||||
|
xmin = new double [sysdim];
|
||||||
|
xmax = new double [sysdim];
|
||||||
|
|
||||||
|
int nave = nit/4;
|
||||||
|
for (int idim=0; idim<sysdim; idim++){
|
||||||
|
alpha_inf[idim] = beta_inf[idim] = 0.;
|
||||||
|
|
||||||
|
for (int i= nit-nave; i<nit; i++){
|
||||||
|
alpha_inf[idim] += alpha[idim][i];
|
||||||
|
beta_inf[idim] += beta[idim][i+1];
|
||||||
|
}
|
||||||
|
|
||||||
|
alpha_inf[idim] /= double(nave);
|
||||||
|
beta_inf[idim] /= double(nave);
|
||||||
|
|
||||||
|
xmin[idim] = alpha_inf[idim] - 2.*beta_inf[idim];
|
||||||
|
xmax[idim] = alpha_inf[idim] + 2.*beta_inf[idim];
|
||||||
|
}
|
||||||
|
|
||||||
|
std::complex<double> Z, z_m_a, r_x, rec_x, rec_x_inv;
|
||||||
|
double sr, si;
|
||||||
|
|
||||||
|
double w = wmin;
|
||||||
|
for (int i=0; i<nw; i++){
|
||||||
|
double a = w*w, ax, bx;
|
||||||
|
Z = std::complex<double>(w*w, epson);
|
||||||
|
|
||||||
|
for (int idim=0; idim<sysdim; idim++){
|
||||||
|
double two_b = 2.*beta_inf[idim]*beta_inf[idim];
|
||||||
|
double rtwob = 1./two_b;
|
||||||
|
|
||||||
|
z_m_a = Z - alpha_inf[idim]*alpha_inf[idim];
|
||||||
|
|
||||||
|
if ( a < xmin[idim] ){
|
||||||
|
r_x = sqrt(-2.*two_b + z_m_a);
|
||||||
|
ax = std::real(r_x) * rtwob;
|
||||||
|
bx = std::imag(r_x) * rtwob;
|
||||||
|
} else if (a > xmax[idim]) {
|
||||||
|
r_x = sqrt(-2.*two_b + z_m_a);
|
||||||
|
ax = -std::real(r_x) * rtwob;
|
||||||
|
bx = -std::imag(r_x) * rtwob;
|
||||||
|
} else {
|
||||||
|
r_x = sqrt(2.*two_b - z_m_a);
|
||||||
|
ax = std::imag(r_x) * rtwob;
|
||||||
|
bx = -std::real(r_x) * rtwob;
|
||||||
|
}
|
||||||
|
|
||||||
|
sr = (a - alpha_inf[idim])*rtwob + ax;
|
||||||
|
si = epson * rtwob + bx;
|
||||||
|
rec_x = std::complex<double> (sr, si);
|
||||||
|
|
||||||
|
for (int j=0; j<nit; j++){
|
||||||
|
rec_x_inv = Z - alpha[idim][nit-j-1] - beta[idim][nit-j]*beta[idim][nit-j]*rec_x;
|
||||||
|
rec_x = 1./rec_x_inv;
|
||||||
|
}
|
||||||
|
ldos[i][idim] = std::imag(rec_x)*w;
|
||||||
|
}
|
||||||
|
w += dw;
|
||||||
|
}
|
||||||
|
delete []alpha_inf;
|
||||||
|
delete []beta_inf;
|
||||||
|
delete []xmin;
|
||||||
|
delete []xmax;
|
||||||
|
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
|
||||||
|
/*------------------------------------------------------------------------------
|
||||||
|
* Private method to compute the LDOS via the recusive method for system with
|
||||||
|
* a few atoms (less than NMAX)
|
||||||
|
*----------------------------------------------------------------------------*/
|
||||||
|
void Green::recursion()
|
||||||
|
{
|
||||||
|
// local variables
|
||||||
|
std::complex<double> Z, rec_x, rec_x_inv;
|
||||||
|
std::complex<double> cunit = std::complex<double>(0.,1.);
|
||||||
|
|
||||||
|
double w = wmin;
|
||||||
|
|
||||||
|
for (int i=0; i<nw; i++){
|
||||||
|
Z = std::complex<double>(w*w, epson);
|
||||||
|
|
||||||
|
for (int idim=0; idim<sysdim; idim++){
|
||||||
|
rec_x = std::complex<double>(0.,0.);
|
||||||
|
|
||||||
|
for (int j=0; j<nit; j++){
|
||||||
|
rec_x_inv = Z - alpha[idim][nit-j-1] - beta[idim][nit-j]*beta[idim][nit-j]*rec_x;
|
||||||
|
rec_x = 1./rec_x_inv;
|
||||||
|
}
|
||||||
|
ldos[i][idim] = std::imag(rec_x)*w;
|
||||||
|
}
|
||||||
|
w += dw;
|
||||||
|
}
|
||||||
|
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
|
||||||
|
/*------------------------------------------------------------------------------*/
|
||||||
25
tools/phonon/green.h
Normal file
25
tools/phonon/green.h
Normal file
@ -0,0 +1,25 @@
|
|||||||
|
#ifndef GREEN_H
|
||||||
|
#define GREEN_H
|
||||||
|
|
||||||
|
#include "memory.h"
|
||||||
|
|
||||||
|
class Green{
|
||||||
|
public:
|
||||||
|
Green(const int, const int, const int, const double, const double,
|
||||||
|
const int, const double, double **, const int, double **);
|
||||||
|
~Green();
|
||||||
|
|
||||||
|
private:
|
||||||
|
void Lanczos();
|
||||||
|
void Recursion();
|
||||||
|
void recursion();
|
||||||
|
|
||||||
|
int ndos;
|
||||||
|
double **ldos;
|
||||||
|
|
||||||
|
int natom, iatom, sysdim, nit, nw, ndim;
|
||||||
|
double dw, wmin, wmax, epson;
|
||||||
|
double **alpha, **beta, **H;
|
||||||
|
Memory *memory;
|
||||||
|
};
|
||||||
|
#endif
|
||||||
310
tools/phonon/interpolate.cpp
Normal file
310
tools/phonon/interpolate.cpp
Normal file
@ -0,0 +1,310 @@
|
|||||||
|
#include "interpolate.h"
|
||||||
|
#include "math.h"
|
||||||
|
|
||||||
|
#define MAXLINE 256
|
||||||
|
#define MIN(a,b) ((a)>(b)?(b):(a))
|
||||||
|
#define MAX(a,b) ((a)>(b)?(a):(b))
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------------
|
||||||
|
* 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; ii<Nx; ii++)
|
||||||
|
for (int jj=0; jj<Ny; jj++)
|
||||||
|
for (int kk=0; kk<Nz; kk++){
|
||||||
|
|
||||||
|
int ip = (ii+1)%Nx, jp = (jj+1)%Ny, kp = (kk+1)%Nz;
|
||||||
|
int im = (ii-1+Nx)%Nx, jm = (jj-1+Ny)%Ny, km = (kk-1+Nz)%Nz;
|
||||||
|
|
||||||
|
int p100 = (ip*Ny+jj)*Nz+kk;
|
||||||
|
int p010 = (ii*Ny+jp)*Nz+kk;
|
||||||
|
int p001 = (ii*Ny+jj)*Nz+kp;
|
||||||
|
int p110 = (ip*Ny+jp)*Nz+kk;
|
||||||
|
int p101 = (ip*Ny+jj)*Nz+kp;
|
||||||
|
int p011 = (ii*Ny+jp)*Nz+kp;
|
||||||
|
int pm00 = (im*Ny+jj)*Nz+kk;
|
||||||
|
int p0m0 = (ii*Ny+jm)*Nz+kk;
|
||||||
|
int p00m = (ii*Ny+jj)*Nz+km;
|
||||||
|
int pmm0 = (im*Ny+jm)*Nz+kk;
|
||||||
|
int pm0m = (im*Ny+jj)*Nz+km;
|
||||||
|
int p0mm = (ii*Ny+jm)*Nz+km;
|
||||||
|
int p1m0 = (ip*Ny+jm)*Nz+kk;
|
||||||
|
int p10m = (ip*Ny+jj)*Nz+km;
|
||||||
|
int p01m = (ii*Ny+jp)*Nz+km;
|
||||||
|
int pm10 = (im*Ny+jp)*Nz+kk;
|
||||||
|
int pm01 = (im*Ny+jj)*Nz+kp;
|
||||||
|
int p0m1 = (ii*Ny+jm)*Nz+kp;
|
||||||
|
int p111 = (ip*Ny+jp)*Nz+kp;
|
||||||
|
int pm11 = (im*Ny+jp)*Nz+kp;
|
||||||
|
int p1m1 = (ip*Ny+jm)*Nz+kp;
|
||||||
|
int p11m = (ip*Ny+jp)*Nz+km;
|
||||||
|
int pm1m = (im*Ny+jp)*Nz+km;
|
||||||
|
int p1mm = (ip*Ny+jm)*Nz+km;
|
||||||
|
int pmm1 = (im*Ny+jm)*Nz+kp;
|
||||||
|
int pmmm = (im*Ny+jm)*Nz+km;
|
||||||
|
|
||||||
|
for (int idim=0; idim<ndim; idim++){
|
||||||
|
Dfdx[n][idim].r = (data[p100][idim].r - data[pm00][idim].r) * half;
|
||||||
|
Dfdx[n][idim].i = (data[p100][idim].i - data[pm00][idim].i) * half;
|
||||||
|
Dfdy[n][idim].r = (data[p010][idim].r - data[p0m0][idim].r) * half;
|
||||||
|
Dfdy[n][idim].i = (data[p010][idim].i - data[p0m0][idim].i) * half;
|
||||||
|
Dfdz[n][idim].r = (data[p001][idim].r - data[p00m][idim].r) * half;
|
||||||
|
Dfdz[n][idim].i = (data[p001][idim].i - data[p00m][idim].i) * half;
|
||||||
|
D2fdxdy[n][idim].r = (data[p110][idim].r - data[p1m0][idim].r - data[pm10][idim].r + data[pmm0][idim].r) * one4;
|
||||||
|
D2fdxdy[n][idim].i = (data[p110][idim].i - data[p1m0][idim].i - data[pm10][idim].i + data[pmm0][idim].i) * one4;
|
||||||
|
D2fdxdz[n][idim].r = (data[p101][idim].r - data[p10m][idim].r - data[pm01][idim].r + data[pm0m][idim].r) * one4;
|
||||||
|
D2fdxdz[n][idim].i = (data[p101][idim].i - data[p10m][idim].i - data[pm01][idim].i + data[pm0m][idim].i) * one4;
|
||||||
|
D2fdydz[n][idim].r = (data[p011][idim].r - data[p01m][idim].r - data[p0m1][idim].r + data[p0mm][idim].r) * one4;
|
||||||
|
D2fdydz[n][idim].i = (data[p011][idim].i - data[p01m][idim].i - data[p0m1][idim].i + data[p0mm][idim].i) * one4;
|
||||||
|
D3fdxdydz[n][idim].r = (data[p111][idim].r-data[pm11][idim].r - data[p1m1][idim].r - data[p11m][idim].r +
|
||||||
|
data[p1mm][idim].r+data[pm1m][idim].r + data[pmm1][idim].r - data[pmmm][idim].r) * one8;
|
||||||
|
D3fdxdydz[n][idim].i = (data[p111][idim].i-data[pm11][idim].i - data[p1m1][idim].i - data[p11m][idim].i +
|
||||||
|
data[p1mm][idim].i+data[pm1m][idim].i + data[pmm1][idim].i - data[pmmm][idim].i) * one8;
|
||||||
|
}
|
||||||
|
n++;
|
||||||
|
}
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------------
|
||||||
|
* Deconstructor used to free memory
|
||||||
|
* ---------------------------------------------------------------------------- */
|
||||||
|
Interpolate::~Interpolate()
|
||||||
|
{
|
||||||
|
data = NULL;
|
||||||
|
memory->destroy(Dfdx);
|
||||||
|
memory->destroy(Dfdy);
|
||||||
|
memory->destroy(Dfdz);
|
||||||
|
memory->destroy(D2fdxdy);
|
||||||
|
memory->destroy(D2fdxdz);
|
||||||
|
memory->destroy(D2fdydz);
|
||||||
|
memory->destroy(D3fdxdydz);
|
||||||
|
delete memory;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------------
|
||||||
|
* Tricubic interpolation, by calling the tricubic library
|
||||||
|
* ---------------------------------------------------------------------------- */
|
||||||
|
void Interpolate::tricubic(double *qin, doublecomplex *DMq)
|
||||||
|
{
|
||||||
|
// qin should be in unit of 2*pi/L
|
||||||
|
double q[3];
|
||||||
|
for (int i=0; i<3; i++) q[i] = qin[i];
|
||||||
|
for (int i=0; i<3; i++){
|
||||||
|
while (q[i] < 0.) q[i] += 1.;
|
||||||
|
while (q[i] >= 1.) q[i] -= 1.;
|
||||||
|
}
|
||||||
|
|
||||||
|
int ix = int(q[0]*double(Nx));
|
||||||
|
int iy = int(q[1]*double(Ny));
|
||||||
|
int iz = int(q[2]*double(Nz));
|
||||||
|
double x = q[0]*double(Nx)-double(ix);
|
||||||
|
double y = q[1]*double(Ny)-double(iy);
|
||||||
|
double z = q[2]*double(Nz)-double(iz);
|
||||||
|
int ixp = (ix+1)%Nx, iyp = (iy+1)%Ny, izp = (iz+1)%Nz;
|
||||||
|
vidx[0] = (ix*Ny+iy)*Nz+iz;
|
||||||
|
vidx[1] = (ixp*Ny+iy)*Nz+iz;
|
||||||
|
vidx[2] = (ix*Ny+iyp)*Nz+iz;
|
||||||
|
vidx[3] = (ixp*Ny+iyp)*Nz+iz;
|
||||||
|
vidx[4] = (ix*Ny+iy)*Nz+izp;
|
||||||
|
vidx[5] = (ixp*Ny+iy)*Nz+izp;
|
||||||
|
vidx[6] = (ix*Ny+iyp)*Nz+izp;
|
||||||
|
vidx[7] = (ixp*Ny+iyp)*Nz+izp;
|
||||||
|
for (int i=0; i<8; i++) if (vidx[i] == 0) UseGamma = 1;
|
||||||
|
|
||||||
|
for (int idim=0; idim<ndim; idim++){
|
||||||
|
for (int i=0; i<8; i++){
|
||||||
|
f[i] = data[vidx[i]][idim].r;
|
||||||
|
dfdx[i] = Dfdx[vidx[i]][idim].r;
|
||||||
|
dfdy[i] = Dfdy[vidx[i]][idim].r;
|
||||||
|
dfdz[i] = Dfdz[vidx[i]][idim].r;
|
||||||
|
d2fdxdy[i] = D2fdxdy[vidx[i]][idim].r;
|
||||||
|
d2fdxdz[i] = D2fdxdz[vidx[i]][idim].r;
|
||||||
|
d2fdydz[i] = D2fdydz[vidx[i]][idim].r;
|
||||||
|
d3fdxdydz[i] = D3fdxdydz[vidx[i]][idim].r;
|
||||||
|
}
|
||||||
|
tricubic_get_coeff(&a[0],&f[0],&dfdx[0],&dfdy[0],&dfdz[0],&d2fdxdy[0],&d2fdxdz[0],&d2fdydz[0],&d3fdxdydz[0]);
|
||||||
|
DMq[idim].r = tricubic_eval(&a[0],x,y,z);
|
||||||
|
|
||||||
|
for (int i=0; i<8; i++){
|
||||||
|
f[i] = data[vidx[i]][idim].i;
|
||||||
|
dfdx[i] = Dfdx[vidx[i]][idim].i;
|
||||||
|
dfdy[i] = Dfdy[vidx[i]][idim].i;
|
||||||
|
dfdz[i] = Dfdz[vidx[i]][idim].i;
|
||||||
|
d2fdxdy[i] = D2fdxdy[vidx[i]][idim].i;
|
||||||
|
d2fdxdz[i] = D2fdxdz[vidx[i]][idim].i;
|
||||||
|
d2fdydz[i] = D2fdydz[vidx[i]][idim].i;
|
||||||
|
d3fdxdydz[i] = D3fdxdydz[vidx[i]][idim].i;
|
||||||
|
}
|
||||||
|
tricubic_get_coeff(&a[0],&f[0],&dfdx[0],&dfdy[0],&dfdz[0],&d2fdxdy[0],&d2fdxdz[0],&d2fdydz[0],&d3fdxdydz[0]);
|
||||||
|
DMq[idim].i = tricubic_eval(&a[0],x,y,z);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------------
|
||||||
|
* method to interpolate the DM at an arbitrary q point;
|
||||||
|
* the input q should be a vector in unit of (2pi/a 2pi/b 2pi/c).
|
||||||
|
* All q components will be rescaled into [0 1).
|
||||||
|
* ---------------------------------------------------------------------------- */
|
||||||
|
void Interpolate::trilinear(double *qin, doublecomplex *DMq)
|
||||||
|
{
|
||||||
|
// rescale q[i] into [0 1)
|
||||||
|
double q[3];
|
||||||
|
for (int i=0; i<3; i++) q[i] = qin[i];
|
||||||
|
for (int i=0; i<3; i++){
|
||||||
|
while (q[i] < 0.) q[i] += 1.;
|
||||||
|
while (q[i] >= 1.) q[i] -= 1.;
|
||||||
|
}
|
||||||
|
|
||||||
|
// find the index of the eight vertice
|
||||||
|
int ix, iy, iz, ixp, iyp, izp;
|
||||||
|
double x, y, z;
|
||||||
|
q[0] *= double(Nx);
|
||||||
|
q[1] *= double(Ny);
|
||||||
|
q[2] *= double(Nz);
|
||||||
|
|
||||||
|
ix = int(q[0])%Nx;
|
||||||
|
iy = int(q[1])%Ny;
|
||||||
|
iz = int(q[2])%Nz;
|
||||||
|
ixp = (ix+1)%Nx;
|
||||||
|
iyp = (iy+1)%Ny;
|
||||||
|
izp = (iz+1)%Nz;
|
||||||
|
x = q[0] - double(ix);
|
||||||
|
y = q[1] - double(iy);
|
||||||
|
z = q[2] - double(iz);
|
||||||
|
|
||||||
|
//--------------------------------------
|
||||||
|
vidx[0] = ((ix*Ny)+iy)*Nz + iz;
|
||||||
|
vidx[1] = ((ixp*Ny)+iy)*Nz + iz;
|
||||||
|
vidx[2] = ((ix*Ny)+iyp)*Nz + iz;
|
||||||
|
vidx[3] = ((ix*Ny)+iy)*Nz + izp;
|
||||||
|
vidx[4] = ((ixp*Ny)+iy)*Nz + izp;
|
||||||
|
vidx[5] = ((ix*Ny)+iyp)*Nz + izp;
|
||||||
|
vidx[6] = ((ixp*Ny)+iyp)*Nz + iz;
|
||||||
|
vidx[7] = ((ixp*Ny)+iyp)*Nz + izp;
|
||||||
|
for (int i=0; i<8; i++) if (vidx[i] == 0) UseGamma = 1;
|
||||||
|
|
||||||
|
double fac[8];
|
||||||
|
fac[0] = (1.-x)*(1.-y)*(1.-z);
|
||||||
|
fac[1] = x*(1.-y)*(1.-z);
|
||||||
|
fac[2] = (1.-x)*y*(1.-z);
|
||||||
|
fac[3] = (1.-x)*(1.-y)*z;
|
||||||
|
fac[4] = x*(1.-y)*z;
|
||||||
|
fac[5] = (1.-x)*y*z;
|
||||||
|
fac[6] = x*y*(1.-z);
|
||||||
|
fac[7] = x*y*z;
|
||||||
|
|
||||||
|
// now to do the interpolation
|
||||||
|
for (int idim=0; idim<ndim; idim++){
|
||||||
|
DMq[idim].r = 0.;
|
||||||
|
DMq[idim].i = 0.;
|
||||||
|
for (int i=0; i<8; i++){
|
||||||
|
DMq[idim].r += data[vidx[i]][idim].r*fac[i];
|
||||||
|
DMq[idim].i += data[vidx[i]][idim].i*fac[i];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------------
|
||||||
|
* To invoke the interpolation
|
||||||
|
* ---------------------------------------------------------------------------- */
|
||||||
|
void Interpolate::execute(double *qin, doublecomplex *DMq)
|
||||||
|
{
|
||||||
|
UseGamma = 0;
|
||||||
|
if (which == 1) // 1: tricubic
|
||||||
|
tricubic(qin, DMq);
|
||||||
|
else // otherwise: trilinear
|
||||||
|
trilinear(qin, DMq);
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------------
|
||||||
|
* Public method, to set/reset the interpolation method
|
||||||
|
* ---------------------------------------------------------------------------- */
|
||||||
|
void Interpolate::set_method()
|
||||||
|
{
|
||||||
|
char str[MAXLINE];
|
||||||
|
int im = 1;
|
||||||
|
printf("\n");for(int i=0; i<80; i++) printf("=");
|
||||||
|
printf("\nWhich interpolation method would you like to use?\n");
|
||||||
|
printf(" 1. Tricubic;\n 2. Trilinear;\n");
|
||||||
|
printf("Your choice[1]: ");
|
||||||
|
fgets(str,MAXLINE,stdin);
|
||||||
|
char *ptr = strtok(str," \t\n\r\f");
|
||||||
|
if (ptr) im = atoi(ptr);
|
||||||
|
|
||||||
|
which =2-im%2;
|
||||||
|
printf("Your selection: %d\n", which);
|
||||||
|
for(int i=0; i<80; i++) printf("="); printf("\n\n");
|
||||||
|
|
||||||
|
if (which == 1) tricubic_init();
|
||||||
|
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------------
|
||||||
|
* Public method, to reset gamma point data; in this case, the gamma point data
|
||||||
|
* will be meaningless. should only be called once.
|
||||||
|
* ---------------------------------------------------------------------------- */
|
||||||
|
void Interpolate::reset_gamma()
|
||||||
|
{
|
||||||
|
if (flag_reset_gamma) return;
|
||||||
|
flag_reset_gamma = 1;
|
||||||
|
|
||||||
|
int p1 = 1%Nx, p2 = 2%Nx;
|
||||||
|
int m1 = (Nx-1), m2 = (Nx-2+Nx)%Nx;
|
||||||
|
|
||||||
|
int ip1 = p1*Ny*Nz, ip2 = p2*Ny*Nz;
|
||||||
|
int im1 = m1*Ny*Nz, im2 = m2*Ny*Nz;
|
||||||
|
|
||||||
|
double const one6 = -1./6., two3 = 2./3.;
|
||||||
|
|
||||||
|
for (int idim=0; idim<ndim; idim++){
|
||||||
|
data[0][idim].i = 0.;
|
||||||
|
data[0][idim].r = (data[im2][idim].r + data[ip2][idim].r) * one6
|
||||||
|
+ (data[im1][idim].r + data[ip1][idim].r) * two3;
|
||||||
|
}
|
||||||
|
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
/* ---------------------------------------------------------------------------- */
|
||||||
43
tools/phonon/interpolate.h
Normal file
43
tools/phonon/interpolate.h
Normal file
@ -0,0 +1,43 @@
|
|||||||
|
#ifndef INTERPOLATION_H
|
||||||
|
#define INTERPOLATION_H
|
||||||
|
|
||||||
|
#include "stdio.h"
|
||||||
|
#include "stdlib.h"
|
||||||
|
#include "string.h"
|
||||||
|
#include "memory.h"
|
||||||
|
#include <tricubic.h>
|
||||||
|
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
|
||||||
18
tools/phonon/main.cpp
Normal file
18
tools/phonon/main.cpp
Normal file
@ -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;
|
||||||
|
}
|
||||||
51
tools/phonon/memory.cpp
Normal file
51
tools/phonon/memory.cpp
Normal file
@ -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);
|
||||||
|
}
|
||||||
430
tools/phonon/memory.h
Normal file
430
tools/phonon/memory.h
Normal file
@ -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 <typename TYPE>
|
||||||
|
TYPE *create(TYPE *&array, int n, const char *name)
|
||||||
|
{
|
||||||
|
bigint nbytes = sizeof(TYPE) * n;
|
||||||
|
array = (TYPE *) smalloc(nbytes,name);
|
||||||
|
return array;
|
||||||
|
};
|
||||||
|
|
||||||
|
template <typename TYPE>
|
||||||
|
TYPE **create(TYPE **&array, int n, const char *name) {fail(name);}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
grow or shrink 1d array
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
template <typename TYPE>
|
||||||
|
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 <typename TYPE>
|
||||||
|
TYPE **grow(TYPE **&array, int n, const char *name) {fail(name);}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
destroy a 1d array
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
template <typename TYPE>
|
||||||
|
void destroy(TYPE *array)
|
||||||
|
{
|
||||||
|
sfree(array);
|
||||||
|
};
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
create a 1d array with index from nlo to nhi inclusive
|
||||||
|
cannot grow it
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
template <typename TYPE>
|
||||||
|
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 <typename TYPE>
|
||||||
|
TYPE **create1d_offset(TYPE **&array, int nlo, int nhi, const char *name)
|
||||||
|
{fail(name);}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
destroy a 1d array with index offset
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
template <typename TYPE>
|
||||||
|
void destroy1d_offset(TYPE *array, int offset)
|
||||||
|
{
|
||||||
|
if (array) sfree(&array[offset]);
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
create a 2d array
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
template <typename TYPE>
|
||||||
|
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 <typename TYPE>
|
||||||
|
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 <typename TYPE>
|
||||||
|
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 <typename TYPE>
|
||||||
|
TYPE ***grow(TYPE ***&array, int n1, int n2, const char *name)
|
||||||
|
{fail(name);}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
destroy a 2d array
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
template <typename TYPE>
|
||||||
|
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 <typename TYPE>
|
||||||
|
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 <typename TYPE>
|
||||||
|
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 <typename TYPE>
|
||||||
|
void destroy2d_offset(TYPE **array, int offset)
|
||||||
|
{
|
||||||
|
if (array == NULL) return;
|
||||||
|
sfree(&array[0][offset]);
|
||||||
|
sfree(array);
|
||||||
|
}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
create a 3d array
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
template <typename TYPE>
|
||||||
|
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 <typename TYPE>
|
||||||
|
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 <typename TYPE>
|
||||||
|
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 <typename TYPE>
|
||||||
|
TYPE ****grow(TYPE ****&array, int n1, int n2, int n3, const char *name)
|
||||||
|
{fail(name);}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
destroy a 3d array
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
template <typename TYPE>
|
||||||
|
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 <typename TYPE>
|
||||||
|
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 <typename TYPE>
|
||||||
|
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 <typename TYPE>
|
||||||
|
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 <typename TYPE>
|
||||||
|
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 <typename TYPE>
|
||||||
|
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 <typename TYPE>
|
||||||
|
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 <typename TYPE>
|
||||||
|
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 <typename TYPE>
|
||||||
|
TYPE *****create(TYPE *****&array, int n1, int n2, int n3, int n4,
|
||||||
|
const char *name)
|
||||||
|
{fail(name);}
|
||||||
|
|
||||||
|
/* ----------------------------------------------------------------------
|
||||||
|
destroy a 4d array
|
||||||
|
------------------------------------------------------------------------- */
|
||||||
|
|
||||||
|
template <typename TYPE>
|
||||||
|
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 <typename TYPE>
|
||||||
|
bigint usage(TYPE *array, int n)
|
||||||
|
{
|
||||||
|
bigint bytes = sizeof(TYPE) * n;
|
||||||
|
return bytes;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename TYPE>
|
||||||
|
bigint usage(TYPE **array, int n1, int n2)
|
||||||
|
{
|
||||||
|
bigint bytes = sizeof(TYPE) * n1*n2;
|
||||||
|
bytes += sizeof(TYPE *) * n1;
|
||||||
|
return bytes;
|
||||||
|
}
|
||||||
|
|
||||||
|
template <typename TYPE>
|
||||||
|
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 <typename TYPE>
|
||||||
|
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
|
||||||
1266
tools/phonon/phonon.cpp
Normal file
1266
tools/phonon/phonon.cpp
Normal file
File diff suppressed because it is too large
Load Diff
59
tools/phonon/phonon.h
Normal file
59
tools/phonon/phonon.h
Normal file
@ -0,0 +1,59 @@
|
|||||||
|
#ifndef PHONON_H
|
||||||
|
#define PHONON_H
|
||||||
|
|
||||||
|
#include "stdio.h"
|
||||||
|
#include "stdlib.h"
|
||||||
|
#include <complex>
|
||||||
|
#include "dynmat.h"
|
||||||
|
#include "memory.h"
|
||||||
|
|
||||||
|
using namespace std;
|
||||||
|
|
||||||
|
class Phonon{
|
||||||
|
public:
|
||||||
|
Phonon(DynMat *);
|
||||||
|
~Phonon();
|
||||||
|
|
||||||
|
DynMat *dynmat;
|
||||||
|
|
||||||
|
private:
|
||||||
|
int nq, ndim, sysdim;
|
||||||
|
double **qpts, *wt;
|
||||||
|
double **eigs;
|
||||||
|
|
||||||
|
int ndos, nlocal, *locals;
|
||||||
|
double *dos, fmin, fmax, df, rdf;
|
||||||
|
double ***ldos;
|
||||||
|
|
||||||
|
Memory *memory;
|
||||||
|
|
||||||
|
void QMesh();
|
||||||
|
void ComputeAll();
|
||||||
|
|
||||||
|
void pdos();
|
||||||
|
void pdisp();
|
||||||
|
void therm();
|
||||||
|
|
||||||
|
void ldos_egv();
|
||||||
|
void ldos_rsgf();
|
||||||
|
void local_therm();
|
||||||
|
|
||||||
|
void dmanyq();
|
||||||
|
void vfanyq();
|
||||||
|
void DMdisp();
|
||||||
|
void vecanyq();
|
||||||
|
|
||||||
|
void 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
|
||||||
41
tools/phonon/timer.cpp
Normal file
41
tools/phonon/timer.cpp
Normal file
@ -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;
|
||||||
|
}
|
||||||
24
tools/phonon/timer.h
Normal file
24
tools/phonon/timer.h
Normal file
@ -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
|
||||||
Reference in New Issue
Block a user