update phana tool for USER-PHONON from github

This commit is contained in:
Axel Kohlmeyer 2020-07-29 10:48:04 -04:00
parent 29c50671da
commit 7df8a63045
No known key found for this signature in database
GPG Key ID: D9B44E93BF0C375A
19 changed files with 5398 additions and 4895 deletions

View File

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

View File

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

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

View File

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

View File

@ -38,36 +38,36 @@
Green::Green(const int ntm, const int sdim, const int niter, const double min, const double max,
const int ndos, const double eps, double **Hessian, const int itm, double **lpdos)
{
const double tpi = 8.*atan(1.);
natom = ntm; sysdim = sdim; nit = niter; epson = eps;
wmin = min*tpi; wmax = max*tpi; nw = ndos + (ndos+1)%2;
H = Hessian; iatom = itm;
ldos = lpdos;
memory = new Memory();
if (natom < 1 || iatom < 0 || iatom >= natom){
printf("\nError: Wrong number of total atoms or wrong index of interested atom!\n");
return;
}
ndim = natom * sysdim;
if (nit < 1){printf("\nError: Wrong input of maximum iterations!\n"); return;}
if (nit > ndim){printf("\nError: # Lanczos iterations is not expected to exceed the degree of freedom!\n"); return;}
if (nw < 1){printf("\nError: Wrong input of points in LDOS!\n"); return;}
// initialize variables and allocate local memories
dw = (wmax - wmin)/double(nw-1);
memory->create(alpha, sysdim,nit, "Green_Green:alpha");
memory->create(beta, sysdim,nit+1,"Green_Green:beta");
//memory->create(ldos, nw,sysdim, "Green_Green:ldos");
// use Lanczos algorithm to diagonalize the Hessian
Lanczos();
// Get the inverser of the treated hessian by continued fractional method
Recursion();
return;
const double tpi = 8.*atan(1.);
natom = ntm; sysdim = sdim; nit = niter; epson = eps;
wmin = min*tpi; wmax = max*tpi; nw = ndos + (ndos+1)%2;
H = Hessian; iatom = itm;
ldos = lpdos;
memory = new Memory();
if (natom < 1 || iatom < 0 || iatom >= natom){
printf("\nError: Wrong number of total atoms or wrong index of interested atom!\n");
return;
}
ndim = natom * sysdim;
if (nit < 1){printf("\nError: Wrong input of maximum iterations!\n"); return;}
if (nit > ndim){printf("\nError: # Lanczos iterations is not expected to exceed the degree of freedom!\n"); return;}
if (nw < 1){printf("\nError: Wrong input of points in LDOS!\n"); return;}
// initialize variables and allocate local memories
dw = (wmax - wmin)/double(nw-1);
memory->create(alpha, sysdim,nit, "Green_Green:alpha");
memory->create(beta, sysdim,nit+1,"Green_Green:beta");
//memory->create(ldos, nw,sysdim, "Green_Green:ldos");
// use Lanczos algorithm to diagonalize the Hessian
Lanczos();
// Get the inverser of the treated hessian by continued fractional method
Recursion();
return;
}
/*------------------------------------------------------------------------------
@ -75,14 +75,14 @@ return;
*----------------------------------------------------------------------------*/
Green::~Green()
{
H = NULL;
ldos = NULL;
memory->destroy(alpha);
memory->destroy(beta);
delete memory;
return;
H = NULL;
ldos = NULL;
memory->destroy(alpha);
memory->destroy(beta);
delete memory;
return;
}
/*------------------------------------------------------------------------------
@ -90,159 +90,159 @@ return;
*----------------------------------------------------------------------------*/
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];
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;
}
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;
}
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
*----------------------------------------------------------------------------*/
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;
// 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];
}
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;
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;
}
ldos[i][idim] = std::imag(rec_x)*w;
}
w += dw;
}
delete []alpha_inf;
delete []beta_inf;
delete []xmin;
delete []xmax;
return;
w += dw;
}
delete []alpha_inf;
delete []beta_inf;
delete []xmin;
delete []xmax;
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)
*----------------------------------------------------------------------------*/
void Green::recursion()
{
// local variables
std::complex<double> Z, rec_x, rec_x_inv;
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;
// 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;
}
ldos[i][idim] = std::imag(rec_x)*w;
}
w += dw;
}
return;
w += dw;
}
return;
}
/*------------------------------------------------------------------------------*/

View File

@ -5,21 +5,21 @@
class Green{
public:
Green(const int, const int, const int, const double, const double,
const int, const double, double **, const int, double **);
~Green();
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;
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

View File

@ -1,144 +1,26 @@
#include "interpolate.h"
#include <math.h>
#include "math.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
* ---------------------------------------------------------------------------- */
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;
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;
}
/* ----------------------------------------------------------------------------
@ -146,77 +28,76 @@ return;
* ---------------------------------------------------------------------------- */
void Interpolate::tricubic_init()
{
// prepare necessary data for tricubic
if (flag_allocated_dfs == 0){
memory->create(Dfdx, Npt, ndim, "Interpolate_Interpolate:Dfdx");
memory->create(Dfdy, Npt, ndim, "Interpolate_Interpolate:Dfdy");
memory->create(Dfdz, Npt, ndim, "Interpolate_Interpolate:Dfdz");
memory->create(D2fdxdy, Npt, ndim, "Interpolate_Interpolate:D2fdxdy");
memory->create(D2fdxdz, Npt, ndim, "Interpolate_Interpolate:D2fdxdz");
memory->create(D2fdydz, Npt, ndim, "Interpolate_Interpolate:D2fdydz");
memory->create(D3fdxdydz, Npt, ndim, "Interpolate_Interpolate:D2fdxdydz");
// prepare necessary data for tricubic
if (flag_allocated_dfs == 0){
memory->create(Dfdx, Npt, ndim, "Interpolate_Interpolate:Dfdx");
memory->create(Dfdy, Npt, ndim, "Interpolate_Interpolate:Dfdy");
memory->create(Dfdz, Npt, ndim, "Interpolate_Interpolate:Dfdz");
memory->create(D2fdxdy, Npt, ndim, "Interpolate_Interpolate:D2fdxdy");
memory->create(D2fdxdz, Npt, ndim, "Interpolate_Interpolate:D2fdxdz");
memory->create(D2fdydz, Npt, ndim, "Interpolate_Interpolate:D2fdydz");
memory->create(D3fdxdydz, Npt, ndim, "Interpolate_Interpolate:D2fdxdydz");
flag_allocated_dfs = 1;
}
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;
// 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;
}
/* ----------------------------------------------------------------------------
@ -224,15 +105,17 @@ return;
* ---------------------------------------------------------------------------- */
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;
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;
return;
}
/* ----------------------------------------------------------------------------
@ -240,60 +123,59 @@ Interpolate::~Interpolate()
* ---------------------------------------------------------------------------- */
void Interpolate::tricubic(double *qin, doublecomplex *DMq)
{
// qin should be in unit of 2*pi/L
double q[3];
for (int i = 0; i < 3; ++i) q[i] = qin[i];
for (int i = 0; i < 3; ++i){
while (q[i] < 0.) q[i] += 1.;
while (q[i] >= 1.) q[i] -= 1.;
}
int ix = int(q[0]*double(Nx));
int iy = int(q[1]*double(Ny));
int iz = int(q[2]*double(Nz));
double x = q[0]*double(Nx)-double(ix);
double y = q[1]*double(Ny)-double(iy);
double z = q[2]*double(Nz)-double(iz);
int ixp = (ix+1)%Nx, iyp = (iy+1)%Ny, izp = (iz+1)%Nz;
vidx[0] = (ix*Ny+iy)*Nz+iz;
vidx[1] = (ixp*Ny+iy)*Nz+iz;
vidx[2] = (ix*Ny+iyp)*Nz+iz;
vidx[3] = (ixp*Ny+iyp)*Nz+iz;
vidx[4] = (ix*Ny+iy)*Nz+izp;
vidx[5] = (ixp*Ny+iy)*Nz+izp;
vidx[6] = (ix*Ny+iyp)*Nz+izp;
vidx[7] = (ixp*Ny+iyp)*Nz+izp;
for (int i=0; i<8; i++) if (vidx[i] == 0) UseGamma = 1;
for (int idim = 0; idim < ndim; ++idim){
for (int i = 0; i < 8; ++i){
f[i] = data[vidx[i]][idim].r;
dfdx[i] = Dfdx[vidx[i]][idim].r;
dfdy[i] = Dfdy[vidx[i]][idim].r;
dfdz[i] = Dfdz[vidx[i]][idim].r;
d2fdxdy[i] = D2fdxdy[vidx[i]][idim].r;
d2fdxdz[i] = D2fdxdz[vidx[i]][idim].r;
d2fdydz[i] = D2fdydz[vidx[i]][idim].r;
d3fdxdydz[i] = D3fdxdydz[vidx[i]][idim].r;
}
tricubic_get_coeff(&a[0],&f[0],&dfdx[0],&dfdy[0],&dfdz[0],&d2fdxdy[0],&d2fdxdz[0],&d2fdydz[0],&d3fdxdydz[0]);
DMq[idim].r = tricubic_eval(&a[0],x,y,z);
for (int i = 0; i < 8; ++i){
f[i] = data[vidx[i]][idim].i;
dfdx[i] = Dfdx[vidx[i]][idim].i;
dfdy[i] = Dfdy[vidx[i]][idim].i;
dfdz[i] = Dfdz[vidx[i]][idim].i;
d2fdxdy[i] = D2fdxdy[vidx[i]][idim].i;
d2fdxdz[i] = D2fdxdz[vidx[i]][idim].i;
d2fdydz[i] = D2fdydz[vidx[i]][idim].i;
d3fdxdydz[i] = D3fdxdydz[vidx[i]][idim].i;
}
tricubic_get_coeff(&a[0],&f[0],&dfdx[0],&dfdy[0],&dfdz[0],&d2fdxdy[0],&d2fdxdz[0],&d2fdydz[0],&d3fdxdydz[0]);
DMq[idim].i = tricubic_eval(&a[0],x,y,z);
}
return;
// qin should be in unit of 2*pi/L
double q[3];
for (int i = 0; i < 3; ++i) q[i] = qin[i];
for (int i = 0; i < 3; ++i){
while (q[i] < 0.) q[i] += 1.;
while (q[i] >= 1.) q[i] -= 1.;
}
int ix = int(q[0]*double(Nx));
int iy = int(q[1]*double(Ny));
int iz = int(q[2]*double(Nz));
double x = q[0]*double(Nx)-double(ix);
double y = q[1]*double(Ny)-double(iy);
double z = q[2]*double(Nz)-double(iz);
int ixp = (ix+1)%Nx, iyp = (iy+1)%Ny, izp = (iz+1)%Nz;
vidx[0] = (ix*Ny+iy)*Nz+iz;
vidx[1] = (ixp*Ny+iy)*Nz+iz;
vidx[2] = (ix*Ny+iyp)*Nz+iz;
vidx[3] = (ixp*Ny+iyp)*Nz+iz;
vidx[4] = (ix*Ny+iy)*Nz+izp;
vidx[5] = (ixp*Ny+iy)*Nz+izp;
vidx[6] = (ix*Ny+iyp)*Nz+izp;
vidx[7] = (ixp*Ny+iyp)*Nz+izp;
for (int i=0; i<8; i++) if (vidx[i] == 0) UseGamma = 1;
for (int idim = 0; idim < ndim; ++idim){
for (int i = 0; i < 8; ++i){
f[i] = data[vidx[i]][idim].r;
dfdx[i] = Dfdx[vidx[i]][idim].r;
dfdy[i] = Dfdy[vidx[i]][idim].r;
dfdz[i] = Dfdz[vidx[i]][idim].r;
d2fdxdy[i] = D2fdxdy[vidx[i]][idim].r;
d2fdxdz[i] = D2fdxdz[vidx[i]][idim].r;
d2fdydz[i] = D2fdydz[vidx[i]][idim].r;
d3fdxdydz[i] = D3fdxdydz[vidx[i]][idim].r;
}
tricubic_get_coeff(&a[0],&f[0],&dfdx[0],&dfdy[0],&dfdz[0],&d2fdxdy[0],&d2fdxdz[0],&d2fdydz[0],&d3fdxdydz[0]);
DMq[idim].r = tricubic_eval(&a[0],x,y,z);
for (int i = 0; i < 8; ++i){
f[i] = data[vidx[i]][idim].i;
dfdx[i] = Dfdx[vidx[i]][idim].i;
dfdy[i] = Dfdy[vidx[i]][idim].i;
dfdz[i] = Dfdz[vidx[i]][idim].i;
d2fdxdy[i] = D2fdxdy[vidx[i]][idim].i;
d2fdxdz[i] = D2fdxdz[vidx[i]][idim].i;
d2fdydz[i] = D2fdydz[vidx[i]][idim].i;
d3fdxdydz[i] = D3fdxdydz[vidx[i]][idim].i;
}
tricubic_get_coeff(&a[0],&f[0],&dfdx[0],&dfdy[0],&dfdz[0],&d2fdxdy[0],&d2fdxdz[0],&d2fdydz[0],&d3fdxdydz[0]);
DMq[idim].i = tricubic_eval(&a[0],x,y,z);
}
return;
}
/* ----------------------------------------------------------------------------
@ -303,63 +185,63 @@ return;
* ---------------------------------------------------------------------------- */
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;
// 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;
}
/* ----------------------------------------------------------------------------
@ -367,12 +249,13 @@ return;
* ---------------------------------------------------------------------------- */
void Interpolate::execute(double *qin, doublecomplex *DMq)
{
UseGamma = 0;
if (which == 1) // 1: tricubic
tricubic(qin, DMq);
else // otherwise: trilinear
trilinear(qin, DMq);
return;
UseGamma = 0;
if (which == 1) // 1: tricubic
tricubic(qin, DMq);
else // otherwise: trilinear
trilinear(qin, DMq);
return;
}
/* ----------------------------------------------------------------------------
@ -380,24 +263,23 @@ return;
* ---------------------------------------------------------------------------- */
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;
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;
}
/* ----------------------------------------------------------------------------
@ -406,22 +288,23 @@ return;
* ---------------------------------------------------------------------------- */
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;
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;
}
/* ---------------------------------------------------------------------------- */

View File

@ -5,36 +5,38 @@
#include "stdlib.h"
#include "string.h"
#include "memory.h"
#include "tricubic.h"
extern "C" typedef struct { double r, i; } doublecomplex;
using namespace std;
extern "C"{
#include "f2c.h"
#include "clapack.h"
}
class Interpolate{
public:
Interpolate(int, int, int, int, doublecomplex **);
~Interpolate();
void set_method();
void execute(double *, doublecomplex *);
void reset_gamma();
int UseGamma;
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];
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

2776
tools/phonon/kpath.cpp Normal 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

File diff suppressed because it is too large Load Diff

View File

@ -11,50 +11,50 @@ using namespace std;
class Phonon{
public:
Phonon(DynMat *);
~Phonon();
DynMat *dynmat;
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 ShowCell();
void smooth(double *, const int);
void writeDOS();
void writeLDOS();
void Normalize();
int count_words(const char *);
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 ShowCell();
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;
int num_atom, *attyp;
double latvec[3][3], **atpos;
#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

View File

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

View File

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