ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'master' into rcb-tiled-tri

Browse Source
This commit is contained in:
Axel Kohlmeyer
2020-07-29 20:37:28 -04:00
parent cb98fa00da b00bfc82c8
commit dd001e80d2
20 changed files with 5412 additions and 4908 deletions
Download Patch File Download Diff File Expand all files Collapse all files

27
src/comm_tiled.cpp
View File

@ -916,6 +916,7 @@ void CommTiled::borders()
// swap atoms with other procs using pack_border(), unpack_border() // swap atoms with other procs using pack_border(), unpack_border()
// use Waitall() instead of Waitany() because calls to unpack_border() // use Waitall() instead of Waitany() because calls to unpack_border()
// must increment per-atom arrays in ascending order // must increment per-atom arrays in ascending order
// For the same reason, sendself unpacks must occur after recvother unpacks
if (ghost_velocity) { if (ghost_velocity) {
if (recvother[iswap]) { if (recvother[iswap]) {
@ -931,13 +932,6 @@ void CommTiled::borders()
MPI_Send(buf_send,n,MPI_DOUBLE,sendproc[iswap][m],0,world); MPI_Send(buf_send,n,MPI_DOUBLE,sendproc[iswap][m],0,world);
} }
} }
if (sendself[iswap]) {
avec->pack_border_vel(sendnum[iswap][nsend],sendlist[iswap][nsend],
buf_send,pbc_flag[iswap][nsend],
pbc[iswap][nsend]);
avec->unpack_border_vel(recvnum[iswap][nrecv],firstrecv[iswap][nrecv],
buf_send);
}
if (recvother[iswap]) { if (recvother[iswap]) {
MPI_Waitall(nrecv,requests,MPI_STATUS_IGNORE); MPI_Waitall(nrecv,requests,MPI_STATUS_IGNORE);
for (m = 0; m < nrecv; m++) for (m = 0; m < nrecv; m++)
@ -945,6 +939,13 @@ void CommTiled::borders()
&buf_recv[size_border* &buf_recv[size_border*
forward_recv_offset[iswap][m]]); forward_recv_offset[iswap][m]]);
} }
if (sendself[iswap]) {
avec->pack_border_vel(sendnum[iswap][nsend],sendlist[iswap][nsend],
buf_send,pbc_flag[iswap][nsend],
pbc[iswap][nsend]);
avec->unpack_border_vel(recvnum[iswap][nrecv],firstrecv[iswap][nrecv],
buf_send);
}
} else { } else {
if (recvother[iswap]) { if (recvother[iswap]) {
@ -960,12 +961,6 @@ void CommTiled::borders()
MPI_Send(buf_send,n,MPI_DOUBLE,sendproc[iswap][m],0,world); MPI_Send(buf_send,n,MPI_DOUBLE,sendproc[iswap][m],0,world);
} }
} }
if (sendself[iswap]) {
avec->pack_border(sendnum[iswap][nsend],sendlist[iswap][nsend],
buf_send,pbc_flag[iswap][nsend],pbc[iswap][nsend]);
avec->unpack_border(recvnum[iswap][nsend],firstrecv[iswap][nsend],
buf_send);
}
if (recvother[iswap]) { if (recvother[iswap]) {
MPI_Waitall(nrecv,requests,MPI_STATUS_IGNORE); MPI_Waitall(nrecv,requests,MPI_STATUS_IGNORE);
for (m = 0; m < nrecv; m++) for (m = 0; m < nrecv; m++)
@ -973,6 +968,12 @@ void CommTiled::borders()
&buf_recv[size_border* &buf_recv[size_border*
forward_recv_offset[iswap][m]]); forward_recv_offset[iswap][m]]);
} }
if (sendself[iswap]) {
avec->pack_border(sendnum[iswap][nsend],sendlist[iswap][nsend],
buf_send,pbc_flag[iswap][nsend],pbc[iswap][nsend]);
avec->unpack_border(recvnum[iswap][nrecv],firstrecv[iswap][nrecv],
buf_send);
}
} }
// increment ghost atoms // increment ghost atoms

58
tools/phonon/Makefile
View File

@ -1,30 +1,37 @@
.SUFFIXES : .o .cpp .SUFFIXES : .o .cpp
# compiler and flags # compiler and flags
CC = g++ -Wall CC = g++ -Wno-unused-result
LINK = $(CC) LINK = $(CC) -static
CFLAGS = -O3 $(DEBUG) $(UFLAG) CFLAGS = -O3 $(DEBUG) $(UFLAG)
#
OFLAGS = -O3 $(DEBUG) OFLAGS = -O3 $(DEBUG)
INC = $(LPKINC) $(TCINC) $(SPGINC) INC = $(LPKINC) $(TCINC) $(SPGINC) $(FFTINC)
LIB = $(LPKLIB) $(TCLIB) $(SPGLIB) LIB = $(LPKLIB) $(TCLIB) $(SPGLIB) $(FFTLIB)
#
# cLapack library needed # cLapack library needed
LPKINC = LPKINC = -I/opt/clapack/3.2.1/include
LPKLIB =-llapack LPKLIB = -L/opt/clapack/3.2.1/lib -lclapack -lblas -lf2c -lm
#
#
# spglib 1.8.2, used to get the irreducible q-points
# if UFLAG is not set, spglib won't be used.
# UFLAG = -DUseSPG # Tricubic library needed
# SPGINC = -I/opt/libs/spglib/1.8.2/include TCINC = -I/opt/tricubic/1.0/include
# SPGLIB = -L/opt/libs/spglib/1.8.2/lib -lsymspg TCLIB = -L/opt/tricubic/1.0/lib -ltricubic
# if spglib other than version 1.8.2 is used, please # spglib, used to get the irreducible q-points
# modify file phonon.cpp, instruction can be found by searching 1.8.2 # if SFLAG is not set, spglib won't be used.
SFLAG = -DUseSPG
SPGINC = -I/opt/spglib/1.9.7/include/spglib
SPGLIB = -L/opt/spglib/1.9.7/lib -lsymspg
# FFTW 3， used to deduce the force constants in real space
# if FFLAG is not set, fftw won't be used.
FFLAG = -DFFTW3
FFTINC = -I/opt/fftw/3.3.7/include
FFTLIB = -L/opt/fftw/3.3.7/lib -lfftw3
# Debug flags # Debug flags
#DEBUG = -g -DDEBUG # DEBUG = -g -DDEBUG
UFLAG = $(SFLAG) $(FFLAG)
#==================================================================== #====================================================================
ROOT = phana ROOT = phana
# executable name # executable name
@ -44,7 +51,7 @@ clean:
rm -f *.o *~ *.mod ${EXE} rm -f *.o *~ *.mod ${EXE}
tar: tar:
rm -f ${ROOT}.tar; tar -czvf ${ROOT}.tar.gz *.cpp *.h Makefile README rm -f ${ROOT}.tar.gz; tar -czvf ${ROOT}.tar.gz *.cpp *.h Makefile README
ver: ver:
@echo "#define VERSION `git log|grep '^commit'|wc -l`" > version.h @echo "#define VERSION `git log|grep '^commit'|wc -l`" > version.h
@ -58,16 +65,3 @@ ver:
$(CC) $(CFLAGS) -c $< $(CC) $(CFLAGS) -c $<
.cpp.o: .cpp.o:
$(CC) $(CFLAGS) $(INC) -c $< $(CC) $(CFLAGS) $(INC) -c $<
#====================================================================
# dependencies
disp.o: disp.cpp phonon.h dynmat.h memory.h interpolate.h green.h timer.h \
global.h
dynmat.o: dynmat.cpp dynmat.h memory.h interpolate.h version.h global.h
green.o: green.cpp green.h memory.h global.h
interpolate.o: interpolate.cpp interpolate.h memory.h global.h
main.o: main.cpp dynmat.h memory.h interpolate.h phonon.h
memory.o: memory.cpp memory.h
phonon.o: phonon.cpp phonon.h dynmat.h memory.h interpolate.h green.h \
timer.h global.h
timer.o: timer.cpp timer.h

27
tools/phonon/README
View File

@ -5,9 +5,13 @@
analyse the phonon related information. analyse the phonon related information.
#------------------------------------------------------------------------------- #-------------------------------------------------------------------------------
1. Dependencies 1. Dependencies
The LAPACK library is needed to solve the eigen problems. The clapack library is needed to solve the eigen problems,
http://www.netlib.org/lapack/ which could be downloaded from:
Intel MKL can be used as well. http://www.netlib.org/clapack/
The tricubic library is also needed to do tricubic interpolations,
which could now be obtained from:
https://github.com/nbigaouette/libtricubic/
The spglib is optionally needed, enabling one to evaluate the The spglib is optionally needed, enabling one to evaluate the
phonon density of states or vibrational thermal properties phonon density of states or vibrational thermal properties
@ -16,6 +20,12 @@
automatic mode. Currently, the 1.8.3 version of spglib is used. automatic mode. Currently, the 1.8.3 version of spglib is used.
It can be obtained from: It can be obtained from:
http://spglib.sourceforge.net/ http://spglib.sourceforge.net/
FFTW 3 might also be needed if you would like to interface with
phonopy: necessary input files for phonopy will be prepared so
that you can make use of the functions provided by phonopy.
FFTW 3 can be downloaded from:
http://www.fftw.org
2. Compilation 2. Compilation
To compile the code, one needs therefore to install the above To compile the code, one needs therefore to install the above
@ -26,17 +36,16 @@
3. Unit system 3. Unit system
The units of the output frequencies by this code is THz for The units of the output frequencies by this code is THz for
LAMMPS units "real", "si", "metal", and "cgs"; in these cases, LAMMPS units "real", "si", "metal", "cgs", "micro", "nano";
the frequencies are $\nu$ instead of $\omega$. in these cases, the frequencies are $\nu$ instead of $\omega$.
4. Updates 4. Updates
For updates of phana, please check: For updates of phana, please check:
http://nes.sjtu.edu.cn/english/research/software.htm https://github.com/lingtikong/phana.git
or
https://github.com/lingtikong/phana.git
5. Bug report 5. Bug report
If any bug found, please drop a line to: konglt(at)sjtu.edu.cn If any bug found, please drop a line to: konglt(at)sjtu.edu.cn
#------------------------------------------------------------------------------- #-------------------------------------------------------------------------------
Author: Ling-Ti Kong, konglt(at)sjtu.edu.cn Author: Ling-Ti Kong, konglt(at)sjtu.edu.cn
Oct 2015 May 2020

2969
tools/phonon/disp.cpp
View File

File diff suppressed because it is too large Load Diff

1001
tools/phonon/dynmat.cpp
View File

File diff suppressed because it is too large Load Diff

10
tools/phonon/dynmat.h
View File

@ -7,8 +7,6 @@
#include "memory.h" #include "memory.h"
#include "interpolate.h" #include "interpolate.h"
using namespace std;
class DynMat { class DynMat {
public: public:
@ -17,7 +15,7 @@ public:
int nx, ny, nz, nucell; int nx, ny, nz, nucell;
int sysdim, fftdim; int sysdim, fftdim;
double eml2f; double eml2f, eml2fc;
char *funit; char *funit;
void getDMq(double *); void getDMq(double *);
@ -30,7 +28,9 @@ public:
doublecomplex **DM_q; doublecomplex **DM_q;
int flag_latinfo; int flag_latinfo;
int npt, fftdim2;
double Tmeasure, basevec[9], ibasevec[9]; double Tmeasure, basevec[9], ibasevec[9];
double *M_inv_sqrt;
double **basis; double **basis;
int *attyp; int *attyp;
@ -40,14 +40,12 @@ private:
Interpolate *interpolate; Interpolate *interpolate;
Memory *memory; Memory *memory;
int npt, fftdim2;
int nasr; int nasr;
void EnforceASR(); void EnforceASR();
char *binfile, *dmfile; char *binfile, *dmfile;
double boltz, q[3]; double boltz, q[3];
double *M_inv_sqrt;
doublecomplex **DM_all; doublecomplex **DM_all;
@ -56,6 +54,8 @@ private:
void GaussJordan(int, double *); void GaussJordan(int, double *);
void help(); void help();
void ShowInfo();
void ShowVersion(); void ShowVersion();
void Define_Conversion_Factor();
}; };
#endif #endif

346
tools/phonon/green.cpp
View File

@ -38,36 +38,36 @@
Green::Green(const int ntm, const int sdim, const int niter, const double min, const double max, 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 int ndos, const double eps, double **Hessian, const int itm, double **lpdos)
{ {
const double tpi = 8.*atan(1.); const double tpi = 8.*atan(1.);
natom = ntm; sysdim = sdim; nit = niter; epson = eps; natom = ntm; sysdim = sdim; nit = niter; epson = eps;
wmin = min*tpi; wmax = max*tpi; nw = ndos + (ndos+1)%2; wmin = min*tpi; wmax = max*tpi; nw = ndos + (ndos+1)%2;
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;
} }
ndim = natom * sysdim; ndim = natom * sysdim;
if (nit < 1){printf("\nError: Wrong input of maximum iterations!\n"); return;} 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 (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;} if (nw < 1){printf("\nError: Wrong input of points in LDOS!\n"); return;}
// initialize variables and allocate local memories // initialize variables and allocate local memories
dw = (wmax - wmin)/double(nw-1); dw = (wmax - wmin)/double(nw-1);
memory->create(alpha, sysdim,nit, "Green_Green:alpha"); memory->create(alpha, sysdim,nit, "Green_Green:alpha");
memory->create(beta, sysdim,nit+1,"Green_Green:beta"); memory->create(beta, sysdim,nit+1,"Green_Green:beta");
//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();
// Get the inverser of the treated hessian by continued fractional method // Get the inverser of the treated hessian by continued fractional method
Recursion(); Recursion();
return; return;
} }
/*------------------------------------------------------------------------------ /*------------------------------------------------------------------------------
@ -75,14 +75,14 @@ return;
*----------------------------------------------------------------------------*/ *----------------------------------------------------------------------------*/
Green::~Green() Green::~Green()
{ {
H = NULL; H = NULL;
ldos = NULL; ldos = NULL;
memory->destroy(alpha); memory->destroy(alpha);
memory->destroy(beta); memory->destroy(beta);
delete memory; delete memory;
return; return;
} }
/*------------------------------------------------------------------------------ /*------------------------------------------------------------------------------
@ -90,159 +90,159 @@ return;
*----------------------------------------------------------------------------*/ *----------------------------------------------------------------------------*/
void Green::Lanczos() void Green::Lanczos()
{ {
double *vp, *v, *w, *ptr; double *vp, *v, *w, *ptr;
vp = new double [ndim]; vp = new double [ndim];
v = new double [ndim]; v = new double [ndim];
w = new double [ndim]; w = new double [ndim];
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;
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;
} }
alpha[idim][i] = sum_a; }
for (int k = 0; k < ndim; ++k) w[k] -= alpha[idim][i]*v[k]; ptr = NULL;
delete []vp;
double gamma = 0.; delete []v;
for (int k = 0; k < ndim; ++k) gamma += w[k]*v[k]; delete []w;
for (int k = 0; k < ndim; ++k) w[k] -= gamma*v[k];
return;
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 recursive method for system with * Private method to compute the LDOS via the recusive method for system with
* many atoms * many atoms
*----------------------------------------------------------------------------*/ *----------------------------------------------------------------------------*/
void Green::Recursion() void Green::Recursion()
{ {
// local variables // local variables
double *alpha_inf, *beta_inf, *xmin, *xmax; double *alpha_inf, *beta_inf, *xmin, *xmax;
alpha_inf = new double [sysdim]; alpha_inf = new double [sysdim];
beta_inf = new double [sysdim]; beta_inf = new double [sysdim];
xmin = new double [sysdim]; xmin = new double [sysdim];
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];
}
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; alpha_inf[idim] /= double(nave);
si = epson * rtwob + bx; beta_inf[idim] /= double(nave);
rec_x = std::complex<double> (sr, si);
xmin[idim] = alpha_inf[idim] - 2.*beta_inf[idim];
for (int j = 0; j < nit; ++j){ xmax[idim] = alpha_inf[idim] + 2.*beta_inf[idim];
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;
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;
} }
ldos[i][idim] = std::imag(rec_x)*w; w += dw;
} }
w += dw; delete []alpha_inf;
} delete []beta_inf;
delete []alpha_inf; delete []xmin;
delete []beta_inf; delete []xmax;
delete []xmin;
delete []xmax; return;
return;
} }
/*------------------------------------------------------------------------------ /*------------------------------------------------------------------------------
* Private method to compute the LDOS via the recursive method for system with * Private method to compute the LDOS via the recusive method for system with
* a few atoms (less than NMAX) * a few atoms (less than NMAX)
*----------------------------------------------------------------------------*/ *----------------------------------------------------------------------------*/
void Green::recursion() void Green::recursion()
{ {
// local variables // local variables
std::complex<double> Z, rec_x, rec_x_inv; std::complex<double> Z, rec_x, rec_x_inv;
std::complex<double> cunit = std::complex<double>(0.,1.);
double w = wmin;
double w = wmin;
for (int i = 0; i < nw; ++i){
Z = std::complex<double>(w*w, epson); 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 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; for (int j = 0; j < nit; ++j){
rec_x = 1./rec_x_inv; 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;
} }
ldos[i][idim] = std::imag(rec_x)*w; w += dw;
} }
w += dw; return;
}
return;
} }
/*------------------------------------------------------------------------------*/ /*------------------------------------------------------------------------------*/

28
tools/phonon/green.h
View File

@ -5,21 +5,21 @@
class Green{ class Green{
public: public:
Green(const int, const int, const int, const double, const double, Green(const int, const int, const int, const double, const double,
const int, const double, double **, const int, double **); const int, const double, double **, const int, double **);
~Green(); ~Green();
private: private:
void Lanczos(); void Lanczos();
void Recursion(); void Recursion();
void recursion(); void recursion();
int ndos; int ndos;
double **ldos; double **ldos;
int natom, iatom, sysdim, nit, nw, ndim; int natom, iatom, sysdim, nit, nw, ndim;
double dw, wmin, wmax, epson; double dw, wmin, wmax, epson;
double **alpha, **beta, **H; double **alpha, **beta, **H;
Memory *memory; Memory *memory;
}; };
#endif #endif

613
tools/phonon/interpolate.cpp
View File

@ -1,144 +1,26 @@
#include "interpolate.h" #include "interpolate.h"
#include <math.h> #include "math.h"
#include "global.h" #include "global.h"
///////////////////////
// tricubic library code
static int A[64][64] = {
{ 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{-3, 3, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{ 2,-2, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{ 9,-9,-9, 9, 0, 0, 0, 0, 6, 3,-6,-3, 0, 0, 0, 0, 6,-6, 3,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 2, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{-6, 6, 6,-6, 0, 0, 0, 0,-3,-3, 3, 3, 0, 0, 0, 0,-4, 4,-2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-2,-1,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{ 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0, 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{-6, 6, 6,-6, 0, 0, 0, 0,-4,-2, 4, 2, 0, 0, 0, 0,-3, 3,-3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-1,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{ 4,-4,-4, 4, 0, 0, 0, 0, 2, 2,-2,-2, 0, 0, 0, 0, 2,-2, 2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0, 0, 0, 0, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9,-9,-9, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 3,-6,-3, 0, 0, 0, 0, 6,-6, 3,-3, 0, 0, 0, 0, 4, 2, 2, 1, 0, 0, 0, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-6, 6, 6,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3,-3, 3, 3, 0, 0, 0, 0,-4, 4,-2, 2, 0, 0, 0, 0,-2,-2,-1,-1, 0, 0, 0, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-6, 6, 6,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4,-2, 4, 2, 0, 0, 0, 0,-3, 3,-3, 3, 0, 0, 0, 0,-2,-1,-2,-1, 0, 0, 0, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4,-4,-4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2,-2,-2, 0, 0, 0, 0, 2,-2, 2,-2, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0},
{-3, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{ 9,-9, 0, 0,-9, 9, 0, 0, 6, 3, 0, 0,-6,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6,-6, 0, 0, 3,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 2, 0, 0, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{-6, 6, 0, 0, 6,-6, 0, 0,-3,-3, 0, 0, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4, 4, 0, 0,-2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-2, 0, 0,-1,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0, 0, 0,-1, 0, 0, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9,-9, 0, 0,-9, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 3, 0, 0,-6,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6,-6, 0, 0, 3,-3, 0, 0, 4, 2, 0, 0, 2, 1, 0, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-6, 6, 0, 0, 6,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3,-3, 0, 0, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4, 4, 0, 0,-2, 2, 0, 0,-2,-2, 0, 0,-1,-1, 0, 0},
{ 9, 0,-9, 0,-9, 0, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 3, 0,-6, 0,-3, 0, 6, 0,-6, 0, 3, 0,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0, 2, 0, 2, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0, 9, 0,-9, 0,-9, 0, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 3, 0,-6, 0,-3, 0, 6, 0,-6, 0, 3, 0,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0, 2, 0, 2, 0, 1, 0},
{-27,27,27,-27,27,-27,-27,27,-18,-9,18, 9,18, 9,-18,-9,-18,18,-9, 9,18,-18, 9,-9,-18,18,18,-18,-9, 9, 9,-9,-12,-6,-6,-3,12, 6, 6, 3,-12,-6,12, 6,-6,-3, 6, 3,-12,12,-6, 6,-6, 6,-3, 3,-8,-4,-4,-2,-4,-2,-2,-1},
{18,-18,-18,18,-18,18,18,-18, 9, 9,-9,-9,-9,-9, 9, 9,12,-12, 6,-6,-12,12,-6, 6,12,-12,-12,12, 6,-6,-6, 6, 6, 6, 3, 3,-6,-6,-3,-3, 6, 6,-6,-6, 3, 3,-3,-3, 8,-8, 4,-4, 4,-4, 2,-2, 4, 4, 2, 2, 2, 2, 1, 1},
{-6, 0, 6, 0, 6, 0,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0,-3, 0, 3, 0, 3, 0,-4, 0, 4, 0,-2, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-2, 0,-1, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0,-6, 0, 6, 0, 6, 0,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0,-3, 0, 3, 0, 3, 0,-4, 0, 4, 0,-2, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-2, 0,-1, 0,-1, 0},
{18,-18,-18,18,-18,18,18,-18,12, 6,-12,-6,-12,-6,12, 6, 9,-9, 9,-9,-9, 9,-9, 9,12,-12,-12,12, 6,-6,-6, 6, 6, 3, 6, 3,-6,-3,-6,-3, 8, 4,-8,-4, 4, 2,-4,-2, 6,-6, 6,-6, 3,-3, 3,-3, 4, 2, 4, 2, 2, 1, 2, 1},
{-12,12,12,-12,12,-12,-12,12,-6,-6, 6, 6, 6, 6,-6,-6,-6, 6,-6, 6, 6,-6, 6,-6,-8, 8, 8,-8,-4, 4, 4,-4,-3,-3,-3,-3, 3, 3, 3, 3,-4,-4, 4, 4,-2,-2, 2, 2,-4, 4,-4, 4,-2, 2,-2, 2,-2,-2,-2,-2,-1,-1,-1,-1},
{ 2, 0, 0, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{-6, 6, 0, 0, 6,-6, 0, 0,-4,-2, 0, 0, 4, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0,-3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{ 4,-4, 0, 0,-4, 4, 0, 0, 2, 2, 0, 0,-2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-6, 6, 0, 0, 6,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4,-2, 0, 0, 4, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0,-3, 3, 0, 0,-2,-1, 0, 0,-2,-1, 0, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4,-4, 0, 0,-4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 0, 0,-2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 2,-2, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0},
{-6, 0, 6, 0, 6, 0,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4, 0,-2, 0, 4, 0, 2, 0,-3, 0, 3, 0,-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0,-2, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0,-6, 0, 6, 0, 6, 0,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4, 0,-2, 0, 4, 0, 2, 0,-3, 0, 3, 0,-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0,-2, 0,-1, 0},
{18,-18,-18,18,-18,18,18,-18,12, 6,-12,-6,-12,-6,12, 6,12,-12, 6,-6,-12,12,-6, 6, 9,-9,-9, 9, 9,-9,-9, 9, 8, 4, 4, 2,-8,-4,-4,-2, 6, 3,-6,-3, 6, 3,-6,-3, 6,-6, 3,-3, 6,-6, 3,-3, 4, 2, 2, 1, 4, 2, 2, 1},
{-12,12,12,-12,12,-12,-12,12,-6,-6, 6, 6, 6, 6,-6,-6,-8, 8,-4, 4, 8,-8, 4,-4,-6, 6, 6,-6,-6, 6, 6,-6,-4,-4,-2,-2, 4, 4, 2, 2,-3,-3, 3, 3,-3,-3, 3, 3,-4, 4,-2, 2,-4, 4,-2, 2,-2,-2,-1,-1,-2,-2,-1,-1},
{ 4, 0,-4, 0,-4, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 2, 0,-2, 0,-2, 0, 2, 0,-2, 0, 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{ 0, 0, 0, 0, 0, 0, 0, 0, 4, 0,-4, 0,-4, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 2, 0,-2, 0,-2, 0, 2, 0,-2, 0, 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0},
{-12,12,12,-12,12,-12,-12,12,-8,-4, 8, 4, 8, 4,-8,-4,-6, 6,-6, 6, 6,-6, 6,-6,-6, 6, 6,-6,-6, 6, 6,-6,-4,-2,-4,-2, 4, 2, 4, 2,-4,-2, 4, 2,-4,-2, 4, 2,-3, 3,-3, 3,-3, 3,-3, 3,-2,-1,-2,-1,-2,-1,-2,-1},
{ 8,-8,-8, 8,-8, 8, 8,-8, 4, 4,-4,-4,-4,-4, 4, 4, 4,-4, 4,-4,-4, 4,-4, 4, 4,-4,-4, 4, 4,-4,-4, 4, 2, 2, 2, 2,-2,-2,-2,-2, 2, 2,-2,-2, 2, 2,-2,-2, 2,-2, 2,-2, 2,-2, 2,-2, 1, 1, 1, 1, 1, 1, 1, 1}};
static int ijk2n(int i, int j, int k) {
return(i+4*j+16*k);
}
/* ---------------------------------------------------------------------------- */
static void tricubic_get_coeff_stacked(double a[64], double x[64]) {
int i,j;
for (i=0;i<64;i++) {
a[i]=(double)(0.0);
for (j=0;j<64;j++) {
a[i]+=A[i][j]*x[j];
}
}
}
static void tricubic_get_coeff(double a[64], double f[8], double dfdx[8], double dfdy[8], double dfdz[8], double d2fdxdy[8], double d2fdxdz[8], double d2fdydz[8], double d3fdxdydz[8]) {
int i;
double x[64];
for (i=0;i<8;i++) {
x[0+i]=f[i];
x[8+i]=dfdx[i];
x[16+i]=dfdy[i];
x[24+i]=dfdz[i];
x[32+i]=d2fdxdy[i];
x[40+i]=d2fdxdz[i];
x[48+i]=d2fdydz[i];
x[56+i]=d3fdxdydz[i];
}
tricubic_get_coeff_stacked(a,x);
}
static double tricubic_eval(double a[64], double x, double y, double z) {
int i,j,k;
double ret=(double)(0.0);
/* TRICUBIC EVAL
This is the short version of tricubic_eval. It is used to compute
the value of the function at a given point (x,y,z). To compute
partial derivatives of f, use the full version with the extra args.
*/
for (i=0;i<4;i++) {
for (j=0;j<4;j++) {
for (k=0;k<4;k++) {
ret+=a[ijk2n(i,j,k)]*pow(x,i)*pow(y,j)*pow(z,k);
}
}
}
return(ret);
}
/* ---------------------------------------------------------------------------- /* ----------------------------------------------------------------------------
* Constructor used to get info from caller, and prepare other necessary data * Constructor used to get info from caller, and prepare other necessary data
* ---------------------------------------------------------------------------- */ * ---------------------------------------------------------------------------- */
Interpolate::Interpolate(int nx, int ny, int nz, int ndm, doublecomplex **DM) Interpolate::Interpolate(int nx, int ny, int nz, int ndm, doublecomplex **DM)
{ {
Nx = nx; Nx = nx;
Ny = ny; Ny = ny;
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;
data = DM; data = DM;
Dfdx = Dfdy = Dfdz = D2fdxdy = D2fdxdz = D2fdydz = D3fdxdydz = NULL; Dfdx = Dfdy = Dfdz = D2fdxdy = D2fdxdz = D2fdydz = D3fdxdydz = NULL;
flag_reset_gamma = flag_allocated_dfs = 0; flag_reset_gamma = flag_allocated_dfs = 0;
return; return;
} }
/* ---------------------------------------------------------------------------- /* ----------------------------------------------------------------------------
@ -146,77 +28,76 @@ return;
* ---------------------------------------------------------------------------- */ * ---------------------------------------------------------------------------- */
void Interpolate::tricubic_init() void Interpolate::tricubic_init()
{ {
// prepare necessary data for tricubic // prepare necessary data for tricubic
if (flag_allocated_dfs == 0){ if (flag_allocated_dfs == 0){
memory->create(Dfdx, Npt, ndim, "Interpolate_Interpolate:Dfdx"); memory->create(Dfdx, Npt, ndim, "Interpolate_Interpolate:Dfdx");
memory->create(Dfdy, Npt, ndim, "Interpolate_Interpolate:Dfdy"); memory->create(Dfdy, Npt, ndim, "Interpolate_Interpolate:Dfdy");
memory->create(Dfdz, Npt, ndim, "Interpolate_Interpolate:Dfdz"); memory->create(Dfdz, Npt, ndim, "Interpolate_Interpolate:Dfdz");
memory->create(D2fdxdy, Npt, ndim, "Interpolate_Interpolate:D2fdxdy"); memory->create(D2fdxdy, Npt, ndim, "Interpolate_Interpolate:D2fdxdy");
memory->create(D2fdxdz, Npt, ndim, "Interpolate_Interpolate:D2fdxdz"); memory->create(D2fdxdz, Npt, ndim, "Interpolate_Interpolate:D2fdxdz");
memory->create(D2fdydz, Npt, ndim, "Interpolate_Interpolate:D2fdydz"); memory->create(D2fdydz, Npt, ndim, "Interpolate_Interpolate:D2fdydz");
memory->create(D3fdxdydz, Npt, ndim, "Interpolate_Interpolate:D2fdxdydz"); memory->create(D3fdxdydz, Npt, ndim, "Interpolate_Interpolate:D2fdxdydz");
flag_allocated_dfs = 1;
}
flag_allocated_dfs = 1; // get the derivatives
} int n=0;
const double half = 0.5, one4 = 0.25, one8 = 0.125;
// get the derivatives for (int ii = 0; ii < Nx; ++ii)
int n=0; for (int jj = 0; jj < Ny; ++jj)
const double half = 0.5, one4 = 0.25, one8 = 0.125; for (int kk = 0; kk < Nz; ++kk){
for (int ii = 0; ii < Nx; ++ii) int ip = (ii+1)%Nx, jp = (jj+1)%Ny, kp = (kk+1)%Nz;
for (int jj = 0; jj < Ny; ++jj) int im = (ii-1+Nx)%Nx, jm = (jj-1+Ny)%Ny, km = (kk-1+Nz)%Nz;
for (int kk = 0; kk < Nz; ++kk){
int p100 = (ip*Ny+jj)*Nz+kk;
int ip = (ii+1)%Nx, jp = (jj+1)%Ny, kp = (kk+1)%Nz; int p010 = (ii*Ny+jp)*Nz+kk;
int im = (ii-1+Nx)%Nx, jm = (jj-1+Ny)%Ny, km = (kk-1+Nz)%Nz; int p001 = (ii*Ny+jj)*Nz+kp;
int p110 = (ip*Ny+jp)*Nz+kk;
int p100 = (ip*Ny+jj)*Nz+kk; int p101 = (ip*Ny+jj)*Nz+kp;
int p010 = (ii*Ny+jp)*Nz+kk; int p011 = (ii*Ny+jp)*Nz+kp;
int p001 = (ii*Ny+jj)*Nz+kp; int pm00 = (im*Ny+jj)*Nz+kk;
int p110 = (ip*Ny+jp)*Nz+kk; int p0m0 = (ii*Ny+jm)*Nz+kk;
int p101 = (ip*Ny+jj)*Nz+kp; int p00m = (ii*Ny+jj)*Nz+km;
int p011 = (ii*Ny+jp)*Nz+kp; int pmm0 = (im*Ny+jm)*Nz+kk;
int pm00 = (im*Ny+jj)*Nz+kk; int pm0m = (im*Ny+jj)*Nz+km;
int p0m0 = (ii*Ny+jm)*Nz+kk; int p0mm = (ii*Ny+jm)*Nz+km;
int p00m = (ii*Ny+jj)*Nz+km; int p1m0 = (ip*Ny+jm)*Nz+kk;
int pmm0 = (im*Ny+jm)*Nz+kk; int p10m = (ip*Ny+jj)*Nz+km;
int pm0m = (im*Ny+jj)*Nz+km; int p01m = (ii*Ny+jp)*Nz+km;
int p0mm = (ii*Ny+jm)*Nz+km; int pm10 = (im*Ny+jp)*Nz+kk;
int p1m0 = (ip*Ny+jm)*Nz+kk; int pm01 = (im*Ny+jj)*Nz+kp;
int p10m = (ip*Ny+jj)*Nz+km; int p0m1 = (ii*Ny+jm)*Nz+kp;
int p01m = (ii*Ny+jp)*Nz+km; int p111 = (ip*Ny+jp)*Nz+kp;
int pm10 = (im*Ny+jp)*Nz+kk; int pm11 = (im*Ny+jp)*Nz+kp;
int pm01 = (im*Ny+jj)*Nz+kp; int p1m1 = (ip*Ny+jm)*Nz+kp;
int p0m1 = (ii*Ny+jm)*Nz+kp; int p11m = (ip*Ny+jp)*Nz+km;
int p111 = (ip*Ny+jp)*Nz+kp; int pm1m = (im*Ny+jp)*Nz+km;
int pm11 = (im*Ny+jp)*Nz+kp; int p1mm = (ip*Ny+jm)*Nz+km;
int p1m1 = (ip*Ny+jm)*Nz+kp; int pmm1 = (im*Ny+jm)*Nz+kp;
int p11m = (ip*Ny+jp)*Nz+km; int pmmm = (im*Ny+jm)*Nz+km;
int pm1m = (im*Ny+jp)*Nz+km;
int p1mm = (ip*Ny+jm)*Nz+km; for (int idim=0; idim<ndim; idim++){
int pmm1 = (im*Ny+jm)*Nz+kp; Dfdx[n][idim].r = (data[p100][idim].r - data[pm00][idim].r) * half;
int pmmm = (im*Ny+jm)*Nz+km; 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;
for (int idim=0; idim<ndim; idim++){ Dfdy[n][idim].i = (data[p010][idim].i - data[p0m0][idim].i) * half;
Dfdx[n][idim].r = (data[p100][idim].r - data[pm00][idim].r) * half; Dfdz[n][idim].r = (data[p001][idim].r - data[p00m][idim].r) * half;
Dfdx[n][idim].i = (data[p100][idim].i - data[pm00][idim].i) * half; Dfdz[n][idim].i = (data[p001][idim].i - data[p00m][idim].i) * half;
Dfdy[n][idim].r = (data[p010][idim].r - data[p0m0][idim].r) * half; D2fdxdy[n][idim].r = (data[p110][idim].r - data[p1m0][idim].r - data[pm10][idim].r + data[pmm0][idim].r) * one4;
Dfdy[n][idim].i = (data[p010][idim].i - data[p0m0][idim].i) * half; D2fdxdy[n][idim].i = (data[p110][idim].i - data[p1m0][idim].i - data[pm10][idim].i + data[pmm0][idim].i) * one4;
Dfdz[n][idim].r = (data[p001][idim].r - data[p00m][idim].r) * half; D2fdxdz[n][idim].r = (data[p101][idim].r - data[p10m][idim].r - data[pm01][idim].r + data[pm0m][idim].r) * one4;
Dfdz[n][idim].i = (data[p001][idim].i - data[p00m][idim].i) * half; D2fdxdz[n][idim].i = (data[p101][idim].i - data[p10m][idim].i - data[pm01][idim].i + data[pm0m][idim].i) * one4;
D2fdxdy[n][idim].r = (data[p110][idim].r - data[p1m0][idim].r - data[pm10][idim].r + data[pmm0][idim].r) * one4; D2fdydz[n][idim].r = (data[p011][idim].r - data[p01m][idim].r - data[p0m1][idim].r + data[p0mm][idim].r) * one4;
D2fdxdy[n][idim].i = (data[p110][idim].i - data[p1m0][idim].i - data[pm10][idim].i + data[pmm0][idim].i) * one4; D2fdydz[n][idim].i = (data[p011][idim].i - data[p01m][idim].i - data[p0m1][idim].i + data[p0mm][idim].i) * one4;
D2fdxdz[n][idim].r = (data[p101][idim].r - data[p10m][idim].r - data[pm01][idim].r + data[pm0m][idim].r) * one4; D3fdxdydz[n][idim].r = (data[p111][idim].r-data[pm11][idim].r - data[p1m1][idim].r - data[p11m][idim].r +
D2fdxdz[n][idim].i = (data[p101][idim].i - data[p10m][idim].i - data[pm01][idim].i + data[pm0m][idim].i) * one4; data[p1mm][idim].r+data[pm1m][idim].r + data[pmm1][idim].r - data[pmmm][idim].r) * one8;
D2fdydz[n][idim].r = (data[p011][idim].r - data[p01m][idim].r - data[p0m1][idim].r + data[p0mm][idim].r) * one4; D3fdxdydz[n][idim].i = (data[p111][idim].i-data[pm11][idim].i - data[p1m1][idim].i - data[p11m][idim].i +
D2fdydz[n][idim].i = (data[p011][idim].i - data[p01m][idim].i - data[p0m1][idim].i + data[p0mm][idim].i) * one4; data[p1mm][idim].i+data[pm1m][idim].i + data[pmm1][idim].i - data[pmmm][idim].i) * one8;
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; n++;
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; return;
}
n++;
}
return;
} }
/* ---------------------------------------------------------------------------- /* ----------------------------------------------------------------------------
@ -224,15 +105,17 @@ return;
* ---------------------------------------------------------------------------- */ * ---------------------------------------------------------------------------- */
Interpolate::~Interpolate() Interpolate::~Interpolate()
{ {
data = NULL; data = NULL;
memory->destroy(Dfdx); memory->destroy(Dfdx);
memory->destroy(Dfdy); memory->destroy(Dfdy);
memory->destroy(Dfdz); memory->destroy(Dfdz);
memory->destroy(D2fdxdy); memory->destroy(D2fdxdy);
memory->destroy(D2fdxdz); memory->destroy(D2fdxdz);
memory->destroy(D2fdydz); memory->destroy(D2fdydz);
memory->destroy(D3fdxdydz); memory->destroy(D3fdxdydz);
delete memory; delete memory;
return;
} }
/* ---------------------------------------------------------------------------- /* ----------------------------------------------------------------------------
@ -240,60 +123,59 @@ Interpolate::~Interpolate()
* ---------------------------------------------------------------------------- */ * ---------------------------------------------------------------------------- */
void Interpolate::tricubic(double *qin, doublecomplex *DMq) 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.;
} }
int ix = int(q[0]*double(Nx)); int ix = int(q[0]*double(Nx));
int iy = int(q[1]*double(Ny)); int iy = int(q[1]*double(Ny));
int iz = int(q[2]*double(Nz)); int iz = int(q[2]*double(Nz));
double x = q[0]*double(Nx)-double(ix); double x = q[0]*double(Nx)-double(ix);
double y = q[1]*double(Ny)-double(iy); double y = q[1]*double(Ny)-double(iy);
double z = q[2]*double(Nz)-double(iz); double z = q[2]*double(Nz)-double(iz);
int ixp = (ix+1)%Nx, iyp = (iy+1)%Ny, izp = (iz+1)%Nz; int ixp = (ix+1)%Nx, iyp = (iy+1)%Ny, izp = (iz+1)%Nz;
vidx[0] = (ix*Ny+iy)*Nz+iz; vidx[0] = (ix*Ny+iy)*Nz+iz;
vidx[1] = (ixp*Ny+iy)*Nz+iz; vidx[1] = (ixp*Ny+iy)*Nz+iz;
vidx[2] = (ix*Ny+iyp)*Nz+iz; vidx[2] = (ix*Ny+iyp)*Nz+iz;
vidx[3] = (ixp*Ny+iyp)*Nz+iz; vidx[3] = (ixp*Ny+iyp)*Nz+iz;
vidx[4] = (ix*Ny+iy)*Nz+izp; vidx[4] = (ix*Ny+iy)*Nz+izp;
vidx[5] = (ixp*Ny+iy)*Nz+izp; vidx[5] = (ixp*Ny+iy)*Nz+izp;
vidx[6] = (ix*Ny+iyp)*Nz+izp; vidx[6] = (ix*Ny+iyp)*Nz+izp;
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;
dfdz[i] = Dfdz[vidx[i]][idim].r; dfdz[i] = Dfdz[vidx[i]][idim].r;
d2fdxdy[i] = D2fdxdy[vidx[i]][idim].r; d2fdxdy[i] = D2fdxdy[vidx[i]][idim].r;
d2fdxdz[i] = D2fdxdz[vidx[i]][idim].r; d2fdxdz[i] = D2fdxdz[vidx[i]][idim].r;
d2fdydz[i] = D2fdydz[vidx[i]][idim].r; d2fdydz[i] = D2fdydz[vidx[i]][idim].r;
d3fdxdydz[i] = D3fdxdydz[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]); 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;
dfdz[i] = Dfdz[vidx[i]][idim].i; dfdz[i] = Dfdz[vidx[i]][idim].i;
d2fdxdy[i] = D2fdxdy[vidx[i]][idim].i; d2fdxdy[i] = D2fdxdy[vidx[i]][idim].i;
d2fdxdz[i] = D2fdxdz[vidx[i]][idim].i; d2fdxdz[i] = D2fdxdz[vidx[i]][idim].i;
d2fdydz[i] = D2fdydz[vidx[i]][idim].i; d2fdydz[i] = D2fdydz[vidx[i]][idim].i;
d3fdxdydz[i] = D3fdxdydz[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]); 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;
return;
} }
/* ---------------------------------------------------------------------------- /* ----------------------------------------------------------------------------
@ -303,63 +185,63 @@ return;
* ---------------------------------------------------------------------------- */ * ---------------------------------------------------------------------------- */
void Interpolate::trilinear(double *qin, doublecomplex *DMq) 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.;
} }
// find the index of the eight vertice // find the index of the eight vertice
int ix, iy, iz, ixp, iyp, izp; int ix, iy, iz, ixp, iyp, izp;
double x, y, z; double x, y, z;
q[0] *= double(Nx); q[0] *= double(Nx);
q[1] *= double(Ny); q[1] *= double(Ny);
q[2] *= double(Nz); q[2] *= double(Nz);
ix = int(q[0])%Nx; ix = int(q[0])%Nx;
iy = int(q[1])%Ny; iy = int(q[1])%Ny;
iz = int(q[2])%Nz; iz = int(q[2])%Nz;
ixp = (ix+1)%Nx; ixp = (ix+1)%Nx;
iyp = (iy+1)%Ny; iyp = (iy+1)%Ny;
izp = (iz+1)%Nz; izp = (iz+1)%Nz;
x = q[0] - double(ix); x = q[0] - double(ix);
y = q[1] - double(iy); y = q[1] - double(iy);
z = q[2] - double(iz); z = q[2] - double(iz);
//-------------------------------------- //--------------------------------------
vidx[0] = ((ix*Ny)+iy)*Nz + iz; vidx[0] = ((ix*Ny)+iy)*Nz + iz;
vidx[1] = ((ixp*Ny)+iy)*Nz + iz; vidx[1] = ((ixp*Ny)+iy)*Nz + iz;
vidx[2] = ((ix*Ny)+iyp)*Nz + iz; vidx[2] = ((ix*Ny)+iyp)*Nz + iz;
vidx[3] = ((ix*Ny)+iy)*Nz + izp; vidx[3] = ((ix*Ny)+iy)*Nz + izp;
vidx[4] = ((ixp*Ny)+iy)*Nz + izp; vidx[4] = ((ixp*Ny)+iy)*Nz + izp;
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);
fac[1] = x*(1.-y)*(1.-z); fac[1] = x*(1.-y)*(1.-z);
fac[2] = (1.-x)*y*(1.-z); fac[2] = (1.-x)*y*(1.-z);
fac[3] = (1.-x)*(1.-y)*z; fac[3] = (1.-x)*(1.-y)*z;
fac[4] = x*(1.-y)*z; fac[4] = x*(1.-y)*z;
fac[5] = (1.-x)*y*z; fac[5] = (1.-x)*y*z;
fac[6] = x*y*(1.-z); fac[6] = x*y*(1.-z);
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];
} }
} }
return; return;
} }
/* ---------------------------------------------------------------------------- /* ----------------------------------------------------------------------------
@ -367,12 +249,13 @@ return;
* ---------------------------------------------------------------------------- */ * ---------------------------------------------------------------------------- */
void Interpolate::execute(double *qin, doublecomplex *DMq) void Interpolate::execute(double *qin, doublecomplex *DMq)
{ {
UseGamma = 0; UseGamma = 0;
if (which == 1) // 1: tricubic if (which == 1) // 1: tricubic
tricubic(qin, DMq); tricubic(qin, DMq);
else // otherwise: trilinear else // otherwise: trilinear
trilinear(qin, DMq); trilinear(qin, DMq);
return;
return;
} }
/* ---------------------------------------------------------------------------- /* ----------------------------------------------------------------------------
@ -380,24 +263,23 @@ return;
* ---------------------------------------------------------------------------- */ * ---------------------------------------------------------------------------- */
void Interpolate::set_method() void Interpolate::set_method()
{ {
char str[MAXLINE]; char str[MAXLINE];
int im = 1; int im = 1;
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);
which =2-im%2; which =2-im%2;
printf("Your selection: %d\n", which); printf("Your selection: %d\n", which);
for(int i=0; i<80; i++) printf("="); for(int i=0; i<80; i++) printf("="); printf("\n\n");
printf("\n\n");
if (which == 1) tricubic_init();
if (which == 1) tricubic_init();
return;
return;
} }
/* ---------------------------------------------------------------------------- /* ----------------------------------------------------------------------------
@ -406,22 +288,23 @@ return;
* ---------------------------------------------------------------------------- */ * ---------------------------------------------------------------------------- */
void Interpolate::reset_gamma() void Interpolate::reset_gamma()
{ {
if (flag_reset_gamma) return; if (flag_reset_gamma) return;
flag_reset_gamma = 1; flag_reset_gamma = 1;
int p1 = 1%Nx, p2 = 2%Nx; int p1 = 1%Nx, p2 = 2%Nx;
int m1 = (Nx-1), m2 = (Nx-2+Nx)%Nx; int m1 = (Nx-1), m2 = (Nx-2+Nx)%Nx;
int ip1 = p1*Ny*Nz, ip2 = p2*Ny*Nz; int ip1 = p1*Ny*Nz, ip2 = p2*Ny*Nz;
int im1 = m1*Ny*Nz, im2 = m2*Ny*Nz; int im1 = m1*Ny*Nz, im2 = m2*Ny*Nz;
double const one6 = -1./6., two3 = 2./3.; double const one6 = -1./6., two3 = 2./3.;
for (int idim=0; idim<ndim; idim++){ for (int idim=0; idim < ndim; idim++){
data[0][idim].i = 0.; data[0][idim].i = 0.;
data[0][idim].r = (data[im2][idim].r + data[ip2][idim].r) * one6 data[0][idim].r = (data[im2][idim].r + data[ip2][idim].r) * one6
+ (data[im1][idim].r + data[ip1][idim].r) * two3; + (data[im1][idim].r + data[ip1][idim].r) * two3;
} }
return; return;
} }
/* ---------------------------------------------------------------------------- */

50
tools/phonon/interpolate.h
View File

@ -5,36 +5,38 @@
#include "stdlib.h" #include "stdlib.h"
#include "string.h" #include "string.h"
#include "memory.h" #include "memory.h"
#include "tricubic.h"
extern "C" typedef struct { double r, i; } doublecomplex; extern "C"{
#include "f2c.h"
using namespace std; #include "clapack.h"
}
class Interpolate{ class Interpolate{
public: public:
Interpolate(int, int, int, int, doublecomplex **); Interpolate(int, int, int, int, doublecomplex **);
~Interpolate(); ~Interpolate();
void set_method(); void set_method();
void execute(double *, doublecomplex *); void execute(double *, doublecomplex *);
void reset_gamma(); void reset_gamma();
int UseGamma; int UseGamma;
private: private:
void tricubic_init(); void tricubic_init();
void tricubic(double *, doublecomplex *); void tricubic(double *, doublecomplex *);
void trilinear(double *, doublecomplex *); void trilinear(double *, doublecomplex *);
Memory *memory; Memory *memory;
int which; int which;
int Nx, Ny, Nz, Npt, ndim; int Nx, Ny, Nz, Npt, ndim;
int flag_reset_gamma, flag_allocated_dfs; int flag_reset_gamma, flag_allocated_dfs;
doublecomplex **data; doublecomplex **data;
doublecomplex **Dfdx, **Dfdy, **Dfdz, **D2fdxdy, **D2fdxdz, **D2fdydz, **D3fdxdydz; 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]; double a[64], f[8], dfdx[8], dfdy[8], dfdz[8], d2fdxdy[8], d2fdxdz[8], d2fdydz[8], d3fdxdydz[8];
int vidx[8]; int vidx[8];
}; };
#endif #endif

2776
tools/phonon/kpath.cpp Normal file
View File

File diff suppressed because it is too large Load Diff

32
tools/phonon/kpath.h Normal file
View File

@ -0,0 +1,32 @@
#ifndef UseSPG
#define KPATH_H
#endif
#ifndef KPATH_H
#define KPATH_H
#include "qnodes.h"
#include "dynmat.h"
#include "memory.h"
class kPath{
public:
kPath(DynMat *, QNodes *);
~kPath();
void kpath();
void show_path();
void show_info();
private:
Memory *memory;
DynMat *dynmat;
QNodes *q;
char symbol[11];
int spgnum, sysdim, fftdim, num_atom, *attyp;
double latvec[3][3], **atpos;
};
#endif

1858
tools/phonon/phonon.cpp
View File

File diff suppressed because it is too large Load Diff

84
tools/phonon/phonon.h
View File

@ -11,50 +11,50 @@ using namespace std;
class Phonon{ class Phonon{
public: public:
Phonon(DynMat *); Phonon(DynMat *);
~Phonon(); ~Phonon();
DynMat *dynmat; DynMat *dynmat;
private: private:
int nq, ndim, sysdim; int nq, ndim, sysdim;
double **qpts, *wt; double **qpts, *wt;
double **eigs; double **eigs;
int ndos, nlocal, *locals; int ndos, nlocal, *locals;
double *dos, fmin, fmax, df, rdf; double *dos, fmin, fmax, df, rdf;
double ***ldos; double ***ldos;
Memory *memory; Memory *memory;
void QMesh(); void QMesh();
void ComputeAll(); void ComputeAll();
void pdos(); void pdos();
void pdisp(); void pdisp();
void therm(); void therm();
void ldos_egv(); void ldos_egv();
void ldos_rsgf(); void ldos_rsgf();
void local_therm(); void local_therm();
void dmanyq(); void dmanyq();
void vfanyq(); void vfanyq();
void DMdisp(); void DMdisp();
void vecanyq(); void vecanyq();
void ShowCell(); void ShowCell();
void smooth(double *, const int); void smooth(double *, const int);
void writeDOS(); void writeDOS();
void writeLDOS(); void writeLDOS();
void Normalize(); void Normalize();
int count_words(const char *); int count_words(const char *);
#ifdef UseSPG #ifdef UseSPG
int num_atom, *attyp; int num_atom, *attyp;
double latvec[3][3], **atpos; double latvec[3][3], **atpos;
#endif #endif
}; };

320
tools/phonon/phonopy.cpp Normal file
View File

@ -0,0 +1,320 @@
#ifdef FFTW3
#include <map>
#include "phonopy.h"
#include "math.h"
#include "kpath.h"
#include "fftw3.h"
/* ----------------------------------------------------------------------------
* Class Phonopy is designed to interface with phonopy.
* ---------------------------------------------------------------------------- */
Phonopy::Phonopy(DynMat *dynmat)
{
dm = dynmat;
memory = new Memory();
sysdim = dm->sysdim;
fftdim = dm->fftdim;
fftdim2 = fftdim * fftdim;
nucell = dm->nucell;
nx = ny = nz = 5;
write(1);
char str[MAXLINE];
if (count_words(fgets(str,MAXLINE,stdin)) >= 3){
nx = atoi(strtok(str," \t\n\r\f"));
ny = atoi(strtok(NULL," \t\n\r\f"));
nz = atoi(strtok(NULL," \t\n\r\f"));
}
if (dm->nx == 1) nx = 1;
if (dm->ny == 1) ny = 1;
if (dm->nz == 1) nz = 1;
if (nx < 1) nx = 1;
if (ny < 1) ny = 1;
if (nz < 1) nz = 1;
npt = nx * ny * nz;
write(2);
memory->create(mass, nucell, "Phonopy:mass");
for (int i = 0; i < nucell; ++i){
double m = 1./dm->M_inv_sqrt[i];
mass[i] = m * m;
}
memory->create(FC_all, npt, fftdim2, "phonopy:FC_all");
get_my_FC();
phonopy();
return;
}
/* ----------------------------------------------------------------------------
* Deconstructor, free the memory used.
* ---------------------------------------------------------------------------- */
Phonopy::~Phonopy()
{
memory->destroy(mass);
memory->destroy(FC_all);
delete memory;
dm = NULL;
return;
}
/* ----------------------------------------------------------------------------
* Subroutine to write the outputs to screen.
* ---------------------------------------------------------------------------- */
void Phonopy::write(int flag)
{
if (flag == 1){ // basic information
for (int ii = 0; ii < 80; ++ii) printf("="); printf("\n");
printf("Now to prepare the input files for phonopy.\n");
printf("The dimension of your present supercell is : %d x %d x %d.\n", dm->nx, dm->ny, dm->nz);
printf("The size of the force constant matrix will be: %d x %d.\n", dm->npt*3, dm->npt*3);
printf("Please define your desired dimension [5 5 5] : ");
} else if (flag == 2){
printf("\nThe new dimension of the supercell will be : %d x %d x %d.\n", nx, ny, nz);
printf("The size of the force constant matrix is then: %d x %d.\n", npt*3, npt*3);
} else if (flag == 3){
printf("\nNow to prepare the phonopy FORCE_CONSTANTS ..."); fflush(stdout);
} else if (flag == 4){
printf("Done!\nThe force constants information is extracted and written to FORCE_CONSTANTS,\n");
printf("the primitive cell is written to POSCAR.primitive, and the input file for\n");
printf("phonopy band evaluation is written to band.conf.\n");
printf("One should be able to obtain the phonon band structure after correcting\n");
printf("the element names in POSCAR.primitive and band.conf by running\n");
printf("`phonopy --readfc -c POSCAR.primitive -p band.conf`.\n");
for (int ii = 0; ii < 80; ++ii) printf("-");
printf("\n*** Remember to change the element names. ***\n");
#ifdef UseSPG
for (int ii = 0; ii < 80; ++ii) printf("-");
#endif
} else if (flag == 5){
for (int ii = 0; ii < 80; ++ii) printf("="); printf("\n");
}
return;
}
/* ----------------------------------------------------------------------------
* Driver to obtain the force constant matrix
* ---------------------------------------------------------------------------- */
void Phonopy::get_my_FC()
{
double q[3];
int ipt = 0;
for (int ix = 0; ix < nx; ++ix)
for (int iy = 0; iy < ny; ++iy)
for (int iz = 0; iz < nz; ++iz){
q[0] = double(ix)/double(nx);
q[1] = double(iy)/double(ny);
q[2] = double(iz)/double(nz);
dm->getDMq(q);
int ndim = 0;
for (int idim = 0; idim < fftdim; ++idim)
for (int jdim = 0; jdim < fftdim; ++jdim){
FC_all[ipt][ndim] = dm->DM_q[idim][jdim];
double m = sqrt(mass[idim/sysdim] * mass[jdim/sysdim]);
FC_all[ipt][ndim].r *= m;
FC_all[ipt][ndim].i *= m;
++ndim;
}
++ipt;
}
return;
}
/* ----------------------------------------------------------------------------
* Method to write out the force constants and related files for
* postprocessing with phonopy.
* ---------------------------------------------------------------------------- */
void Phonopy::phonopy()
{
// output info
write(3);
fftw_complex *in, *out;
double **fc;
memory->create(in, npt, "phonopy:in");
memory->create(out, npt, "phonopy:in");
memory->create(fc, npt, fftdim2, "phonopy:in");
fftw_plan plan = fftw_plan_dft_3d(nx, ny, nz, in, out, -1, FFTW_ESTIMATE);
double factor = dm->eml2fc / double(npt);
for (int idim = 0; idim < fftdim2; ++idim){
for (int i = 0; i < npt; ++i){
in[i][0] = FC_all[i][idim].r;
in[i][1] = FC_all[i][idim].i;
}
fftw_execute(plan);
for (int i = 0; i < npt; ++i) fc[i][idim] = out[i][0] * factor;
}
fftw_destroy_plan(plan);
memory->destroy(in);
memory->destroy(out);
// in POSCAR, atoms are sorted/aggregated by type, while for LAMMPS there is no such requirment
int type_id[nucell], num_type[nucell], ntype = 0;
for (int i = 0; i < nucell; ++i) num_type[i] = 0;
for (int i = 0; i < nucell; ++i){
int ip = ntype;
for (int j = 0; j < ntype; ++j){
if (dm->attyp[i] == type_id[j]) ip = j;
}
if (ip == ntype) type_id[ntype++] = dm->attyp[i];
num_type[ip]++;
}
std::map<int, int> iu_by_type;
iu_by_type.clear();
int id_new = 0;
for (int i = 0; i < ntype; ++i){
for (int j = 0; j < nucell; ++j){
if (dm->attyp[j] == type_id[i]) iu_by_type[id_new++] = j;
}
}
// write the FORCE_CONSTANTS file
FILE *fp = fopen("FORCE_CONSTANTS", "w");
int natom = npt * nucell;
fprintf(fp, "%d %d\n", natom, natom);
for (int i = 0; i < natom; ++i){
int iu = i / npt;
int iz = (i % npt) / (nx * ny);
int iy = (i % (nx *ny) ) / nx;
int ix = i % nx;
iu = iu_by_type[iu];
for (int j = 0; j < natom; ++j){
int ju = j / npt;
int jz = (j % npt) / (nx * ny);
int jy = (j % (nx *ny) ) / nx;
int jx = j % nx;
int dx = abs(ix - jx);
int dy = abs(iy - jy);
int dz = abs(iz - jz);
int id = (dx * ny + dy) * nz + dz;
ju = iu_by_type[ju];
fprintf(fp, "%d %d\n", i+1, j+1);
for (int idim = iu * sysdim; idim < (iu+1)*sysdim; ++idim){
for (int jdim = ju * sysdim; jdim < (ju+1)*sysdim; ++jdim){
int dd = idim * fftdim + jdim;
fprintf(fp, " %lg", fc[id][dd]);
}
fprintf(fp, "\n");
}
}
}
fclose(fp);
iu_by_type.clear();
memory->destroy(fc);
// write the primitive cell in POSCAR format
fp = fopen("POSCAR.primitive", "w");
fprintf(fp, "Fix-phonon unit cell");
for (int ip = 0; ip < ntype; ++ip) fprintf(fp, ", Elem-%d: %lg", type_id[ip], mass[ip]);
fprintf(fp, "\n1.\n");
int ndim = 0;
for (int idim = 0; idim < 3; ++idim){
for (int jdim = 0; jdim < 3; ++jdim) fprintf(fp, "%lg ", dm->basevec[ndim++]);
fprintf(fp, "\n");
}
for (int ip = 0; ip < ntype; ++ip) fprintf(fp, "Elem-%d ", type_id[ip]); fprintf(fp, "\n");
for (int ip = 0; ip < ntype; ++ip) fprintf(fp, "%d ", num_type[ip]);
fprintf(fp, "\nDirect\n");
for (int ip = 0; ip < ntype; ++ip){
for (int i = 0; i < nucell; ++i){
if (dm->attyp[i] == type_id[ip]){
fprintf(fp, "%lg %lg %lg\n", dm->basis[i][0], dm->basis[i][1], dm->basis[i][2]);
}
}
}
fclose(fp);
#ifdef UseSPG
// Get high symmetry k-path
QNodes *q = new QNodes();
kPath *kp = new kPath(dm, q);
kp->kpath();
#endif
// band.conf
fp = fopen("band.conf", "w");
fprintf(fp, "# From Fix-phonon");
for (int ip = 0; ip < ntype; ++ip) fprintf(fp, ", Elem-%d: %lg", type_id[ip], mass[ip]);
fprintf(fp, "\n\nATOM_NAME = ");
for (int ip = 0; ip < ntype; ++ip) fprintf(fp, "Elem-%d ", type_id[ip]);
fprintf(fp, "\nDIM = %d %d %d\nBAND = ", nx, ny, nz);
#ifdef UseSPG
int nsect = q->qs.size();
for (int i = 0; i < nsect; ++i){
fprintf(fp, " %lg %lg %lg", q->qs[i][0], q->qs[i][1], q->qs[i][2]);
if (i+1 < nsect){
double dq = 0.;
for (int j = 0; j < 3; ++j) dq += (q->qe[i][j] - q->qs[i+1][j]) * (q->qe[i][j] - q->qs[i+1][j]);
if (dq > ZERO) {
fprintf(fp, " %lg %lg %lg,", q->qe[i][0], q->qe[i][1], q->qe[i][2]);
}
} else if (i+1 == nsect){
fprintf(fp, " %lg %lg %lg\n", q->qe[i][0], q->qe[i][1], q->qe[i][2]);
}
}
#endif
fprintf(fp, "\nBAND_POINTS = 21\nBAND_LABELS =");
#ifdef UseSPG
for (int i = 0; i < q->ndstr.size(); ++i){
std::size_t found = q->ndstr[i].find("{/Symbol G}");
if (found != std::string::npos) q->ndstr[i].replace(found, found+11, "$\\Gamma$");
found = q->ndstr[i].find("/");
if (found != std::string::npos) q->ndstr[i].replace(found, found, " ");
fprintf(fp, " %s", q->ndstr[i].c_str());
}
#endif
fprintf(fp, "\nFORCE_CONSTANTS = READ\nBAND_CONNECTION = .TRUE.\n");
// output info
write(4);
#ifdef UseSPG
kp->show_path();
delete kp;
delete q;
#endif
write(5);
return;
}
/*------------------------------------------------------------------------------
* Method to count # of words in a string, without destroying the string
*----------------------------------------------------------------------------*/
int Phonopy::count_words(const char *line)
{
int n = strlen(line) + 1;
char *copy;
memory->create(copy, n, "count_words:copy");
strcpy(copy,line);
char *ptr;
if (ptr = strchr(copy,'#')) *ptr = '\0';
if (strtok(copy," \t\n\r\f") == NULL) {
memory->destroy(copy);
return 0;
}
n = 1;
while (strtok(NULL," \t\n\r\f")) n++;
memory->destroy(copy);
return n;
}
/*----------------------------------------------------------------------------*/
#endif

36
tools/phonon/phonopy.h Normal file
View File

@ -0,0 +1,36 @@
#ifndef FFTW3
#define PHONOPY_H
#endif
#ifndef PHONOPY_H
#define PHONOPY_H
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "memory.h"
#include "qnodes.h"
#include "dynmat.h"
#include "global.h"
class Phonopy {
public:
Phonopy(DynMat *);
~Phonopy();
private:
Memory *memory;
char str[MAXLINE];
int npt, fftdim2; // local variables
int nx, ny, nz, nucell; // local variables
int sysdim, fftdim; // local variables
double *mass;
doublecomplex **FC_all;
DynMat *dm;
void write(int);
void get_my_FC();
void phonopy();
int count_words(const char *line);
};
#endif

30
tools/phonon/qnodes.cpp Normal file
View File

@ -0,0 +1,30 @@
#include "qnodes.h"
/* ----------------------------------------------------------------------------
* Class QNodes stores the high symmetry k-path nodes for a given lattice.
* The constructor and the deconstructor simply empties the data.
* ---------------------------------------------------------------------------- */
QNodes::QNodes()
{
nodes.clear();
ndstr.clear();
qs.clear();
qe.clear();
nqbin.clear();
return;
}
/* ----------------------------------------------------------------------------
* The constructor and the deconstructor simply empties the data.
* ---------------------------------------------------------------------------- */
QNodes::~QNodes()
{
nodes.clear();
ndstr.clear();
qs.clear();
qe.clear();
nqbin.clear();
return;
}

17
tools/phonon/qnodes.h Normal file
View File

@ -0,0 +1,17 @@
#ifndef QNODES_H
#define QNODES_H
#include <vector>
#include <string>
class QNodes {
public:
QNodes();
~QNodes();
std::vector<double> nodes;
std::vector<std::string> ndstr;
std::vector<double *> qs, qe;
std::vector<int> nqbin;
};
#endif

36
tools/phonon/timer.cpp
View File

@ -4,9 +4,10 @@
* -------------------------------------------------------------------------- */ * -------------------------------------------------------------------------- */
Timer::Timer() Timer::Timer()
{ {
flag = 0; flag = 0;
start(); start();
return;
return;
} }
/* ----------------------------------------------------------------------------- /* -----------------------------------------------------------------------------
@ -14,10 +15,9 @@ return;
* -------------------------------------------------------------------------- */ * -------------------------------------------------------------------------- */
void Timer::start() void Timer::start()
{ {
t1 = clock(); t1 = clock();
flag |= 1; flag |= 1;
return;
return;
} }
/* ----------------------------------------------------------------------------- /* -----------------------------------------------------------------------------
@ -25,11 +25,11 @@ return;
* -------------------------------------------------------------------------- */ * -------------------------------------------------------------------------- */
void Timer::stop() void Timer::stop()
{ {
if ( flag&1 ) { if ( flag&1 ) {
t2 = clock(); t2 = clock();
flag |= 2; flag |= 2;
} }
return; return;
} }
/* ----------------------------------------------------------------------------- /* -----------------------------------------------------------------------------
@ -37,12 +37,12 @@ return;
* -------------------------------------------------------------------------- */ * -------------------------------------------------------------------------- */
void Timer::print() void Timer::print()
{ {
if ( (flag&3) != 3) return; if ( (flag&3) != 3) return;
cpu_time_used = ((double) (t2 - t1)) / CLOCKS_PER_SEC; cpu_time_used = ((double) (t2 - t1)) / CLOCKS_PER_SEC;
printf("Total CPU time used: %g seconds.\n", cpu_time_used); printf("Total CPU time used: %g seconds.\n", cpu_time_used);
return; return;
} }
/* ----------------------------------------------------------------------------- /* -----------------------------------------------------------------------------
@ -50,6 +50,6 @@ return;
* -------------------------------------------------------------------------- */ * -------------------------------------------------------------------------- */
double Timer::elapse() double Timer::elapse()
{ {
if ( (flag&3) != 3) return 0.; if ( (flag&3) != 3) return 0.;
else return ((double) (t2 - t1)) / CLOCKS_PER_SEC; else return ((double) (t2 - t1)) / CLOCKS_PER_SEC;
} }

2
tools/phonon/version.h
View File

@ -1 +1 @@
#define VERSION 8 #define VERSION 21