From d252e4429d8a7a2b492ab400c85b18d860d7404d Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Sun, 14 Jul 2024 12:51:13 +0900 Subject: [PATCH 01/15] Isolate recprocal space grid code --- CMakeLists.txt | 65 +++++++-- c/_phono3py.cpp | 8 -- c/_recgrid.cpp | 179 +++++++++++++++++++++++++ c/bzgrid.c | 121 +++++------------ c/bzgrid.h | 65 +-------- c/gridsys.c | 20 +-- c/gridsys.h | 2 +- c/interaction.c | 25 +++- c/interaction.h | 15 ++- c/phono3py.c | 158 +++------------------- c/pp_collision.c | 37 ++--- c/pp_collision.h | 24 ++-- c/recgrid.c | 248 ++++++++++++++++++++++++++++++++++ c/recgrid.h | 289 ++++++++++++++++++++++++++++++++++++++++ c/triplet.c | 10 +- c/triplet.h | 10 +- c/triplet_grid.c | 133 +++++++++--------- c/triplet_grid.h | 2 +- c/triplet_iw.c | 29 ++-- c/triplet_iw.h | 10 +- phono3py/phonon/grid.py | 32 ++--- 21 files changed, 1017 insertions(+), 465 deletions(-) create mode 100644 c/_recgrid.cpp create mode 100644 c/recgrid.c create mode 100644 c/recgrid.h diff --git a/CMakeLists.txt b/CMakeLists.txt index fc372955..2cc7bf39 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,19 +126,54 @@ 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 @@ -145,7 +181,6 @@ if(BUILD_PHPHCALC_LIB OR BUILD_NANOBIND_MODULE) ${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 @@ -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/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/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/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..cdd6d268 100644 --- a/c/pp_collision.c +++ b/c/pp_collision.c @@ -43,6 +43,7 @@ #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/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..d9370fce 100644 --- a/c/triplet_grid.c +++ b/c/triplet_grid.c @@ -46,7 +46,7 @@ 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 +56,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 +84,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 +94,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 +124,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 +135,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 +150,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 +209,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 +246,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 +261,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 +271,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 +332,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 +347,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 +409,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 +425,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 +449,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 +462,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 +475,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..3df801f8 100644 --- a/c/triplet_iw.c +++ b/c/triplet_iw.c @@ -51,23 +51,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]; @@ -204,7 +204,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 +276,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 +291,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 +304,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 +318,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 +327,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/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, From 78c1e1d70bcc523836cacf6492667a82466265f0 Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Sun, 14 Jul 2024 15:13:11 +0900 Subject: [PATCH 02/15] Rename phonoc_utils -> funcs --- CMakeLists.txt | 4 ++-- c/collision_matrix.c | 4 ++-- c/{phonoc_utils.c => funcs.c} | 8 ++++---- c/{phonoc_utils.h => funcs.h} | 10 +++++----- c/imag_self_energy_with_g.c | 6 +++--- c/isotope.c | 10 +++++----- c/pp_collision.c | 2 +- c/real_self_energy.c | 6 +++--- c/triplet_grid.c | 3 --- c/triplet_iw.c | 11 +++++------ 10 files changed, 30 insertions(+), 34 deletions(-) rename c/{phonoc_utils.c => funcs.c} (90%) rename c/{phonoc_utils.h => funcs.h} (87%) diff --git a/CMakeLists.txt b/CMakeLists.txt index 2cc7bf39..a4d236e2 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -176,7 +176,7 @@ if(BUILD_PHPHCALC_LIB OR BUILD_NANOBIND_MODULE) ${PROJECT_SOURCE_DIR}/c/isotope.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 @@ -358,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 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/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/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/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/pp_collision.c b/c/pp_collision.c index cdd6d268..d083a6c8 100644 --- a/c/pp_collision.c +++ b/c/pp_collision.c @@ -37,11 +37,11 @@ #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" 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/triplet_grid.c b/c/triplet_grid.c index d9370fce..e5cd980e 100644 --- a/c/triplet_grid.c +++ b/c/triplet_grid.c @@ -39,10 +39,7 @@ #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], diff --git a/c/triplet_iw.c b/c/triplet_iw.c index 3df801f8..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, @@ -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); } } } From 399b50555d609638b7efc85685c362dae6fd00b6 Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Sun, 14 Jul 2024 15:32:03 +0900 Subject: [PATCH 03/15] Use same dynmat.c as phonopy --- c/dynmat.c | 39 +++++++++++++++++++-------------------- 1 file changed, 19 insertions(+), 20 deletions(-) diff --git a/c/dynmat.c b/c/dynmat.c index e179a1a6..6bf955ad 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, @@ -244,6 +241,7 @@ long dym_get_dynamical_matrices_openmp_over_qpoints( 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 +489,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 +600,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; From b15c337458de45174477bc699f8c07f03e9fbd47 Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Sun, 14 Jul 2024 17:13:07 +0900 Subject: [PATCH 04/15] Use same dynmat.c as phonopy --- c/dynmat.c | 25 ------------------------- c/dynmat.h | 5 ----- test/phonon/test_velocity_operator.py | 2 +- 3 files changed, 1 insertion(+), 31 deletions(-) diff --git a/c/dynmat.c b/c/dynmat.c index 6bf955ad..6be2ba1e 100644 --- a/c/dynmat.c +++ b/c/dynmat.c @@ -193,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, @@ -218,29 +216,6 @@ 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, 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/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) From 75a04ae23b03a61a626fc0171a7e5573785008eb Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Mon, 15 Jul 2024 08:34:39 +0900 Subject: [PATCH 05/15] Set version 3.3.1 --- doc/changelog.md | 4 ++++ doc/conf.py | 2 +- phono3py/version.py | 2 +- 3 files changed, 6 insertions(+), 2 deletions(-) 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/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" From 3ad372de10696d3e701db3284209f1ff6323e7de Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 15 Jul 2024 23:07:29 +0000 Subject: [PATCH 06/15] [pre-commit.ci] pre-commit autoupdate MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit updates: - [github.com/astral-sh/ruff-pre-commit: v0.5.1 → v0.5.2](https://github.com/astral-sh/ruff-pre-commit/compare/v0.5.1...v0.5.2) --- .pre-commit-config.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index e89949ff..293ee4b8 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.1 + rev: v0.5.2 hooks: - id: ruff args: [ "--fix", "--show-fixes" ] From 682ac9d9c7b9980e02de08ee37f6de567f971bec Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Tue, 16 Jul 2024 18:12:50 +0900 Subject: [PATCH 07/15] Fix showing fc2 filename --- phono3py/cui/load.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/phono3py/cui/load.py b/phono3py/cui/load.py index 14730c85..8c33403e 100644 --- a/phono3py/cui/load.py +++ b/phono3py/cui/load.py @@ -558,7 +558,7 @@ def _get_dataset_phonon_dataset_or_fc2( ph3py.fc2 = fc2 read_fc2 = True if log_level: - print(f'fc2 was read from "{fc2_filename}".') + print(f'fc2 was read from "{_fc2_filename}".') elif ( ph3py_yaml is not None and ph3py_yaml.phonon_dataset is not None From b4019dd312fd2fe54d990543e4eb6a31729bdeae Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Sun, 21 Jul 2024 15:55:34 +0900 Subject: [PATCH 08/15] Fix load dataset wo forces --- phono3py/api_phono3py.py | 17 +++--- phono3py/cui/create_force_constants.py | 72 ++++++++++++++++++++------ phono3py/cui/load.py | 51 +++++++++++++----- phono3py/cui/phono3py_script.py | 5 +- 4 files changed, 108 insertions(+), 37 deletions(-) diff --git a/phono3py/api_phono3py.py b/phono3py/api_phono3py.py index b26cc2a5..51c5b592 100644 --- a/phono3py/api_phono3py.py +++ b/phono3py/api_phono3py.py @@ -846,6 +846,9 @@ def displacements(self): """ dataset = self._dataset + if self._dataset is None: + raise RuntimeError("displacement dataset is not set.") + if "first_atoms" in dataset: num_scells = len(dataset["first_atoms"]) for disp1 in dataset["first_atoms"]: @@ -1218,10 +1221,10 @@ def run_phonon_solver(self, grid_points=None): def generate_displacements( self, - distance=0.03, - cutoff_pair_distance=None, - is_plusminus="auto", - is_diagonal=True, + distance: float = 0.03, + cutoff_pair_distance: Optional[float] = None, + is_plusminus: Union[bool, str] = "auto", + is_diagonal: bool = True, number_of_snapshots: Optional[int] = None, random_seed: Optional[int] = None, is_random_distance: bool = False, @@ -1410,7 +1413,7 @@ def produce_fc3( self, symmetrize_fc3r: bool = False, is_compact_fc: bool = False, - fc_calculator: Optional[str] = None, + fc_calculator: Optional[Union[str, dict]] = None, fc_calculator_options: Optional[Union[str, dict]] = None, ): """Calculate fc3 from displacements and forces. @@ -2150,7 +2153,9 @@ def run_thermal_conductivity( log_level=_log_level, ) - def save(self, filename="phono3py_params.yaml", settings=None): + def save( + self, filename: str = "phono3py_params.yaml", settings: Optional[dict] = None + ): """Save parameters in Phono3py instants into file. Parameters diff --git a/phono3py/cui/create_force_constants.py b/phono3py/cui/create_force_constants.py index 9f6f778a..d0fb4d44 100644 --- a/phono3py/cui/create_force_constants.py +++ b/phono3py/cui/create_force_constants.py @@ -95,8 +95,6 @@ def create_phono3py_force_constants( symmetrize_fc3r = settings.is_symmetrize_fc3_r or settings.fc_symmetry symmetrize_fc2 = settings.is_symmetrize_fc2 or settings.fc_symmetry - (fc_calculator, fc_calculator_options) = get_fc_calculator_params(settings) - if log_level: show_phono3py_force_constants_settings(settings) @@ -113,6 +111,9 @@ def create_phono3py_force_constants( ): pass else: + (fc_calculator, fc_calculator_options) = get_fc_calculator_params( + settings, log_level=(not settings.read_fc3) * 1 + ) if settings.read_fc3: _read_phono3py_fc3(phono3py, symmetrize_fc3r, input_filename, log_level) else: # fc3 from FORCES_FC3 or ph3py_yaml @@ -292,10 +293,17 @@ def parse_forces( # Type-1 FORCES_FC*. # dataset comes either from disp_fc*.yaml or phono3py*.yaml. if not forces_in_dataset(dataset): - if fc_type == "phonon_fc2": - parse_FORCES_FC2(dataset, filename=force_filename) - else: - parse_FORCES_FC3(dataset, filename=force_filename) + if force_filename is not None: + if fc_type == "phonon_fc2": + parse_FORCES_FC2(dataset, filename=force_filename) + else: + parse_FORCES_FC3(dataset, filename=force_filename) + + if log_level: + print( + f'Sets of supercell forces were read from "{force_filename}".', + flush=True, + ) # Unit of displacements is already converted. # Therefore, only unit of forces is converted. @@ -305,15 +313,13 @@ def parse_forces( force_to_eVperA=physical_units["force_to_eVperA"], ) - if log_level: - print('Sets of supercell forces were read from "%s".' % force_filename) - sys.stdout.flush() - return dataset def forces_in_dataset(dataset: dict) -> bool: """Return whether forces in dataset or not.""" + if dataset is None: + return False return "forces" in dataset or ( "first_atoms" in dataset and "forces" in dataset["first_atoms"][0] ) @@ -326,7 +332,7 @@ def displacements_in_dataset(dataset: Optional[dict]) -> bool: return "displacements" in dataset or "first_atoms" in dataset -def get_fc_calculator_params(settings): +def get_fc_calculator_params(settings, log_level=0): """Return fc_calculator and fc_calculator_params from settings.""" fc_calculator = None fc_calculator_list = [] @@ -339,13 +345,49 @@ def get_fc_calculator_params(settings): if fc_calculator_list: fc_calculator = "|".join(fc_calculator_list) - fc_calculator_options = None - if settings.fc_calculator_options is not None: - fc_calculator_options = settings.fc_calculator_options + fc_calculator_options = settings.fc_calculator_options + if settings.cutoff_pair_distance: + if fc_calculator_list and fc_calculator_list[-1] in ("alm", "symfc"): + if fc_calculator_list[-1] == "alm": + cutoff_str = f"-1 {settings.cutoff_pair_distance}" + if fc_calculator_list[-1] == "symfc": + cutoff_str = f"{settings.cutoff_pair_distance}" + fc_calculator_options = _set_cutoff_in_fc_calculator_options( + fc_calculator_options, + cutoff_str, + log_level, + ) return fc_calculator, fc_calculator_options +def _set_cutoff_in_fc_calculator_options( + fc_calculator_options: Optional[str], + cutoff_str: str, + log_level: int, +): + str_appended = f"cutoff={cutoff_str}" + calc_opts = fc_calculator_options + if calc_opts is None: + calc_opts = "|" + if "|" in calc_opts: + calc_opts_fc2, calc_opts_fc3 = [v.strip() for v in calc_opts.split("|")][:2] + else: + calc_opts_fc2 = calc_opts + calc_opts_fc3 = calc_opts + + if calc_opts_fc3 == "": + calc_opts_fc3 += f"{str_appended}" + if log_level: + print(f'Set "{str_appended}" to fc_calculator_options for fc3.') + elif "cutoff" not in calc_opts_fc3: + calc_opts_fc3 += f", {str_appended}" + if log_level: + print(f'Appended "{str_appended}" to fc_calculator_options for fc3.') + + return f"{calc_opts_fc2}|{calc_opts_fc3}" + + def _read_phono3py_fc3(phono3py: Phono3py, symmetrize_fc3r, input_filename, log_level): if input_filename is None: filename = "fc3.hdf5" @@ -417,7 +459,7 @@ def _read_phono3py_fc2(phono3py, symmetrize_fc2, input_filename, log_level): def read_type2_dataset(natom, filename="FORCES_FC3", log_level=0) -> Optional[dict]: """Read type-2 FORCES_FC3.""" - if not pathlib.Path(filename).exists(): + if filename is None or not pathlib.Path(filename).exists(): return None with open(filename, "r") as f: diff --git a/phono3py/cui/load.py b/phono3py/cui/load.py index 8c33403e..12fb9cbf 100644 --- a/phono3py/cui/load.py +++ b/phono3py/cui/load.py @@ -453,7 +453,7 @@ def compute_force_constants_from_datasets( fc3_calculator = extract_fc2_fc3_calculators(fc_calculator, 3) fc2_calculator = extract_fc2_fc3_calculators(fc_calculator, 2) if not read_fc["fc3"] and (ph3py.dataset or ph3py.mlp_dataset): - if ph3py.mlp_dataset and use_pypolymlp: + if use_pypolymlp and forces_in_dataset(ph3py.mlp_dataset): run_pypolymlp_to_compute_forces( ph3py, mlp_params=mlp_params, @@ -462,17 +462,22 @@ def compute_force_constants_from_datasets( random_seed=random_seed, log_level=log_level, ) - ph3py.produce_fc3( - symmetrize_fc3r=symmetrize_fc, - is_compact_fc=is_compact_fc, - fc_calculator=fc3_calculator, - fc_calculator_options=extract_fc2_fc3_calculators(fc_calculator_options, 3), - ) + if forces_in_dataset(ph3py.dataset): + ph3py.produce_fc3( + symmetrize_fc3r=symmetrize_fc, + is_compact_fc=is_compact_fc, + fc_calculator=fc3_calculator, + fc_calculator_options=extract_fc2_fc3_calculators( + fc_calculator_options, 3 + ), + ) - if log_level and symmetrize_fc and fc_calculator is None: - print("fc3 was symmetrized.") + if log_level and symmetrize_fc and fc_calculator is None: + print("fc3 was symmetrized.") - if not read_fc["fc2"] and (ph3py.dataset or ph3py.phonon_dataset): + if not read_fc["fc2"] and ( + forces_in_dataset(ph3py.dataset) or forces_in_dataset(ph3py.phonon_dataset) + ): ph3py.produce_fc2( symmetrize_fc2=symmetrize_fc, is_compact_fc=is_compact_fc, @@ -495,6 +500,7 @@ def _get_dataset_or_fc3( p2s_map = ph3py.primitive.p2s_map read_fc3 = False dataset = None + if fc3_filename is not None or pathlib.Path("fc3.hdf5").exists(): if fc3_filename is None: _fc3_filename = "fc3.hdf5" @@ -532,8 +538,17 @@ def _get_dataset_or_fc3( cutoff_pair_distance, log_level, ) - if not forces_in_dataset(dataset): - dataset = None + elif ph3py_yaml is not None and ph3py_yaml.dataset is not None: + # not forces_in_dataset(ph3py_yaml.dataset) + # but want to read displacement dataset. + dataset = _get_dataset_for_fc3( + ph3py, + ph3py_yaml, + None, + phono3py_yaml_filename, + cutoff_pair_distance, + log_level, + ) return read_fc3, dataset @@ -585,8 +600,16 @@ def _get_dataset_phonon_dataset_or_fc2( "phonon_fc2", log_level, ) - if not forces_in_dataset(phonon_dataset): - phonon_dataset = None + elif ph3py_yaml is not None and ph3py_yaml.phonon_dataset is not None: + # not forces_in_dataset(ph3py_yaml.dataset) + # but want to read displacement dataset. + phonon_dataset = _get_dataset_for_fc2( + ph3py, + ph3py_yaml, + None, + "phonon_fc2", + log_level, + ) return read_fc2, phonon_dataset diff --git a/phono3py/cui/phono3py_script.py b/phono3py/cui/phono3py_script.py index a1fc34be..3bf870bf 100644 --- a/phono3py/cui/phono3py_script.py +++ b/phono3py/cui/phono3py_script.py @@ -528,8 +528,6 @@ def store_force_constants( if log_level: print("-" * 29 + " Force constants " + "-" * 30) - (fc_calculator, fc_calculator_options) = get_fc_calculator_params(settings) - read_fc = set_dataset_and_force_constants( phono3py, ph3py_yaml=ph3py_yaml, @@ -538,6 +536,9 @@ def store_force_constants( use_pypolymlp=settings.use_pypolymlp, log_level=log_level, ) + (fc_calculator, fc_calculator_options) = get_fc_calculator_params( + settings, log_level=(not read_fc["fc3"]) * 1 + ) try: compute_force_constants_from_datasets( phono3py, From f96dfab5f7c702796a1f276736419b667336a953 Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Mon, 22 Jul 2024 13:32:30 +0900 Subject: [PATCH 09/15] Add test of phono3py.load that reads displacements --- phono3py/api_phono3py.py | 9 +++++++-- phono3py/cui/create_force_constants.py | 17 +---------------- phono3py/cui/load.py | 2 +- phono3py/phonon3/dataset.py | 11 +++++++++++ test/conftest.py | 16 ++++++++++++++++ test/cui/test_phono3py_load.py | 18 ++++++++++++++++++ test/cui/test_phono3py_load_script.py | 2 +- 7 files changed, 55 insertions(+), 20 deletions(-) create mode 100644 test/cui/test_phono3py_load.py diff --git a/phono3py/api_phono3py.py b/phono3py/api_phono3py.py index 51c5b592..294ac2dd 100644 --- a/phono3py/api_phono3py.py +++ b/phono3py/api_phono3py.py @@ -81,7 +81,7 @@ from phono3py.interface.fc_calculator import get_fc3 from phono3py.interface.phono3py_yaml import Phono3pyYaml from phono3py.phonon.grid import BZGrid -from phono3py.phonon3.dataset import get_displacements_and_forces_fc3 +from phono3py.phonon3.dataset import forces_in_dataset, get_displacements_and_forces_fc3 from phono3py.phonon3.displacement_fc3 import ( direction_to_displacement, get_third_order_displacements, @@ -866,7 +866,7 @@ def displacements(self): for disp2 in disp1["second_atoms"]: displacements[i, disp2["number"]] = disp2["displacement"] i += 1 - elif "forces" in dataset or "displacements" in dataset: + elif "displacements" in dataset: displacements = dataset["displacements"] else: raise RuntimeError("displacement dataset has wrong format.") @@ -2546,6 +2546,11 @@ def _get_forces_energies( Return None if tagert data is not found rather than raising exception. """ + if self._dataset is None: + return None + if not forces_in_dataset(self._dataset): + return None + if target in self._dataset: # type-2 return self._dataset[target] elif "first_atoms" in self._dataset: # type-1 diff --git a/phono3py/cui/create_force_constants.py b/phono3py/cui/create_force_constants.py index d0fb4d44..fbe9a5e3 100644 --- a/phono3py/cui/create_force_constants.py +++ b/phono3py/cui/create_force_constants.py @@ -68,6 +68,7 @@ ) from phono3py.interface.fc_calculator import extract_fc2_fc3_calculators from phono3py.interface.phono3py_yaml import Phono3pyYaml +from phono3py.phonon3.dataset import forces_in_dataset from phono3py.phonon3.fc3 import ( set_permutation_symmetry_fc3, set_translational_invariance_fc3, @@ -316,22 +317,6 @@ def parse_forces( return dataset -def forces_in_dataset(dataset: dict) -> bool: - """Return whether forces in dataset or not.""" - if dataset is None: - return False - return "forces" in dataset or ( - "first_atoms" in dataset and "forces" in dataset["first_atoms"][0] - ) - - -def displacements_in_dataset(dataset: Optional[dict]) -> bool: - """Return whether displacements in dataset or not.""" - if dataset is None: - return False - return "displacements" in dataset or "first_atoms" in dataset - - def get_fc_calculator_params(settings, log_level=0): """Return fc_calculator and fc_calculator_params from settings.""" fc_calculator = None diff --git a/phono3py/cui/load.py b/phono3py/cui/load.py index 12fb9cbf..a5acd263 100644 --- a/phono3py/cui/load.py +++ b/phono3py/cui/load.py @@ -49,13 +49,13 @@ from phono3py import Phono3py from phono3py.cui.create_force_constants import ( - forces_in_dataset, parse_forces, run_pypolymlp_to_compute_forces, ) from phono3py.file_IO import read_fc2_from_hdf5, read_fc3_from_hdf5 from phono3py.interface.fc_calculator import extract_fc2_fc3_calculators from phono3py.interface.phono3py_yaml import Phono3pyYaml +from phono3py.phonon3.dataset import forces_in_dataset from phono3py.phonon3.fc3 import show_drift_fc3 diff --git a/phono3py/phonon3/dataset.py b/phono3py/phonon3/dataset.py index 74fdea86..777910f9 100644 --- a/phono3py/phonon3/dataset.py +++ b/phono3py/phonon3/dataset.py @@ -34,6 +34,8 @@ # ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # POSSIBILITY OF SUCH DAMAGE. +from typing import Optional + import numpy as np @@ -94,3 +96,12 @@ def get_displacements_and_forces_fc3(disp_dataset): return disp_dataset["displacements"], disp_dataset["forces"] else: raise RuntimeError("disp_dataset doesn't contain correct information.") + + +def forces_in_dataset(dataset: Optional[dict]) -> bool: + """Return whether forces in dataset or not.""" + if dataset is None: + return False + return "forces" in dataset or ( + "first_atoms" in dataset and "forces" in dataset["first_atoms"][0] + ) diff --git a/test/conftest.py b/test/conftest.py index 61c2408d..340eda03 100644 --- a/test/conftest.py +++ b/test/conftest.py @@ -69,6 +69,22 @@ def si_pbesol(request) -> Phono3py: ) +@pytest.fixture(scope="session") +def si_pbesol_without_forcesets(request) -> Phono3py: + """Return Phono3py instance of Si 2x2x2 without force sets. + + * with symmetry + + """ + yaml_filename = cwd / "phono3py_si_pbesol.yaml" + enable_v2 = request.config.getoption("--v2") + return phono3py.load( + yaml_filename, + make_r0_average=not enable_v2, + log_level=1, + ) + + @pytest.fixture(scope="session") def si_pbesol_grg(request) -> Phono3py: """Return Phono3py instance of Si 2x2x2. diff --git a/test/cui/test_phono3py_load.py b/test/cui/test_phono3py_load.py new file mode 100644 index 00000000..4f7df4aa --- /dev/null +++ b/test/cui/test_phono3py_load.py @@ -0,0 +1,18 @@ +"""Tests of Phono3py load.""" + +from __future__ import annotations + +from phono3py import Phono3py + + +def test_phono3py_load(si_pbesol_without_forcesets: Phono3py): + """Test phono3py.load. + + Check phono3py.load can read displacements from phono3py_disp.yaml like file that + doesn't contain forces. + + """ + ph3 = si_pbesol_without_forcesets + assert ph3.dataset is not None + assert ph3.displacements.shape == (111, 64, 3) + assert ph3.forces is None diff --git a/test/cui/test_phono3py_load_script.py b/test/cui/test_phono3py_load_script.py index e6987194..4437e06d 100644 --- a/test/cui/test_phono3py_load_script.py +++ b/test/cui/test_phono3py_load_script.py @@ -1,4 +1,4 @@ -"""Tests of Phono3py API.""" +"""Tests of phono3py-load script.""" from __future__ import annotations From d1dbb3a60ff7b2641f8db71aefd0421601d77e05 Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Mon, 22 Jul 2024 13:43:57 +0900 Subject: [PATCH 10/15] Drop py3.8 tests in github workflow --- .github/workflows/phono3py-pytest-conda.yml | 2 +- .github/workflows/publish-to-test-pypi.yml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/phono3py-pytest-conda.yml b/.github/workflows/phono3py-pytest-conda.yml index e9f0f1e2..b384ea41 100644 --- a/.github/workflows/phono3py-pytest-conda.yml +++ b/.github/workflows/phono3py-pytest-conda.yml @@ -12,7 +12,7 @@ jobs: shell: bash -l {0} strategy: matrix: - python-version: ["3.8", "3.9", "3.10", "3.11", "3.12"] + python-version: ["3.9", "3.10", "3.11", "3.12"] steps: - uses: actions/checkout@v4 diff --git a/.github/workflows/publish-to-test-pypi.yml b/.github/workflows/publish-to-test-pypi.yml index 65606fbe..f0d3a455 100644 --- a/.github/workflows/publish-to-test-pypi.yml +++ b/.github/workflows/publish-to-test-pypi.yml @@ -11,7 +11,7 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - python-version: [3.12, ] + python-version: ["3.12"] steps: - uses: actions/checkout@v4 From a77b73a5d2ff05978f44bbce0c10ab9bfd75bf12 Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Mon, 22 Jul 2024 13:52:44 +0900 Subject: [PATCH 11/15] Set version 3.3.2 --- doc/changelog.md | 5 +++++ doc/conf.py | 2 +- phono3py/version.py | 2 +- test/cui/test_phono3py_load.py | 4 ++-- 4 files changed, 9 insertions(+), 4 deletions(-) diff --git a/doc/changelog.md b/doc/changelog.md index fe894aa2..480a9b91 100644 --- a/doc/changelog.md +++ b/doc/changelog.md @@ -2,6 +2,11 @@ # Change Log +## Jul-22-2024: Version 3.3.2 + +- Minor fix of `phono3py.load` function for reading displacements from + `phono3py_disp.yaml` like file that doesn't contain forces. + ## Jul-8-2024: Version 3.3.1 - Major refactoring to isolate reciprocal space grid code. diff --git a/doc/conf.py b/doc/conf.py index b7536d4a..bf1a79a5 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.1" +release = "3.3.2" # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. diff --git a/phono3py/version.py b/phono3py/version.py index dcc4a40d..a30ed251 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.1" +__version__ = "3.3.2" diff --git a/test/cui/test_phono3py_load.py b/test/cui/test_phono3py_load.py index 4f7df4aa..5a651bd9 100644 --- a/test/cui/test_phono3py_load.py +++ b/test/cui/test_phono3py_load.py @@ -8,8 +8,8 @@ def test_phono3py_load(si_pbesol_without_forcesets: Phono3py): """Test phono3py.load. - Check phono3py.load can read displacements from phono3py_disp.yaml like file that - doesn't contain forces. + Check phono3py.load can read displacements from phono3py_disp.yaml like file + that doesn't contain forces. """ ph3 = si_pbesol_without_forcesets From 31603648df3ab4fa94c25299a6e1bb42444e3eef Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 22 Jul 2024 23:15:49 +0000 Subject: [PATCH 12/15] [pre-commit.ci] pre-commit autoupdate MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit updates: - [github.com/astral-sh/ruff-pre-commit: v0.5.2 → v0.5.4](https://github.com/astral-sh/ruff-pre-commit/compare/v0.5.2...v0.5.4) --- .pre-commit-config.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 293ee4b8..1cc6f884 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.2 + rev: v0.5.4 hooks: - id: ruff args: [ "--fix", "--show-fixes" ] From 1cd20edc33648acbc74867dd8539bde73a6fab67 Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Tue, 23 Jul 2024 14:36:45 +0900 Subject: [PATCH 13/15] Refactor workflows for running tests --- .../workflows/phono3py-pytest-conda-mkl-phphmtblas.yml | 9 +++------ .github/workflows/phono3py-pytest-conda-mkl-v2.yml | 9 +++------ .github/workflows/phono3py-pytest-conda-mkl.yml | 1 + .github/workflows/phono3py-pytest-conda-numpy2.yml | 7 ++----- .github/workflows/phono3py-pytest-conda-phphmtblas.yml | 9 +++------ .github/workflows/phono3py-pytest-conda.yml | 9 +++------ 6 files changed, 15 insertions(+), 29 deletions(-) diff --git a/.github/workflows/phono3py-pytest-conda-mkl-phphmtblas.yml b/.github/workflows/phono3py-pytest-conda-mkl-phphmtblas.yml index f8c54bed..f32a9b58 100644 --- a/.github/workflows/phono3py-pytest-conda-mkl-phphmtblas.yml +++ b/.github/workflows/phono3py-pytest-conda-mkl-phphmtblas.yml @@ -24,7 +24,7 @@ jobs: run: | conda activate test conda install --yes python=${{ matrix.python-version }} - conda install --yes matplotlib-base pyyaml "libblas=*=*mkl" mkl-include h5py scipy pytest codecov pytest-cov spglib alm cmake c-compiler cxx-compiler + conda install --yes matplotlib-base pyyaml "libblas=*=*mkl" mkl-include h5py scipy pytest spglib alm cmake c-compiler cxx-compiler - name: Install symfc develop branch run: | conda activate test @@ -45,8 +45,5 @@ jobs: PHPHCALC_USE_MTBLAS=ON pip install -e . -vvv - name: Run pytest run: | - pytest -v --cov=./ --cov-report=xml test - - name: Upload coverage to Codecov - uses: codecov/codecov-action@v3 - with: - verbose: true + conda activate test + pytest -v test diff --git a/.github/workflows/phono3py-pytest-conda-mkl-v2.yml b/.github/workflows/phono3py-pytest-conda-mkl-v2.yml index 91c32ee9..7a285923 100644 --- a/.github/workflows/phono3py-pytest-conda-mkl-v2.yml +++ b/.github/workflows/phono3py-pytest-conda-mkl-v2.yml @@ -24,7 +24,7 @@ jobs: run: | conda activate test conda install --yes python=${{ matrix.python-version }} - conda install --yes matplotlib-base pyyaml "libblas=*=*mkl" mkl-include h5py scipy pytest codecov pytest-cov spglib alm cmake c-compiler cxx-compiler + conda install --yes matplotlib-base pyyaml "libblas=*=*mkl" mkl-include h5py scipy pytest spglib alm cmake c-compiler cxx-compiler - name: Install symfc develop branch run: | conda activate test @@ -45,8 +45,5 @@ jobs: pip install -e . -vvv - name: Run pytest run: | - pytest --v2 -v --cov=./ --cov-report=xml test - - name: Upload coverage to Codecov - uses: codecov/codecov-action@v3 - with: - verbose: true + conda activate test + pytest --v2 -v test diff --git a/.github/workflows/phono3py-pytest-conda-mkl.yml b/.github/workflows/phono3py-pytest-conda-mkl.yml index 72616e31..a9be299a 100644 --- a/.github/workflows/phono3py-pytest-conda-mkl.yml +++ b/.github/workflows/phono3py-pytest-conda-mkl.yml @@ -45,6 +45,7 @@ jobs: pip install -e . -vvv - name: Run pytest run: | + conda activate test pytest -v --cov=./ --cov-report=xml test - name: Upload coverage to Codecov uses: codecov/codecov-action@v3 diff --git a/.github/workflows/phono3py-pytest-conda-numpy2.yml b/.github/workflows/phono3py-pytest-conda-numpy2.yml index 92f010cd..fdb106b2 100644 --- a/.github/workflows/phono3py-pytest-conda-numpy2.yml +++ b/.github/workflows/phono3py-pytest-conda-numpy2.yml @@ -24,7 +24,7 @@ jobs: run: | conda activate test conda install --yes python=${{ matrix.python-version }} - conda install --yes matplotlib-base pyyaml "libblas=*=*openblas" openblas h5py "numpy=2" scipy pytest codecov pytest-cov cmake c-compiler cxx-compiler + conda install --yes matplotlib-base pyyaml "libblas=*=*openblas" openblas h5py "numpy=2" scipy pytest cmake c-compiler cxx-compiler - name: Install spglib develop branch run: | conda activate test @@ -52,8 +52,5 @@ jobs: pip install -e . -vvv - name: Run pytest run: | + conda activate test pytest -v test - - name: Upload coverage to Codecov - uses: codecov/codecov-action@v3 - with: - verbose: true diff --git a/.github/workflows/phono3py-pytest-conda-phphmtblas.yml b/.github/workflows/phono3py-pytest-conda-phphmtblas.yml index 37473838..59f41d31 100644 --- a/.github/workflows/phono3py-pytest-conda-phphmtblas.yml +++ b/.github/workflows/phono3py-pytest-conda-phphmtblas.yml @@ -24,7 +24,7 @@ jobs: run: | conda activate test conda install --yes python=${{ matrix.python-version }} - conda install --yes matplotlib-base pyyaml "libblas=*=*openblas" openblas h5py scipy pytest codecov pytest-cov spglib alm cmake c-compiler cxx-compiler + conda install --yes matplotlib-base pyyaml "libblas=*=*openblas" openblas h5py scipy pytest spglib alm cmake c-compiler cxx-compiler - name: Install symfc develop branch run: | conda activate test @@ -45,8 +45,5 @@ jobs: PHPHCALC_USE_MTBLAS=ON pip install -e . -vvv - name: Run pytest run: | - pytest -v --cov=./ --cov-report=xml test - - name: Upload coverage to Codecov - uses: codecov/codecov-action@v3 - with: - verbose: true + conda activate test + pytest -v test diff --git a/.github/workflows/phono3py-pytest-conda.yml b/.github/workflows/phono3py-pytest-conda.yml index b384ea41..ad620453 100644 --- a/.github/workflows/phono3py-pytest-conda.yml +++ b/.github/workflows/phono3py-pytest-conda.yml @@ -24,7 +24,7 @@ jobs: run: | conda activate test conda install --yes python=${{ matrix.python-version }} - conda install --yes matplotlib-base pyyaml "libblas=*=*openblas" openblas h5py scipy pytest codecov pytest-cov spglib alm cmake c-compiler cxx-compiler + conda install --yes matplotlib-base pyyaml "libblas=*=*openblas" openblas h5py scipy pytest spglib alm cmake c-compiler cxx-compiler - name: Install phonopy develop branch run: | conda activate test @@ -38,8 +38,5 @@ jobs: pip install -e . -vvv - name: Run pytest run: | - pytest -v --cov=./ --cov-report=xml test - - name: Upload coverage to Codecov - uses: codecov/codecov-action@v3 - with: - verbose: true + conda activate test + pytest -v test From 323e2bfc9ac4a2823bb16bb5fac3caf7173de251 Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Thu, 25 Jul 2024 12:13:18 +0900 Subject: [PATCH 14/15] Update write_FORCES_FC3 for writing random displacements --- phono3py/cui/create_force_sets.py | 9 ++++-- phono3py/file_IO.py | 49 ++++++++++++++++++++++++------- 2 files changed, 45 insertions(+), 13 deletions(-) diff --git a/phono3py/cui/create_force_sets.py b/phono3py/cui/create_force_sets.py index 8916aa17..38bfe428 100644 --- a/phono3py/cui/create_force_sets.py +++ b/phono3py/cui/create_force_sets.py @@ -107,6 +107,8 @@ def create_FORCES_FC3_and_FORCES_FC2( print("%s could not be created." % "FORCES_FC2") print_error() sys.exit(1) + else: + calc_dataset_fc2 = None if settings.save_params: fc3_yaml_filename = "phono3py_params.yaml" @@ -349,7 +351,10 @@ def _get_force_sets_fc3( def _set_forces_and_nac_params( - ph3py_yaml: Phono3pyYaml, settings, calc_dataset_fc3: dict, calc_dataset_fc2: dict + ph3py_yaml: Phono3pyYaml, + settings, + calc_dataset_fc3: dict, + calc_dataset_fc2: Optional[dict], ): if "first_atoms" in ph3py_yaml.dataset: count = len(ph3py_yaml.dataset["first_atoms"]) @@ -375,7 +380,7 @@ def _set_forces_and_nac_params( calc_dataset_fc3["supercell_energies"], dtype="double" ) - if settings.create_forces_fc2: + if settings.create_forces_fc2 and calc_dataset_fc2: if "first_atoms" in ph3py_yaml.phonon_dataset: for i, d in enumerate(ph3py_yaml.phonon_dataset["first_atoms"]): d["forces"] = calc_dataset_fc2["forces"][i] diff --git a/phono3py/file_IO.py b/phono3py/file_IO.py index e43dfa7f..58a94ffc 100644 --- a/phono3py/file_IO.py +++ b/phono3py/file_IO.py @@ -38,7 +38,7 @@ import os import warnings from collections.abc import Sequence -from typing import Optional, Union +from typing import Optional, TextIO, Union import h5py import numpy as np @@ -46,7 +46,11 @@ # This import is deactivated for a while. # from phonopy.file_IO import write_force_constants_to_hdf5 -from phonopy.file_IO import check_force_constants_indices, get_cell_from_disp_yaml +from phonopy.file_IO import ( + check_force_constants_indices, + get_cell_from_disp_yaml, + write_FORCE_SETS, +) from phono3py.version import __version__ @@ -211,18 +215,46 @@ def write_FORCES_FC2(disp_dataset, forces_fc2=None, fp=None, filename="FORCES_FC w.close() -def write_FORCES_FC3(disp_dataset, forces_fc3=None, fp=None, filename="FORCES_FC3"): +def write_FORCES_FC3( + disp_dataset: dict, + forces_fc3: Optional[Sequence] = None, + fp: Optional[TextIO] = None, + filename: str = "FORCES_FC3", +): """Write FORCES_FC3. fp : IO object, optional, default=None When this is given, FORCES_FC3 content is written into this IO object. """ - if fp is None: - w = open(filename, "w") + if "first_atoms" in disp_dataset: + if fp is None: + with open(filename, "w") as w: + _write_FORCES_FC3_typeI(disp_dataset, w, forces_fc3=forces_fc3) + else: + _write_FORCES_FC3_typeI(disp_dataset, fp, forces_fc3=forces_fc3) else: - w = fp + if "forces" in disp_dataset: + write_FORCE_SETS(disp_dataset, filename="FORCES_FC3") + else: + if forces_fc3 is None: + raise RuntimeError("No forces are found.") + dataset = disp_dataset.copy() + dataset["forces"] = forces_fc3 + write_FORCE_SETS(dataset, filename="FORCES_FC3") + + +def _write_FORCES_FC3_typeI( + disp_dataset: dict, + w: Optional[TextIO], + forces_fc3: Optional[Sequence] = None, +): + """Write FORCES_FC3. + + w : TextIO object, optional, default=None + When this is given, FORCES_FC3 content is written into this IO object. + """ natom = disp_dataset["natom"] num_disp1 = len(disp_dataset["first_atoms"]) count = num_disp1 @@ -253,15 +285,10 @@ def write_FORCES_FC3(disp_dataset, forces_fc3=None, fp=None, filename="FORCES_FC w.write("%15.10f %15.10f %15.10f\n" % tuple(force)) file_count += 1 else: - # for forces in forces_fc3[i]: - # w.write("%15.10f %15.10f %15.10f\n" % (tuple(forces))) for _ in range(natom): w.write("%15.10f %15.10f %15.10f\n" % (0, 0, 0)) count += 1 - if fp is None: - w.close() - def write_fc3_to_hdf5(fc3, filename="fc3.hdf5", p2s_map=None, compression="gzip"): """Write fc3 in fc3.hdf5. From 10d66419829332e09060fba91986657900652312 Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Sat, 27 Jul 2024 12:13:00 +0900 Subject: [PATCH 15/15] Update installation doc --- doc/install.md | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/doc/install.md b/doc/install.md index f3103668..04a10c1f 100644 --- a/doc/install.md +++ b/doc/install.md @@ -136,12 +136,13 @@ wrong python libraries can be imported. % cd phonopy % pip install -e . -vvv % cd ../phono3py - % pip install -e . -vvv + % pip install . -vvv ``` - The conda packages dependency can often change and this recipe may not work - properly. So if you find this instruction doesn't work, it is very - appreciated if letting us know it in the phonopy mailing list. + The editable install (`pip install -e`) may not work depending on the + computer environment. The conda packages dependency can often change and this + recipe may not work properly. So if you find this instruction doesn't work, + it is very appreciated if letting us know it in the phonopy mailing list. (install_lapacke)= ## Installation of LAPACKE