diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 4e7888ba..e89949ff 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -11,7 +11,7 @@ repos: - id: check-added-large-files - repo: https://github.com/astral-sh/ruff-pre-commit - rev: v0.5.0 + rev: v0.5.1 hooks: - id: ruff args: [ "--fix", "--show-fixes" ] diff --git a/CMakeLists.txt b/CMakeLists.txt index fc372955..a4d236e2 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,6 +2,7 @@ cmake_minimum_required(VERSION 3.20) option(BUILD_PHPHCALC_LIB "Option to build phph calculation module" OFF) option(BUILD_PHONONCALC_LIB "Option to build phonon calculation module" OFF) +option(BUILD_RECGRID_LIB "Option to build reciprocal space grid module" OFF) option(BUILD_GRIDSYS_LIB "Option to build gridsys module" OFF) option(PHONO3PY_WITH_Fortran "enable fortran interface" OFF) option(PHONO3PY_USE_OMP "Option to search OpenMP library" ON) @@ -125,27 +126,61 @@ if(BUILD_PHPHCALC_LIB endif() # ################################################################################### -# Build phphcalc module # +# Reciprocal space grid library # +# ################################################################################### +if(BUILD_RECGRID_LIB + OR BUILD_PHPHCALC_LIB + OR BUILD_NANOBIND_MODULE) + # Source code + set(SOURCES_PHPHCALC + ${PROJECT_SOURCE_DIR}/c/bzgrid.c ${PROJECT_SOURCE_DIR}/c/grgrid.c + ${PROJECT_SOURCE_DIR}/c/lagrid.c ${PROJECT_SOURCE_DIR}/c/snf3x3.c + ${PROJECT_SOURCE_DIR}/c/recgrid.c) + + if(BUILD_SHARED_LIBS) + # Shared library + add_library(recgrid_lib SHARED ${SOURCES_PHPHCALC}) + else() + # Static link library + add_library(recgrid_lib STATIC ${SOURCES_PHPHCALC}) + endif() + + if(NOT BUILD_NANOBIND_MODULE) + if(BUILD_SHARED_LIBS) + set_property(TARGET recgrid_lib PROPERTY VERSION ${SERIAL}) + set_property(TARGET recgrid_lib PROPERTY SOVERSION ${SOSERIAL}) + install(TARGETS recgrid_lib LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}) + else() + set_property(TARGET recgrid_lib PROPERTY VERSION ${SERIAL}) + set_property(TARGET recgrid_lib PROPERTY SOVERSION ${SOSERIAL}) + set_property(TARGET recgrid_lib PROPERTY OUTPUT_NAME recgrid_lib) + install(TARGETS recgrid_lib ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR}) + endif() + + # Header file + install(FILES ${PROJECT_SOURCE_DIR}/c/recgrid.h + DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}) + endif() +endif() + +# ################################################################################### +# Ph-ph calculation library # # ################################################################################### if(BUILD_PHPHCALC_LIB OR BUILD_NANOBIND_MODULE) # Source code set(SOURCES_PHPHCALC - ${PROJECT_SOURCE_DIR}/c/bzgrid.c ${PROJECT_SOURCE_DIR}/c/collision_matrix.c ${PROJECT_SOURCE_DIR}/c/fc3.c - ${PROJECT_SOURCE_DIR}/c/grgrid.c ${PROJECT_SOURCE_DIR}/c/imag_self_energy_with_g.c ${PROJECT_SOURCE_DIR}/c/interaction.c ${PROJECT_SOURCE_DIR}/c/isotope.c - ${PROJECT_SOURCE_DIR}/c/lagrid.c ${PROJECT_SOURCE_DIR}/c/lapack_wrapper.c ${PROJECT_SOURCE_DIR}/c/phono3py.c - ${PROJECT_SOURCE_DIR}/c/phonoc_utils.c + ${PROJECT_SOURCE_DIR}/c/funcs.c ${PROJECT_SOURCE_DIR}/c/pp_collision.c ${PROJECT_SOURCE_DIR}/c/real_self_energy.c ${PROJECT_SOURCE_DIR}/c/real_to_reciprocal.c ${PROJECT_SOURCE_DIR}/c/reciprocal_to_normal.c - ${PROJECT_SOURCE_DIR}/c/snf3x3.c ${PROJECT_SOURCE_DIR}/c/tetrahedron_method.c ${PROJECT_SOURCE_DIR}/c/triplet.c ${PROJECT_SOURCE_DIR}/c/triplet_grid.c @@ -156,10 +191,12 @@ if(BUILD_PHPHCALC_LIB OR BUILD_NANOBIND_MODULE) add_library(phphcalc_lib SHARED ${SOURCES_PHPHCALC}) if(OpenMP_FOUND) - target_link_libraries(phphcalc_lib PRIVATE BLAS::BLAS LAPACK::LAPACK - OpenMP::OpenMP_C) + target_link_libraries( + phphcalc_lib PRIVATE recgrid_lib BLAS::BLAS LAPACK::LAPACK + OpenMP::OpenMP_C) else() - target_link_libraries(phphcalc_lib PRIVATE BLAS::BLAS LAPACK::LAPACK) + target_link_libraries(phphcalc_lib PRIVATE recgrid_lib BLAS::BLAS + LAPACK::LAPACK) endif() target_include_directories(phphcalc_lib PRIVATE ${MY_INCLUDES}) @@ -188,10 +225,10 @@ if(BUILD_PHPHCALC_LIB OR BUILD_NANOBIND_MODULE) add_library(phphcalc_lib STATIC ${SOURCES_PHPHCALC}) if(OpenMP_FOUND) - target_link_libraries(phphcalc_lib BLAS::BLAS LAPACK::LAPACK + target_link_libraries(phphcalc_lib recgrid_lib BLAS::BLAS LAPACK::LAPACK OpenMP::OpenMP_C) else() - target_link_libraries(phphcalc_lib BLAS::BLAS LAPACK::LAPACK) + target_link_libraries(phphcalc_lib recgrid_lib BLAS::BLAS LAPACK::LAPACK) endif() target_include_directories(phphcalc_lib PRIVATE ${MY_INCLUDES}) @@ -236,7 +273,7 @@ if(BUILD_PHPHCALC_LIB OR BUILD_NANOBIND_MODULE) endif() # ################################################################################### -# phononcalc # +# Phonon calculation library # # ################################################################################### if(BUILD_PHONONCALC_LIB OR BUILD_NANOBIND_MODULE) # Source code @@ -311,7 +348,7 @@ if(BUILD_PHONONCALC_LIB OR BUILD_NANOBIND_MODULE) endif() # ################################################################################### -# grid # +# gridsys # # ################################################################################### if(BUILD_GRIDSYS_LIB) # Source code @@ -321,7 +358,7 @@ if(BUILD_GRIDSYS_LIB) ${PROJECT_SOURCE_DIR}/c/gridsys.c ${PROJECT_SOURCE_DIR}/c/lagrid.c ${PROJECT_SOURCE_DIR}/c/niggli.c - ${PROJECT_SOURCE_DIR}/c/phonoc_utils.c + ${PROJECT_SOURCE_DIR}/c/funcs.c ${PROJECT_SOURCE_DIR}/c/snf3x3.c ${PROJECT_SOURCE_DIR}/c/tetrahedron_method.c ${PROJECT_SOURCE_DIR}/c/triplet.c @@ -388,11 +425,15 @@ if(BUILD_NANOBIND_MODULE) ${PROJECT_SOURCE_DIR}/c/_phono3py.cpp) nanobind_add_module(_phononcalc STABLE_ABI ${PROJECT_SOURCE_DIR}/c/phononcalc.h ${PROJECT_SOURCE_DIR}/c/_phononcalc.cpp) + nanobind_add_module(_recgrid STABLE_ABI ${PROJECT_SOURCE_DIR}/c/recgrid.h + ${PROJECT_SOURCE_DIR}/c/_recgrid.cpp) target_link_libraries(_phono3py PRIVATE phphcalc_lib) target_link_libraries(_phononcalc PRIVATE phononcalc_lib) + target_link_libraries(_recgrid PRIVATE recgrid_lib) target_compile_definitions(_phono3py PRIVATE THM_EPSILON=1e-10) install(TARGETS _phono3py LIBRARY DESTINATION ${SKBUILD_PROJECT_NAME}) install(TARGETS _phononcalc LIBRARY DESTINATION ${SKBUILD_PROJECT_NAME}) + install(TARGETS _recgrid LIBRARY DESTINATION ${SKBUILD_PROJECT_NAME}) endif() diff --git a/c/_phono3py.cpp b/c/_phono3py.cpp index 5befc1cd..061631a5 100644 --- a/c/_phono3py.cpp +++ b/c/_phono3py.cpp @@ -1179,14 +1179,6 @@ NB_MODULE(_phono3py, m) { m.def("triplets_integration_weights", &py_get_triplets_integration_weights); m.def("triplets_integration_weights_with_sigma", &py_get_triplets_integration_weights_with_sigma); - m.def("grid_index_from_address", &py_get_grid_index_from_address); - m.def("ir_grid_map", &py_get_ir_grid_map); - m.def("gr_grid_addresses", &py_get_gr_grid_addresses); - m.def("reciprocal_rotations", &py_get_reciprocal_rotations); - m.def("transform_rotations", &py_transform_rotations); - m.def("snf3x3", &py_get_snf3x3); - m.def("bz_grid_addresses", &py_get_bz_grid_addresses); - m.def("rotate_bz_grid_index", &py_rotate_bz_grid_addresses); m.def("diagonalize_collision_matrix", &py_diagonalize_collision_matrix); m.def("pinv_from_eigensolution", &py_pinv_from_eigensolution); m.def("default_colmat_solver", &py_get_default_colmat_solver); diff --git a/c/_recgrid.cpp b/c/_recgrid.cpp new file mode 100644 index 00000000..bff4b7d1 --- /dev/null +++ b/c/_recgrid.cpp @@ -0,0 +1,179 @@ +#include +#include + +#include "recgrid.h" + +namespace nb = nanobind; + +long py_get_grid_index_from_address(nb::ndarray<> py_address, + nb::ndarray<> py_D_diag) { + long *address; + long *D_diag; + long gp; + + address = (long *)py_address.data(); + D_diag = (long *)py_D_diag.data(); + + gp = recgrid_get_grid_index_from_address(address, D_diag); + + return gp; +} + +long py_get_ir_grid_map(nb::ndarray<> py_grid_mapping_table, + nb::ndarray<> py_D_diag, nb::ndarray<> py_is_shift, + nb::ndarray<> py_rotations) { + long *D_diag; + long *is_shift; + long(*rot)[3][3]; + long num_rot; + + long *grid_mapping_table; + long num_ir; + + D_diag = (long *)py_D_diag.data(); + is_shift = (long *)py_is_shift.data(); + rot = (long(*)[3][3])py_rotations.data(); + num_rot = (long)py_rotations.shape(0); + grid_mapping_table = (long *)py_grid_mapping_table.data(); + + num_ir = recgrid_get_ir_grid_map(grid_mapping_table, rot, num_rot, D_diag, + is_shift); + return num_ir; +} + +void py_get_gr_grid_addresses(nb::ndarray<> py_gr_grid_addresses, + nb::ndarray<> py_D_diag) { + long(*gr_grid_addresses)[3]; + long *D_diag; + + gr_grid_addresses = (long(*)[3])py_gr_grid_addresses.data(); + D_diag = (long *)py_D_diag.data(); + + recgrid_get_all_grid_addresses(gr_grid_addresses, D_diag); +} + +long py_get_reciprocal_rotations(nb::ndarray<> py_rec_rotations, + nb::ndarray<> py_rotations, + long is_time_reversal) { + long(*rec_rotations)[3][3]; + long(*rotations)[3][3]; + long num_rot, num_rec_rot; + + rec_rotations = (long(*)[3][3])py_rec_rotations.data(); + rotations = (long(*)[3][3])py_rotations.data(); + num_rot = (long)py_rotations.shape(0); + + num_rec_rot = recgrid_get_reciprocal_point_group( + rec_rotations, rotations, num_rot, is_time_reversal, 1); + + return num_rec_rot; +} + +bool py_transform_rotations(nb::ndarray<> py_transformed_rotations, + nb::ndarray<> py_rotations, nb::ndarray<> py_D_diag, + nb::ndarray<> py_Q) { + long(*transformed_rotations)[3][3]; + long(*rotations)[3][3]; + long *D_diag; + long(*Q)[3]; + long num_rot, succeeded; + + transformed_rotations = (long(*)[3][3])py_transformed_rotations.data(); + rotations = (long(*)[3][3])py_rotations.data(); + D_diag = (long *)py_D_diag.data(); + Q = (long(*)[3])py_Q.data(); + num_rot = (long)py_transformed_rotations.shape(0); + + succeeded = recgrid_transform_rotations(transformed_rotations, rotations, + num_rot, D_diag, Q); + if (succeeded) { + return true; + } else { + return false; + } +} + +bool py_get_snf3x3(nb::ndarray<> py_D_diag, nb::ndarray<> py_P, + nb::ndarray<> py_Q, nb::ndarray<> py_A) { + long *D_diag; + long(*P)[3]; + long(*Q)[3]; + long(*A)[3]; + long succeeded; + + D_diag = (long *)py_D_diag.data(); + P = (long(*)[3])py_P.data(); + Q = (long(*)[3])py_Q.data(); + A = (long(*)[3])py_A.data(); + + succeeded = recgrid_get_snf3x3(D_diag, P, Q, A); + if (succeeded) { + return true; + } else { + return false; + } +} + +long py_get_bz_grid_addresses(nb::ndarray<> py_bz_grid_addresses, + nb::ndarray<> py_bz_map, nb::ndarray<> py_bzg2grg, + nb::ndarray<> py_D_diag, nb::ndarray<> py_Q, + nb::ndarray<> py_PS, + nb::ndarray<> py_reciprocal_lattice, long type) { + long(*bz_grid_addresses)[3]; + long *bz_map; + long *bzg2grg; + long *D_diag; + long(*Q)[3]; + long *PS; + double(*reciprocal_lattice)[3]; + long num_total_gp; + + bz_grid_addresses = (long(*)[3])py_bz_grid_addresses.data(); + bz_map = (long *)py_bz_map.data(); + bzg2grg = (long *)py_bzg2grg.data(); + D_diag = (long *)py_D_diag.data(); + Q = (long(*)[3])py_Q.data(); + PS = (long *)py_PS.data(); + reciprocal_lattice = (double(*)[3])py_reciprocal_lattice.data(); + + num_total_gp = + recgrid_get_bz_grid_addresses(bz_grid_addresses, bz_map, bzg2grg, + D_diag, Q, PS, reciprocal_lattice, type); + + return num_total_gp; +} + +long py_rotate_bz_grid_addresses(long bz_grid_index, nb::ndarray<> py_rotation, + nb::ndarray<> py_bz_grid_addresses, + nb::ndarray<> py_bz_map, + nb::ndarray<> py_D_diag, nb::ndarray<> py_PS, + long type) { + long(*bz_grid_addresses)[3]; + long(*rotation)[3]; + long *bz_map; + long *D_diag; + long *PS; + long ret_bz_gp; + + bz_grid_addresses = (long(*)[3])py_bz_grid_addresses.data(); + rotation = (long(*)[3])py_rotation.data(); + bz_map = (long *)py_bz_map.data(); + D_diag = (long *)py_D_diag.data(); + PS = (long *)py_PS.data(); + + ret_bz_gp = recgrid_rotate_bz_grid_index( + bz_grid_index, rotation, bz_grid_addresses, bz_map, D_diag, PS, type); + + return ret_bz_gp; +} + +NB_MODULE(_recgrid, m) { + m.def("grid_index_from_address", &py_get_grid_index_from_address); + m.def("ir_grid_map", &py_get_ir_grid_map); + m.def("gr_grid_addresses", &py_get_gr_grid_addresses); + m.def("reciprocal_rotations", &py_get_reciprocal_rotations); + m.def("transform_rotations", &py_transform_rotations); + m.def("snf3x3", &py_get_snf3x3); + m.def("bz_grid_addresses", &py_get_bz_grid_addresses); + m.def("rotate_bz_grid_index", &py_rotate_bz_grid_addresses); +} diff --git a/c/bzgrid.c b/c/bzgrid.c index fb938e1a..d70c64c2 100644 --- a/c/bzgrid.c +++ b/c/bzgrid.c @@ -40,9 +40,9 @@ #include "grgrid.h" #include "lagrid.h" +#include "recgrid.h" #define BZG_NUM_BZ_SEARCH_SPACE 125 -#define GRID_TOLERANCE_FACTOR 0.01 static long bz_search_space[BZG_NUM_BZ_SEARCH_SPACE][3] = { {0, 0, 0}, {0, 0, 1}, {0, 0, 2}, {0, 0, -2}, {0, 0, -1}, {0, 1, 0}, {0, 1, 1}, {0, 1, 2}, {0, 1, -2}, {0, 1, -1}, @@ -70,21 +70,26 @@ static long bz_search_space[BZG_NUM_BZ_SEARCH_SPACE][3] = { {-1, -2, 0}, {-1, -2, 1}, {-1, -2, 2}, {-1, -2, -2}, {-1, -2, -1}, {-1, -1, 0}, {-1, -1, 1}, {-1, -1, 2}, {-1, -1, -2}, {-1, -1, -1}}; -static void get_bz_grid_addresses_type1(BZGrid *bzgrid, const long Qinv[3][3]); -static void get_bz_grid_addresses_type2(BZGrid *bzgrid, const long Qinv[3][3]); +static void get_bz_grid_addresses_type1(RecgridBZGrid *bzgrid, + const long Qinv[3][3]); +static void get_bz_grid_addresses_type2(RecgridBZGrid *bzgrid, + const long Qinv[3][3]); static void set_bz_address(long address[3], const long bz_index, const long grid_address[3], const long D_diag[3], const long nint[3], const long Qinv[3][3]); static double get_bz_distances(long nint[3], double distances[], - const BZGrid *bzgrid, const long grid_address[3], + const RecgridBZGrid *bzgrid, + const long grid_address[3], const double tolerance); static void multiply_matrix_vector_d3(double v[3], const double a[3][3], const double b[3]); +static void multiply_matrix_vector_ld3(double v[3], const long a[3][3], + const double b[3]); static long get_inverse_unimodular_matrix_l3(long m[3][3], const long a[3][3]); static double norm_squared_d3(const double a[3]); long bzg_rotate_grid_index(const long bz_grid_index, const long rotation[3][3], - const ConstBZGrid *bzgrid) { + const RecgridConstBZGrid *bzgrid) { long i, gp, num_bzgp, num_grgp; long dadrs[3], dadrs_rot[3], adrs_rot[3]; @@ -124,7 +129,7 @@ long bzg_rotate_grid_index(const long bz_grid_index, const long rotation[3][3], return bzgrid->gp_map[gp]; } -long bzg_get_bz_grid_addresses(BZGrid *bzgrid) { +long bzg_get_bz_grid_addresses(RecgridBZGrid *bzgrid) { long det; long Qinv[3][3]; @@ -142,91 +147,15 @@ long bzg_get_bz_grid_addresses(BZGrid *bzgrid) { return 1; } -/* Note: Tolerance in squared distance. */ -double bzg_get_tolerance_for_BZ_reduction(const BZGrid *bzgrid) { - long i, j; - double tolerance; - double length[3]; - double reclatQ[3][3]; - - for (i = 0; i < 3; i++) { - for (j = 0; j < 3; j++) { - reclatQ[i][j] = bzgrid->reclat[i][0] * bzgrid->Q[0][j] + - bzgrid->reclat[i][1] * bzgrid->Q[1][j] + - bzgrid->reclat[i][2] * bzgrid->Q[2][j]; - } - } - - for (i = 0; i < 3; i++) { - length[i] = 0; - for (j = 0; j < 3; j++) { - length[i] += reclatQ[j][i] * reclatQ[j][i]; - } - length[i] /= bzgrid->D_diag[i] * bzgrid->D_diag[i]; - } - tolerance = length[0]; - for (i = 1; i < 3; i++) { - if (tolerance < length[i]) { - tolerance = length[i]; - } - } - tolerance *= GRID_TOLERANCE_FACTOR; - - return tolerance; -} - -RotMats *bzg_alloc_RotMats(const long size) { - RotMats *rotmats; - - rotmats = NULL; - - if ((rotmats = (RotMats *)malloc(sizeof(RotMats))) == NULL) { - warning_print("Memory could not be allocated."); - return NULL; - } - - rotmats->size = size; - if (size > 0) { - if ((rotmats->mat = (long(*)[3][3])malloc(sizeof(long[3][3]) * size)) == - NULL) { - warning_print("Memory could not be allocated "); - warning_print("(RotMats, line %d, %s).\n", __LINE__, __FILE__); - free(rotmats); - rotmats = NULL; - return NULL; - } - } - return rotmats; -} - -void bzg_free_RotMats(RotMats *rotmats) { - if (rotmats->size > 0) { - free(rotmats->mat); - rotmats->mat = NULL; - } - free(rotmats); -} - -void bzg_multiply_matrix_vector_ld3(double v[3], const long a[3][3], - const double b[3]) { - long i; - double c[3]; - for (i = 0; i < 3; i++) { - c[i] = a[i][0] * b[0] + a[i][1] * b[1] + a[i][2] * b[2]; - } - for (i = 0; i < 3; i++) { - v[i] = c[i]; - } -} - -static void get_bz_grid_addresses_type1(BZGrid *bzgrid, const long Qinv[3][3]) { +static void get_bz_grid_addresses_type1(RecgridBZGrid *bzgrid, + const long Qinv[3][3]) { double tolerance, min_distance; double distances[BZG_NUM_BZ_SEARCH_SPACE]; long bzmesh[3], bz_address_double[3], nint[3], gr_adrs[3]; long i, j, k, boundary_num_gp, total_num_gp, bzgp, gp, num_bzmesh; long count, id_shift; - tolerance = bzg_get_tolerance_for_BZ_reduction(bzgrid); + tolerance = recgrid_get_tolerance_for_BZ_reduction(bzgrid); for (j = 0; j < 3; j++) { bzmesh[j] = bzgrid->D_diag[j] * 2; } @@ -278,13 +207,14 @@ static void get_bz_grid_addresses_type1(BZGrid *bzgrid, const long Qinv[3][3]) { bzgrid->size = boundary_num_gp + total_num_gp; } -static void get_bz_grid_addresses_type2(BZGrid *bzgrid, const long Qinv[3][3]) { +static void get_bz_grid_addresses_type2(RecgridBZGrid *bzgrid, + const long Qinv[3][3]) { double tolerance, min_distance; double distances[BZG_NUM_BZ_SEARCH_SPACE]; long nint[3], gr_adrs[3]; long i, j, num_gp; - tolerance = bzg_get_tolerance_for_BZ_reduction(bzgrid); + tolerance = recgrid_get_tolerance_for_BZ_reduction(bzgrid); num_gp = 0; /* The first element of gp_map is always 0. */ bzgrid->gp_map[0] = 0; @@ -324,7 +254,8 @@ static void set_bz_address(long address[3], const long bz_index, } static double get_bz_distances(long nint[3], double distances[], - const BZGrid *bzgrid, const long grid_address[3], + const RecgridBZGrid *bzgrid, + const long grid_address[3], const double tolerance) { long i, j; long dadrs[3]; @@ -336,7 +267,7 @@ static double get_bz_distances(long nint[3], double distances[], for (i = 0; i < 3; i++) { q_red[i] = dadrs[i] / (2.0 * bzgrid->D_diag[i]); } - bzg_multiply_matrix_vector_ld3(q_red, bzgrid->Q, q_red); + multiply_matrix_vector_ld3(q_red, bzgrid->Q, q_red); for (i = 0; i < 3; i++) { nint[i] = lagmat_Nint(q_red[i]); q_red[i] -= nint[i]; @@ -376,6 +307,18 @@ static void multiply_matrix_vector_d3(double v[3], const double a[3][3], } } +static void multiply_matrix_vector_ld3(double v[3], const long a[3][3], + const double b[3]) { + long i; + double c[3]; + for (i = 0; i < 3; i++) { + c[i] = a[i][0] * b[0] + a[i][1] * b[1] + a[i][2] * b[2]; + } + for (i = 0; i < 3; i++) { + v[i] = c[i]; + } +} + static long get_inverse_unimodular_matrix_l3(long m[3][3], const long a[3][3]) { long det; long c[3][3]; diff --git a/c/bzgrid.h b/c/bzgrid.h index b1e700de..c0fe2db1 100644 --- a/c/bzgrid.h +++ b/c/bzgrid.h @@ -35,69 +35,10 @@ #ifndef __bzgrid_H__ #define __bzgrid_H__ -typedef struct { - long size; - long (*mat)[3][3]; -} RotMats; - -/* Data structure of Brillouin zone grid - * - * size : long - * Number of grid points in Brillouin zone including its surface. - * D_diag : long array - * Diagonal part of matrix D of SNF. - * shape=(3, ) - * PS : long array - * Matrix P of SNF multiplied by shift. - * shape=(3, ) - * gp_map : long array - * Type1 : Twice enlarged grid space along basis vectors. - * Grid index is recovered in the same way as regular grid. - * shape=(prod(mesh * 2), ) - * Type2 : In the last index, multiplicity and array index of - * each address of the grid point are stored. Here, - * multiplicity means the number of translationally - * equivalent grid points in BZ. - * shape=(prod(mesh), 2) -> flattened. - * addresses : long array - * Grid point addresses. - * shape=(size, 3) - * reclat : double array - * Reciprocal basis vectors given as column vectors. - * shape=(3, 3) - * type : long - * 1 or 2. */ -typedef struct { - long size; - long D_diag[3]; - long Q[3][3]; - long PS[3]; - long *gp_map; - long *bzg2grg; - long (*addresses)[3]; - double reclat[3][3]; - long type; -} BZGrid; - -typedef struct { - long size; - long D_diag[3]; - long Q[3][3]; - long PS[3]; - const long *gp_map; - const long *bzg2grg; - const long (*addresses)[3]; - double reclat[3][3]; - long type; -} ConstBZGrid; +#include "recgrid.h" long bzg_rotate_grid_index(const long grid_index, const long rotation[3][3], - const ConstBZGrid *bzgrid); -long bzg_get_bz_grid_addresses(BZGrid *bzgrid); -double bzg_get_tolerance_for_BZ_reduction(const BZGrid *bzgrid); -RotMats *bzg_alloc_RotMats(const long size); -void bzg_free_RotMats(RotMats *rotmats); -void bzg_multiply_matrix_vector_ld3(double v[3], const long a[3][3], - const double b[3]); + const RecgridConstBZGrid *bzgrid); +long bzg_get_bz_grid_addresses(RecgridBZGrid *bzgrid); #endif diff --git a/c/collision_matrix.c b/c/collision_matrix.c index 8ae7a3c0..75d17849 100644 --- a/c/collision_matrix.c +++ b/c/collision_matrix.c @@ -38,8 +38,8 @@ #include #include +#include "funcs.h" #include "phonoc_array.h" -#include "phonoc_utils.h" static void get_collision_matrix( double *collision_matrix, const double *fc3_normal_squared, @@ -255,7 +255,7 @@ static long get_inv_sinh(double *inv_sinh, const long gp, for (i = 0; i < num_band; i++) { f = frequencies[gp2 * num_band + i]; if (f > cutoff_frequency) { - inv_sinh[i] = phonoc_inv_sinh_occupation(f, temperature); + inv_sinh[i] = funcs_inv_sinh_occupation(f, temperature); } else { inv_sinh[i] = 0; } diff --git a/c/dynmat.c b/c/dynmat.c index e179a1a6..6be2ba1e 100644 --- a/c/dynmat.c +++ b/c/dynmat.c @@ -89,9 +89,9 @@ static void get_q_cart(double q_cart[3], const double q[3], const double reciprocal_lattice[3][3]); /// @brief Calculate dynamical matrices with openmp over q-points -/// use_Wang_NAC gives NAC by Wang et al. -/// dd_0 and NULL gives NAC by Gonze and Lee -/// Otherwise no-NAC. +/// use_Wang_NAC: Wang et al. +/// !use_Wang_NAC and dd_0 is NULL: no-NAC +/// !use_Wang_NAC and dd_0 is not NULL: NAC by Gonze and Lee. /// @param reciprocal_lattice in column vectors /// @param q_direction in Crystallographic coordinates. long dym_dynamical_matrices_with_dd_openmp_over_qpoints( @@ -126,28 +126,25 @@ long dym_dynamical_matrices_with_dd_openmp_over_qpoints( #pragma omp parallel for private(charge_sum, q_cart, q_norm) #endif for (i = 0; i < n_qpoints; i++) { - charge_sum = (double(*)[3][3])malloc(sizeof(double[3][3]) * - num_patom * num_patom); get_q_cart(q_cart, qpoints[i], reciprocal_lattice); q_norm = sqrt(q_cart[0] * q_cart[0] + q_cart[1] * q_cart[1] + q_cart[2] * q_cart[2]); - if (q_norm < q_zero_tolerance) { - if (q_direction) { - dym_get_charge_sum( - charge_sum, num_patom, - nac_factor / n / - get_dielectric_part(q_direction, dielectric), - q_dir_cart, born); - } else { - free(charge_sum); - charge_sum = NULL; - } - } else { + if (q_norm < q_zero_tolerance && q_direction) { + charge_sum = (double(*)[3][3])malloc(sizeof(double[3][3]) * + num_patom * num_patom); + dym_get_charge_sum( + charge_sum, num_patom, + nac_factor / n / + get_dielectric_part(q_direction, dielectric), + q_dir_cart, born); + } else if (!(q_norm < q_zero_tolerance)) { + charge_sum = (double(*)[3][3])malloc(sizeof(double[3][3]) * + num_patom * num_patom); dym_get_charge_sum( charge_sum, num_patom, nac_factor / n / get_dielectric_part(q_cart, dielectric), q_cart, born); - } + } // (q_norm < q_zero_tolerance && !q_direction) -> no NAC dym_get_dynamical_matrix_at_q(dynamical_matrices + adrs_shift * i, num_patom, num_satom, fc, qpoints[i], svecs, multi, masses, s2p_map, @@ -166,7 +163,7 @@ long dym_dynamical_matrices_with_dd_openmp_over_qpoints( num_patom, num_satom, fc, qpoints[i], svecs, multi, masses, s2p_map, p2s_map, charge_sum, 0); - if (dd_q0) { + if (dd_q0) { // NAC by Gonze and Lee if dd_in is not NULL add_dynmat_dd_at_q(dynamical_matrices + adrs_shift * i, qpoints[i], fc, positions, num_patom, masses, born, dielectric, reciprocal_lattice, @@ -196,8 +193,6 @@ static void add_dynmat_dd_at_q( double q_cart[3]; double mm; - q_dir_cart = NULL; - dd = (double(*)[2])malloc(sizeof(double[2]) * num_patom * num_patom * 9); get_q_cart(q_cart, q, reciprocal_lattice); dym_get_recip_dipole_dipole(dd, dd_q0, G_list, num_G_points, num_patom, @@ -221,29 +216,7 @@ static void add_dynmat_dd_at_q( dd = NULL; } -long dym_get_dynamical_matrices_openmp_over_qpoints( - double (*dynamical_matrices)[2], // [q-points, num_band, num_band, - // (real, imag)] - const long num_patom, const long num_satom, const double *fc, - const double (*qpoints)[3], const long n_qpoints, const double (*svecs)[3], - const long (*multi)[2], const double *mass, const long *s2p_map, - const long *p2s_map, const double (*charge_sum)[3][3]) { - long i, adrs_shift; - - adrs_shift = num_patom * num_patom * 9; - -#ifdef _OPENMP -#pragma omp parallel for -#endif - for (i = 0; i < n_qpoints; i++) { - dym_get_dynamical_matrix_at_q( - dynamical_matrices + adrs_shift * i, num_patom, num_satom, fc, - qpoints[i], svecs, multi, mass, s2p_map, p2s_map, charge_sum, 0); - } - - return 0; -} - +/// @brief charge_sum is NULL if G-L NAC or no-NAC. long dym_get_dynamical_matrix_at_q(double (*dynamical_matrix)[2], const long num_patom, const long num_satom, const double *fc, const double q[3], @@ -491,6 +464,7 @@ void dym_transform_dynmat_to_fc(double *fc, const double (*dm)[2], } } +/// @brief charge_sum is NULL if G-L NAC or no-NAC. static void get_dynmat_ij(double (*dynamical_matrix)[2], const long num_patom, const long num_satom, const double *fc, const double q[3], const double (*svecs)[3], @@ -601,7 +575,7 @@ static void get_dd(double (*dd_part)[2], /* [natom, 3, natom, 3, (real,imag)] */ /* q_direction is NULL for summation over G */ #ifdef _OPENMP #pragma omp parallel for private(i, j, q_K, norm, \ - dielectric_part) if (use_openmp) + dielectric_part) if (use_openmp) #endif for (g = 0; g < num_G; g++) { norm = 0; diff --git a/c/dynmat.h b/c/dynmat.h index 0e982fc3..491f768a 100644 --- a/c/dynmat.h +++ b/c/dynmat.h @@ -45,11 +45,6 @@ long dym_dynamical_matrices_with_dd_openmp_over_qpoints( const double *q_direction, const double nac_factor, const double (*dd_q0)[2], const double (*G_list)[3], const long num_G_points, const double lambda, const long use_Wang_NAC); -long dym_get_dynamical_matrices_openmp_over_qpoints( - double (*dynamical_matrices)[2], const long num_patom, const long num_satom, - const double *fc, const double (*qpoints)[3], const long n_qpoints, - const double (*svecs)[3], const long (*multi)[2], const double *mass, - const long *s2p_map, const long *p2s_map, const double (*charge_sum)[3][3]); long dym_get_dynamical_matrix_at_q(double (*dynamical_matrix)[2], const long num_patom, const long num_satom, const double *fc, const double q[3], diff --git a/c/phonoc_utils.c b/c/funcs.c similarity index 90% rename from c/phonoc_utils.c rename to c/funcs.c index f86236bd..5ffab41e 100644 --- a/c/phonoc_utils.c +++ b/c/funcs.c @@ -32,7 +32,7 @@ /* ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE */ /* POSSIBILITY OF SUCH DAMAGE. */ -#include "phonoc_utils.h" +#include "funcs.h" #include @@ -41,14 +41,14 @@ #define THZTOEVPARKB 47.992398658977166 #define INVSQRT2PI 0.3989422804014327 -double phonoc_bose_einstein(const double x, const double t) { +double funcs_bose_einstein(const double x, const double t) { return 1.0 / (exp(THZTOEVPARKB * x / t) - 1); } -double phonoc_gaussian(const double x, const double sigma) { +double funcs_gaussian(const double x, const double sigma) { return INVSQRT2PI / sigma * exp(-x * x / 2 / sigma / sigma); } -double phonoc_inv_sinh_occupation(const double x, const double t) { +double funcs_inv_sinh_occupation(const double x, const double t) { return 1.0 / sinh(x * THZTOEVPARKB / 2 / t); } diff --git a/c/phonoc_utils.h b/c/funcs.h similarity index 87% rename from c/phonoc_utils.h rename to c/funcs.h index 0334e424..952e4d1c 100644 --- a/c/phonoc_utils.h +++ b/c/funcs.h @@ -32,11 +32,11 @@ /* ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE */ /* POSSIBILITY OF SUCH DAMAGE. */ -#ifndef __phonoc_utils_H__ -#define __phonoc_utils_H__ +#ifndef __funcs_H__ +#define __funcs_H__ -double phonoc_bose_einstein(const double x, const double t); -double phonoc_gaussian(const double x, const double sigma); -double phonoc_inv_sinh_occupation(const double x, const double t); +double funcs_bose_einstein(const double x, const double t); +double funcs_gaussian(const double x, const double sigma); +double funcs_inv_sinh_occupation(const double x, const double t); #endif diff --git a/c/gridsys.c b/c/gridsys.c index 4c8aa645..4a7e9458 100644 --- a/c/gridsys.c +++ b/c/gridsys.c @@ -41,6 +41,7 @@ #include "grgrid.h" #include "lagrid.h" #include "niggli.h" +#include "recgrid.h" #include "tetrahedron_method.h" #include "triplet.h" #include "triplet_iw.h" @@ -111,13 +112,13 @@ long gridsys_get_bz_grid_addresses(long (*bz_grid_addresses)[3], long *bz_map, const long Q[3][3], const long PS[3], const double rec_lattice[3][3], const long bz_grid_type) { - BZGrid *bzgrid; + RecgridBZGrid *bzgrid; long i, j, size; long inv_Mpr_int[3][3]; double inv_Lr[3][3], inv_Mpr[3][3]; double niggli_lattice[9]; - if ((bzgrid = (BZGrid *)malloc(sizeof(BZGrid))) == NULL) { + if ((bzgrid = (RecgridBZGrid *)malloc(sizeof(RecgridBZGrid))) == NULL) { warning_print("Memory could not be allocated."); return 0; } @@ -171,10 +172,11 @@ long gridsys_rotate_bz_grid_index(const long bz_grid_index, const long (*bz_grid_addresses)[3], const long *bz_map, const long D_diag[3], const long PS[3], const long bz_grid_type) { - ConstBZGrid *bzgrid; + RecgridConstBZGrid *bzgrid; long i, rot_bz_gp; - if ((bzgrid = (ConstBZGrid *)malloc(sizeof(ConstBZGrid))) == NULL) { + if ((bzgrid = (RecgridConstBZGrid *)malloc(sizeof(RecgridConstBZGrid))) == + NULL) { warning_print("Memory could not be allocated."); return 0; } @@ -212,10 +214,11 @@ long gridsys_get_bz_triplets_at_q(long (*ir_triplets)[3], const long num_map_triplets, const long D_diag[3], const long Q[3][3], const long bz_grid_type) { - ConstBZGrid *bzgrid; + RecgridConstBZGrid *bzgrid; long i, j, num_ir; - if ((bzgrid = (ConstBZGrid *)malloc(sizeof(ConstBZGrid))) == NULL) { + if ((bzgrid = (RecgridConstBZGrid *)malloc(sizeof(RecgridConstBZGrid))) == + NULL) { warning_print("Memory could not be allocated."); return 0; } @@ -267,10 +270,11 @@ long gridsys_get_integration_weight( const long bz_grid_type, const double *frequencies1, const long num_band1, const double *frequencies2, const long num_band2, const long tp_type, const long openmp_per_triplets) { - ConstBZGrid *bzgrid; + RecgridConstBZGrid *bzgrid; long i; - if ((bzgrid = (ConstBZGrid *)malloc(sizeof(ConstBZGrid))) == NULL) { + if ((bzgrid = (RecgridConstBZGrid *)malloc(sizeof(RecgridConstBZGrid))) == + NULL) { warning_print("Memory could not be allocated."); return 0; } diff --git a/c/gridsys.h b/c/gridsys.h index a7aa8431..622b79b7 100644 --- a/c/gridsys.h +++ b/c/gridsys.h @@ -177,7 +177,7 @@ void gridsys_get_ir_grid_map(long *ir_grid_map, const long (*rotations)[3][3], /** * @brief Find shortest grid points from Gamma considering periodicity of - * reciprocal lattice. See the details in docstring of BZGrid + * reciprocal lattice. See the details in docstring of RecgridBZGrid * * @param bz_grid_addresses Grid point addresses of shortest grid points * @param bz_map List of accumulated numbers of BZ grid points from the diff --git a/c/imag_self_energy_with_g.c b/c/imag_self_energy_with_g.c index d9aac798..b4f1e019 100644 --- a/c/imag_self_energy_with_g.c +++ b/c/imag_self_energy_with_g.c @@ -38,9 +38,9 @@ #include #include +#include "funcs.h" #include "lagrid.h" #include "phonoc_array.h" -#include "phonoc_utils.h" #include "triplet.h" static long set_g_pos_frequency_point(long (*g_pos)[4], const long num_band0, @@ -455,12 +455,12 @@ static void set_occupations(double *n1, double *n2, const long num_band, f1 = frequencies[triplet[1] * num_band + j]; f2 = frequencies[triplet[2] * num_band + j]; if (f1 > cutoff_frequency) { - n1[j] = phonoc_bose_einstein(f1, temperature); + n1[j] = funcs_bose_einstein(f1, temperature); } else { n1[j] = -1; } if (f2 > cutoff_frequency) { - n2[j] = phonoc_bose_einstein(f2, temperature); + n2[j] = funcs_bose_einstein(f2, temperature); } else { n2[j] = -1; } diff --git a/c/interaction.c b/c/interaction.c index cc7ba929..7ec03798 100644 --- a/c/interaction.c +++ b/c/interaction.c @@ -42,6 +42,7 @@ #include "lapack_wrapper.h" #include "phonoc_array.h" #include "real_to_reciprocal.h" +#include "recgrid.h" #include "reciprocal_to_normal.h" static const long index_exchange[6][3] = {{0, 1, 2}, {2, 0, 1}, {1, 2, 0}, @@ -73,9 +74,10 @@ static void real_to_normal_sym_q( void itr_get_interaction( Darray *fc3_normal_squared, const char *g_zero, const Darray *frequencies, const lapack_complex_double *eigenvectors, const long (*triplets)[3], - const long num_triplets, const ConstBZGrid *bzgrid, const double *fc3, - const long is_compact_fc3, const AtomTriplets *atom_triplets, - const double *masses, const long *band_indices, const long symmetrize_fc3_q, + const long num_triplets, const RecgridConstBZGrid *bzgrid, + const double *fc3, const long is_compact_fc3, + const AtomTriplets *atom_triplets, const double *masses, + const long *band_indices, const long symmetrize_fc3_q, const double cutoff_frequency, const long openmp_per_triplets) { long(*g_pos)[4]; long i; @@ -112,9 +114,9 @@ void itr_get_interaction_at_triplet( double *fc3_normal_squared, const long num_band0, const long num_band, const long (*g_pos)[4], const long num_g_pos, const double *frequencies, const lapack_complex_double *eigenvectors, const long triplet[3], - const ConstBZGrid *bzgrid, const double *fc3, const long is_compact_fc3, - const AtomTriplets *atom_triplets, const double *masses, - const long *band_indices, const long symmetrize_fc3_q, + const RecgridConstBZGrid *bzgrid, const double *fc3, + const long is_compact_fc3, const AtomTriplets *atom_triplets, + const double *masses, const long *band_indices, const long symmetrize_fc3_q, const double cutoff_frequency, const long triplet_index, /* only for print */ const long num_triplets, /* only for print */ @@ -123,13 +125,22 @@ void itr_get_interaction_at_triplet( double *freqs[3]; lapack_complex_double *eigvecs[3]; double q_vecs[3][3]; + double vec_tmp[3]; for (j = 0; j < 3; j++) { for (k = 0; k < 3; k++) { q_vecs[j][k] = ((double)bzgrid->addresses[triplet[j]][k]) / bzgrid->D_diag[k]; } - bzg_multiply_matrix_vector_ld3(q_vecs[j], bzgrid->Q, q_vecs[j]); + + for (k = 0; k < 3; k++) { + vec_tmp[k] = bzgrid->Q[k][0] * q_vecs[j][0] + + bzgrid->Q[k][1] * q_vecs[j][1] + + bzgrid->Q[k][2] * q_vecs[j][2]; + } + for (k = 0; k < 3; k++) { + q_vecs[j][k] = vec_tmp[k]; + } } if (symmetrize_fc3_q) { diff --git a/c/interaction.h b/c/interaction.h index c345f742..a3c27c71 100644 --- a/c/interaction.h +++ b/c/interaction.h @@ -35,25 +35,26 @@ #ifndef __interaction_H__ #define __interaction_H__ -#include "bzgrid.h" #include "lapack_wrapper.h" #include "phonoc_array.h" #include "real_to_reciprocal.h" +#include "recgrid.h" void itr_get_interaction( Darray *fc3_normal_squared, const char *g_zero, const Darray *frequencies, const lapack_complex_double *eigenvectors, const long (*triplets)[3], - const long num_triplets, const ConstBZGrid *bzgrid, const double *fc3, - const long is_compact_fc3, const AtomTriplets *atom_triplets, - const double *masses, const long *band_indices, const long symmetrize_fc3_q, + const long num_triplets, const RecgridConstBZGrid *bzgrid, + const double *fc3, const long is_compact_fc3, + const AtomTriplets *atom_triplets, const double *masses, + const long *band_indices, const long symmetrize_fc3_q, const double cutoff_frequency, const long openmp_per_triplets); void itr_get_interaction_at_triplet( double *fc3_normal_squared, const long num_band0, const long num_band, const long (*g_pos)[4], const long num_g_pos, const double *frequencies, const lapack_complex_double *eigenvectors, const long triplet[3], - const ConstBZGrid *bzgrid, const double *fc3, const long is_compact_fc3, - const AtomTriplets *atom_triplets, const double *masses, - const long *band_indices, const long symmetrize_fc3_q, + const RecgridConstBZGrid *bzgrid, const double *fc3, + const long is_compact_fc3, const AtomTriplets *atom_triplets, + const double *masses, const long *band_indices, const long symmetrize_fc3_q, const double cutoff_frequency, const long triplet_index, /* only for print */ const long num_triplets, /* only for print */ diff --git a/c/isotope.c b/c/isotope.c index 8e8cd5c4..6a146bc2 100644 --- a/c/isotope.c +++ b/c/isotope.c @@ -36,9 +36,9 @@ #include +#include "funcs.h" #include "lapack_wrapper.h" #include "phonoc_const.h" -#include "phonoc_utils.h" void iso_get_isotope_scattering_strength( double *gamma, const long grid_point, const long *ir_grid_points, @@ -75,8 +75,8 @@ void iso_get_isotope_scattering_strength( } sum_g = 0; #ifdef _OPENMP -#pragma omp parallel for private(gp, k, l, m, f, e1_r, e1_i, a, b, dist, sum_g_k) reduction(+ \ - : sum_g) +#pragma omp parallel for private(gp, k, l, m, f, e1_r, e1_i, a, b, dist, \ + sum_g_k) reduction(+ : sum_g) #endif for (j = 0; j < num_grid_points; j++) { gp = ir_grid_points[j]; @@ -86,7 +86,7 @@ void iso_get_isotope_scattering_strength( if (f < cutoff_frequency) { continue; } - dist = phonoc_gaussian(f - f0[i], sigma); + dist = funcs_gaussian(f - f0[i], sigma); for (l = 0; l < num_band / 3; l++) { /* elements */ a = 0; b = 0; @@ -162,7 +162,7 @@ void iso_get_thm_isotope_scattering_strength( #ifdef _OPENMP #pragma omp parallel for private(j, k, l, m, f, gp, e1_r, e1_i, a, b, dist, \ - sum_g_k) + sum_g_k) #endif for (i = 0; i < num_grid_points; i++) { gp = ir_grid_points[i]; diff --git a/c/phono3py.c b/c/phono3py.c index e5cee178..4f7f625b 100644 --- a/c/phono3py.c +++ b/c/phono3py.c @@ -37,10 +37,8 @@ #include #include -#include "bzgrid.h" #include "collision_matrix.h" #include "fc3.h" -#include "grgrid.h" #include "imag_self_energy_with_g.h" #include "interaction.h" #include "isotope.h" @@ -50,6 +48,7 @@ #include "pp_collision.h" #include "real_self_energy.h" #include "real_to_reciprocal.h" +#include "recgrid.h" #include "tetrahedron_method.h" #include "triplet.h" #include "triplet_iw.h" @@ -69,11 +68,12 @@ long ph3py_get_interaction( const long *band_indices, const long symmetrize_fc3_q, const long make_r0_average, const char *all_shortest, const double cutoff_frequency, const long openmp_per_triplets) { - ConstBZGrid *bzgrid; + RecgridConstBZGrid *bzgrid; AtomTriplets *atom_triplets; long i, j; - if ((bzgrid = (ConstBZGrid *)malloc(sizeof(ConstBZGrid))) == NULL) { + if ((bzgrid = (RecgridConstBZGrid *)malloc(sizeof(RecgridConstBZGrid))) == + NULL) { warning_print("Memory could not be allocated."); return 0; } @@ -132,11 +132,12 @@ long ph3py_get_pp_collision( const long symmetrize_fc3_q, const long make_r0_average, const char *all_shortest, const double cutoff_frequency, const long openmp_per_triplets) { - ConstBZGrid *bzgrid; + RecgridConstBZGrid *bzgrid; AtomTriplets *atom_triplets; long i, j; - if ((bzgrid = (ConstBZGrid *)malloc(sizeof(ConstBZGrid))) == NULL) { + if ((bzgrid = (RecgridConstBZGrid *)malloc(sizeof(RecgridConstBZGrid))) == + NULL) { warning_print("Memory could not be allocated."); return 0; } @@ -196,11 +197,12 @@ long ph3py_get_pp_collision_with_sigma( const long symmetrize_fc3_q, const long make_r0_average, const char *all_shortest, const double cutoff_frequency, const long openmp_per_triplets) { - ConstBZGrid *bzgrid; + RecgridConstBZGrid *bzgrid; AtomTriplets *atom_triplets; long i, j; - if ((bzgrid = (ConstBZGrid *)malloc(sizeof(ConstBZGrid))) == NULL) { + if ((bzgrid = (RecgridConstBZGrid *)malloc(sizeof(RecgridConstBZGrid))) == + NULL) { warning_print("Memory could not be allocated."); return 0; } @@ -397,10 +399,11 @@ long ph3py_get_BZ_triplets_at_q(long (*triplets)[3], const long grid_point, const long num_map_triplets, const long D_diag[3], const long Q[3][3], const long bz_grid_type) { - ConstBZGrid *bzgrid; + RecgridConstBZGrid *bzgrid; long i, j, num_ir; - if ((bzgrid = (ConstBZGrid *)malloc(sizeof(ConstBZGrid))) == NULL) { + if ((bzgrid = (RecgridConstBZGrid *)malloc(sizeof(RecgridConstBZGrid))) == + NULL) { warning_print("Memory could not be allocated."); return 0; } @@ -436,10 +439,11 @@ long ph3py_get_integration_weight( const long bz_grid_type, const double *frequencies1, const long num_band1, const double *frequencies2, const long num_band2, const long tp_type, const long openmp_per_triplets) { - ConstBZGrid *bzgrid; + RecgridConstBZGrid *bzgrid; long i; - if ((bzgrid = (ConstBZGrid *)malloc(sizeof(ConstBZGrid))) == NULL) { + if ((bzgrid = (RecgridConstBZGrid *)malloc(sizeof(RecgridConstBZGrid))) == + NULL) { warning_print("Memory could not be allocated."); return 0; } @@ -471,126 +475,6 @@ void ph3py_get_integration_weight_with_sigma( num_triplets, frequencies, num_band, tp_type); } -/* From single address to grid index */ -long ph3py_get_grid_index_from_address(const long address[3], - const long D_diag[3]) { - return grg_get_grid_index(address, D_diag); -} - -void ph3py_get_gr_grid_addresses(long gr_grid_addresses[][3], - const long D_diag[3]) { - grg_get_all_grid_addresses(gr_grid_addresses, D_diag); -} - -long ph3py_get_reciprocal_rotations(long rec_rotations[48][3][3], - const long (*rotations)[3][3], - const long num_rot, - const long is_time_reversal) { - return grg_get_reciprocal_point_group(rec_rotations, rotations, num_rot, - is_time_reversal, 1); -} - -/* Rotation matrices with respect to reciprocal basis vectors are - * transformed to those for GRGrid. This set of the rotations are - * used always in GRGrid handling. */ -long ph3py_transform_rotations(long (*transformed_rots)[3][3], - const long (*rotations)[3][3], - const long num_rot, const long D_diag[3], - const long Q[3][3]) { - return grg_transform_rotations(transformed_rots, rotations, num_rot, D_diag, - Q); -} - -long ph3py_get_snf3x3(long D_diag[3], long P[3][3], long Q[3][3], - const long A[3][3]) { - return grg_get_snf3x3(D_diag, P, Q, A); -} - -/* The rotations are those after proper transformation in GRGrid. */ -long ph3py_get_ir_grid_map(long *ir_grid_map, const long D_diag[3], - const long PS[3], const long (*grg_rotations)[3][3], - const long num_rot) { - long num_ir, i; - - grg_get_ir_grid_map(ir_grid_map, grg_rotations, num_rot, D_diag, PS); - - num_ir = 0; - for (i = 0; i < D_diag[0] * D_diag[1] * D_diag[2]; i++) { - if (ir_grid_map[i] == i) { - num_ir++; - } - } - - return num_ir; -} - -long ph3py_get_bz_grid_addresses(long (*bz_grid_addresses)[3], long *bz_map, - long *bzg2grg, const long D_diag[3], - const long Q[3][3], const long PS[3], - const double rec_lattice[3][3], - const long type) { - BZGrid *bzgrid; - long i, j, size; - - if ((bzgrid = (BZGrid *)malloc(sizeof(BZGrid))) == NULL) { - warning_print("Memory could not be allocated."); - return 0; - } - - bzgrid->addresses = bz_grid_addresses; - bzgrid->gp_map = bz_map; - bzgrid->bzg2grg = bzg2grg; - bzgrid->type = type; - for (i = 0; i < 3; i++) { - bzgrid->D_diag[i] = D_diag[i]; - bzgrid->PS[i] = PS[i]; - for (j = 0; j < 3; j++) { - bzgrid->Q[i][j] = Q[i][j]; - bzgrid->reclat[i][j] = rec_lattice[i][j]; - } - } - - if (bzg_get_bz_grid_addresses(bzgrid)) { - size = bzgrid->size; - } else { - size = 0; - } - - free(bzgrid); - bzgrid = NULL; - - return size; -} - -long ph3py_rotate_bz_grid_index(const long bz_grid_index, - const long rotation[3][3], - const long (*bz_grid_addresses)[3], - const long *bz_map, const long D_diag[3], - const long PS[3], const long bz_grid_type) { - ConstBZGrid *bzgrid; - long i, rot_bz_gp; - - if ((bzgrid = (ConstBZGrid *)malloc(sizeof(ConstBZGrid))) == NULL) { - warning_print("Memory could not be allocated."); - return 0; - } - - bzgrid->addresses = bz_grid_addresses; - bzgrid->gp_map = bz_map; - bzgrid->type = bz_grid_type; - for (i = 0; i < 3; i++) { - bzgrid->D_diag[i] = D_diag[i]; - bzgrid->PS[i] = PS[i]; - } - - rot_bz_gp = bzg_rotate_grid_index(bz_grid_index, rotation, bzgrid); - - free(bzgrid); - bzgrid = NULL; - - return rot_bz_gp; -} - void ph3py_symmetrize_collision_matrix(double *collision_matrix, const long num_column, const long num_temp, @@ -721,9 +605,10 @@ long ph3py_get_neighboring_gird_points( const long bz_grid_type, const long num_grid_points, const long num_relative_grid_address) { long i; - ConstBZGrid *bzgrid; + RecgridConstBZGrid *bzgrid; - if ((bzgrid = (ConstBZGrid *)malloc(sizeof(ConstBZGrid))) == NULL) { + if ((bzgrid = (RecgridConstBZGrid *)malloc(sizeof(RecgridConstBZGrid))) == + NULL) { warning_print("Memory could not be allocated."); return 0; } @@ -766,9 +651,10 @@ long ph3py_get_thm_integration_weights_at_grid_points( long i, j, k, bi; long vertices[24][4]; double freq_vertices[24][4]; - ConstBZGrid *bzgrid; + RecgridConstBZGrid *bzgrid; - if ((bzgrid = (ConstBZGrid *)malloc(sizeof(ConstBZGrid))) == NULL) { + if ((bzgrid = (RecgridConstBZGrid *)malloc(sizeof(RecgridConstBZGrid))) == + NULL) { warning_print("Memory could not be allocated."); return 0; } diff --git a/c/pp_collision.c b/c/pp_collision.c index 9cd6b9e8..d083a6c8 100644 --- a/c/pp_collision.c +++ b/c/pp_collision.c @@ -37,12 +37,13 @@ #include #include +#include "funcs.h" #include "imag_self_energy_with_g.h" #include "interaction.h" #include "lapack_wrapper.h" #include "phonoc_array.h" -#include "phonoc_utils.h" #include "real_to_reciprocal.h" +#include "recgrid.h" #include "triplet.h" #include "triplet_iw.h" @@ -51,9 +52,10 @@ static void get_collision( const long num_temps, const double *temperatures, const double *g, const char *g_zero, const double *frequencies, const lapack_complex_double *eigenvectors, const long triplet[3], - const long triplet_weight, const ConstBZGrid *bzgrid, const double *fc3, - const long is_compact_fc3, const AtomTriplets *atom_triplets, - const double *masses, const long *band_indices, const long symmetrize_fc3_q, + const long triplet_weight, const RecgridConstBZGrid *bzgrid, + const double *fc3, const long is_compact_fc3, + const AtomTriplets *atom_triplets, const double *masses, + const long *band_indices, const long symmetrize_fc3_q, const double cutoff_frequency, const long openmp_per_triplets); static void finalize_ise(double *imag_self_energy, const double *ise, const long (*bz_grid_address)[3], @@ -66,11 +68,12 @@ void ppc_get_pp_collision( const long relative_grid_address[24][4][3], /* thm */ const double *frequencies, const lapack_complex_double *eigenvectors, const long (*triplets)[3], const long num_triplets, - const long *triplet_weights, const ConstBZGrid *bzgrid, const double *fc3, - const long is_compact_fc3, const AtomTriplets *atom_triplets, - const double *masses, const Larray *band_indices, - const Darray *temperatures, const long is_NU, const long symmetrize_fc3_q, - const double cutoff_frequency, const long openmp_per_triplets) { + const long *triplet_weights, const RecgridConstBZGrid *bzgrid, + const double *fc3, const long is_compact_fc3, + const AtomTriplets *atom_triplets, const double *masses, + const Larray *band_indices, const Darray *temperatures, const long is_NU, + const long symmetrize_fc3_q, const double cutoff_frequency, + const long openmp_per_triplets) { long i; long num_band, num_band0, num_band_prod, num_temps; double *ise, *freqs_at_gp, *g; @@ -136,11 +139,12 @@ void ppc_get_pp_collision_with_sigma( double *imag_self_energy, const double sigma, const double sigma_cutoff, const double *frequencies, const lapack_complex_double *eigenvectors, const long (*triplets)[3], const long num_triplets, - const long *triplet_weights, const ConstBZGrid *bzgrid, const double *fc3, - const long is_compact_fc3, const AtomTriplets *atom_triplets, - const double *masses, const Larray *band_indices, - const Darray *temperatures, const long is_NU, const long symmetrize_fc3_q, - const double cutoff_frequency, const long openmp_per_triplets) { + const long *triplet_weights, const RecgridConstBZGrid *bzgrid, + const double *fc3, const long is_compact_fc3, + const AtomTriplets *atom_triplets, const double *masses, + const Larray *band_indices, const Darray *temperatures, const long is_NU, + const long symmetrize_fc3_q, const double cutoff_frequency, + const long openmp_per_triplets) { long i; long num_band, num_band0, num_band_prod, num_temps; long const_adrs_shift; @@ -207,9 +211,10 @@ static void get_collision( const long num_temps, const double *temperatures, const double *g, const char *g_zero, const double *frequencies, const lapack_complex_double *eigenvectors, const long triplet[3], - const long triplet_weight, const ConstBZGrid *bzgrid, const double *fc3, - const long is_compact_fc3, const AtomTriplets *atom_triplets, - const double *masses, const long *band_indices, const long symmetrize_fc3_q, + const long triplet_weight, const RecgridConstBZGrid *bzgrid, + const double *fc3, const long is_compact_fc3, + const AtomTriplets *atom_triplets, const double *masses, + const long *band_indices, const long symmetrize_fc3_q, const double cutoff_frequency, const long openmp_per_triplets) { long i; long num_band_prod, num_g_pos; diff --git a/c/pp_collision.h b/c/pp_collision.h index b9ccbace..e5fb0d5e 100644 --- a/c/pp_collision.h +++ b/c/pp_collision.h @@ -35,29 +35,31 @@ #ifndef __pp_collision_H__ #define __pp_collision_H__ -#include "bzgrid.h" #include "lapack_wrapper.h" #include "phonoc_array.h" #include "real_to_reciprocal.h" +#include "recgrid.h" void ppc_get_pp_collision( double *imag_self_energy, const long relative_grid_address[24][4][3], /* thm */ const double *frequencies, const lapack_complex_double *eigenvectors, const long (*triplets)[3], const long num_triplets, - const long *triplet_weights, const ConstBZGrid *bzgrid, const double *fc3, - const long is_compact_fc3, const AtomTriplets *atom_triplets, - const double *masses, const Larray *band_indices, - const Darray *temperatures, const long is_NU, const long symmetrize_fc3_q, - const double cutoff_frequency, const long openmp_per_triplets); + const long *triplet_weights, const RecgridConstBZGrid *bzgrid, + const double *fc3, const long is_compact_fc3, + const AtomTriplets *atom_triplets, const double *masses, + const Larray *band_indices, const Darray *temperatures, const long is_NU, + const long symmetrize_fc3_q, const double cutoff_frequency, + const long openmp_per_triplets); void ppc_get_pp_collision_with_sigma( double *imag_self_energy, const double sigma, const double sigma_cutoff, const double *frequencies, const lapack_complex_double *eigenvectors, const long (*triplets)[3], const long num_triplets, - const long *triplet_weights, const ConstBZGrid *bzgrid, const double *fc3, - const long is_compact_fc3, const AtomTriplets *atom_triplets, - const double *masses, const Larray *band_indices, - const Darray *temperatures, const long is_NU, const long symmetrize_fc3_q, - const double cutoff_frequency, const long openmp_per_triplets); + const long *triplet_weights, const RecgridConstBZGrid *bzgrid, + const double *fc3, const long is_compact_fc3, + const AtomTriplets *atom_triplets, const double *masses, + const Larray *band_indices, const Darray *temperatures, const long is_NU, + const long symmetrize_fc3_q, const double cutoff_frequency, + const long openmp_per_triplets); #endif diff --git a/c/real_self_energy.c b/c/real_self_energy.c index c723def8..f5d9875f 100644 --- a/c/real_self_energy.c +++ b/c/real_self_energy.c @@ -37,8 +37,8 @@ #include #include +#include "funcs.h" #include "phonoc_array.h" -#include "phonoc_utils.h" #include "real_to_reciprocal.h" static double get_real_self_energy_at_band( @@ -160,10 +160,10 @@ static double sum_real_self_energy_at_band( shift = 0; for (i = 0; i < num_band; i++) { if (freqs1[i] > cutoff_frequency) { - n1 = phonoc_bose_einstein(freqs1[i], temperature); + n1 = funcs_bose_einstein(freqs1[i], temperature); for (j = 0; j < num_band; j++) { if (freqs2[j] > cutoff_frequency) { - n2 = phonoc_bose_einstein(freqs2[j], temperature); + n2 = funcs_bose_einstein(freqs2[j], temperature); f1 = fpoint + freqs1[i] + freqs2[j]; f2 = fpoint - freqs1[i] - freqs2[j]; f3 = fpoint - freqs1[i] + freqs2[j]; diff --git a/c/recgrid.c b/c/recgrid.c new file mode 100644 index 00000000..79c6001e --- /dev/null +++ b/c/recgrid.c @@ -0,0 +1,248 @@ +/* Copyright (C) 2021 Atsushi Togo */ +/* All rights reserved. */ + +/* This file is part of phonopy. */ + +/* Redistribution and use in source and binary forms, with or without */ +/* modification, are permitted provided that the following conditions */ +/* are met: */ + +/* * Redistributions of source code must retain the above copyright */ +/* notice, this list of conditions and the following disclaimer. */ + +/* * Redistributions in binary form must reproduce the above copyright */ +/* notice, this list of conditions and the following disclaimer in */ +/* the documentation and/or other materials provided with the */ +/* distribution. */ + +/* * Neither the name of the phonopy project nor the names of its */ +/* contributors may be used to endorse or promote products derived */ +/* from this software without specific prior written permission. */ + +/* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS */ +/* "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT */ +/* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS */ +/* FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE */ +/* COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, */ +/* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, */ +/* BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; */ +/* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER */ +/* CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT */ +/* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN */ +/* ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE */ +/* POSSIBILITY OF SUCH DAMAGE. */ + +#include "recgrid.h" + +#include +#include + +#include "bzgrid.h" +#include "grgrid.h" +#include "gridsys.h" +#include "lagrid.h" + +#define GRID_TOLERANCE_FACTOR 0.01 + +void recgrid_get_all_grid_addresses(long (*gr_grid_addresses)[3], + const long D_diag[3]) { + grg_get_all_grid_addresses(gr_grid_addresses, D_diag); +} + +void recgrid_get_double_grid_address(long address_double[3], + const long address[3], const long PS[3]) { + grg_get_double_grid_address(address_double, address, PS); +} + +void recgrid_get_grid_address_from_index(long address[3], const long grid_index, + const long D_diag[3]) { + grg_get_grid_address_from_index(address, grid_index, D_diag); +} + +long recgrid_get_double_grid_index(const long address_double[3], + const long D_diag[3], const long PS[3]) { + return grg_get_double_grid_index(address_double, D_diag, PS); +} + +long recgrid_get_grid_index_from_address(const long address[3], + const long D_diag[3]) { + return grg_get_grid_index(address, D_diag); +} + +long recgrid_rotate_grid_index(const long grid_index, const long rotation[3][3], + const long D_diag[3], const long PS[3]) { + return grg_rotate_grid_index(grid_index, rotation, D_diag, PS); +} + +long recgrid_get_reciprocal_point_group(long rec_rotations[48][3][3], + const long (*rotations)[3][3], + const long num_rot, + const long is_time_reversal, + const long is_transpose) { + return grg_get_reciprocal_point_group(rec_rotations, rotations, num_rot, + is_time_reversal, is_transpose); +} + +long recgrid_get_snf3x3(long D_diag[3], long P[3][3], long Q[3][3], + const long A[3][3]) { + return grg_get_snf3x3(D_diag, P, Q, A); +} + +long recgrid_transform_rotations(long (*transformed_rots)[3][3], + const long (*rotations)[3][3], + const long num_rot, const long D_diag[3], + const long Q[3][3]) { + long succeeded; + succeeded = grg_transform_rotations(transformed_rots, rotations, num_rot, + D_diag, Q); + return succeeded; +} + +long recgrid_get_ir_grid_map(long *ir_grid_map, const long (*rotations)[3][3], + const long num_rot, const long D_diag[3], + const long PS[3]) { + long num_ir, i; + + grg_get_ir_grid_map(ir_grid_map, rotations, num_rot, D_diag, PS); + + num_ir = 0; + for (i = 0; i < D_diag[0] * D_diag[1] * D_diag[2]; i++) { + if (ir_grid_map[i] == i) { + num_ir++; + } + } + return num_ir; +} + +long recgrid_get_bz_grid_addresses(long (*bz_grid_addresses)[3], long *bz_map, + long *bzg2grg, const long D_diag[3], + const long Q[3][3], const long PS[3], + const double rec_lattice[3][3], + const long bz_grid_type) { + RecgridBZGrid *bzgrid; + long i, j, size; + + if ((bzgrid = (RecgridBZGrid *)malloc(sizeof(RecgridBZGrid))) == NULL) { + warning_print("Memory could not be allocated."); + return 0; + } + + bzgrid->addresses = bz_grid_addresses; + bzgrid->gp_map = bz_map; + bzgrid->bzg2grg = bzg2grg; + bzgrid->type = bz_grid_type; + for (i = 0; i < 3; i++) { + bzgrid->D_diag[i] = D_diag[i]; + bzgrid->PS[i] = PS[i]; + for (j = 0; j < 3; j++) { + bzgrid->Q[i][j] = Q[i][j]; + bzgrid->reclat[i][j] = rec_lattice[i][j]; + } + } + + if (bzg_get_bz_grid_addresses(bzgrid)) { + size = bzgrid->size; + } else { + size = 0; + } + + free(bzgrid); + bzgrid = NULL; + + return size; +} + +long recgrid_rotate_bz_grid_index(const long bz_grid_index, + const long rotation[3][3], + const long (*bz_grid_addresses)[3], + const long *bz_map, const long D_diag[3], + const long PS[3], const long bz_grid_type) { + RecgridConstBZGrid *bzgrid; + long i, rot_bz_gp; + + if ((bzgrid = (RecgridConstBZGrid *)malloc(sizeof(RecgridConstBZGrid))) == + NULL) { + warning_print("Memory could not be allocated."); + return 0; + } + + bzgrid->addresses = bz_grid_addresses; + bzgrid->gp_map = bz_map; + bzgrid->type = bz_grid_type; + for (i = 0; i < 3; i++) { + bzgrid->D_diag[i] = D_diag[i]; + bzgrid->PS[i] = PS[i]; + } + + rot_bz_gp = bzg_rotate_grid_index(bz_grid_index, rotation, bzgrid); + + free(bzgrid); + bzgrid = NULL; + + return rot_bz_gp; +} + +/* Note: Tolerance in squared distance. */ +double recgrid_get_tolerance_for_BZ_reduction(const RecgridBZGrid *bzgrid) { + long i, j; + double tolerance; + double length[3]; + double reclatQ[3][3]; + + for (i = 0; i < 3; i++) { + for (j = 0; j < 3; j++) { + reclatQ[i][j] = bzgrid->reclat[i][0] * bzgrid->Q[0][j] + + bzgrid->reclat[i][1] * bzgrid->Q[1][j] + + bzgrid->reclat[i][2] * bzgrid->Q[2][j]; + } + } + + for (i = 0; i < 3; i++) { + length[i] = 0; + for (j = 0; j < 3; j++) { + length[i] += reclatQ[j][i] * reclatQ[j][i]; + } + length[i] /= bzgrid->D_diag[i] * bzgrid->D_diag[i]; + } + tolerance = length[0]; + for (i = 1; i < 3; i++) { + if (tolerance < length[i]) { + tolerance = length[i]; + } + } + tolerance *= GRID_TOLERANCE_FACTOR; + + return tolerance; +} + +RecgridMats *recgrid_alloc_RotMats(const long size) { + RecgridMats *rotmats; + + rotmats = NULL; + + if ((rotmats = (RecgridMats *)malloc(sizeof(RecgridMats))) == NULL) { + warning_print("Memory could not be allocated."); + return NULL; + } + + rotmats->size = size; + if (size > 0) { + if ((rotmats->mat = (long(*)[3][3])malloc(sizeof(long[3][3]) * size)) == + NULL) { + warning_print("Memory could not be allocated "); + warning_print("(RecgridMats, line %d, %s).\n", __LINE__, __FILE__); + free(rotmats); + rotmats = NULL; + return NULL; + } + } + return rotmats; +} + +void recgrid_free_RotMats(RecgridMats *rotmats) { + if (rotmats->size > 0) { + free(rotmats->mat); + rotmats->mat = NULL; + } + free(rotmats); +} diff --git a/c/recgrid.h b/c/recgrid.h new file mode 100644 index 00000000..ac0f988e --- /dev/null +++ b/c/recgrid.h @@ -0,0 +1,289 @@ +/* Copyright (C) 2021 Atsushi Togo */ +/* All rights reserved. */ + +/* This file is part of phonopy. */ + +/* Redistribution and use in source and binary forms, with or without */ +/* modification, are permitted provided that the following conditions */ +/* are met: */ + +/* * Redistributions of source code must retain the above copyright */ +/* notice, this list of conditions and the following disclaimer. */ + +/* * Redistributions in binary form must reproduce the above copyright */ +/* notice, this list of conditions and the following disclaimer in */ +/* the documentation and/or other materials provided with the */ +/* distribution. */ + +/* * Neither the name of the phonopy project nor the names of its */ +/* contributors may be used to endorse or promote products derived */ +/* from this software without specific prior written permission. */ + +/* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS */ +/* "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT */ +/* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS */ +/* FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE */ +/* COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, */ +/* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, */ +/* BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; */ +/* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER */ +/* CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT */ +/* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN */ +/* ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE */ +/* POSSIBILITY OF SUCH DAMAGE. */ + +#ifndef __recgrid_H__ +#define __recgrid_H__ + +#ifdef __cplusplus +extern "C" { +#endif + +typedef struct { + long size; + long (*mat)[3][3]; +} RecgridMats; + +/* Data structure of Brillouin zone grid + * + * size : long + * Number of grid points in Brillouin zone including its surface. + * D_diag : long array + * Diagonal part of matrix D of SNF. + * shape=(3, ) + * PS : long array + * Matrix P of SNF multiplied by shift. + * shape=(3, ) + * gp_map : long array + * Type1 : Twice enlarged grid space along basis vectors. + * Grid index is recovered in the same way as regular grid. + * shape=(prod(mesh * 2), ) + * Type2 : In the last index, multiplicity and array index of + * each address of the grid point are stored. Here, + * multiplicity means the number of translationally + * equivalent grid points in BZ. + * shape=(prod(mesh), 2) -> flattened. + * addresses : long array + * Grid point addresses. + * shape=(size, 3) + * reclat : double array + * Reciprocal basis vectors given as column vectors. + * shape=(3, 3) + * type : long + * 1 or 2. */ +typedef struct { + long size; + long D_diag[3]; + long Q[3][3]; + long PS[3]; + long *gp_map; + long *bzg2grg; + long (*addresses)[3]; + double reclat[3][3]; + long type; +} RecgridBZGrid; + +typedef struct { + long size; + long D_diag[3]; + long Q[3][3]; + long PS[3]; + const long *gp_map; + const long *bzg2grg; + const long (*addresses)[3]; + double reclat[3][3]; + long type; +} RecgridConstBZGrid; + +/* Generalized regular (GR) grid + +Integer grid matrix M_g is unimodular transformed to integer diagnonal matrix D +by D = P M_g Q, which can be achieved by Smith normal form like transformation. +P and Q used are integer unimodular matrices with determinant=1. + +S in PS is doubled shift with respect to microzone basis vectors, i.e., +half-grid shift along an axis corresponds to 1. + +*/ + +/** + * @brief Return all GR-grid addresses with respect to n_1, n_2, n_3 + * + * @param gr_grid_addresses All GR-grid addresses + * @param D_diag Numbers of divisions along a, b, c directions of GR-grid + * @return void + */ +void recgrid_get_all_grid_addresses(long (*gr_grid_addresses)[3], + const long D_diag[3]); + +/** + * @brief Return double grid address in GR-grid + * + * @param address_double Double grid address, i.e., possibly with shift in + * GR-grid + * @param address Single grid address in GR-grid + * @param PS Shift in GR-grid + * @return void + */ +void recgrid_get_double_grid_address(long address_double[3], + const long address[3], const long PS[3]); + +/** + * @brief Return single grid address in GR-grid with given grid point index. + * + * @param address Single grid address in GR-grid + * @param grid_index Grid point index in GR-grid + * @param D_diag Numbers of divisions along a, b, c directions of GR-grid + * @return void + */ +void recgrid_get_grid_address_from_index(long address[3], const long grid_index, + const long D_diag[3]); + +/** + * @brief Return grid point index of double grid address in GR-grid + * + * @param address_double Double grid address, i.e., possibly with shift in + * GR-grid + * @param D_diag Numbers of divisions along a, b, c directions of GR-grid + * @param PS Shift in GR-grid + * @return long + */ +long recgrid_get_double_grid_index(const long address_double[3], + const long D_diag[3], const long PS[3]); + +/** + * @brief Return grid point index of single grid address in GR-grid + * + * @param address Single grid address in GR-grid + * @param D_diag Numbers of divisions along a, b, c directions of GR-grid + * @return long + */ +long recgrid_get_grid_index_from_address(const long address[3], + const long D_diag[3]); + +/** + * @brief Return grid point index of rotated address of given grid point index. + * + * @param grid_index Grid point index in GR-grid + * @param rotation Transformed rotation in reciprocal space tilde-R^T + * @param D_diag Numbers of divisions along a, b, c directions of GR-grid + * @param PS Shift in GR-grid + * @return long + */ +long recgrid_rotate_grid_index(const long grid_index, const long rotation[3][3], + const long D_diag[3], const long PS[3]); + +/** + * @brief Return {R^T} of crystallographic point group {R} with and without time + * reversal symmetry. + * + * @param rec_rotations Rotations in reciprocal space {R^T} + * @param rotations Rotations in direct space {R} + * @param num_rot Number of given rotations |{R}| + * @param is_time_reversal With (1) or without (0) time reversal symmetry + * @param is_transpose With (1) or without (0) transpose of rotation matrices + * @return long + */ +long recgrid_get_reciprocal_point_group(long rec_rotations[48][3][3], + const long (*rotations)[3][3], + const long num_rot, + const long is_time_reversal, + const long is_transpose); +/** + * @brief Return D, P, Q of Smith normal form of A. + * + * @param D_diag Diagonal elements of diagnoal matrix D + * @param P Unimodular matrix P + * @param Q Unimodular matrix Q + * @param A Integer matrix + * @return long + */ +long recgrid_get_snf3x3(long D_diag[3], long P[3][3], long Q[3][3], + const long A[3][3]); + +/** + * @brief Transform {R^T} to {R^T} with respect to transformed microzone basis + * vectors in GR-grid + * + * @param transformed_rots Transformed rotation matrices in reciprocal space + * {tilde-R^T} + * @param rotations Original rotations matrices in reciprocal space {R^T} + * @param num_rot Number of rotation matrices + * @param D_diag Diagonal elements of diagnoal matrix D of Smith normal form + * @param Q Unimodular matrix Q of Smith normal form + * @return long + */ +long recgrid_transform_rotations(long (*transformed_rots)[3][3], + const long (*rotations)[3][3], + const long num_rot, const long D_diag[3], + const long Q[3][3]); + +/** + * @brief Return mapping table from GR-grid points to GR-ir-grid points + * + * @param ir_grid_map Grid point index mapping to ir-grid point indices with + * array size of prod(D_diag) + * @param rotations Transformed rotation matrices in reciprocal space + * @param num_rot Number of rotation matrices + * @param D_diag Diagonal elements of diagnoal matrix D of Smith normal form + * @param PS Shift in GR-grid + * @return long Number of ir_grid_points. + */ +long recgrid_get_ir_grid_map(long *ir_grid_map, const long (*rotations)[3][3], + const long num_rot, const long D_diag[3], + const long PS[3]); + +/** + * @brief Find shortest grid points from Gamma considering periodicity of + * reciprocal lattice. See the details in docstring of RecgridBZGrid + * + * @param bz_grid_addresses Grid point addresses of shortest grid points + * @param bz_map List of accumulated numbers of BZ grid points from the + * first GR grid point to the last grid point. In type-II, [0, 1, 3, 4, ...] + * means multiplicities of [1, 2, 1, ...], with len(bz_map)=product(D_diag) + 1. + * @param bzg2grg Mapping table of bz_grid_addresses to gr_grid_addresses. In + * type-II, len(bzg2grg) == len(bz_grid_addresses) <= (D_diag[0] + 1) * + * (D_diag[1] + 1) * (D_diag[2] + 1). + * @param D_diag Diagonal elements of diagnoal matrix D of Smith normal form + * @param Q Unimodular matrix Q of Smith normal form + * @param PS Shift in GR-grid + * @param rec_lattice Reduced reciprocal basis vectors in column vectors + * @param bz_grid_type Data structure type I (old and sparse) or II (new and + * dense, recommended) of bz_map + * @return long Number of bz_grid_addresses stored. + */ +long recgrid_get_bz_grid_addresses(long (*bz_grid_addresses)[3], long *bz_map, + long *bzg2grg, const long D_diag[3], + const long Q[3][3], const long PS[3], + const double rec_lattice[3][3], + const long bz_grid_type); + +/** + * @brief Return index of rotated bz grid point + * + * @param bz_grid_index BZ grid point index + * @param rotation Transformed rotation in reciprocal space tilde-R^T + * @param bz_grid_addresses BZ grid point adddresses + * @param bz_map List of accumulated numbers of BZ grid points from the + * first GR grid point to the last grid point. In type-II, [0, 1, 3, 4, ...] + * means multiplicities of [1, 2, 1, ...], with len(bz_map)=product(D_diag) + 1. + * @param D_diag Numbers of divisions along a, b, c directions of GR-grid + * @param PS Shift in GR-grid + * @param bz_grid_type Data structure type I (old and sparse) or II (new and + * dense, recommended) of bz_map + * @return long + */ +long recgrid_rotate_bz_grid_index(const long bz_grid_index, + const long rotation[3][3], + const long (*bz_grid_addresses)[3], + const long *bz_map, const long D_diag[3], + const long PS[3], const long bz_grid_type); +double recgrid_get_tolerance_for_BZ_reduction(const RecgridBZGrid *bzgrid); +RecgridMats *recgrid_alloc_RotMats(const long size); +void recgrid_free_RotMats(RecgridMats *rotmats); + +#ifdef __cplusplus +} +#endif + +#endif diff --git a/c/triplet.c b/c/triplet.c index 27453fc2..b943a243 100644 --- a/c/triplet.c +++ b/c/triplet.c @@ -36,12 +36,12 @@ #include "triplet.h" -#include "bzgrid.h" +#include "recgrid.h" #include "triplet_grid.h" #include "triplet_iw.h" long tpl_get_BZ_triplets_at_q(long (*triplets)[3], const long grid_point, - const ConstBZGrid *bzgrid, + const RecgridConstBZGrid *bzgrid, const long *map_triplets) { return tpk_get_BZ_triplets_at_q(triplets, grid_point, bzgrid, map_triplets); } @@ -62,9 +62,9 @@ void tpl_get_integration_weight( double *iw, char *iw_zero, const double *frequency_points, const long num_band0, const long relative_grid_address[24][4][3], const long (*triplets)[3], const long num_triplets, - const ConstBZGrid *bzgrid, const double *frequencies1, const long num_band1, - const double *frequencies2, const long num_band2, const long tp_type, - const long openmp_per_triplets) { + const RecgridConstBZGrid *bzgrid, const double *frequencies1, + const long num_band1, const double *frequencies2, const long num_band2, + const long tp_type, const long openmp_per_triplets) { long i, num_band_prod; long tp_relative_grid_address[2][24][4][3]; diff --git a/c/triplet.h b/c/triplet.h index 73c2ad51..5d5a478b 100644 --- a/c/triplet.h +++ b/c/triplet.h @@ -39,7 +39,7 @@ #include -#include "bzgrid.h" +#include "recgrid.h" /* Irreducible triplets of k-points are searched under conservation of */ /* :math:``\mathbf{k}_1 + \mathbf{k}_2 + \mathbf{k}_3 = \mathbf{G}``. */ @@ -58,15 +58,15 @@ long tpl_get_triplets_reciprocal_mesh_at_q( /* triplets[num_ir_triplets][3] = number of non-zero triplets weights*/ /* Number of ir-triplets is returned. */ long tpl_get_BZ_triplets_at_q(long (*triplets)[3], const long grid_point, - const ConstBZGrid *bzgrid, + const RecgridConstBZGrid *bzgrid, const long *map_triplets); void tpl_get_integration_weight( double *iw, char *iw_zero, const double *frequency_points, const long num_band0, const long relative_grid_address[24][4][3], const long (*triplets)[3], const long num_triplets, - const ConstBZGrid *bzgrid, const double *frequencies1, const long num_band1, - const double *frequencies2, const long num_band2, const long tp_type, - const long openmp_per_triplets); + const RecgridConstBZGrid *bzgrid, const double *frequencies1, + const long num_band1, const double *frequencies2, const long num_band2, + const long tp_type, const long openmp_per_triplets); void tpl_get_integration_weight_with_sigma( double *iw, char *iw_zero, const double sigma, const double sigma_cutoff, const double *frequency_points, const long num_band0, diff --git a/c/triplet_grid.c b/c/triplet_grid.c index 112e1dd6..e5cd980e 100644 --- a/c/triplet_grid.c +++ b/c/triplet_grid.c @@ -39,14 +39,11 @@ #include #include -#include "bzgrid.h" #include "grgrid.h" -#include "lagrid.h" -#include "triplet.h" static long get_ir_triplets_at_q(long *map_triplets, long *map_q, const long grid_point, const long D_diag[3], - const RotMats *rot_reciprocal, + const RecgridMats *rot_reciprocal, const long swappable); static long get_ir_triplets_at_q_perm_q1q2(long *map_triplets, const long *map_q, @@ -56,27 +53,27 @@ static long get_ir_triplets_at_q_noperm(long *map_triplets, const long *map_q, const long grid_point, const long D_diag[3]); static long get_BZ_triplets_at_q(long (*triplets)[3], const long grid_point, - const ConstBZGrid *bzgrid, + const RecgridConstBZGrid *bzgrid, const long *map_triplets); static void get_BZ_triplets_at_q_type1(long (*triplets)[3], const long grid_point, - const ConstBZGrid *bzgrid, + const RecgridConstBZGrid *bzgrid, const long *ir_q1_gps, const long num_ir); static void get_BZ_triplets_at_q_type2(long (*triplets)[3], const long grid_point, - const ConstBZGrid *bzgrid, + const RecgridConstBZGrid *bzgrid, const long *ir_q1_gps, const long num_ir); static double get_squared_distance(const long G[3], const double LQD_inv[3][3]); -static void get_LQD_inv(double LQD_inv[3][3], const ConstBZGrid *bzgrid); -static RotMats *get_reciprocal_point_group_with_q(const RotMats *rot_reciprocal, - const long D_diag[3], - const long grid_point); -static RotMats *get_reciprocal_point_group(const long (*rec_rotations_in)[3][3], - const long num_rot, - const long is_time_reversal, - const long is_transpose); +static void get_LQD_inv(double LQD_inv[3][3], const RecgridConstBZGrid *bzgrid); +static RecgridMats *get_reciprocal_point_group_with_q( + const RecgridMats *rot_reciprocal, const long D_diag[3], + const long grid_point); +static RecgridMats *get_reciprocal_point_group( + const long (*rec_rotations_in)[3][3], const long num_rot, + const long is_time_reversal, const long is_transpose); +static void copy_matrix_l3(long a[3][3], const long b[3][3]); long tpk_get_ir_triplets_at_q(long *map_triplets, long *map_q, const long grid_point, const long D_diag[3], @@ -84,7 +81,7 @@ long tpk_get_ir_triplets_at_q(long *map_triplets, long *map_q, const long (*rec_rotations_in)[3][3], const long num_rot, const long swappable) { long num_ir; - RotMats *rotations; + RecgridMats *rotations; rotations = get_reciprocal_point_group(rec_rotations_in, num_rot, is_time_reversal, 0); @@ -94,25 +91,25 @@ long tpk_get_ir_triplets_at_q(long *map_triplets, long *map_q, num_ir = get_ir_triplets_at_q(map_triplets, map_q, grid_point, D_diag, rotations, swappable); - bzg_free_RotMats(rotations); + recgrid_free_RotMats(rotations); rotations = NULL; return num_ir; } long tpk_get_BZ_triplets_at_q(long (*triplets)[3], const long grid_point, - const ConstBZGrid *bzgrid, + const RecgridConstBZGrid *bzgrid, const long *map_triplets) { return get_BZ_triplets_at_q(triplets, grid_point, bzgrid, map_triplets); } static long get_ir_triplets_at_q(long *map_triplets, long *map_q, const long grid_point, const long D_diag[3], - const RotMats *rot_reciprocal, + const RecgridMats *rot_reciprocal, const long swappable) { long i, num_ir_q, num_ir_triplets; long PS[3]; - RotMats *rot_reciprocal_q; + RecgridMats *rot_reciprocal_q; rot_reciprocal_q = NULL; @@ -124,14 +121,8 @@ static long get_ir_triplets_at_q(long *map_triplets, long *map_q, rot_reciprocal_q = get_reciprocal_point_group_with_q(rot_reciprocal, D_diag, grid_point); - grg_get_ir_grid_map(map_q, rot_reciprocal_q->mat, rot_reciprocal_q->size, - D_diag, PS); - num_ir_q = 0; - for (i = 0; i < D_diag[0] * D_diag[1] * D_diag[2]; i++) { - if (map_q[i] == i) { - num_ir_q++; - } - } + num_ir_q = recgrid_get_ir_grid_map(map_q, rot_reciprocal_q->mat, + rot_reciprocal_q->size, D_diag, PS); if (swappable) { num_ir_triplets = get_ir_triplets_at_q_perm_q1q2(map_triplets, map_q, @@ -141,7 +132,7 @@ static long get_ir_triplets_at_q(long *map_triplets, long *map_q, grid_point, D_diag); } - bzg_free_RotMats(rot_reciprocal_q); + recgrid_free_RotMats(rot_reciprocal_q); rot_reciprocal_q = NULL; return num_ir_triplets; @@ -156,21 +147,21 @@ static long get_ir_triplets_at_q_perm_q1q2(long *map_triplets, num_ir_triplets = 0; num_grid = D_diag[0] * D_diag[1] * D_diag[2]; - grg_get_grid_address_from_index(adrs0, grid_point, D_diag); + recgrid_get_grid_address_from_index(adrs0, grid_point, D_diag); // #ifdef _OPENMP // #pragma omp parallel for private(j, gp2, adrs1, adrs2) // #endif for (gp1 = 0; gp1 < num_grid; gp1++) { if (map_q[gp1] == gp1) { - grg_get_grid_address_from_index(adrs1, gp1, D_diag); + recgrid_get_grid_address_from_index(adrs1, gp1, D_diag); for (j = 0; j < 3; j++) { adrs2[j] = -adrs0[j] - adrs1[j]; } /* If map_q[gp2] is smaller than current gp1, map_q[gp2] should */ /* equal to a previous gp1 for which map_triplets is already */ /* filled. So the counter is not incremented. */ - gp2 = grg_get_grid_index(adrs2, D_diag); + gp2 = recgrid_get_grid_index_from_address(adrs2, D_diag); if (map_q[gp2] < gp1) { map_triplets[gp1] = map_q[gp2]; } else { @@ -215,7 +206,7 @@ static long get_ir_triplets_at_q_noperm(long *map_triplets, const long *map_q, } static long get_BZ_triplets_at_q(long (*triplets)[3], const long grid_point, - const ConstBZGrid *bzgrid, + const RecgridConstBZGrid *bzgrid, const long *map_triplets) { long gp1, num_ir; long *ir_q1_gps; @@ -252,7 +243,7 @@ static long get_BZ_triplets_at_q(long (*triplets)[3], const long grid_point, static void get_BZ_triplets_at_q_type1(long (*triplets)[3], const long grid_point, - const ConstBZGrid *bzgrid, + const RecgridConstBZGrid *bzgrid, const long *ir_q1_gps, const long num_ir) { long i, j, gp2, num_gp, num_bzgp, bz0, bz1, bz2; @@ -267,7 +258,7 @@ static void get_BZ_triplets_at_q_type1(long (*triplets)[3], bz_adrs = bzgrid->addresses; get_LQD_inv(LQD_inv, bzgrid); /* This tolerance is used to be consistent to BZ reduction in bzgrid. */ - tolerance = bzg_get_tolerance_for_BZ_reduction((BZGrid *)bzgrid); + tolerance = recgrid_get_tolerance_for_BZ_reduction((RecgridBZGrid *)bzgrid); for (i = 0; i < 3; i++) { bz_adrs0[i] = bz_adrs[grid_point][i]; @@ -277,14 +268,14 @@ static void get_BZ_triplets_at_q_type1(long (*triplets)[3], #ifdef _OPENMP #pragma omp parallel for private(j, gp2, bzgp, G, bz_adrs1, bz_adrs2, d2, \ - min_d2, bz0, bz1, bz2) + min_d2, bz0, bz1, bz2) #endif for (i = 0; i < num_ir; i++) { for (j = 0; j < 3; j++) { bz_adrs1[j] = bz_adrs[ir_q1_gps[i]][j]; bz_adrs2[j] = -bz_adrs0[j] - bz_adrs1[j]; } - gp2 = grg_get_grid_index(bz_adrs2, bzgrid->D_diag); + gp2 = recgrid_get_grid_index_from_address(bz_adrs2, bzgrid->D_diag); /* Negative value is the signal to initialize min_d2 later. */ min_d2 = -1; for (bz0 = 0; bz0 < gp_map[num_bzgp + grid_point + 1] - @@ -338,7 +329,7 @@ static void get_BZ_triplets_at_q_type1(long (*triplets)[3], static void get_BZ_triplets_at_q_type2(long (*triplets)[3], const long grid_point, - const ConstBZGrid *bzgrid, + const RecgridConstBZGrid *bzgrid, const long *ir_q1_gps, const long num_ir) { long i, j, gp0, gp2; @@ -353,23 +344,23 @@ static void get_BZ_triplets_at_q_type2(long (*triplets)[3], bz_adrs = bzgrid->addresses; get_LQD_inv(LQD_inv, bzgrid); /* This tolerance is used to be consistent to BZ reduction in bzgrid. */ - tolerance = bzg_get_tolerance_for_BZ_reduction((BZGrid *)bzgrid); + tolerance = recgrid_get_tolerance_for_BZ_reduction((RecgridBZGrid *)bzgrid); for (i = 0; i < 3; i++) { bz_adrs0[i] = bz_adrs[grid_point][i]; } - gp0 = grg_get_grid_index(bz_adrs0, bzgrid->D_diag); + gp0 = recgrid_get_grid_index_from_address(bz_adrs0, bzgrid->D_diag); #ifdef _OPENMP #pragma omp parallel for private(j, gp2, bzgp, G, bz_adrs1, bz_adrs2, d2, \ - min_d2) + min_d2) #endif for (i = 0; i < num_ir; i++) { for (j = 0; j < 3; j++) { bz_adrs1[j] = bz_adrs[gp_map[ir_q1_gps[i]]][j]; bz_adrs2[j] = -bz_adrs0[j] - bz_adrs1[j]; } - gp2 = grg_get_grid_index(bz_adrs2, bzgrid->D_diag); + gp2 = recgrid_get_grid_index_from_address(bz_adrs2, bzgrid->D_diag); /* Negative value is the signal to initialize min_d2 later. */ min_d2 = -1; for (bzgp[0] = gp_map[gp0]; bzgp[0] < gp_map[gp0 + 1]; bzgp[0]++) { @@ -415,7 +406,8 @@ static double get_squared_distance(const long G[3], return d2; } -static void get_LQD_inv(double LQD_inv[3][3], const ConstBZGrid *bzgrid) { +static void get_LQD_inv(double LQD_inv[3][3], + const RecgridConstBZGrid *bzgrid) { long i, j, k; /* LQD^-1 */ @@ -430,19 +422,19 @@ static void get_LQD_inv(double LQD_inv[3][3], const ConstBZGrid *bzgrid) { } /* Return NULL if failed */ -static RotMats *get_reciprocal_point_group_with_q(const RotMats *rot_reciprocal, - const long D_diag[3], - const long grid_point) { - long i, num_rot, gp_rot; +static RecgridMats *get_reciprocal_point_group_with_q( + const RecgridMats *rot_reciprocal, const long D_diag[3], + const long grid_point) { + long i, j, num_rot, gp_rot; long *ir_rot; long adrs[3], adrs_rot[3]; - RotMats *rot_reciprocal_q; + RecgridMats *rot_reciprocal_q; ir_rot = NULL; rot_reciprocal_q = NULL; num_rot = 0; - grg_get_grid_address_from_index(adrs, grid_point, D_diag); + recgrid_get_grid_address_from_index(adrs, grid_point, D_diag); if ((ir_rot = (long *)malloc(sizeof(long) * rot_reciprocal->size)) == NULL) { @@ -454,9 +446,12 @@ static RotMats *get_reciprocal_point_group_with_q(const RotMats *rot_reciprocal, ir_rot[i] = -1; } for (i = 0; i < rot_reciprocal->size; i++) { - lagmat_multiply_matrix_vector_l3(adrs_rot, rot_reciprocal->mat[i], - adrs); - gp_rot = grg_get_grid_index(adrs_rot, D_diag); + for (j = 0; j < 3; j++) { + adrs_rot[j] = rot_reciprocal->mat[i][j][0] * adrs[0] + + rot_reciprocal->mat[i][j][1] * adrs[1] + + rot_reciprocal->mat[i][j][2] * adrs[2]; + } + gp_rot = recgrid_get_grid_index_from_address(adrs_rot, D_diag); if (gp_rot == grid_point) { ir_rot[num_rot] = i; @@ -464,10 +459,10 @@ static RotMats *get_reciprocal_point_group_with_q(const RotMats *rot_reciprocal, } } - if ((rot_reciprocal_q = bzg_alloc_RotMats(num_rot)) != NULL) { + if ((rot_reciprocal_q = recgrid_alloc_RotMats(num_rot)) != NULL) { for (i = 0; i < num_rot; i++) { - lagmat_copy_matrix_l3(rot_reciprocal_q->mat[i], - rot_reciprocal->mat[ir_rot[i]]); + copy_matrix_l3(rot_reciprocal_q->mat[i], + rot_reciprocal->mat[ir_rot[i]]); } } @@ -477,25 +472,36 @@ static RotMats *get_reciprocal_point_group_with_q(const RotMats *rot_reciprocal, return rot_reciprocal_q; } -static RotMats *get_reciprocal_point_group(const long (*rec_rotations_in)[3][3], - const long num_rot, - const long is_time_reversal, - const long is_transpose) { +static RecgridMats *get_reciprocal_point_group( + const long (*rec_rotations_in)[3][3], const long num_rot, + const long is_time_reversal, const long is_transpose) { long i, num_rot_out; long rec_rotations_out[48][3][3]; - RotMats *rec_rotations; + RecgridMats *rec_rotations; - num_rot_out = - grg_get_reciprocal_point_group(rec_rotations_out, rec_rotations_in, - num_rot, is_time_reversal, is_transpose); + num_rot_out = recgrid_get_reciprocal_point_group( + rec_rotations_out, rec_rotations_in, num_rot, is_time_reversal, + is_transpose); if (num_rot_out == 0) { return NULL; } - rec_rotations = bzg_alloc_RotMats(num_rot_out); + rec_rotations = recgrid_alloc_RotMats(num_rot_out); for (i = 0; i < num_rot_out; i++) { - lagmat_copy_matrix_l3(rec_rotations->mat[i], rec_rotations_out[i]); + copy_matrix_l3(rec_rotations->mat[i], rec_rotations_out[i]); } return rec_rotations; } + +static void copy_matrix_l3(long a[3][3], const long b[3][3]) { + a[0][0] = b[0][0]; + a[0][1] = b[0][1]; + a[0][2] = b[0][2]; + a[1][0] = b[1][0]; + a[1][1] = b[1][1]; + a[1][2] = b[1][2]; + a[2][0] = b[2][0]; + a[2][1] = b[2][1]; + a[2][2] = b[2][2]; +} diff --git a/c/triplet_grid.h b/c/triplet_grid.h index 4b22af2c..77a186e6 100644 --- a/c/triplet_grid.h +++ b/c/triplet_grid.h @@ -46,7 +46,7 @@ long tpk_get_ir_triplets_at_q(long *map_triplets, long *map_q, const long (*rec_rotations_in)[3][3], const long num_rot, const long swappable); long tpk_get_BZ_triplets_at_q(long (*triplets)[3], const long grid_point, - const ConstBZGrid *bzgrid, + const RecgridConstBZGrid *bzgrid, const long *map_triplets); #endif diff --git a/c/triplet_iw.c b/c/triplet_iw.c index b49b4e5d..163c6bd9 100644 --- a/c/triplet_iw.c +++ b/c/triplet_iw.c @@ -36,10 +36,9 @@ #include +#include "funcs.h" #include "grgrid.h" -#include "phonoc_utils.h" #include "tetrahedron_method.h" -#include "triplet.h" static void set_freq_vertices(double freq_vertices[3][24][4], const double *frequencies1, @@ -51,23 +50,23 @@ static long set_g(double g[3], const double f0, const double freq_vertices[3][24][4], const long max_i); static void get_triplet_tetrahedra_vertices( long vertices[2][24][4], const long tp_relative_grid_address[2][24][4][3], - const long triplet[3], const ConstBZGrid *bzgrid); + const long triplet[3], const RecgridConstBZGrid *bzgrid); static void get_neighboring_grid_points_type1( long *neighboring_grid_points, const long grid_point, const long (*relative_grid_address)[3], - const long num_relative_grid_address, const ConstBZGrid *bzgrid); + const long num_relative_grid_address, const RecgridConstBZGrid *bzgrid); static void get_neighboring_grid_points_type2( long *neighboring_grid_points, const long grid_point, const long (*relative_grid_address)[3], - const long num_relative_grid_address, const ConstBZGrid *bzgrid); + const long num_relative_grid_address, const RecgridConstBZGrid *bzgrid); void tpi_get_integration_weight( double *iw, char *iw_zero, const double *frequency_points, const long num_band0, const long tp_relative_grid_address[2][24][4][3], - const long triplets[3], const long num_triplets, const ConstBZGrid *bzgrid, - const double *frequencies1, const long num_band1, - const double *frequencies2, const long num_band2, const long tp_type, - const long openmp_per_triplets) { + const long triplets[3], const long num_triplets, + const RecgridConstBZGrid *bzgrid, const double *frequencies1, + const long num_band1, const double *frequencies2, const long num_band2, + const long tp_type, const long openmp_per_triplets) { long max_i, j, b1, b2, b12, num_band_prod, adrs_shift; long vertices[2][24][4]; double g[3]; @@ -159,9 +158,9 @@ void tpi_get_integration_weight_with_sigma( g2 = 0; } else { iw_zero[adrs_shift] = 0; - g0 = phonoc_gaussian(f0 + f1 - f2, sigma); - g1 = phonoc_gaussian(f0 - f1 + f2, sigma); - g2 = phonoc_gaussian(f0 - f1 - f2, sigma); + g0 = funcs_gaussian(f0 + f1 - f2, sigma); + g1 = funcs_gaussian(f0 - f1 + f2, sigma); + g2 = funcs_gaussian(f0 - f1 - f2, sigma); } if (tp_type == 2) { iw[adrs_shift] = g2; @@ -182,7 +181,7 @@ void tpi_get_integration_weight_with_sigma( iw[adrs_shift] = 0; } else { iw_zero[adrs_shift] = 0; - iw[adrs_shift] = phonoc_gaussian(f0 + f1 - f2, sigma); + iw[adrs_shift] = funcs_gaussian(f0 + f1 - f2, sigma); } } } @@ -204,7 +203,7 @@ void tpi_get_neighboring_grid_points(long *neighboring_grid_points, const long grid_point, const long (*relative_grid_address)[3], const long num_relative_grid_address, - const ConstBZGrid *bzgrid) { + const RecgridConstBZGrid *bzgrid) { if (bzgrid->type == 1) { get_neighboring_grid_points_type1(neighboring_grid_points, grid_point, relative_grid_address, @@ -276,7 +275,7 @@ static long set_g(double g[3], const double f0, static void get_triplet_tetrahedra_vertices( long vertices[2][24][4], const long tp_relative_grid_address[2][24][4][3], - const long triplet[3], const ConstBZGrid *bzgrid) { + const long triplet[3], const RecgridConstBZGrid *bzgrid) { long i, j; for (i = 0; i < 2; i++) { @@ -291,7 +290,7 @@ static void get_triplet_tetrahedra_vertices( static void get_neighboring_grid_points_type1( long *neighboring_grid_points, const long grid_point, const long (*relative_grid_address)[3], - const long num_relative_grid_address, const ConstBZGrid *bzgrid) { + const long num_relative_grid_address, const RecgridConstBZGrid *bzgrid) { long bzmesh[3], bz_address[3]; long i, j, bz_gp, prod_bz_mesh; @@ -304,10 +303,11 @@ static void get_neighboring_grid_points_type1( bz_address[j] = bzgrid->addresses[grid_point][j] + relative_grid_address[i][j]; } - bz_gp = bzgrid->gp_map[grg_get_grid_index(bz_address, bzmesh)]; + bz_gp = bzgrid->gp_map[recgrid_get_grid_index_from_address(bz_address, + bzmesh)]; if (bz_gp == prod_bz_mesh) { neighboring_grid_points[i] = - grg_get_grid_index(bz_address, bzgrid->D_diag); + recgrid_get_grid_index_from_address(bz_address, bzgrid->D_diag); } else { neighboring_grid_points[i] = bz_gp; } @@ -317,7 +317,7 @@ static void get_neighboring_grid_points_type1( static void get_neighboring_grid_points_type2( long *neighboring_grid_points, const long grid_point, const long (*relative_grid_address)[3], - const long num_relative_grid_address, const ConstBZGrid *bzgrid) { + const long num_relative_grid_address, const RecgridConstBZGrid *bzgrid) { long bz_address[3]; long i, j, gp; @@ -326,7 +326,7 @@ static void get_neighboring_grid_points_type2( bz_address[j] = bzgrid->addresses[grid_point][j] + relative_grid_address[i][j]; } - gp = grg_get_grid_index(bz_address, bzgrid->D_diag); + gp = recgrid_get_grid_index_from_address(bz_address, bzgrid->D_diag); neighboring_grid_points[i] = bzgrid->gp_map[gp]; if (bzgrid->gp_map[gp + 1] - bzgrid->gp_map[gp] > 1) { for (j = bzgrid->gp_map[gp]; j < bzgrid->gp_map[gp + 1]; j++) { diff --git a/c/triplet_iw.h b/c/triplet_iw.h index aacdae48..75468766 100644 --- a/c/triplet_iw.h +++ b/c/triplet_iw.h @@ -40,10 +40,10 @@ void tpi_get_integration_weight( double *iw, char *iw_zero, const double *frequency_points, const long num_band0, const long tp_relative_grid_address[2][24][4][3], - const long triplets[3], const long num_triplets, const ConstBZGrid *bzgrid, - const double *frequencies1, const long num_band1, - const double *frequencies2, const long num_band2, const long tp_type, - const long openmp_per_triplets); + const long triplets[3], const long num_triplets, + const RecgridConstBZGrid *bzgrid, const double *frequencies1, + const long num_band1, const double *frequencies2, const long num_band2, + const long tp_type, const long openmp_per_triplets); void tpi_get_integration_weight_with_sigma( double *iw, char *iw_zero, const double sigma, const double cutoff, const double *frequency_points, const long num_band0, const long triplet[3], @@ -53,6 +53,6 @@ void tpi_get_neighboring_grid_points(long *neighboring_grid_points, const long grid_point, const long (*relative_grid_address)[3], const long num_relative_grid_address, - const ConstBZGrid *bzgrid); + const RecgridConstBZGrid *bzgrid); #endif diff --git a/doc/changelog.md b/doc/changelog.md index 77c32354..fe894aa2 100644 --- a/doc/changelog.md +++ b/doc/changelog.md @@ -2,6 +2,10 @@ # Change Log +## Jul-8-2024: Version 3.3.1 + +- Major refactoring to isolate reciprocal space grid code. + ## Jul-8-2024: Version 3.3.0 - Build system of phono3py was renewed. Now nanobind, cmake, and diff --git a/doc/conf.py b/doc/conf.py index ab73c002..b7536d4a 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -60,7 +60,7 @@ # The short X.Y version. version = "3.3" # The full version, including alpha/beta/rc tags. -release = "3.3.0" +release = "3.3.1" # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. diff --git a/phono3py/api_phono3py.py b/phono3py/api_phono3py.py index 2be91174..b26cc2a5 100644 --- a/phono3py/api_phono3py.py +++ b/phono3py/api_phono3py.py @@ -1908,35 +1908,36 @@ def run_spectral_function( def run_thermal_conductivity( self, - is_LBTE=False, - temperatures=None, - is_isotope=False, - mass_variances=None, - grid_points=None, - boundary_mfp=None, # in micrometre - solve_collective_phonon=False, - use_ave_pp=False, - is_reducible_collision_matrix=False, - is_kappa_star=True, - gv_delta_q=None, # for group velocity - is_full_pp=False, - pinv_cutoff=1.0e-8, # for pseudo-inversion of collision matrix - pinv_method=0, # for pseudo-inversion of collision matrix - pinv_solver=0, # solver of pseudo-inversion of collision matrix - write_gamma=False, - read_gamma=False, - is_N_U=False, - conductivity_type=None, - write_kappa=False, - write_gamma_detail=False, - write_collision=False, - read_collision=False, - write_pp=False, - read_pp=False, - write_LBTE_solution=False, - compression="gzip", - input_filename=None, - output_filename=None, + is_LBTE: bool = False, + temperatures: Optional[Sequence] = None, + is_isotope: bool = False, + mass_variances: Optional[Sequence] = None, + grid_points: Optional[Sequence[int]] = None, + boundary_mfp: Optional[float] = None, # in micrometre + solve_collective_phonon: bool = False, + use_ave_pp: bool = False, + is_reducible_collision_matrix: bool = False, + is_kappa_star: bool = True, + gv_delta_q: Optional[float] = None, # for group velocity + is_full_pp: bool = False, + pinv_cutoff: float = 1.0e-8, # for pseudo-inversion of collision matrix + pinv_method: int = 0, # for pseudo-inversion of collision matrix + pinv_solver: int = 0, # solver of pseudo-inversion of collision matrix + write_gamma: bool = False, + read_gamma: bool = False, + is_N_U: bool = False, + conductivity_type: Optional[str] = None, + write_kappa: bool = False, + write_gamma_detail: bool = False, + write_collision: bool = False, + read_collision: bool = False, + write_pp: bool = False, + read_pp: bool = False, + write_LBTE_solution: bool = False, + compression: str = "gzip", + input_filename: Optional[str] = None, + output_filename: Optional[str] = None, + log_level: Optional[int] = None, ): """Run thermal conductivity calculation. @@ -2076,6 +2077,11 @@ def run_thermal_conductivity( ) raise RuntimeError(msg) + if log_level is None: + _log_level = self._log_level + else: + _log_level = log_level + if is_LBTE: if temperatures is None: _temperatures = [ @@ -2110,7 +2116,7 @@ def run_thermal_conductivity( compression=compression, input_filename=input_filename, output_filename=output_filename, - log_level=self._log_level, + log_level=_log_level, ) else: if temperatures is None: @@ -2141,7 +2147,7 @@ def run_thermal_conductivity( compression=compression, input_filename=input_filename, output_filename=output_filename, - log_level=self._log_level, + log_level=_log_level, ) def save(self, filename="phono3py_params.yaml", settings=None): diff --git a/phono3py/conductivity/direct_solution.py b/phono3py/conductivity/direct_solution.py index 0b292431..8477e089 100644 --- a/phono3py/conductivity/direct_solution.py +++ b/phono3py/conductivity/direct_solution.py @@ -1329,8 +1329,7 @@ def _set_kappa_at_sigmas(self, weights): text += "for sigma=%s -----------" % sigma else: text += "with tetrahedron method -----------" - print(text) - sys.stdout.flush() + print(text, flush=True) for k, t in enumerate(self._temperatures): if t > 0: @@ -1360,11 +1359,10 @@ def _set_kappa_at_sigmas(self, weights): (" %6s " + " %10.3f" * 6) % (("(RTA)",) + tuple(self._kappa_RTA[j, k])) ) - print("-" * 76) - sys.stdout.flush() + print("-" * 76, flush=True) if self._log_level: - print("") + print("", flush=True) def _set_kappa(self, i_sigma, i_temp, weights): if self._is_reducible_collision_matrix: @@ -1532,8 +1530,7 @@ def _set_kappa_at_sigmas(self, weights): text += "for sigma=%s -----------" % sigma else: text += "with tetrahedron method -----------" - print(text) - sys.stdout.flush() + print(text, flush=True) for k, t in enumerate(self._temperatures): if t > 0: @@ -1587,11 +1584,10 @@ def _set_kappa_at_sigmas(self, weights): + tuple(self._kappa_C[j, k] + self._kappa_P_exact[j, k]) ) ) - print("-" * 76) - sys.stdout.flush() + print("-" * 76, flush=True) if self._log_level: - print("") + print("", flush=True) def _set_kappa(self, i_sigma, i_temp, weights): if self._is_reducible_collision_matrix: diff --git a/phono3py/conductivity/rta.py b/phono3py/conductivity/rta.py index 44a786e8..ac26f057 100644 --- a/phono3py/conductivity/rta.py +++ b/phono3py/conductivity/rta.py @@ -34,6 +34,7 @@ # ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # POSSIBILITY OF SUCH DAMAGE. +import warnings from abc import abstractmethod from typing import Type, Union, cast @@ -141,6 +142,16 @@ def get_gamma_detail_at_q(self): return self._gamma_detail_at_q def get_number_of_ignored_phonon_modes(self): + """Return number of ignored phonon modes.""" + warnings.warn( + "Use attribute, number_of_ignored_phonon_modes", + DeprecationWarning, + stacklevel=2, + ) + return self.number_of_ignored_phonon_modes + + @property + def number_of_ignored_phonon_modes(self): """Return number of ignored phonon modes.""" return self._num_ignored_phonon_modes diff --git a/phono3py/conductivity/utils.py b/phono3py/conductivity/utils.py index dcc4c63b..ea5797bc 100644 --- a/phono3py/conductivity/utils.py +++ b/phono3py/conductivity/utils.py @@ -841,7 +841,7 @@ def kappa_RTA(br: "ConductivityRTA", log_level): temperatures = br.temperatures sigmas = br.sigmas kappa = br.kappa - num_ignored_phonon_modes = br.get_number_of_ignored_phonon_modes() + num_ignored_phonon_modes = br.number_of_ignored_phonon_modes num_band = br.frequencies.shape[1] num_phonon_modes = br.get_number_of_sampling_grid_points() * num_band for i, sigma in enumerate(sigmas): @@ -872,7 +872,7 @@ def kappa_RTA(br: "ConductivityRTA", log_level): ) for t, k in zip(temperatures, kappa[i]): print(("%7.1f " + " %10.3f" * 6) % ((t,) + tuple(k))) - print("") + print("", flush=True) @staticmethod def kappa_Wigner_RTA(br: "ConductivityWignerRTA", log_level): @@ -882,7 +882,7 @@ def kappa_Wigner_RTA(br: "ConductivityWignerRTA", log_level): kappa_TOT_RTA = br.kappa_TOT_RTA kappa_P_RTA = br.kappa_P_RTA kappa_C = br.kappa_C - num_ignored_phonon_modes = br.get_number_of_ignored_phonon_modes() + num_ignored_phonon_modes = br.number_of_ignored_phonon_modes num_band = br.frequencies.shape[1] num_phonon_modes = br.get_number_of_sampling_grid_points() * num_band for i, sigma in enumerate(sigmas): @@ -943,7 +943,7 @@ def kappa_Wigner_RTA(br: "ConductivityWignerRTA", log_level): print(" ") for t, k in zip(temperatures, kappa_TOT_RTA[i]): print("K_T\t" + ("%7.1f " + " %10.3f" * 6) % ((t,) + tuple(k))) - print("") + print("", flush=True) def write_pp( diff --git a/phono3py/file_IO.py b/phono3py/file_IO.py index 8875fb35..e43dfa7f 100644 --- a/phono3py/file_IO.py +++ b/phono3py/file_IO.py @@ -308,7 +308,7 @@ def read_fc3_from_hdf5(filename="fc3.hdf5", p2s_map=None): check_force_constants_indices( fc3.shape[:2], p2s_map_in_file, p2s_map, filename ) - if fc3.dtype == np.double and fc3.flags.c_contiguous: + if fc3.dtype == np.dtype("double") and fc3.flags.c_contiguous: return fc3 else: msg = ( @@ -316,7 +316,6 @@ def read_fc3_from_hdf5(filename="fc3.hdf5", p2s_map=None): "dtype='double' and c_contiguous." % filename ) raise TypeError(msg) - return None def write_fc2_to_hdf5( diff --git a/phono3py/phonon/grid.py b/phono3py/phonon/grid.py index b42112ea..a3407b96 100644 --- a/phono3py/phonon/grid.py +++ b/phono3py/phonon/grid.py @@ -466,7 +466,7 @@ def _set_rotations(self): Terminate when symmetry of grid is broken. """ - import phono3py._phono3py as phono3c + import phono3py._recgrid as recgrid if self._symmetry_dataset is None: direct_rotations = np.eye(3, dtype="int_", order="C").reshape(1, 3, 3) @@ -475,7 +475,7 @@ def _set_rotations(self): self._symmetry_dataset["rotations"], dtype="int_", order="C" ) rec_rotations = np.zeros((48, 3, 3), dtype="int_", order="C") - num_rec_rot = phono3c.reciprocal_rotations( + num_rec_rot = recgrid.reciprocal_rotations( rec_rotations, direct_rotations, self._is_time_reversal ) self._reciprocal_operations = np.array( @@ -492,7 +492,7 @@ def _set_rotations(self): self._rotations = np.zeros( self._reciprocal_operations.shape, dtype="int_", order="C" ) - if not phono3c.transform_rotations( + if not recgrid.transform_rotations( self._rotations, self._reciprocal_operations, self._D_diag, self._Q ): msg = "Grid symmetry is broken. Use generalized regular grid." @@ -778,9 +778,9 @@ def _set_GRG_mesh( if (np.diag(gm_diag) == _grid_matrix).all() and not force_SNF: self._D_diag = np.array(gm_diag, dtype="int_") else: - import phono3py._phono3py as phono3c + import phono3py._recgrid as recgrid - if not phono3c.snf3x3(self._D_diag, self._P, self._Q, _grid_matrix): + if not recgrid.snf3x3(self._D_diag, self._P, self._Q, _grid_matrix): msg = "SNF3x3 failed." raise RuntimeError(msg) @@ -879,17 +879,17 @@ def get_grid_point_from_address(address, D_diag): shape=(n, ), dtype='int_' """ - import phono3py._phono3py as phono3c + import phono3py._recgrid as recgrid adrs_array = np.array(address, dtype="int_", order="C") mesh_array = np.array(D_diag, dtype="int_") if adrs_array.ndim == 1: - return phono3c.grid_index_from_address(adrs_array, mesh_array) + return recgrid.grid_index_from_address(adrs_array, mesh_array) gps = np.zeros(adrs_array.shape[0], dtype="int_") for i, adrs in enumerate(adrs_array): - gps[i] = phono3c.grid_index_from_address(adrs, mesh_array) + gps[i] = recgrid.grid_index_from_address(adrs, mesh_array) return gps @@ -977,11 +977,11 @@ def _get_grid_points_by_bz_rotations(bz_gp, bz_grid: BZGrid, rotations, lang="C" def _get_grid_points_by_bz_rotations_c(bz_gp, bz_grid: BZGrid, rotations): - import phono3py._phono3py as phono3c + import phono3py._recgrid as recgrid bzgps = np.zeros(len(rotations), dtype="int_") for i, r in enumerate(rotations): - bzgps[i] = phono3c.rotate_bz_grid_index( + bzgps[i] = recgrid.rotate_bz_grid_index( bz_gp, r, bz_grid.addresses, @@ -1053,10 +1053,10 @@ def _get_grid_address(D_diag): shape=(prod(D_diag), 3), dtype='int_' """ - import phono3py._phono3py as phono3c + import phono3py._recgrid as recgrid gr_grid_addresses = np.zeros((np.prod(D_diag), 3), dtype="int_") - phono3c.gr_grid_addresses(gr_grid_addresses, np.array(D_diag, dtype="int_")) + recgrid.gr_grid_addresses(gr_grid_addresses, np.array(D_diag, dtype="int_")) return gr_grid_addresses @@ -1122,7 +1122,7 @@ def _relocate_BZ_grid_address( shape=(prod(mesh) + 1, ) """ - import phono3py._phono3py as phono3c + import phono3py._recgrid as recgrid if PS is None: _PS = np.zeros(3, dtype="int_") @@ -1143,7 +1143,7 @@ def _relocate_BZ_grid_address( tmat_inv_int = np.rint(tmat_inv).astype("int_") assert (np.abs(tmat_inv - tmat_inv_int) < 1e-5).all() - num_gp = phono3c.bz_grid_addresses( + num_gp = recgrid.bz_grid_addresses( bz_grid_addresses, bz_map, bzg2grg, @@ -1183,7 +1183,7 @@ def _get_ir_grid_map(D_diag, grg_rotations, PS=None): dtype='int_', shape=(prod(mesh),) """ - import phono3py._phono3py as phono3c + import phono3py._recgrid as recgrid ir_grid_map = np.zeros(np.prod(D_diag), dtype="int_") if PS is None: @@ -1191,7 +1191,7 @@ def _get_ir_grid_map(D_diag, grg_rotations, PS=None): else: _PS = np.array(PS, dtype="int_") - num_ir = phono3c.ir_grid_map( + num_ir = recgrid.ir_grid_map( ir_grid_map, np.array(D_diag, dtype="int_"), _PS, diff --git a/phono3py/phonon/solver.py b/phono3py/phonon/solver.py index ff2a76db..a694c5bd 100644 --- a/phono3py/phonon/solver.py +++ b/phono3py/phonon/solver.py @@ -118,9 +118,9 @@ def run_phonon_solver_c( is_nac_q_zero = True _nac_q_direction = np.array(nac_q_direction, dtype="double") - assert grid_points.dtype == "int_" + assert grid_points.dtype == np.dtype("int_") assert grid_points.flags.c_contiguous - assert QDinv.dtype == "double" + assert QDinv.dtype == np.dtype("double") assert QDinv.flags.c_contiguous assert lapack_zheev_uplo in ("L", "U") diff --git a/phono3py/version.py b/phono3py/version.py index 87c2d2f1..dcc4a40d 100644 --- a/phono3py/version.py +++ b/phono3py/version.py @@ -34,4 +34,4 @@ # ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # POSSIBILITY OF SUCH DAMAGE. -__version__ = "3.3.0" +__version__ = "3.3.1" diff --git a/test/phonon/test_velocity_operator.py b/test/phonon/test_velocity_operator.py index 320f2dd7..d77de804 100644 --- a/test/phonon/test_velocity_operator.py +++ b/test/phonon/test_velocity_operator.py @@ -401,7 +401,7 @@ def test_gv_operator_si(ph_si: Phonopy): ph_si.dynamical_matrix, symmetry=ph_si.primitive_symmetry ) - ph_si.dynamical_matrix.run([[0.1, 0.22, 0.33]]) + ph_si.dynamical_matrix.run([0.1, 0.22, 0.33]) dm = ph_si.dynamical_matrix.dynamical_matrix eigvals, eigvecs = np.linalg.eigh(dm)