From 478058119b051ea256a519f843abecc8c9a7bb43 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 29 Apr 2023 02:59:04 -0400 Subject: [PATCH] integrate CMake build procedure for tools/phonon --- cmake/CMakeLists.txt | 2 +- cmake/Modules/Tools.cmake | 2 + tools/phonon/CMakeLists.spglib | 54 +++++ tools/phonon/CMakeLists.txt | 99 ++++++++ tools/phonon/Makefile | 67 ------ tools/phonon/README | 56 +++-- tools/phonon/disp.cpp | 43 ++-- tools/phonon/dynmat.cpp | 72 ++++-- tools/phonon/dynmat.h | 20 +- tools/phonon/global.h | 2 +- tools/phonon/green.cpp | 13 +- tools/phonon/green.h | 5 +- tools/phonon/input.cpp | 35 +++ tools/phonon/input.h | 17 ++ tools/phonon/interpolate.cpp | 21 +- tools/phonon/interpolate.h | 15 +- tools/phonon/kpath.cpp | 39 ++- tools/phonon/kpath.h | 17 +- tools/phonon/main.cpp | 4 - tools/phonon/memory.cpp | 6 +- tools/phonon/memory.h | 308 +++++------------------- tools/phonon/phonon.cpp | 176 +++++++++----- tools/phonon/phonon.h | 16 +- tools/phonon/phonopy.cpp | 141 ++++++----- tools/phonon/phonopy.h | 16 +- tools/phonon/timer.cpp | 3 + tools/phonon/timer.h | 4 +- tools/phonon/tricubic/CMakeLists.txt | 18 ++ tools/phonon/tricubic/LICENSE | 340 +++++++++++++++++++++++++++ tools/phonon/tricubic/README.md | 28 +++ tools/phonon/tricubic/tricubic.cpp | 185 +++++++++++++++ tools/phonon/tricubic/tricubic.h | 10 + tools/phonon/version.h | 1 - tools/phonon/zheevd.h | 16 ++ 34 files changed, 1260 insertions(+), 591 deletions(-) create mode 100644 tools/phonon/CMakeLists.spglib create mode 100644 tools/phonon/CMakeLists.txt delete mode 100644 tools/phonon/Makefile create mode 100644 tools/phonon/input.cpp create mode 100644 tools/phonon/input.h create mode 100644 tools/phonon/tricubic/CMakeLists.txt create mode 100644 tools/phonon/tricubic/LICENSE create mode 100644 tools/phonon/tricubic/README.md create mode 100644 tools/phonon/tricubic/tricubic.cpp create mode 100644 tools/phonon/tricubic/tricubic.h delete mode 100644 tools/phonon/version.h create mode 100644 tools/phonon/zheevd.h diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt index a71347c2c4..88b7ec422e 100644 --- a/cmake/CMakeLists.txt +++ b/cmake/CMakeLists.txt @@ -440,7 +440,7 @@ if(BUILD_OMP) target_link_libraries(lmp PRIVATE OpenMP::OpenMP_CXX) endif() -if(PKG_MSCG OR PKG_ATC OR PKG_AWPMD OR PKG_ML-QUIP OR PKG_ML-POD OR PKG_ELECTRODE) +if(PKG_MSCG OR PKG_ATC OR PKG_AWPMD OR PKG_ML-QUIP OR PKG_ML-POD OR PKG_ELECTRODE OR BUILD_TOOLS) enable_language(C) if (NOT USE_INTERNAL_LINALG) find_package(LAPACK) diff --git a/cmake/Modules/Tools.cmake b/cmake/Modules/Tools.cmake index c4c33cc40c..ba91a557bc 100644 --- a/cmake/Modules/Tools.cmake +++ b/cmake/Modules/Tools.cmake @@ -63,6 +63,8 @@ if(BUILD_LAMMPS_SHELL) install(TARGETS lammps-shell EXPORT LAMMPS_Targets DESTINATION ${CMAKE_INSTALL_BINDIR}) install(DIRECTORY ${LAMMPS_TOOLS_DIR}/lammps-shell/icons DESTINATION ${CMAKE_INSTALL_DATAROOTDIR}/) install(FILES ${LAMMPS_TOOLS_DIR}/lammps-shell/lammps-shell.desktop DESTINATION ${CMAKE_INSTALL_DATAROOTDIR}/applications/) + + add_subdirectory(${LAMMPS_TOOLS_DIR}/phonon ${CMAKE_BINARY_DIR}/phana_build) endif() diff --git a/tools/phonon/CMakeLists.spglib b/tools/phonon/CMakeLists.spglib new file mode 100644 index 0000000000..9e22445b80 --- /dev/null +++ b/tools/phonon/CMakeLists.spglib @@ -0,0 +1,54 @@ +cmake_minimum_required(VERSION 3.10) +project(spglib C) +set(CMAKE_MACOSX_RPATH 1) +set(CMAKE_C_FLAGS_RELEASE "-Wall -O2") +set(CMAKE_C_FLAGS_DEBUG "-g -DSPGDEBUG -DSPGWARNING") +if(NOT CMAKE_BUILD_TYPE) + set(CMAKE_BUILD_TYPE Release) +endif(NOT CMAKE_BUILD_TYPE) +message(STATUS "Build type: ${CMAKE_BUILD_TYPE}") + +set(CMAKE_POSITION_INDEPENDENT_CODE ON) + + +# Version numbers +file(READ ${PROJECT_SOURCE_DIR}/src/version.h version_file) +string(REGEX MATCH "SPGLIB_MAJOR_VERSION ([0-9]+)" spglib_major_version ${version_file}) +set(spglib_major_version ${CMAKE_MATCH_1}) +string(REGEX MATCH "SPGLIB_MINOR_VERSION ([0-9]+)" spglib_minor_version ${version_file}) +set(spglib_minor_version ${CMAKE_MATCH_1}) +string(REGEX MATCH "SPGLIB_MICRO_VERSION ([0-9]+)" spglib_micro_version ${version_file}) +set(spglib_micro_version ${CMAKE_MATCH_1}) +set(serial "${spglib_major_version}.${spglib_minor_version}.${spglib_micro_version}") +set(soserial "1") + +# Source code +include_directories("${PROJECT_SOURCE_DIR}/src") +set(SOURCES ${PROJECT_SOURCE_DIR}/src/arithmetic.c + ${PROJECT_SOURCE_DIR}/src/cell.c + ${PROJECT_SOURCE_DIR}/src/debug.c + ${PROJECT_SOURCE_DIR}/src/delaunay.c + ${PROJECT_SOURCE_DIR}/src/determination.c + ${PROJECT_SOURCE_DIR}/src/hall_symbol.c + ${PROJECT_SOURCE_DIR}/src/kgrid.c + ${PROJECT_SOURCE_DIR}/src/kpoint.c + ${PROJECT_SOURCE_DIR}/src/mathfunc.c + ${PROJECT_SOURCE_DIR}/src/niggli.c + ${PROJECT_SOURCE_DIR}/src/overlap.c + ${PROJECT_SOURCE_DIR}/src/pointgroup.c + ${PROJECT_SOURCE_DIR}/src/primitive.c + ${PROJECT_SOURCE_DIR}/src/refinement.c + ${PROJECT_SOURCE_DIR}/src/site_symmetry.c + ${PROJECT_SOURCE_DIR}/src/sitesym_database.c + ${PROJECT_SOURCE_DIR}/src/spacegroup.c + ${PROJECT_SOURCE_DIR}/src/spg_database.c + ${PROJECT_SOURCE_DIR}/src/spglib.c + ${PROJECT_SOURCE_DIR}/src/spin.c + ${PROJECT_SOURCE_DIR}/src/symmetry.c) + +add_library(symspg STATIC ${SOURCES}) +install(TARGETS symspg LIBRARY) + +# Header file +install(FILES ${PROJECT_SOURCE_DIR}/src/spglib.h DESTINATION ${CMAKE_INSTALL_PREFIX}/include) + diff --git a/tools/phonon/CMakeLists.txt b/tools/phonon/CMakeLists.txt new file mode 100644 index 0000000000..283b995590 --- /dev/null +++ b/tools/phonon/CMakeLists.txt @@ -0,0 +1,99 @@ + +# Support Linux from Ubuntu 20.04LTS onward, CentOS 7.x (with EPEL), +# macOS, MSVC 2019 (=Version 16) +cmake_minimum_required(VERSION 3.16) + +# set timestamp of downloaded files to that of archive +if(POLICY CMP0135) + cmake_policy(SET CMP0135 NEW) +endif() + +# set up project +set(PHANA_MINOR_VERSION 48) +project(phonon VERSION ${PHANA_MINOR_VERSION} + DESCRIPTION "Fix phonon post-processor" + LANGUAGES CXX C) +set(CMAKE_POSITION_INDEPENDENT_CODE ON) +if(NOT CMAKE_BUILD_TYPE) + set(CMAKE_BUILD_TYPE RelWithDebInfo) +endif() + +# hacks for MSVC to prevent lots of pointless warnings about "unsafe" functions, +# padding and Spectre mitigation +if(MSVC) + add_compile_options(/wd4244) + add_compile_options(/wd4267) + add_compile_options(/wd4711) + add_compile_options(/wd4820) + add_compile_options(/wd5045) + add_compile_definitions(_CRT_SECURE_NO_WARNINGS) +endif() + +set(CMAKE_PROJECT_VERSION ${PHANA_MINOR_VERSION}) +configure_file(version.h.in version.h @ONLY) +add_executable(phana + main.cpp + disp.cpp + dynmat.cpp + green.cpp + input.cpp + interpolate.cpp + kpath.cpp + memory.cpp + phonon.cpp + phonopy.cpp + qnodes.cpp + timer.cpp +) +target_include_directories(phana PUBLIC $) + +find_package(FFTW3) +if(FFTW3_FOUND) + target_compile_definitions(phana PRIVATE FFTW3) + target_link_libraries(phana PRIVATE FFTW3::FFTW3) +endif() + +# build bundeled libraries +add_subdirectory(tricubic) + +option(USE_SPGLIB "Download and use spglib for phonon DOS and other optional properties" ON) +if(USE_SPGLIB) + set(SPGLIB_URL "https://github.com/spglib/spglib/archive/refs/tags/v1.11.2.1.tar.gz" CACHE STRING "URL for spglib v1.x tarball") + set(SPGLIB_MD5 "3089782bc85b5034dd4765a18ee70bc7" CACHE STRING "MD5 checksum for spglib tarball") + mark_as_advanced(SPGLIB_URL) + mark_as_advanced(SPGLIB_MD5) + GetFallbackURL(SPGLIB_URL SPGLIB_FALLBACK) + + include(ExternalProject) + ExternalProject_Add(spglib_build + URL ${SPGLIB_URL} ${SPGLIB_FALLBACK} + URL_MD5 ${SPGLIB_MD5} + PREFIX ${CMAKE_CURRENT_BINARY_DIR}/spglib_build_ext + CMAKE_ARGS -DCMAKE_INSTALL_PREFIX=${CMAKE_CURRENT_BINARY_DIR}/spglib_build_ext + -DCMAKE_C_COMPILER=${CMAKE_C_COMPILER} + -DCMAKE_BUILD_TYPE=${CMAKE_BUILD_TYPE} + -DCMAKE_MAKE_PROGRAM=${CMAKE_MAKE_PROGRAM} + -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TOOLCHAIN_FILE} + -DCMAKE_POSITION_INDEPENDENT_CODE=ON + UPDATE_COMMAND ${CMAKE_COMMAND} -E copy_if_different ${CMAKE_CURRENT_SOURCE_DIR}/CMakeLists.spglib ${CMAKE_CURRENT_BINARY_DIR}/spglib_build_ext/src/spglib_build/CMakeLists.txt + INSTALL_COMMAND ${CMAKE_COMMAND} --build ${CMAKE_CURRENT_BINARY_DIR}/spglib_build_ext/src/spglib_build-build --target install + BUILD_BYPRODUCTS "${CMAKE_CURRENT_BINARY_DIR}/spglib_build_ext/lib/${CMAKE_STATIC_LIBRARY_PREFIX}symspg${CMAKE_STATIC_LIBRARY_SUFFIX}" + ) + + # workaround for older CMake versions + file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/spglib_build_ext/lib) + file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/spglib_build_ext/include) + + add_library(SPGLIB::SYMSPG UNKNOWN IMPORTED) + add_dependencies(SPGLIB::SYMSPG spglib_build) + set_target_properties(SPGLIB::SYMSPG PROPERTIES + IMPORTED_LOCATION "${CMAKE_CURRENT_BINARY_DIR}/spglib_build_ext/lib/${CMAKE_STATIC_LIBRARY_PREFIX}symspg${CMAKE_STATIC_LIBRARY_SUFFIX}" + INTERFACE_INCLUDE_DIRECTORIES ${CMAKE_CURRENT_BINARY_DIR}/spglib_build_ext/include + ) + + target_compile_definitions(phana PRIVATE UseSPG) + target_link_libraries(phana PRIVATE SPGLIB::SYMSPG) +endif() + +target_link_libraries(phana PRIVATE tricubic ${LAPACK_LIBRARIES}) +install(TARGETS phana EXPORT LAMMPS_Targets DESTINATION ${CMAKE_INSTALL_BINDIR}) diff --git a/tools/phonon/Makefile b/tools/phonon/Makefile deleted file mode 100644 index 6afae5f312..0000000000 --- a/tools/phonon/Makefile +++ /dev/null @@ -1,67 +0,0 @@ -.SUFFIXES : .o .cpp -# compiler and flags -CC = g++ -Wno-unused-result -LINK = $(CC) -static -CFLAGS = -O3 $(DEBUG) $(UFLAG) - -OFLAGS = -O3 $(DEBUG) -INC = $(LPKINC) $(TCINC) $(SPGINC) $(FFTINC) -LIB = $(LPKLIB) $(TCLIB) $(SPGLIB) $(FFTLIB) - -# cLapack library needed -LPKINC = -I/opt/clapack/3.2.1/include -LPKLIB = -L/opt/clapack/3.2.1/lib -lclapack -lblas -lf2c -lm - -# Tricubic library needed -TCINC = -I/opt/tricubic/1.0/include -TCLIB = -L/opt/tricubic/1.0/lib -ltricubic - -# spglib, 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 -UFLAG = $(SFLAG) $(FFLAG) - -#==================================================================== -ROOT = phana -# executable name -EXE = $(ROOT) -#==================================================================== -# source and rules -SRC = $(wildcard *.cpp) -OBJ = $(SRC:.cpp=.o) - -#==================================================================== -all: ${EXE} - -${EXE}: $(OBJ) - $(LINK) $(OFLAGS) $(OBJ) $(LIB) -o $@ - -clean: - rm -f *.o *~ *.mod ${EXE} - -tar: - rm -f ${ROOT}.tar.gz; tar -czvf ${ROOT}.tar.gz *.cpp *.h Makefile README - -ver: - @echo "#define VERSION `git log|grep '^commit'|wc -l`" > version.h - -#==================================================================== -.f.o: - $(FC) $(FFLAGS) $(FREE) $(MPI) ${INC} -c $< -.f90.o: - $(FC) $(FFLAGS) $(FREE) $(MPI) ${INC} -c $< -.c.o: - $(CC) $(CFLAGS) -c $< -.cpp.o: - $(CC) $(CFLAGS) $(INC) -c $< diff --git a/tools/phonon/README b/tools/phonon/README index db69ac50c8..c1430f5940 100644 --- a/tools/phonon/README +++ b/tools/phonon/README @@ -5,34 +5,38 @@ analyse the phonon related information. #------------------------------------------------------------------------------- 1. Dependencies - The clapack library is needed to solve the eigen problems, - which could be downloaded from: - http://www.netlib.org/clapack/ - - The tricubic library is also needed to do tricubic interpolations, - which could now be obtained from: - https://github.com/nbigaouette/libtricubic/ - + The ZHEEVD LAPACK function is needed to solve the eigen problems. + A C++ compilable version based on CLAPACK is included in the linalg folder + and will be automatically built. + + The tricubic library is also needed to do tricubic interpolations. + A copy is included and will be automatically built. + The spglib is optionally needed, enabling one to evaluate the phonon density of states or vibrational thermal properties using only the irreducible q-points in the first Brillouin zone, as well as to evaluate the phonon dispersion curvers with the - automatic mode. Currently, the 1.8.3 version of spglib is used. - It can be obtained from: - http://spglib.sourceforge.net/ + automatic mode. Currently, version 1.11.2.1 of spglib is used. + It is automatically downloaded and compiled unless the -DUSE_SPGLIB=off + variable is set during CMake configuration. 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 - + It is autodetected and used if available. + + FFTW 3 can be downloaded from: http://www.fftw.org + 2. Compilation - To compile the code, one needs therefore to install the above - libraries and set the paths correctly in the Makefile. - Once this is done, by typing - make - will yield the executable "phana". + To compile the code, one needs to have CMake version 3.16 + or later installed. + + The CMake configuration is done with: + cmake -S . -B build + And compilation then performed with: + cmake --build build + The phana (or phana.exe) executable is then available in + the "build" folder 3. Unit system The units of the output frequencies by this code is THz for @@ -46,6 +50,18 @@ 5. Bug report If any bug found, please drop a line to: konglt(at)sjtu.edu.cn +6. Precompiled executable + The "precompiled" folder contains a precompiled and statically + linked Linux executable for x86_64 CPUs. It should work on *any* + Linux machine with using the x86_64 architecture. It includes + spglib support but not fftw3. + +7. Portability + Build and use of phana has been successfully tested on: + - Fedora Linux 38 using GCC, Clang, and MinGW Linux2Windows cross-compiler + - macOS 12 (Monterey) using Xcode + - Windows 11 using Visual Studio 2022 with MSVC and Clang + #------------------------------------------------------------------------------- Author: Ling-Ti Kong, konglt(at)sjtu.edu.cn -May 2020 +Aug 2021 diff --git a/tools/phonon/disp.cpp b/tools/phonon/disp.cpp index 8a53873383..79a22aeee4 100644 --- a/tools/phonon/disp.cpp +++ b/tools/phonon/disp.cpp @@ -1,10 +1,18 @@ -#include "string.h" -#include "qnodes.h" -#include "global.h" + #include "phonon.h" -#include "green.h" -#include "timer.h" + +#include "dynmat.h" +#include "global.h" +#include "input.h" #include "kpath.h" +#include "qnodes.h" + +#include +#include +#include +#include +#include +#include /*------------------------------------------------------------------------------ * Private method to evaluate the phonon dispersion curves @@ -13,19 +21,22 @@ void Phonon::pdisp() { // ask the output file name and write the header. char str[MAXLINE]; - for (int ii = 0; ii < 80; ++ii) printf("="); printf("\n"); + puts("================================================================================"); + #ifdef UseSPG // ask method to generate q-lines int method = 2; printf("Please select your method to generate the phonon dispersion:\n"); printf(" 1. Manual, should always work;\n"); printf(" 2. Automatic, works only for 3D crystals (CMS49-299).\nYour choice [2]: "); - if (count_words(fgets(str,MAXLINE,stdin)) > 0) method = atoi(strtok(str," \t\n\r\f")); + input->read_stdin(str); + if (count_words(str) > 0) method = atoi(strtok(str," \t\n\r\f")); method = 2 - method%2; printf("Your selection: %d\n", method); #endif printf("\nPlease input the filename to output the dispersion data [pdisp.dat]:"); - if (count_words(fgets(str,MAXLINE,stdin)) < 1) strcpy(str, "pdisp.dat"); + input->read_stdin(str); + if (count_words(str) < 1) strcpy(str, "pdisp.dat"); char *ptr = strtok(str," \t\n\r\f"); char *fname = new char[strlen(ptr)+1]; strcpy(fname,ptr); @@ -45,9 +56,9 @@ void Phonon::pdisp() while (1){ for (int i = 0; i < 3; ++i) qstr[i] = qend[i]; - int quit = 0; printf("\nPlease input the start q-point in unit of B1->B3, q to exit [%g %g %g]: ", qstr[0], qstr[1], qstr[2]); - int n = count_words(fgets(str, MAXLINE, stdin)); + input->read_stdin(str); + int n = count_words(str); ptr = strtok(str, " \t\n\r\f"); if ((n == 1) && (strcmp(ptr,"q") == 0)) break; else if (n >= 3){ @@ -56,14 +67,18 @@ void Phonon::pdisp() qstr[2] = atof(strtok(NULL, " \t\n\r\f")); } - do printf("Please input the end q-point in unit of B1->B3: "); - while (count_words(fgets(str, MAXLINE, stdin)) < 3); + while ( 1 ){ + printf("Please input the end q-point in unit of B1->B3: "); + input->read_stdin(str); + if (count_words(str) >= 3) break; + } qend[0] = atof(strtok(str, " \t\n\r\f")); qend[1] = atof(strtok(NULL, " \t\n\r\f")); qend[2] = atof(strtok(NULL, " \t\n\r\f")); printf("Please input the # of points along the line [%d]: ", nq); - if (count_words(fgets(str, MAXLINE, stdin)) > 0) nq = atoi(strtok(str," \t\n\r\f")); + input->read_stdin(str); + if (count_words(str) > 0) nq = atoi(strtok(str," \t\n\r\f")); nq = MAX(nq,2); double *qtmp = new double [3]; @@ -147,7 +162,7 @@ void Phonon::pdisp() printf("\nPhonon dispersion data are written to: %s, you can visualize the results\n", fname); printf("by invoking: `gnuplot pdisp.gnuplot; gv pdisp.eps`\n"); } - for (int ii = 0; ii < 80; ++ii) printf("="); printf("\n"); + puts("================================================================================"); delete []fname; delete qnodes; diff --git a/tools/phonon/dynmat.cpp b/tools/phonon/dynmat.cpp index a6c4105547..aeb8fe6430 100644 --- a/tools/phonon/dynmat.cpp +++ b/tools/phonon/dynmat.cpp @@ -1,7 +1,16 @@ + #include "dynmat.h" -#include "math.h" -#include "version.h" + #include "global.h" +#include "input.h" +#include "interpolate.h" +#include "memory.h" +#include "version.h" +#include "zheevd.h" + +#include +#include +#include /* ---------------------------------------------------------------------------- * Class DynMat stores the Dynamic Matrix read from the binary file from @@ -9,6 +18,7 @@ * ---------------------------------------------------------------------------- */ DynMat::DynMat(int narg, char **arg) { + input = NULL; attyp = NULL; memory = NULL; M_inv_sqrt = NULL; @@ -19,6 +29,8 @@ DynMat::DynMat(int narg, char **arg) attyp = NULL; basis = NULL; flag_reset_gamma = flag_skip = 0; + symprec = -1.; + int flag_save = 0; // analyze the command line options int iarg = 1; @@ -29,9 +41,16 @@ DynMat::DynMat(int narg, char **arg) } else if (strcmp(arg[iarg], "-r") == 0){ flag_reset_gamma = 1; + } else if (strcmp(arg[iarg], "-p") == 0){ + if (++iarg >= narg) help(); + else symprec = fabs(atof(arg[iarg])); + } else if (strcmp(arg[iarg], "-h") == 0){ help(); + } else if (strcmp(arg[iarg], "-save") == 0){ + flag_save = 1; + } else { if (binfile) delete []binfile; int n = strlen(arg[iarg]) + 1; @@ -43,6 +62,8 @@ DynMat::DynMat(int narg, char **arg) } ShowVersion(); + input = new UserInput(flag_save); + // get the binary file name from user input if not found in command line char str[MAXLINE]; if (binfile == NULL) { @@ -50,7 +71,7 @@ DynMat::DynMat(int narg, char **arg) printf("\n"); do { printf("Please input the binary file name from fix_phonon: "); - fgets(str,MAXLINE,stdin); + input->read_stdin(str); ptr = strtok(str, " \n\t\r\f"); } while (ptr == NULL); @@ -137,17 +158,17 @@ DynMat::DynMat(int narg, char **arg) fclose(fp); exit(3); } - if (fread(basis[0], sizeof(double), fftdim, fp) != fftdim){ + if (fread(basis[0], sizeof(double), fftdim, fp) != (size_t)fftdim){ printf("\nError while reading basis info from file: %s\n", binfile); fclose(fp); exit(3); } - if (fread(&attyp[0], sizeof(int), nucell, fp) != nucell){ + if (fread(&attyp[0], sizeof(int), nucell, fp) != (size_t)nucell){ printf("\nError while reading atom types from file: %s\n", binfile); fclose(fp); exit(3); } - if (fread(&M_inv_sqrt[0], sizeof(double), nucell, fp) != nucell){ + if (fread(&M_inv_sqrt[0], sizeof(double), nucell, fp) != (size_t)nucell){ printf("\nError while reading atomic masses from file: %s\n", binfile); fclose(fp); exit(3); @@ -159,6 +180,7 @@ DynMat::DynMat(int narg, char **arg) // initialize interpolation interpolate = new Interpolate(nx,ny,nz,fftdim2,DM_all); + interpolate->input = input; if (flag_reset_gamma) interpolate->reset_gamma(); // Enforcing Austic Sum Rule @@ -217,7 +239,7 @@ void DynMat::writeDMq(double *q) printf("\n"); while ( 1 ){ printf("Please input the filename to output the DM at selected q: "); - fgets(str,MAXLINE,stdin); + input->read_stdin(str); ptr = strtok(str, " \r\t\n\f"); if (ptr) break; } @@ -264,9 +286,9 @@ void DynMat::writeDMq(double *q, const double qr, FILE *fp) int DynMat::geteigen(double *egv, int flag) { char jobz, uplo; - integer n, lda, lwork, lrwork, *iwork, liwork, info; + int n, lda, lwork, lrwork, *iwork, liwork, info; doublecomplex *work; - doublereal *w = &egv[0], *rwork; + double *w = &egv[0], *rwork; n = fftdim; if (flag) jobz = 'V'; @@ -348,7 +370,7 @@ void DynMat::EnforceASR() char str[MAXLINE]; int nasr = 20; if (nucell <= 1) nasr = 1; - printf("\n"); for (int i = 0; i < 80; ++i) printf("="); + printf("\n================================================================================"); // compute and display eigenvalues of Phi at gamma before ASR if (nucell > 100){ @@ -356,7 +378,7 @@ void DynMat::EnforceASR() fflush(stdout); } - double egvs[fftdim]; + double *egvs = new double[fftdim]; for (int i = 0; i < fftdim; ++i) for (int j = 0; j < fftdim; ++j) DM_q[i][j] = DM_all[0][i*fftdim+j]; geteigen(egvs, 0); @@ -370,11 +392,11 @@ void DynMat::EnforceASR() // ask for iterations to enforce ASR printf("Please input the # of iterations to enforce ASR [%d]: ", nasr); - fgets(str,MAXLINE,stdin); + input->read_stdin(str); char *ptr = strtok(str," \t\n\r\f"); if (ptr) nasr = atoi(ptr); if (nasr < 1){ - for (int i=0; i<80; i++) printf("="); printf("\n"); + puts("================================================================================"); return; } @@ -439,9 +461,8 @@ void DynMat::EnforceASR() if (i%10 == 9) printf("\n"); if (i == 99){ printf("...... (%d more skiped)", fftdim-100); break;} } - printf("\n"); - for (int i = 0; i < 80; ++i) printf("="); printf("\n\n"); - + delete[] egvs; + puts("\n================================================================================\n"); return; } @@ -468,7 +489,7 @@ void DynMat::real2rec() for (int i = 0; i < 9; ++i) ibasevec[i] *= vol; - printf("\n"); for (int i = 0; i < 80; ++i) printf("="); + printf("\n================================================================================"); printf("\nBasis vectors of the unit cell in real space:"); for (int i = 0; i < sysdim; ++i){ printf("\n A%d: ", i+1); @@ -479,8 +500,7 @@ void DynMat::real2rec() printf("\n B%d: ", i+1); for (int j = 0; j < sysdim; ++j) printf("%8.4f ", ibasevec[i*3+j]); } - printf("\n"); for (int i = 0; i < 80; ++i) printf("="); printf("\n"); - + puts("\n================================================================================"); return; } @@ -500,16 +520,17 @@ void DynMat::GaussJordan(int n, double *Mat) indxc = new int[n]; indxr = new int[n]; ipiv = new int[n]; - + + irow = icol = -1; for (i = 0; i < n; ++i) ipiv[i] = 0; for (i = 0; i < n; ++i){ - big = 0.; + big = 0.0; for (j = 0; j < n; ++j){ if (ipiv[j] != 1){ for (k = 0; k < n; ++k){ if (ipiv[k] == 0){ idr = j * n + k; - nmjk = abs(Mat[idr]); + nmjk = fabs(Mat[idr]); if (nmjk >= big){ big = nmjk; irow = j; @@ -602,6 +623,9 @@ void DynMat::help() printf(" will also inform the code to skip all q-points that is in the vicinity\n"); printf(" of the gamma point when evaluating phonon DOS and/or phonon dispersion.\n\n"); printf(" By default, this is not set; and not expected for uncharged systems.\n\n"); + printf(" -p prec To define the precision for symmetry identification with spglib.\n"); + printf(" By default, 1.e-3.\n\n"); + printf(" -save To record user input in `script.inp`, facilitating scripting.\n\n"); printf(" -h To print out this help info.\n\n"); printf(" file To define the filename that carries the binary dynamical matrice generated\n"); printf(" by fix-phonon. If not provided, the code will ask for it.\n"); @@ -680,13 +704,13 @@ void DynMat::Define_Conversion_Factor() * ---------------------------------------------------------------------------- */ void DynMat::ShowInfo() { - printf("\n"); for (int i = 0; i < 80; ++i) printf("="); printf("\n"); + puts("\n================================================================================"); printf("Dynamical matrix is read from file: %s\n", binfile); printf("The system size in three dimension: %d x %d x %d\n", nx, ny, nz); printf("Number of atoms per unit cell : %d\n", nucell); printf("System dimension : %d\n", sysdim); printf("Boltzmann constant in used units : %g\n", boltz); - for (int i = 0; i < 80; ++i) printf("="); printf("\n"); + puts("================================================================================"); return; } /* --------------------------------------------------------------------*/ diff --git a/tools/phonon/dynmat.h b/tools/phonon/dynmat.h index 10caa1ce2a..e11badfde6 100644 --- a/tools/phonon/dynmat.h +++ b/tools/phonon/dynmat.h @@ -1,11 +1,9 @@ #ifndef DYNMAT_H #define DYNMAT_H -#include "stdio.h" -#include "stdlib.h" -#include "string.h" -#include "memory.h" -#include "interpolate.h" +#include "zheevd.h" + +#include class DynMat { public: @@ -15,7 +13,7 @@ public: int nx, ny, nz, nucell; int sysdim, fftdim; - double eml2f, eml2fc; + double eml2f, eml2fc, symprec; char *funit; void getDMq(double *); @@ -34,18 +32,18 @@ public: double **basis; int *attyp; + class UserInput *input; + private: int flag_skip, flag_reset_gamma; - Interpolate *interpolate; - - Memory *memory; + class Interpolate *interpolate; + class Memory *memory; - int nasr; void EnforceASR(); char *binfile, *dmfile; - double boltz, q[3]; + double boltz; doublecomplex **DM_all; diff --git a/tools/phonon/global.h b/tools/phonon/global.h index c20a91f317..bd1608ee07 100644 --- a/tools/phonon/global.h +++ b/tools/phonon/global.h @@ -1,7 +1,7 @@ #ifndef GLOBAL_H #define GLOBAL_H -#define ZERO 1.e-8 +#define ZERO 1.0e-8 #define MAXLINE 512 #define MIN(a,b) ((a)>(b)?(b):(a)) diff --git a/tools/phonon/green.cpp b/tools/phonon/green.cpp index 41f5b14886..0f687fd834 100644 --- a/tools/phonon/green.cpp +++ b/tools/phonon/green.cpp @@ -1,10 +1,10 @@ -#include -#include -#include -#include #include "green.h" + +#include "memory.h" + #include -#include "global.h" +#include +#include /******************************************************************************* * The class of Green is designed to evaluate the LDOS via the Green's Function @@ -59,7 +59,6 @@ Green::Green(const int ntm, const int sdim, const int niter, const double min, c 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(); @@ -224,8 +223,6 @@ void Green::recursion() { // local variables std::complex Z, rec_x, rec_x_inv; - std::complex cunit = std::complex(0.,1.); - double w = wmin; for (int i = 0; i < nw; ++i){ diff --git a/tools/phonon/green.h b/tools/phonon/green.h index 1a137e53ba..21e3b091dc 100644 --- a/tools/phonon/green.h +++ b/tools/phonon/green.h @@ -1,8 +1,6 @@ #ifndef GREEN_H #define GREEN_H -#include "memory.h" - class Green{ public: Green(const int, const int, const int, const double, const double, @@ -14,12 +12,11 @@ private: 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; + class Memory *memory; }; #endif diff --git a/tools/phonon/input.cpp b/tools/phonon/input.cpp new file mode 100644 index 0000000000..c2059043c7 --- /dev/null +++ b/tools/phonon/input.cpp @@ -0,0 +1,35 @@ +#include "input.h" + +#include "global.h" + +/* ------------------------------------------------------------------- + * Constructor. If flag = 1, output user inputs as script.inp + * ---------------------------------------------------------------- */ +UserInput::UserInput(int flag) +{ + fp = NULL; + if (flag) fp = fopen("script.inp", "w"); + + return; +} + +/* ------------------------------------------------------------------- + * Deconstructor. Output user inputs as required and clear workspace. + * ---------------------------------------------------------------- */ +UserInput::~UserInput() +{ + if (fp) fclose(fp); + fp = NULL; +} + +/* ------------------------------------------------------------------- + * Read stdin and keep a record of it. + * ---------------------------------------------------------------- */ +void UserInput::read_stdin(char *str) +{ + fgets(str, MAXLINE, stdin); + if (fp) fprintf(fp, "%s", str); + + return; +} +/* ---------------------------------------------------------------- */ diff --git a/tools/phonon/input.h b/tools/phonon/input.h new file mode 100644 index 0000000000..c931895c4d --- /dev/null +++ b/tools/phonon/input.h @@ -0,0 +1,17 @@ +#ifndef INPUT_H +#define INPUT_H + +#include + +class UserInput { +public: + UserInput(int); + ~UserInput(); + + void read_stdin(char *); + +private: + FILE *fp; + +}; +#endif diff --git a/tools/phonon/interpolate.cpp b/tools/phonon/interpolate.cpp index 09e261c763..8ea551d1a8 100644 --- a/tools/phonon/interpolate.cpp +++ b/tools/phonon/interpolate.cpp @@ -1,6 +1,14 @@ + #include "interpolate.h" -#include "math.h" + #include "global.h" +#include "input.h" +#include "memory.h" +#include "tricubic.h" + +#include +#include +#include /* ---------------------------------------------------------------------------- * Constructor used to get info from caller, and prepare other necessary data @@ -19,6 +27,7 @@ Interpolate::Interpolate(int nx, int ny, int nz, int ndm, doublecomplex **DM) data = DM; Dfdx = Dfdy = Dfdz = D2fdxdy = D2fdxdz = D2fdydz = D3fdxdydz = NULL; flag_reset_gamma = flag_allocated_dfs = 0; + input = NULL; return; } @@ -265,17 +274,19 @@ 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"); + if (input == NULL) input = new UserInput(0); + + puts("\n================================================================================"); + printf("Which interpolation method would you like to use?\n"); printf(" 1. Tricubic;\n 2. Trilinear;\n"); printf("Your choice [1]: "); - fgets(str,MAXLINE,stdin); + input->read_stdin(str); 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"); + puts("================================================================================\n"); if (which == 1) tricubic_init(); diff --git a/tools/phonon/interpolate.h b/tools/phonon/interpolate.h index c650b30908..b9e0242b96 100644 --- a/tools/phonon/interpolate.h +++ b/tools/phonon/interpolate.h @@ -1,16 +1,7 @@ #ifndef INTERPOLATION_H #define INTERPOLATION_H -#include "stdio.h" -#include "stdlib.h" -#include "string.h" -#include "memory.h" -#include "tricubic.h" - -extern "C"{ -#include "f2c.h" -#include "clapack.h" -} +#include "zheevd.h" class Interpolate{ public: @@ -23,11 +14,13 @@ public: int UseGamma; + class UserInput *input; + private: void tricubic_init(); void tricubic(double *, doublecomplex *); void trilinear(double *, doublecomplex *); - Memory *memory; + class Memory *memory; int which; int Nx, Ny, Nz, Npt, ndim; diff --git a/tools/phonon/kpath.cpp b/tools/phonon/kpath.cpp index 842e680611..49730b42b6 100644 --- a/tools/phonon/kpath.cpp +++ b/tools/phonon/kpath.cpp @@ -1,11 +1,19 @@ -#include "global.h" + #include "kpath.h" +#include "global.h" +#include "dynmat.h" +#include "memory.h" +#include "qnodes.h" + #ifdef UseSPG extern "C"{ #include "spglib.h" } -#include "math.h" +#include +#include +#include +#include /* ---------------------------------------------------------------------------- * Class kPath will help to find the high symmetry k-path for a given lattice. @@ -47,11 +55,15 @@ kPath::kPath(DynMat *dm, QNodes *qn) for (int idim = 0; idim < sysdim; ++idim) atpos[i][idim] = dynmat->basis[i][idim]; // get the space group number - double symprec = 1.e-4, pos[num_atom][3]; + double symprec = 1.0e-3; + double **pos; + memory->create(pos,num_atom,3,"kpath:pos"); + if (dynmat->symprec > 0.0) symprec = dynmat->symprec; + for (int i = 0; i < num_atom; ++i) for (int j = 0; j < 3; ++j) pos[i][j] = atpos[i][j]; - spgnum = spg_get_international(symbol, latvec, pos, attyp, num_atom, symprec); - + spgnum = spg_get_international(symbol, latvec, (double (*)[3])pos, attyp, num_atom, symprec); + memory->destroy(pos); return; } @@ -61,7 +73,7 @@ kPath::kPath(DynMat *dm, QNodes *qn) void kPath::show_info() { // display the unit cell info read - for (int ii = 0; ii < 80; ++ii) printf("-"); printf("\n"); + puts("--------------------------------------------------------------------------------"); printf("The basis vectors of the unit cell:\n"); for (int idim = 0; idim < 3; ++idim){ printf(" A%d =", idim+1); @@ -76,12 +88,10 @@ void kPath::show_info() if (num_atom > NUMATOM) printf(" ... (%d atoms omitted.)\n", num_atom-NUMATOM); printf("The space group number of your unit cell is: %d => %s\n", spgnum, symbol); - for (int ii = 0; ii < 80; ++ii) printf("-"); printf("\n"); - + puts("--------------------------------------------------------------------------------"); return; } - /* ---------------------------------------------------------------------------- * Free the memeory used by kPath. * ---------------------------------------------------------------------------- */ @@ -2765,9 +2775,16 @@ void kPath::show_path() if (q == NULL) return; int nbin = q->ndstr.size(); if (nbin > 0){ - printf("\nk-path for the current lattice will be:\n\t%s", q->ndstr[0].c_str()); + puts("\n--------------------------------------------------------------------------------"); + printf("k-path for the current lattice will be:\n %s", q->ndstr[0].c_str()); for (int is = 1; is < nbin; ++is) printf("-%s", q->ndstr[is].c_str()); - printf("\n"); + + printf("\n\nThe fractional coordinates of these paths are:\n"); + for (int is = 0; is < nbin-1; ++is) + printf(" [%6.4f %6.4f %6.4f] --> [%6.4f %6.4f %6.4f] (%s - %s)\n", q->qs[is][0], + q->qs[is][1], q->qs[is][2], q->qe[is][0], q->qe[is][1], q->qe[is][2], + q->ndstr[is].c_str(), q->ndstr[is+1].c_str() ); + puts("--------------------------------------------------------------------------------"); } return; diff --git a/tools/phonon/kpath.h b/tools/phonon/kpath.h index bbcedf15d4..ade1e3a27c 100644 --- a/tools/phonon/kpath.h +++ b/tools/phonon/kpath.h @@ -4,14 +4,9 @@ #ifndef KPATH_H #define KPATH_H -#include "qnodes.h" -#include "dynmat.h" -#include "memory.h" - class kPath{ public: - - kPath(DynMat *, QNodes *); + kPath(class DynMat *, class QNodes *); ~kPath(); void kpath(); @@ -19,13 +14,11 @@ public: void show_info(); private: - - Memory *memory; - - DynMat *dynmat; - QNodes *q; + class Memory *memory; + class DynMat *dynmat; + class QNodes *q; char symbol[11]; - int spgnum, sysdim, fftdim, num_atom, *attyp; + int spgnum, sysdim, num_atom, *attyp; double latvec[3][3], **atpos; }; diff --git a/tools/phonon/main.cpp b/tools/phonon/main.cpp index d7d5baa1cc..d0193a037c 100644 --- a/tools/phonon/main.cpp +++ b/tools/phonon/main.cpp @@ -1,10 +1,6 @@ -#include "stdio.h" -#include "stdlib.h" #include "dynmat.h" #include "phonon.h" -using namespace std; - int main(int argc, char** argv) { diff --git a/tools/phonon/memory.cpp b/tools/phonon/memory.cpp index 4d65e83ef8..18dfa2fa43 100644 --- a/tools/phonon/memory.cpp +++ b/tools/phonon/memory.cpp @@ -1,8 +1,8 @@ -#include "stdio.h" -#include "stdlib.h" -#include "string.h" #include "memory.h" +#include +#include + /* ---------------------------------------------------------------------- safe malloc ------------------------------------------------------------------------- */ diff --git a/tools/phonon/memory.h b/tools/phonon/memory.h index ae2feceba3..13eeca4b14 100644 --- a/tools/phonon/memory.h +++ b/tools/phonon/memory.h @@ -4,18 +4,16 @@ #define __STDC_LIMIT_MACROS #define __STDC_FORMAT_MACROS -#include "stdio.h" -#include "stdlib.h" -#include "limits.h" -#include "stdint.h" -#include "inttypes.h" +#include +#include +#include typedef int64_t bigint; #define BIGINT_FORMAT "%" PRId64 #define ATOBIGINT atoll class Memory { - public: +public: Memory(){}; void *smalloc(bigint n, const char *); @@ -24,11 +22,11 @@ class Memory { void fail(const char *); /* ---------------------------------------------------------------------- - create a 1d array -------------------------------------------------------------------------- */ + create a 1d array + ------------------------------------------------------------------------- */ template - TYPE *create(TYPE *&array, int n, const char *name) + TYPE *create(TYPE *&array, int n, const char *name) { bigint nbytes = sizeof(TYPE) * n; array = (TYPE *) smalloc(nbytes,name); @@ -36,42 +34,42 @@ class Memory { }; template - TYPE **create(TYPE **&array, int n, const char *name) {fail(name);} + TYPE **create(TYPE **&, int, const char *name) {fail(name);} /* ---------------------------------------------------------------------- grow or shrink 1d array -------------------------------------------------------------------------- */ + ------------------------------------------------------------------------- */ template - TYPE *grow(TYPE *&array, int n, const char *name) + TYPE *grow(TYPE *&array, int n, const char *name) { if (array == NULL) return create(array,n,name); - + bigint nbytes = sizeof(TYPE) * n; array = (TYPE *) srealloc(array,nbytes,name); return array; }; template - TYPE **grow(TYPE **&array, int n, const char *name) {fail(name);} + TYPE **grow(TYPE **&, int, const char *name) {fail(name);} /* ---------------------------------------------------------------------- - destroy a 1d array -------------------------------------------------------------------------- */ + destroy a 1d array + ------------------------------------------------------------------------- */ template - void destroy(TYPE *array) + void destroy(TYPE *array) { sfree(array); }; /* ---------------------------------------------------------------------- - create a 1d array with index from nlo to nhi inclusive + create a 1d array with index from nlo to nhi inclusive cannot grow it -------------------------------------------------------------------------- */ + ------------------------------------------------------------------------- */ template - TYPE *create1d_offset(TYPE *&array, int nlo, int nhi, const char *name) + TYPE *create1d_offset(TYPE *&array, int nlo, int nhi, const char *name) { bigint nbytes = sizeof(TYPE) * (nhi-nlo+1); array = (TYPE *) smalloc(nbytes,name); @@ -80,76 +78,76 @@ class Memory { } template - TYPE **create1d_offset(TYPE **&array, int nlo, int nhi, const char *name) + TYPE **create1d_offset(TYPE **&, int, int, const char *name) {fail(name);} /* ---------------------------------------------------------------------- - destroy a 1d array with index offset -------------------------------------------------------------------------- */ + destroy a 1d array with index offset + ------------------------------------------------------------------------- */ template - void destroy1d_offset(TYPE *array, int offset) + void destroy1d_offset(TYPE *array, int offset) { if (array) sfree(&array[offset]); } /* ---------------------------------------------------------------------- - create a 2d array -------------------------------------------------------------------------- */ + create a 2d array + ------------------------------------------------------------------------- */ template - TYPE **create(TYPE **&array, int n1, int n2, const char *name) + TYPE **create(TYPE **&array, int n1, int n2, const char *name) { bigint nbytes = sizeof(TYPE) * n1*n2; TYPE *data = (TYPE *) smalloc(nbytes,name); nbytes = sizeof(TYPE *) * n1; array = (TYPE **) smalloc(nbytes,name); - + int n = 0; for (int i = 0; i < n1; i++) { - array[i] = &data[n]; - n += n2; + array[i] = &data[n]; + n += n2; } return array; } template - TYPE ***create(TYPE ***&array, int n1, int n2, const char *name) + TYPE ***create(TYPE ***&, int, int, const char *name) {fail(name);} /* ---------------------------------------------------------------------- grow or shrink 1st dim of a 2d array last dim must stay the same -------------------------------------------------------------------------- */ + ------------------------------------------------------------------------- */ template - TYPE **grow(TYPE **&array, int n1, int n2, const char *name) + TYPE **grow(TYPE **&array, int n1, int n2, const char *name) { if (array == NULL) return create(array,n1,n2,name); - + bigint nbytes = sizeof(TYPE) * n1*n2; TYPE *data = (TYPE *) srealloc(array[0],nbytes,name); nbytes = sizeof(TYPE *) * n1; array = (TYPE **) srealloc(array,nbytes,name); - + int n = 0; for (int i = 0; i < n1; i++) { - array[i] = &data[n]; - n += n2; + array[i] = &data[n]; + n += n2; } return array; } template - TYPE ***grow(TYPE ***&array, int n1, int n2, const char *name) + TYPE ***grow(TYPE ***&, int, int, const char *name) {fail(name);} /* ---------------------------------------------------------------------- - destroy a 2d array -------------------------------------------------------------------------- */ + destroy a 2d array + ------------------------------------------------------------------------- */ template - void destroy(TYPE **array) + void destroy(TYPE **array) { if (array == NULL) return; sfree(array[0]); @@ -157,42 +155,11 @@ class Memory { } /* ---------------------------------------------------------------------- - create a 2d array with 2nd index from n2lo to n2hi inclusive - cannot grow it -------------------------------------------------------------------------- */ + create a 3d array + ------------------------------------------------------------------------- */ template - TYPE **create2d_offset(TYPE **&array, int n1, int n2lo, int n2hi, - const char *name) - { - int n2 = n2hi - n2lo + 1; - create(array,n1,n2,name); - for (int i = 0; i < n1; i++) array[i] -= n2lo; - return array; - } - - template - TYPE ***create2d_offset(TYPE ***&array, int n1, int n2lo, int n2hi, - const char *name) {fail(name);} - -/* ---------------------------------------------------------------------- - destroy a 2d array with 2nd index offset -------------------------------------------------------------------------- */ - - template - void destroy2d_offset(TYPE **array, int offset) - { - if (array == NULL) return; - sfree(&array[0][offset]); - sfree(array); - } - -/* ---------------------------------------------------------------------- - create a 3d array -------------------------------------------------------------------------- */ - - template - TYPE ***create(TYPE ***&array, int n1, int n2, int n3, const char *name) + TYPE ***create(TYPE ***&array, int n1, int n2, int n3, const char *name) { bigint nbytes = sizeof(TYPE) * n1*n2*n3; TYPE *data = (TYPE *) smalloc(nbytes,name); @@ -200,62 +167,62 @@ class Memory { TYPE **plane = (TYPE **) smalloc(nbytes,name); nbytes = sizeof(TYPE **) * n1; array = (TYPE ***) smalloc(nbytes,name); - + int i,j; int n = 0; for (i = 0; i < n1; i++) { - array[i] = &plane[i*n2]; - for (j = 0; j < n2; j++) { - plane[i*n2+j] = &data[n]; - n += n3; - } + array[i] = &plane[i*n2]; + for (j = 0; j < n2; j++) { + plane[i*n2+j] = &data[n]; + n += n3; + } } return array; } template - TYPE ****create(TYPE ****&array, int n1, int n2, int n3, const char *name) + TYPE ****create(TYPE ****&, int, int, int, const char *name) {fail(name);} /* ---------------------------------------------------------------------- grow or shrink 1st dim of a 3d array last 2 dims must stay the same -------------------------------------------------------------------------- */ + ------------------------------------------------------------------------- */ template - TYPE ***grow(TYPE ***&array, int n1, int n2, int n3, const char *name) + TYPE ***grow(TYPE ***&array, int n1, int n2, int n3, const char *name) { if (array == NULL) return create(array,n1,n2,n3,name); - + bigint nbytes = sizeof(TYPE) * n1*n2*n3; TYPE *data = (TYPE *) srealloc(array[0][0],nbytes,name); nbytes = sizeof(TYPE *) * n1*n2; TYPE **plane = (TYPE **) srealloc(array[0],nbytes,name); nbytes = sizeof(TYPE **) * n1; array = (TYPE ***) srealloc(array,nbytes,name); - + int i,j; int n = 0; for (i = 0; i < n1; i++) { - array[i] = &plane[i*n2]; - for (j = 0; j < n2; j++) { - plane[i*n2+j] = &data[n]; - n += n3; - } + array[i] = &plane[i*n2]; + for (j = 0; j < n2; j++) { + plane[i*n2+j] = &data[n]; + n += n3; + } } return array; } template - TYPE ****grow(TYPE ****&array, int n1, int n2, int n3, const char *name) + TYPE ****grow(TYPE ****&, int, int, int, const char *name) {fail(name);} /* ---------------------------------------------------------------------- - destroy a 3d array -------------------------------------------------------------------------- */ + destroy a 3d array + ------------------------------------------------------------------------- */ template - void destroy(TYPE ***array) + void destroy(TYPE ***array) { if (array == NULL) return; sfree(array[0][0]); @@ -263,168 +230,23 @@ class Memory { sfree(array); } -/* ---------------------------------------------------------------------- - create a 3d array with 1st index from n1lo to n1hi inclusive - cannot grow it -------------------------------------------------------------------------- */ - - template - TYPE ***create3d_offset(TYPE ***&array, int n1lo, int n1hi, - int n2, int n3, const char *name) - { - int n1 = n1hi - n1lo + 1; - create(array,n1,n2,n3,name); - array -= n1lo; - return array; - } - - template - TYPE ****create3d_offset(TYPE ****&array, int n1lo, int n1hi, - int n2, int n3, const char *name) - {fail(name);} - -/* ---------------------------------------------------------------------- - free a 3d array with 1st index offset -------------------------------------------------------------------------- */ - - template - void destroy3d_offset(TYPE ***array, int offset) - { - if (array) destroy(&array[offset]); - } - -/* ---------------------------------------------------------------------- - create a 3d array with - 1st index from n1lo to n1hi inclusive, - 2nd index from n2lo to n2hi inclusive, - 3rd index from n3lo to n3hi inclusive - cannot grow it -------------------------------------------------------------------------- */ - - template - TYPE ***create3d_offset(TYPE ***&array, int n1lo, int n1hi, - int n2lo, int n2hi, int n3lo, int n3hi, - const char *name) - { - int n1 = n1hi - n1lo + 1; - int n2 = n2hi - n2lo + 1; - int n3 = n3hi - n3lo + 1; - create(array,n1,n2,n3,name); - - for (int i = 0; i < n1*n2; i++) array[0][i] -= n3lo; - for (int i = 0; i < n1; i++) array[i] -= n2lo; - array -= n1lo; - return array; - } - - template - TYPE ****create3d_offset(TYPE ****&array, int n1lo, int n1hi, - int n2lo, int n2hi, int n3lo, int n3hi, - const char *name) - {fail(name);} - -/* ---------------------------------------------------------------------- - free a 3d array with all 3 indices offset -------------------------------------------------------------------------- */ - - template - void destroy3d_offset(TYPE ***array, - int n1_offset, int n2_offset, int n3_offset) - { - if (array == NULL) return; - sfree(&array[n1_offset][n2_offset][n3_offset]); - sfree(&array[n1_offset][n2_offset]); - sfree(&array[n1_offset]); - } - -/* ---------------------------------------------------------------------- - create a 4d array -------------------------------------------------------------------------- */ - - template - TYPE ****create(TYPE ****&array, int n1, int n2, int n3, int n4, - const char *name) - { - bigint nbytes = sizeof(TYPE) * n1*n2*n3*n4; - TYPE *data = (double *) smalloc(nbytes,name); - nbytes = sizeof(TYPE *) * n1*n2*n3; - TYPE **cube = (double **) smalloc(nbytes,name); - nbytes = sizeof(TYPE **) * n1*n2; - TYPE ***plane = (double ***) smalloc(nbytes,name); - nbytes = sizeof(TYPE ***) * n1; - array = (double ****) smalloc(nbytes,name); - - int i,j,k; - int n = 0; - for (i = 0; i < n1; i++) { - array[i] = &plane[i*n2]; - for (j = 0; j < n2; j++) { - plane[i*n2+j] = &cube[i*n2*n3+j*n3]; - for (k = 0; k < n3; k++) { - cube[i*n2*n3+j*n3+k] = &data[n]; - n += n4; - } - } - } - return array; - } - - template - TYPE *****create(TYPE *****&array, int n1, int n2, int n3, int n4, - const char *name) - {fail(name);} - -/* ---------------------------------------------------------------------- - destroy a 4d array -------------------------------------------------------------------------- */ - - template - void destroy(TYPE ****array) - { - if (array == NULL) return; - sfree(array[0][0][0]); - sfree(array[0][0]); - sfree(array[0]); - sfree(array); - } - /* ---------------------------------------------------------------------- memory usage of arrays, including pointers -------------------------------------------------------------------------- */ + ------------------------------------------------------------------------- */ template - bigint usage(TYPE *array, int n) + bigint usage(TYPE *, int n) { bigint bytes = sizeof(TYPE) * n; return bytes; } template - bigint usage(TYPE **array, int n1, int n2) + bigint usage(TYPE **, int n1, int n2) { bigint bytes = sizeof(TYPE) * n1*n2; bytes += sizeof(TYPE *) * n1; return bytes; } - - template - bigint usage(TYPE ***array, int n1, int n2, int n3) - { - bigint bytes = sizeof(TYPE) * n1*n2*n3; - bytes += sizeof(TYPE *) * n1*n2; - bytes += sizeof(TYPE **) * n1; - return bytes; - } - - template - bigint usage(TYPE ****array, int n1, int n2, int n3, int n4) - { - bigint bytes = sizeof(TYPE) * n1*n2*n3*n4; - bytes += sizeof(TYPE *) * n1*n2*n3; - bytes += sizeof(TYPE **) * n1*n2; - bytes += sizeof(TYPE ***) * n1; - return bytes; - } }; - #endif diff --git a/tools/phonon/phonon.cpp b/tools/phonon/phonon.cpp index 56fa409e06..06372dcd1b 100644 --- a/tools/phonon/phonon.cpp +++ b/tools/phonon/phonon.cpp @@ -1,9 +1,18 @@ -#include -#include "string.h" + #include "phonon.h" -#include "green.h" -#include "timer.h" + #include "global.h" +#include "dynmat.h" +#include "green.h" +#include "input.h" +#include "memory.h" +#include "timer.h" +#include "zheevd.h" + +#include +#include +#include +#include #ifdef UseSPG extern "C"{ @@ -27,6 +36,7 @@ Phonon::Phonon(DynMat *dm) dynmat = dm; sysdim = dynmat->sysdim; ndim = dynmat->fftdim; + input = dm->input; dos = NULL; ldos = NULL; qpts = NULL; @@ -42,10 +52,7 @@ Phonon::Phonon(DynMat *dm) // display the menu char str[MAXLINE]; while ( 1 ){ - printf("\n"); - for (int i = 0; i < 37; ++i) printf("="); - printf(" Menu "); - for (int i = 0; i < 37; ++i) printf("="); printf("\n"); + puts("\n===================================== Menu ====================================="); printf(" 1. Phonon DOS evaluation;\n"); printf(" 2. Phonon dispersion curves;\n"); printf(" 3. Dynamical matrix at arbitrary q;\n"); @@ -64,9 +71,10 @@ Phonon::Phonon(DynMat *dm) // read user choice int job = 0; printf("Your choice [0]: "); - if (count_words(fgets(str,MAXLINE,stdin)) > 0) job = atoi(strtok(str," \t\n\r\f")); + input->read_stdin(str); + if (count_words(str) > 0) job = atoi(strtok(str," \t\n\r\f")); printf("\nYour selection: %d\n", job); - for (int i = 0; i < 80; ++i) printf("=");printf("\n\n"); + puts("================================================================================\n"); // now to do the job according to user's choice if (job == 1) pdos(); @@ -138,7 +146,8 @@ void Phonon::pdos() // Now to ask for the output frequency range printf("\nThe frequency range of all q-points are: [%g %g]\n", fmin, fmax); printf("Please input the desired range to get DOS [%g %g]: ", fmin, fmax); - if (count_words(fgets(str,MAXLINE,stdin)) >= 2){ + input->read_stdin(str); + if (count_words(str) >= 2){ fmin = atof(strtok(str," \t\n\r\f")); fmax = atof(strtok(NULL," \t\n\r\f")); } @@ -147,7 +156,8 @@ void Phonon::pdos() ndos = 201; printf("Please input the number of intervals [%d]: ", ndos); - if (count_words(fgets(str,MAXLINE,stdin)) > 0) ndos = atoi(strtok(str," \t\n\r\f")); + input->read_stdin(str); + if (count_words(str) > 0) ndos = atoi(strtok(str," \t\n\r\f")); ndos += (ndos+1)%2; ndos = MAX(2,ndos); @@ -170,7 +180,8 @@ void Phonon::pdos() // smooth dos ? printf("Would you like to smooth the phonon dos? (y/n)[n]: "); - if (count_words(fgets(str,MAXLINE,stdin)) > 0){ + input->read_stdin(str); + if (count_words(str) > 0){ char *flag = strtok(str," \t\n\r\f"); if (strcmp(flag,"y") == 0 || strcmp(flag,"Y") == 0) smooth(dos, ndos); } @@ -194,7 +205,8 @@ void Phonon::writeDOS() char str[MAXLINE]; // now to output the phonon DOS printf("\nPlease input the filename to write DOS [pdos.dat]: "); - if (count_words(fgets(str,MAXLINE,stdin)) < 1) strcpy(str, "pdos.dat"); + input->read_stdin(str); + if (count_words(str) < 1) strcpy(str, "pdos.dat"); char *fname = strtok(str," \t\n\r\f"); printf("The total phonon DOS will be written to file: %s\n", fname); @@ -234,7 +246,7 @@ void Phonon::writeLDOS() const double one3 = 1./double(sysdim); char str[MAXLINE]; for (int ilocal = 0; ilocal < nlocal; ++ilocal){ - sprintf(str,"pldos_%d.dat", locals[ilocal]); + snprintf(str, MAXLINE-1, "pldos_%d.dat", locals[ilocal]); char *fname = strtok(str," \t\n\r\f"); FILE *fp = fopen(fname, "w"); fname = NULL; @@ -281,7 +293,7 @@ void Phonon::ldos_rsgf() fmin = fmax = egvs[0]; for (int i = 1; i < ndim; ++i){fmin = MIN(fmin, egvs[i]); fmax = MAX(fmax, egvs[i]);} - delete []egvs; + delete[] egvs; } else { fmin = 0.; fmax = 20.; @@ -297,7 +309,8 @@ void Phonon::ldos_rsgf() printf("\nThere are %d atoms in each unit cell of your lattice.\n", dynmat->nucell); printf("Please input the index/index range/index range and increment of atom(s)\n"); printf("in the unit cell to evaluate LDOS, q to exit [%d]: ", ik); - int nr = count_words( fgets(str,MAXLINE,stdin) ); + input->read_stdin(str); + int nr = count_words(str); if (nr < 1){ istr = iend = ik; iinc = 1; @@ -327,7 +340,8 @@ void Phonon::ldos_rsgf() } printf("Please input the frequency range to evaluate LDOS [%g %g]: ", fmin, fmax); - if (count_words(fgets(str,MAXLINE,stdin)) >= 2){ + input->read_stdin(str); + if (count_words(str) >= 2){ fmin = atof(strtok(str," \t\n\r\f")); fmax = atof(strtok(NULL," \t\n\r\f")); } @@ -335,16 +349,19 @@ void Phonon::ldos_rsgf() printf("The frequency range for your LDOS is [%g %g].\n", fmin, fmax); printf("Please input the desired number of points in LDOS [%d]: ", ndos); - if (count_words(fgets(str,MAXLINE,stdin)) > 0) ndos = atoi(strtok(str," \t\n\r\f")); + input->read_stdin(str); + if (count_words(str) > 0) ndos = atoi(strtok(str," \t\n\r\f")); if (ndos < 2) break; ndos += (ndos+1)%2; printf("Please input the maximum # of Lanczos iterations [%d]: ", nit); - if (count_words(fgets(str,MAXLINE,stdin)) > 0) nit = atoi(strtok(str," \t\n\r\f")); + input->read_stdin(str); + if (count_words(str) > 0) nit = atoi(strtok(str," \t\n\r\f")); if (nit < 1) break; printf("Please input the value of epsilon for delta-function [%g]: ", eps); - if (count_words(fgets(str,MAXLINE,stdin)) > 0) eps = atof(strtok(str," \t\n\r\f")); + input->read_stdin(str); + if (count_words(str) > 0) eps = atof(strtok(str," \t\n\r\f")); if (eps <= 0.) break; // prepare array for local pdos @@ -395,8 +412,11 @@ void Phonon::dmanyq() { char str[MAXLINE]; double q[3]; - do printf("Please input the q-point to output the dynamical matrix:"); - while (count_words(fgets(str,MAXLINE,stdin)) < 3); + while ( 1 ){ + printf("Please input the q-point to output the dynamical matrix: "); + input->read_stdin(str); + if (count_words(str) >= 3) break; + } q[0] = atof(strtok(str," \t\n\r\f")); q[1] = atof(strtok(NULL," \t\n\r\f")); q[2] = atof(strtok(NULL," \t\n\r\f")); @@ -413,11 +433,13 @@ void Phonon::dmanyq() void Phonon::vfanyq() { char str[MAXLINE]; - double q[3], egvs[ndim]; + double q[3]; + double *egvs = new double[ndim]; while ( 1 ){ printf("Please input the q-point to compute the frequencies, q to exit: "); - if (count_words(fgets(str,MAXLINE,stdin)) < 3) break; + input->read_stdin(str); + if (count_words(str) < 3) break; q[0] = atof(strtok(str, " \t\n\r\f")); q[1] = atof(strtok(NULL," \t\n\r\f")); @@ -427,9 +449,11 @@ void Phonon::vfanyq() dynmat->geteigen(egvs, 0); printf("q-point: [%lg %lg %lg], ", q[0], q[1], q[2]); printf("vibrational frequencies at this q-point:\n"); - for (int i = 0; i < ndim; ++i) printf("%lg ", egvs[i]); printf("\n\n"); + for (int i = 0; i < ndim; ++i) printf("%lg ", egvs[i]); + printf("\n\n"); } - + + delete[] egvs; return; } @@ -439,15 +463,18 @@ void Phonon::vfanyq() void Phonon::vecanyq() { char str[MAXLINE]; - double q[3], egvs[ndim]; + double q[3]; + double *egvs = new double[ndim]; doublecomplex **eigvec = dynmat->DM_q; printf("Please input the filename to output the result [eigvec.dat]: "); - if (count_words(fgets(str,MAXLINE,stdin)) < 1) strcpy(str,"eigvec.dat"); + input->read_stdin(str); + if (count_words(str) < 1) strcpy(str,"eigvec.dat"); FILE *fp = fopen(strtok(str," \t\n\r\f"), "w"); while ( 1 ){ printf("Please input the q-point to compute the frequencies, q to exit: "); - if (count_words(fgets(str,MAXLINE,stdin)) < 3) break; + input->read_stdin(str); + if (count_words(str) < 3) break; q[0] = atof(strtok(str, " \t\n\r\f")); q[1] = atof(strtok(NULL," \t\n\r\f")); @@ -475,6 +502,7 @@ void Phonon::vecanyq() fprintf(fp,"\n"); } fclose(fp); + delete[] egvs; eigvec = NULL; return; } @@ -488,7 +516,8 @@ void Phonon::DMdisp() char str[MAXLINE]; printf("Please input the filename to output the DM data [DMDisp.dat]: "); - if (count_words(fgets(str,MAXLINE,stdin)) < 1) strcpy(str, "DMDisp.dat"); + input->read_stdin(str); + if (count_words(str) < 1) strcpy(str, "DMDisp.dat"); char *fname = strtok(str," \t\n\r\f"); FILE *fp = fopen(fname, "w"); fname = NULL; @@ -503,7 +532,8 @@ void Phonon::DMdisp() for (int i = 0; i < 3; ++i) qstr[i] = qend[i]; printf("\nPlease input the start q-point in unit of B1->B3, q to exit [%g %g %g]: ", qstr[0], qstr[1], qstr[2]); - int n = count_words(fgets(str,MAXLINE,stdin)); + input->read_stdin(str); + int n = count_words(str); char *ptr = strtok(str," \t\n\r\f"); if ((n == 1) && (strcmp(ptr,"q") == 0)) break; else if (n >= 3){ @@ -512,14 +542,18 @@ void Phonon::DMdisp() qstr[2] = atof(strtok(NULL," \t\n\r\f")); } - do printf("Please input the end q-point in unit of B1->B3: "); - while (count_words(fgets(str,MAXLINE,stdin)) < 3); + while ( 1 ){ + printf("Please input the end q-point in unit of B1->B3: "); + input->read_stdin(str); + if (count_words(str) >= 3) break; + } qend[0] = atof(strtok(str," \t\n\r\f")); qend[1] = atof(strtok(NULL," \t\n\r\f")); qend[2] = atof(strtok(NULL," \t\n\r\f")); printf("Please input the # of points along the line [%d]: ", nq); - if (count_words(fgets(str,MAXLINE,stdin)) > 0) nq = atoi(strtok(str," \t\n\r\f")); + input->read_stdin(str); + if (count_words(str) > 0) nq = atoi(strtok(str," \t\n\r\f")); nq = MAX(nq,2); for (int i=0; i<3; i++) qinc[i] = (qend[i]-qstr[i])/double(nq-1); @@ -588,7 +622,8 @@ void Phonon::therm() char str[MAXLINE]; printf("\nPlease input the filename to output thermal properties [therm.dat]:"); - if (count_words(fgets(str,MAXLINE,stdin)) < 1) strcpy(str, "therm.dat"); + input->read_stdin(str); + if (count_words(str) < 1) strcpy(str, "therm.dat"); char *fname = strtok(str," \t\n\r\f"); FILE *fp = fopen(fname, "a"); fname = NULL; // header line @@ -630,7 +665,8 @@ void Phonon::therm() fprintf(fp,"%lg %lg %lg %lg %lg %lg\n", T, Uvib, Svib, Fvib, ZPE, Cvib); printf("Please input the desired temperature (K), enter to exit: "); - if (count_words(fgets(str,MAXLINE,stdin)) < 1) break; + input->read_stdin(str); + if (count_words(str) < 1) break; T = atof(strtok(str," \t\n\r\f")); } while (T > 0.); @@ -646,12 +682,14 @@ void Phonon::local_therm() { char str[MAXLINE]; printf("\nWould you like to compute the local thermal properties (y/n)[n]: "); - if (count_words(fgets(str,MAXLINE,stdin)) < 1) return; + input->read_stdin(str); + if (count_words(str) < 1) return; char *ptr = strtok(str," \t\n\r\f"); if (strcmp(ptr,"y") != 0 && strcmp(ptr, "Y") != 0 && strcmp(ptr, "yes") != 0) return; printf("Please input the filename to output vibrational thermal info [localtherm.dat]: "); - if (count_words(fgets(str,MAXLINE,stdin)) < 1) strcpy(str, "localtherm.dat"); + input->read_stdin(str); + if (count_words(str) < 1) strcpy(str, "localtherm.dat"); FILE *fp = fopen(strtok(str," \t\n\r\f"), "w"); fprintf(fp,"# atom Temp U_vib (eV) S_vib (kB) F_vib (eV) C_vib (kB) ZPE (eV)\n"); @@ -672,7 +710,8 @@ void Phonon::local_therm() while ( 1 ){ printf("\nPlease input the temperature at which to evaluate the local vibrational\n"); printf("thermal properties, non-positive number to exit [%g]: ", T); - if (count_words(fgets(str,MAXLINE,stdin)) > 0){ + input->read_stdin(str); + if (count_words(str) > 0){ T = atoi(strtok(str," \t\n\r\f")); if (T <= 0.) break; } @@ -765,7 +804,8 @@ void Phonon::QMesh() printf("\nThe q-mesh size from the read dynamical matrix is: %d x %d x %d\n", nx, ny, nz); printf("A denser mesh can be interpolated, but NOTE a too dense mesh can cause segmentation fault.\n"); printf("Please input your desired q-mesh size [%d %d %d]: ", nx, ny, nz); - if (count_words(fgets(str,MAXLINE,stdin)) >= 3){ + input->read_stdin(str); + if (count_words(str) >= 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")); @@ -780,7 +820,8 @@ void Phonon::QMesh() int method = 2; printf("Please select your method to generate the q-points:\n"); printf(" 1. uniform;\n 2. Monkhost-Pack mesh;\nYour choice [2]: "); - if (count_words(fgets(str,MAXLINE,stdin)) > 0) method = atoi(strtok(str," \t\n\r\f")); + input->read_stdin(str); + if (count_words(str) > 0) method = atoi(strtok(str," \t\n\r\f")); method = 2 - method%2; printf("Your selection: %d\n", method); #endif @@ -831,7 +872,7 @@ void Phonon::QMesh() for (int idim = 0; idim < sysdim; ++idim) atpos[i][idim] = dynmat->basis[i][idim]; // display the unit cell info read - printf("\n");for (int ii = 0; ii < 80; ++ii) printf("="); printf("\n"); + puts("\n================================================================================"); printf("The basis vectors of the unit cell:\n"); for (int idim = 0; idim < 3; ++idim) printf(" A%d = %lg %lg %lg\n", idim+1, latvec[0][idim], latvec[1][idim], latvec[2][idim]); @@ -845,14 +886,21 @@ void Phonon::QMesh() mesh[0] = nx; mesh[1] = ny; mesh[2] = nz; shift[0] = shift[1] = shift[2] = 0; int num_grid = mesh[0]*mesh[1]*mesh[2]; - int grid_point[num_grid][3], map[num_grid]; - double symprec = 1.e-4, pos[num_atom][3]; + int **grid_point; + memory->create(grid_point, num_grid, 3, "phonon:grid_point"); + int *map = new int[num_grid]; + double symprec = 1.0e-3; + double **pos; + memory->create(pos, num_atom, 3, "phonon:pos"); + if (dynmat->symprec > 0.) symprec = dynmat->symprec; for (int i = 0; i < num_atom; ++i) - for (int j = 0; j < 3; ++j) pos[i][j] = atpos[i][j]; + for (int j = 0; j < 3; ++j) + pos[i][j] = atpos[i][j]; // if spglib >= 1.0.3 is used - nq = spg_get_ir_reciprocal_mesh(grid_point, map, mesh, shift, is_time_reversal, latvec, pos, attyp, num_atom, symprec); + nq = spg_get_ir_reciprocal_mesh((int (*)[3])grid_point, map, mesh, shift, is_time_reversal, + latvec, (double (*)[3])pos, attyp, num_atom, symprec); memory->create(wt, nq, "QMesh:wt"); memory->create(qpts, nq,3,"QMesh:qpts"); @@ -873,11 +921,14 @@ void Phonon::QMesh() qpts[numq][2] = double(grid_point[i][2])/double(mesh[2]); numq++; } - wt[iq2idx[iq]] += 1.; + wt[iq2idx[iq]] += 1.0; } - delete []iq2idx; - - double wsum = 0.; + delete[] iq2idx; + delete[] map; + memory->destroy(grid_point); + memory->destroy(pos); + + double wsum = 0.0; for (int iq = 0; iq < nq; ++iq) wsum += wt[iq]; for (int iq = 0; iq < nq; ++iq) wt[iq] /= wsum; @@ -898,7 +949,8 @@ void Phonon::ldos_egv() char str[MAXLINE], *ptr; printf("\nThe # of atoms per cell is: %d, please input the atom IDs to compute\n", dynmat->nucell); printf("local PDOS, IDs begin with 0: "); - int nmax = count_words(fgets(str,MAXLINE,stdin)); + input->read_stdin(str); + int nmax = count_words(str); if (nmax < 1) return; memory->destroy(locals); @@ -920,7 +972,8 @@ void Phonon::ldos_egv() fmin = 0.; fmax = 10.; printf("Please input the freqency (nv, THz) range to compute PDOS [%g %g]: ", fmin, fmax); - if (count_words(fgets(str,MAXLINE,stdin)) >= 2) { + input->read_stdin(str); + if (count_words(str) >= 2) { fmin = atof(strtok(str," \t\n\r\f")); fmax = atof(strtok(NULL," \t\n\r\f")); } @@ -928,7 +981,8 @@ void Phonon::ldos_egv() ndos = 201; printf("Please input your desired # of points in PDOS [%d]: ", ndos); - if (count_words(fgets(str,MAXLINE,stdin)) > 0) ndos = atoi(strtok(str," \t\n\r\f")); + input->read_stdin(str); + if (count_words(str) > 0) ndos = atoi(strtok(str," \t\n\r\f")); if (ndos < 2) return; ndos += (ndos+1)%2; @@ -957,7 +1011,8 @@ void Phonon::ldos_egv() Timer *time = new Timer(); // memory and pointer for eigenvalues and eigenvectors - double egval[ndim], offset=fmin-0.5*df; + double offset=fmin-0.5*df; + double *egval = new double[ndim]; doublecomplex **egvec = dynmat->DM_q; printf("\nNow to compute the phonons and DOSs "); fflush(stdout); @@ -985,6 +1040,7 @@ void Phonon::ldos_egv() } } } + delete[] egval; egvec = NULL; printf("Done!\nNow to normalize the DOSs ..."); fflush(stdout); @@ -1008,10 +1064,7 @@ void Phonon::ldos_egv() * ---------------------------------------------------------------------------- */ void Phonon::ShowCell() { - printf("\n"); - for (int i = 0; i < 30; ++i) printf("="); - printf(" Unit Cell Info "); - for (int i = 0; i < 30; ++i) printf("="); printf("\n"); + puts("============================== Unit Cell Info =============================="); printf("Number of atoms in the unit cell: %d\n", dynmat->nucell); printf("Basis vectors of the unit cell:\n"); printf(" %15.8f %15.8f %15.8f\n", dynmat->basevec[0], dynmat->basevec[1], dynmat->basevec[2]); @@ -1024,8 +1077,7 @@ void Phonon::ShowCell() printf("Atomic type and fractional coordinates:\n"); for (int i = 0; i < dynmat->nucell; ++i) printf("%4d %12.8f %12.8f %12.8f\n", dynmat->attyp[i], dynmat->basis[i][0], dynmat->basis[i][1], dynmat->basis[i][2]); - for (int i = 0; i < 80; ++i) printf("="); - printf("\n"); + puts("================================================================================"); return; } @@ -1101,7 +1153,7 @@ int Phonon::count_words(const char *line) strcpy(copy,line); char *ptr; - if (ptr = strchr(copy,'#')) *ptr = '\0'; + if ((ptr = strchr(copy,'#'))) *ptr = '\0'; if (strtok(copy," \t\n\r\f") == NULL) { memory->destroy(copy); diff --git a/tools/phonon/phonon.h b/tools/phonon/phonon.h index 69a1fe5d50..46ce17f268 100644 --- a/tools/phonon/phonon.h +++ b/tools/phonon/phonon.h @@ -1,22 +1,16 @@ #ifndef PHONON_H #define PHONON_H -#include "stdio.h" -#include "stdlib.h" -#include -#include "dynmat.h" -#include "memory.h" - -using namespace std; - class Phonon{ public: - Phonon(DynMat *); + Phonon(class DynMat *); ~Phonon(); - DynMat *dynmat; + class DynMat *dynmat; private: + class UserInput *input; + int nq, ndim, sysdim; double **qpts, *wt; double **eigs; @@ -25,7 +19,7 @@ private: double *dos, fmin, fmax, df, rdf; double ***ldos; - Memory *memory; + class Memory *memory; void QMesh(); void ComputeAll(); diff --git a/tools/phonon/phonopy.cpp b/tools/phonon/phonopy.cpp index cfe32ab61a..a6b991efac 100644 --- a/tools/phonon/phonopy.cpp +++ b/tools/phonon/phonopy.cpp @@ -1,9 +1,20 @@ + #ifdef FFTW3 -#include + #include "phonopy.h" -#include "math.h" -#include "kpath.h" -#include "fftw3.h" + +#include "global.h" +#include "dynmat.h" +#include "input.h" +#include "memory.h" + +#include + +#include +#include +#include +#include +#include /* ---------------------------------------------------------------------------- * Class Phonopy is designed to interface with phonopy. @@ -14,13 +25,16 @@ Phonopy::Phonopy(DynMat *dynmat) memory = new Memory(); sysdim = dm->sysdim; fftdim = dm->fftdim; + input = dm->input; fftdim2 = fftdim * fftdim; nucell = dm->nucell; nx = ny = nz = 5; write(1); char str[MAXLINE]; - if (count_words(fgets(str,MAXLINE,stdin)) >= 3){ + if (input == NULL) input = new UserInput(0); + input->read_stdin(str); + if (count_words(str) >= 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")); @@ -36,7 +50,7 @@ Phonopy::Phonopy(DynMat *dynmat) memory->create(mass, nucell, "Phonopy:mass"); for (int i = 0; i < nucell; ++i){ - double m = 1./dm->M_inv_sqrt[i]; + double m = 1.0/dm->M_inv_sqrt[i]; mass[i] = m * m; } @@ -68,7 +82,7 @@ return; void Phonopy::write(int flag) { if (flag == 1){ // basic information - for (int ii = 0; ii < 80; ++ii) printf("="); printf("\n"); + puts("================================================================================"); 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); @@ -84,19 +98,18 @@ void Phonopy::write(int flag) } 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 + printf("phonopy band evaluation is written to band.conf.\n\n"); + printf("One should be able to obtain the phonon band structure after\n"); + printf(" 1) Correcting the `element names` in POSCAR.primitive and band.conf;\n"); + printf(" 2) Running `phonopy --readfc -c POSCAR.primitive -p band.conf`.\n\n"); + printf("Or the phonon density of states after\n"); + printf(" 1) Correcting the `element names` in POSCAR.primitive and mesh.conf;\n"); + printf(" 2) Running `phonopy --readfc -c POSCAR.primitive -p mesh.conf`.\n"); + puts("--------------------------------------------------------------------------------"); + printf("*** Remember to modify the `element names`. ***\n"); } else if (flag == 5){ - for (int ii = 0; ii < 80; ++ii) printf("="); printf("\n"); - + puts("================================================================================"); } return; } @@ -162,7 +175,9 @@ void Phonopy::phonopy() 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; + int *type_id = new int[nucell]; + int *num_type = new int[nucell]; + int ntype = 0; for (int i = 0; i < nucell; ++i) num_type[i] = 0; for (int i = 0; i < nucell; ++i){ int ip = ntype; @@ -221,7 +236,14 @@ void Phonopy::phonopy() // 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]); + for (int ip = 0; ip < ntype; ++ip){ + for (int i = 0; i < nucell; ++i){ + if (dm->attyp[i] == type_id[ip]){ + fprintf(fp, ", Elem-%d: %lg", type_id[ip], mass[i]); + break; + } + } + } fprintf(fp, "\n1.\n"); int ndim = 0; for (int idim = 0; idim < 3; ++idim){ @@ -240,57 +262,48 @@ void Phonopy::phonopy() } fclose(fp); -#ifdef UseSPG - // Get high symmetry k-path - QNodes *q = new QNodes(); - kPath *kp = new kPath(dm, q); - kp->kpath(); -#endif - + // mesh.conf + fp = fopen("mesh.conf", "w"); + fprintf(fp, "# From Fix-phonon"); + for (int ip = 0; ip < ntype; ++ip){ + for (int i = 0; i < nucell; ++i){ + if (dm->attyp[i] == type_id[ip]){ + fprintf(fp, ", Elem-%d: %lg", type_id[ip], mass[i]); + break; + } + } + } + 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\n", nx, ny, nz); + fprintf(fp, "MP = 31 31 31\nFORCE_CONSTANTS = READ\n"); + fprintf(fp, "#FC_SYMMETRY = .TRUE.\n#SYMMETRY_TOLERANCE = 0.01\n"); + fclose(fp); + + // 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]); + for (int ip = 0; ip < ntype; ++ip){ + for (int i = 0; i < nucell; ++i){ + if (dm->attyp[i] == type_id[ip]){ + fprintf(fp, ", Elem-%d: %lg", type_id[ip], mass[i]); + break; + } + } + } 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"); + fprintf(fp, "\nDIM = %d %d %d\nBAND = AUTO\n", nx, ny, nz); + fprintf(fp, "BAND_POINTS = 21\nFORCE_CONSTANTS = READ\nBAND_CONNECTION = .TRUE.\n"); + fprintf(fp, "#FC_SYMMETRY = .TRUE.\n#SYMMETRY_TOLERANCE = 0.01\n"); // output info write(4); -#ifdef UseSPG - kp->show_path(); - delete kp; - delete q; -#endif write(5); - -return; + delete[] type_id; + delete[] num_type; + return; } /*------------------------------------------------------------------------------ @@ -304,7 +317,7 @@ int Phonopy::count_words(const char *line) strcpy(copy,line); char *ptr; - if (ptr = strchr(copy,'#')) *ptr = '\0'; + if ((ptr = strchr(copy,'#'))) *ptr = '\0'; if (strtok(copy," \t\n\r\f") == NULL) { memory->destroy(copy); @@ -314,7 +327,7 @@ int Phonopy::count_words(const char *line) while (strtok(NULL," \t\n\r\f")) n++; memory->destroy(copy); -return n; + return n; } /*----------------------------------------------------------------------------*/ #endif diff --git a/tools/phonon/phonopy.h b/tools/phonon/phonopy.h index 2f3006e791..e1200b8242 100644 --- a/tools/phonon/phonopy.h +++ b/tools/phonon/phonopy.h @@ -4,22 +4,16 @@ #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" +#include "zheevd.h" class Phonopy { public: - Phonopy(DynMat *); + Phonopy(class DynMat *); ~Phonopy(); private: - Memory *memory; - char str[MAXLINE]; + class UserInput *input; + class Memory *memory; int npt, fftdim2; // local variables int nx, ny, nz, nucell; // local variables int sysdim, fftdim; // local variables @@ -27,7 +21,7 @@ private: doublecomplex **FC_all; - DynMat *dm; + class DynMat *dm; void write(int); void get_my_FC(); void phonopy(); diff --git a/tools/phonon/timer.cpp b/tools/phonon/timer.cpp index 953257340d..b94bed40b1 100644 --- a/tools/phonon/timer.cpp +++ b/tools/phonon/timer.cpp @@ -1,4 +1,7 @@ #include "timer.h" + +#include + /* ----------------------------------------------------------------------------- * Initialization of time * -------------------------------------------------------------------------- */ diff --git a/tools/phonon/timer.h b/tools/phonon/timer.h index cd04dea56d..993a33f9f3 100644 --- a/tools/phonon/timer.h +++ b/tools/phonon/timer.h @@ -1,9 +1,7 @@ #ifndef TIMER_H #define TIMER_H -#include "stdio.h" -#include "stdlib.h" -#include "time.h" +#include class Timer { public: diff --git a/tools/phonon/tricubic/CMakeLists.txt b/tools/phonon/tricubic/CMakeLists.txt new file mode 100644 index 0000000000..d70b44c2b2 --- /dev/null +++ b/tools/phonon/tricubic/CMakeLists.txt @@ -0,0 +1,18 @@ + +# Support Linux from Ubuntu 20.04LTS onward, CentOS 7.x (with EPEL), +# macOS, MSVC 2019 (=Version 16) +cmake_minimum_required(VERSION 3.16) + +# set up project +project(tricubic VERSION 1.1 DESCRIPTION "Tricubic library" LANGUAGES CXX) +set(CMAKE_POSITION_INDEPENDENT_CODE ON) + +# hacks for MSVC to prevent lots of pointless warnings about "unsafe" functions +if(MSVC) + add_compile_options(/wd4244) + add_compile_options(/wd4267) + add_compile_definitions(_CRT_SECURE_NO_WARNINGS) +endif() + +add_library(tricubic STATIC tricubic.cpp) +target_include_directories(tricubic PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) diff --git a/tools/phonon/tricubic/LICENSE b/tools/phonon/tricubic/LICENSE new file mode 100644 index 0000000000..d60c31a97a --- /dev/null +++ b/tools/phonon/tricubic/LICENSE @@ -0,0 +1,340 @@ + GNU GENERAL PUBLIC LICENSE + Version 2, June 1991 + + Copyright (C) 1989, 1991 Free Software Foundation, Inc. + 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + + Preamble + + The licenses for most software are designed to take away your +freedom to share and change it. By contrast, the GNU General Public +License is intended to guarantee your freedom to share and change free +software--to make sure the software is free for all its users. This +General Public License applies to most of the Free Software +Foundation's software and to any other program whose authors commit to +using it. (Some other Free Software Foundation software is covered by +the GNU Library General Public License instead.) You can apply it to +your programs, too. + + When we speak of free software, we are referring to freedom, not +price. Our General Public Licenses are designed to make sure that you +have the freedom to distribute copies of free software (and charge for +this service if you wish), that you receive source code or can get it +if you want it, that you can change the software or use pieces of it +in new free programs; and that you know you can do these things. + + To protect your rights, we need to make restrictions that forbid +anyone to deny you these rights or to ask you to surrender the rights. +These restrictions translate to certain responsibilities for you if you +distribute copies of the software, or if you modify it. + + For example, if you distribute copies of such a program, whether +gratis or for a fee, you must give the recipients all the rights that +you have. You must make sure that they, too, receive or can get the +source code. And you must show them these terms so they know their +rights. + + We protect your rights with two steps: (1) copyright the software, and +(2) offer you this license which gives you legal permission to copy, +distribute and/or modify the software. + + Also, for each author's protection and ours, we want to make certain +that everyone understands that there is no warranty for this free +software. If the software is modified by someone else and passed on, we +want its recipients to know that what they have is not the original, so +that any problems introduced by others will not reflect on the original +authors' reputations. + + Finally, any free program is threatened constantly by software +patents. We wish to avoid the danger that redistributors of a free +program will individually obtain patent licenses, in effect making the +program proprietary. To prevent this, we have made it clear that any +patent must be licensed for everyone's free use or not licensed at all. + + The precise terms and conditions for copying, distribution and +modification follow. + + GNU GENERAL PUBLIC LICENSE + TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION + + 0. This License applies to any program or other work which contains +a notice placed by the copyright holder saying it may be distributed +under the terms of this General Public License. The "Program", below, +refers to any such program or work, and a "work based on the Program" +means either the Program or any derivative work under copyright law: +that is to say, a work containing the Program or a portion of it, +either verbatim or with modifications and/or translated into another +language. (Hereinafter, translation is included without limitation in +the term "modification".) Each licensee is addressed as "you". + +Activities other than copying, distribution and modification are not +covered by this License; they are outside its scope. The act of +running the Program is not restricted, and the output from the Program +is covered only if its contents constitute a work based on the +Program (independent of having been made by running the Program). +Whether that is true depends on what the Program does. + + 1. You may copy and distribute verbatim copies of the Program's +source code as you receive it, in any medium, provided that you +conspicuously and appropriately publish on each copy an appropriate +copyright notice and disclaimer of warranty; keep intact all the +notices that refer to this License and to the absence of any warranty; +and give any other recipients of the Program a copy of this License +along with the Program. + +You may charge a fee for the physical act of transferring a copy, and +you may at your option offer warranty protection in exchange for a fee. + + 2. You may modify your copy or copies of the Program or any portion +of it, thus forming a work based on the Program, and copy and +distribute such modifications or work under the terms of Section 1 +above, provided that you also meet all of these conditions: + + a) You must cause the modified files to carry prominent notices + stating that you changed the files and the date of any change. + + b) You must cause any work that you distribute or publish, that in + whole or in part contains or is derived from the Program or any + part thereof, to be licensed as a whole at no charge to all third + parties under the terms of this License. + + c) If the modified program normally reads commands interactively + when run, you must cause it, when started running for such + interactive use in the most ordinary way, to print or display an + announcement including an appropriate copyright notice and a + notice that there is no warranty (or else, saying that you provide + a warranty) and that users may redistribute the program under + these conditions, and telling the user how to view a copy of this + License. (Exception: if the Program itself is interactive but + does not normally print such an announcement, your work based on + the Program is not required to print an announcement.) + +These requirements apply to the modified work as a whole. If +identifiable sections of that work are not derived from the Program, +and can be reasonably considered independent and separate works in +themselves, then this License, and its terms, do not apply to those +sections when you distribute them as separate works. But when you +distribute the same sections as part of a whole which is a work based +on the Program, the distribution of the whole must be on the terms of +this License, whose permissions for other licensees extend to the +entire whole, and thus to each and every part regardless of who wrote it. + +Thus, it is not the intent of this section to claim rights or contest +your rights to work written entirely by you; rather, the intent is to +exercise the right to control the distribution of derivative or +collective works based on the Program. + +In addition, mere aggregation of another work not based on the Program +with the Program (or with a work based on the Program) on a volume of +a storage or distribution medium does not bring the other work under +the scope of this License. + + 3. You may copy and distribute the Program (or a work based on it, +under Section 2) in object code or executable form under the terms of +Sections 1 and 2 above provided that you also do one of the following: + + a) Accompany it with the complete corresponding machine-readable + source code, which must be distributed under the terms of Sections + 1 and 2 above on a medium customarily used for software interchange; or, + + b) Accompany it with a written offer, valid for at least three + years, to give any third party, for a charge no more than your + cost of physically performing source distribution, a complete + machine-readable copy of the corresponding source code, to be + distributed under the terms of Sections 1 and 2 above on a medium + customarily used for software interchange; or, + + c) Accompany it with the information you received as to the offer + to distribute corresponding source code. (This alternative is + allowed only for noncommercial distribution and only if you + received the program in object code or executable form with such + an offer, in accord with Subsection b above.) + +The source code for a work means the preferred form of the work for +making modifications to it. For an executable work, complete source +code means all the source code for all modules it contains, plus any +associated interface definition files, plus the scripts used to +control compilation and installation of the executable. However, as a +special exception, the source code distributed need not include +anything that is normally distributed (in either source or binary +form) with the major components (compiler, kernel, and so on) of the +operating system on which the executable runs, unless that component +itself accompanies the executable. + +If distribution of executable or object code is made by offering +access to copy from a designated place, then offering equivalent +access to copy the source code from the same place counts as +distribution of the source code, even though third parties are not +compelled to copy the source along with the object code. + + 4. You may not copy, modify, sublicense, or distribute the Program +except as expressly provided under this License. Any attempt +otherwise to copy, modify, sublicense or distribute the Program is +void, and will automatically terminate your rights under this License. +However, parties who have received copies, or rights, from you under +this License will not have their licenses terminated so long as such +parties remain in full compliance. + + 5. You are not required to accept this License, since you have not +signed it. However, nothing else grants you permission to modify or +distribute the Program or its derivative works. These actions are +prohibited by law if you do not accept this License. Therefore, by +modifying or distributing the Program (or any work based on the +Program), you indicate your acceptance of this License to do so, and +all its terms and conditions for copying, distributing or modifying +the Program or works based on it. + + 6. Each time you redistribute the Program (or any work based on the +Program), the recipient automatically receives a license from the +original licensor to copy, distribute or modify the Program subject to +these terms and conditions. You may not impose any further +restrictions on the recipients' exercise of the rights granted herein. +You are not responsible for enforcing compliance by third parties to +this License. + + 7. If, as a consequence of a court judgment or allegation of patent +infringement or for any other reason (not limited to patent issues), +conditions are imposed on you (whether by court order, agreement or +otherwise) that contradict the conditions of this License, they do not +excuse you from the conditions of this License. If you cannot +distribute so as to satisfy simultaneously your obligations under this +License and any other pertinent obligations, then as a consequence you +may not distribute the Program at all. For example, if a patent +license would not permit royalty-free redistribution of the Program by +all those who receive copies directly or indirectly through you, then +the only way you could satisfy both it and this License would be to +refrain entirely from distribution of the Program. + +If any portion of this section is held invalid or unenforceable under +any particular circumstance, the balance of the section is intended to +apply and the section as a whole is intended to apply in other +circumstances. + +It is not the purpose of this section to induce you to infringe any +patents or other property right claims or to contest validity of any +such claims; this section has the sole purpose of protecting the +integrity of the free software distribution system, which is +implemented by public license practices. Many people have made +generous contributions to the wide range of software distributed +through that system in reliance on consistent application of that +system; it is up to the author/donor to decide if he or she is willing +to distribute software through any other system and a licensee cannot +impose that choice. + +This section is intended to make thoroughly clear what is believed to +be a consequence of the rest of this License. + + 8. If the distribution and/or use of the Program is restricted in +certain countries either by patents or by copyrighted interfaces, the +original copyright holder who places the Program under this License +may add an explicit geographical distribution limitation excluding +those countries, so that distribution is permitted only in or among +countries not thus excluded. In such case, this License incorporates +the limitation as if written in the body of this License. + + 9. The Free Software Foundation may publish revised and/or new versions +of the General Public License from time to time. Such new versions will +be similar in spirit to the present version, but may differ in detail to +address new problems or concerns. + +Each version is given a distinguishing version number. If the Program +specifies a version number of this License which applies to it and "any +later version", you have the option of following the terms and conditions +either of that version or of any later version published by the Free +Software Foundation. If the Program does not specify a version number of +this License, you may choose any version ever published by the Free Software +Foundation. + + 10. If you wish to incorporate parts of the Program into other free +programs whose distribution conditions are different, write to the author +to ask for permission. For software which is copyrighted by the Free +Software Foundation, write to the Free Software Foundation; we sometimes +make exceptions for this. Our decision will be guided by the two goals +of preserving the free status of all derivatives of our free software and +of promoting the sharing and reuse of software generally. + + NO WARRANTY + + 11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY +FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN +OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES +PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED +OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF +MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS +TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE +PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING, +REPAIR OR CORRECTION. + + 12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING +WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR +REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, +INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING +OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED +TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY +YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER +PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE +POSSIBILITY OF SUCH DAMAGES. + + END OF TERMS AND CONDITIONS + + How to Apply These Terms to Your New Programs + + If you develop a new program, and you want it to be of the greatest +possible use to the public, the best way to achieve this is to make it +free software which everyone can redistribute and change under these terms. + + To do so, attach the following notices to the program. It is safest +to attach them to the start of each source file to most effectively +convey the exclusion of warranty; and each file should have at least +the "copyright" line and a pointer to where the full notice is found. + + + Copyright (C) + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + + +Also add information on how to contact you by electronic and paper mail. + +If the program is interactive, make it output a short notice like this +when it starts in an interactive mode: + + Gnomovision version 69, Copyright (C) year name of author + Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'. + This is free software, and you are welcome to redistribute it + under certain conditions; type `show c' for details. + +The hypothetical commands `show w' and `show c' should show the appropriate +parts of the General Public License. Of course, the commands you use may +be called something other than `show w' and `show c'; they could even be +mouse-clicks or menu items--whatever suits your program. + +You should also get your employer (if you work as a programmer) or your +school, if any, to sign a "copyright disclaimer" for the program, if +necessary. Here is a sample; alter the names: + + Yoyodyne, Inc., hereby disclaims all copyright interest in the program + `Gnomovision' (which makes passes at compilers) written by James Hacker. + + , 1 April 1989 + Ty Coon, President of Vice + +This General Public License does not permit incorporating your program into +proprietary programs. If your program is a subroutine library, you may +consider it more useful to permit linking proprietary applications with the +library. If this is what you want to do, use the GNU Library General +Public License instead of this License. diff --git a/tools/phonon/tricubic/README.md b/tools/phonon/tricubic/README.md new file mode 100644 index 0000000000..0e8058339e --- /dev/null +++ b/tools/phonon/tricubic/README.md @@ -0,0 +1,28 @@ +# libtricubic + +This folder contains a slightly refactored version of version 1.0 of +the _libtricubic_ library, developed by François Lekien in 2004. + +The following paper explaining the method was published in 2005: + +> Lekien, F., & Marsden, J. (2005). _Tricubic interpolation in three +> dimensions._ International Journal for Numerical Methods in Engineering, +> 63(3), 455–471. doi:10.1002/nme.1296 + +Some additional notes and the full matrix can be found in the technical notes: + +> Lekien, F., Coulliette, C., & Marsden, J. (2004). _Tricubic Engine - Technical +> Notes and Full Matrix. + +The main article refers to the author's website +(http://gyre.cds.caltech.edu/pub/software/tricubic/) to download the code and +documentation. Unfortunately, this website no longer exists. Even the +archive.org snapshot is useless; [A single snapshot from December 3rd 2009 is +available](https://web.archive.org/web/20091203115835/http://gyre.cds.caltech.edu/pub/software/tricubic) +which seems like an FTP listing. No files are accessible. + +The source code was obtained from https://github.com/nbigaouette/libtricubic/ +Only the sources for the library were retained. No functional changes were +made, but some common programming conventions were applied, source files +merged, and build support for CMake added. + diff --git a/tools/phonon/tricubic/tricubic.cpp b/tools/phonon/tricubic/tricubic.cpp new file mode 100644 index 0000000000..974ed6aa52 --- /dev/null +++ b/tools/phonon/tricubic/tricubic.cpp @@ -0,0 +1,185 @@ + +#include "tricubic.h" + +#include + +static const 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 const char tricubic_version_stored[] = "1.1"; + +static int ijk2n(int i, int j, int k) { + return(i+4*j+16*k); +} + +static void point2xyz(int p, int *x, int *y, int *z) { + switch (p) { + case 0: *x=0; *y=0; *z=0; break; + case 1: *x=1; *y=0; *z=0; break; + case 2: *x=0; *y=1; *z=0; break; + case 3: *x=1; *y=1; *z=0; break; + case 4: *x=0; *y=0; *z=1; break; + case 5: *x=1; *y=0; *z=1; break; + case 6: *x=0; *y=1; *z=1; break; + case 7: *x=1; *y=1; *z=1; break; + default:*x=0; *y=0; *z=0; + } +} + +const char *tricubic_version(void) { + return(tricubic_version_stored); +} + +void tricubic_pointID2xyz(int id, int *x, int *y, int *z) { + point2xyz(id,x,y,z); +} + +void tricubic_pointID2xyz(int id, double *x, double *y, double *z) { + int x2,y2,z2; + point2xyz(id,&x2,&y2,&z2); + *x=(double)(x2); + *y=(double)(y2); + *z=(double)(z2); +} + +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]; + } + } +} + +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); +} + +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); +} + +double tricubic_eval(double a[64], double x, double y, double z, int derx, int dery, int derz) { + int i,j,k; + double ret=(double)(0.0); + double cont; + int w; + /* TRICUBIC_EVAL + The full version takes 3 extra integers args that allows to evaluate + any partial derivative of f at the point + derx=dery=derz=0 => f + derx=2 dery=derz=0 => d2f/dx2 + derx=dery=derz=1 =? d3f/dxdydz + NOTICE that (derx>3)||(dery>3)||(derz>3) => returns 0.0 + this computes \frac{\partial ^{derx+dery+derz} d}{\partial x ^{derx} \partial y ^{dery} \partial z ^{derz}} + */ + for (i=derx;i<4;i++) { + for (j=dery;j<4;j++) { + for (k=derz;k<4;k++) { + cont=a[ijk2n(i,j,k)]*pow(x,i-derx)*pow(y,j-dery)*pow(z,k-derz); + for (w=0;w