diff --git a/.github/workflows/phono3py-pytest-conda-nolapacke.yml b/.github/workflows/phono3py-pytest-conda-nolapacke.yml new file mode 100644 index 00000000..52d42829 --- /dev/null +++ b/.github/workflows/phono3py-pytest-conda-nolapacke.yml @@ -0,0 +1,42 @@ +name: Pytest without linking BLAS and LAPACK in C + +on: + pull_request: + branches: [ develop ] + +jobs: + build-linux: + runs-on: ubuntu-latest + defaults: + run: + shell: bash -l {0} + strategy: + matrix: + python-version: ["3.12"] + + steps: + - uses: actions/checkout@v4 + # Use conda-incubator/setup-miniconda for precise control of conda infrastructure + - uses: conda-incubator/setup-miniconda@v3 + with: + miniforge-version: latest + - name: Install dependent packages + run: | + conda activate test + conda install --yes python=${{ matrix.python-version }} + conda install --yes matplotlib-base pyyaml h5py scipy pytest spglib cmake c-compiler cxx-compiler + - name: Install phonopy develop branch + run: | + conda activate test + git clone https://github.com/phonopy/phonopy.git + cd phonopy + pip install -e . -vvv + cd .. + - name: Install phono3py + run: | + conda activate test + BUILD_WITHOUT_LAPACKE=ON pip install -e . -vvv + - name: Run pytest + run: | + conda activate test + pytest -v test diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 39c7ee27..eebc7b06 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -11,7 +11,7 @@ repos: exclude: ^example/AlN-LDA/ - repo: https://github.com/astral-sh/ruff-pre-commit - rev: v0.7.3 + rev: v0.7.4 hooks: - id: ruff args: [ "--fix", "--show-fixes" ] diff --git a/CMakeLists.txt b/CMakeLists.txt index a4d236e2..11dc20d1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -9,6 +9,7 @@ option(PHONO3PY_USE_OMP "Option to search OpenMP library" ON) option(PHONO3PY_USE_MTBLAS "Use multithread BLAS if it exists" ON) option(PHONO3PY_WITH_TESTS "build unit tests" OFF) option(BUILD_SHARED_LIBS "Option to build shared library" OFF) +option(BUILD_WITHOUT_LAPACKE "Option to build without LAPACKE" ON) if(PHONO3PY_WITH_Fortran) enable_language(Fortran) @@ -89,7 +90,8 @@ endif() if(BUILD_PHPHCALC_LIB OR BUILD_PHONONCALC_LIB - OR BUILD_NANOBIND_MODULE) + OR BUILD_NANOBIND_MODULE + AND (NOT BUILD_WITHOUT_LAPACKE)) find_package(BLAS REQUIRED) # set BLAS_LIBRARIES if(BLAS_FOUND) @@ -105,7 +107,7 @@ if(BUILD_PHPHCALC_LIB endif() if(BLAS_LIBRARIES MATCHES "libmkl") - message(STATUS "MKL detected: Set C-macro MKL_LAPACKE.") + message(STATUS "MKL detected.") if(PHONO3PY_USE_MTBLAS) message( @@ -190,68 +192,95 @@ if(BUILD_PHPHCALC_LIB OR BUILD_NANOBIND_MODULE) # Shared library add_library(phphcalc_lib SHARED ${SOURCES_PHPHCALC}) - if(OpenMP_FOUND) - target_link_libraries( - phphcalc_lib PRIVATE recgrid_lib BLAS::BLAS LAPACK::LAPACK - OpenMP::OpenMP_C) - else() - target_link_libraries(phphcalc_lib PRIVATE recgrid_lib BLAS::BLAS - LAPACK::LAPACK) - endif() - - target_include_directories(phphcalc_lib PRIVATE ${MY_INCLUDES}) - - if(BLAS_LIBRARIES MATCHES "libmkl") - if(PHONO3PY_USE_MTBLAS) - target_compile_definitions( - phphcalc_lib PRIVATE MKL_LAPACKE MULTITHREADED_BLAS - THM_EPSILON=1e-10) + if(BUILD_WITHOUT_LAPACKE) + if(OpenMP_FOUND) + target_link_libraries(phphcalc_lib PRIVATE recgrid_lib + OpenMP::OpenMP_C) else() - target_compile_definitions(phphcalc_lib PRIVATE MKL_LAPACKE - THM_EPSILON=1e-10) + target_link_libraries(phphcalc_lib PRIVATE recgrid_lib) endif() - endif() - if(BLAS_LIBRARIES MATCHES "libopenblas") - if(PHONO3PY_USE_MTBLAS) - target_compile_definitions(phphcalc_lib PRIVATE MULTITHREADED_BLAS - THM_EPSILON=1e-10) + target_compile_definitions(phphcalc_lib PRIVATE NO_INCLUDE_LAPACKE + THM_EPSILON=1e-10) + else() + if(OpenMP_FOUND) + target_link_libraries( + phphcalc_lib PRIVATE recgrid_lib BLAS::BLAS LAPACK::LAPACK + OpenMP::OpenMP_C) else() - target_compile_definitions(phphcalc_lib PRIVATE THM_EPSILON=1e-10) + target_link_libraries(phphcalc_lib PRIVATE recgrid_lib BLAS::BLAS + LAPACK::LAPACK) endif() - endif() - else() - # Static link library - add_library(phphcalc_lib STATIC ${SOURCES_PHPHCALC}) - if(OpenMP_FOUND) - target_link_libraries(phphcalc_lib recgrid_lib BLAS::BLAS LAPACK::LAPACK - OpenMP::OpenMP_C) - else() - target_link_libraries(phphcalc_lib recgrid_lib BLAS::BLAS LAPACK::LAPACK) + if(BLAS_LIBRARIES MATCHES "libmkl") + if(PHONO3PY_USE_MTBLAS) + target_compile_definitions( + phphcalc_lib PRIVATE MKL_BLAS MULTITHREADED_BLAS + THM_EPSILON=1e-10) + else() + target_compile_definitions(phphcalc_lib + PRIVATE MKL_BLAS THM_EPSILON=1e-10) + endif() + endif() + + if(BLAS_LIBRARIES MATCHES "libopenblas") + if(PHONO3PY_USE_MTBLAS) + target_compile_definitions( + phphcalc_lib PRIVATE MULTITHREADED_BLAS THM_EPSILON=1e-10) + else() + target_compile_definitions(phphcalc_lib + PRIVATE THM_EPSILON=1e-10) + endif() + endif() endif() target_include_directories(phphcalc_lib PRIVATE ${MY_INCLUDES}) - if(BLAS_LIBRARIES MATCHES "libmkl") - if(PHONO3PY_USE_MTBLAS) - target_compile_definitions( - phphcalc_lib PRIVATE MKL_LAPACKE MULTITHREADED_BLAS - THM_EPSILON=1e-10) + else() + # Static link library + add_library(phphcalc_lib STATIC ${SOURCES_PHPHCALC}) + + if(BUILD_WITHOUT_LAPACKE) + if(OpenMP_FOUND) + target_link_libraries(phphcalc_lib recgrid_lib OpenMP::OpenMP_C) else() - target_compile_definitions(phphcalc_lib PRIVATE MKL_LAPACKE - THM_EPSILON=1e-10) + target_link_libraries(phphcalc_lib recgrid_lib) endif() - endif() - if(BLAS_LIBRARIES MATCHES "libopenblas") - if(PHONO3PY_USE_MTBLAS) - target_compile_definitions(phphcalc_lib PRIVATE MULTITHREADED_BLAS - THM_EPSILON=1e-10) + target_compile_definitions(phphcalc_lib PRIVATE NO_INCLUDE_LAPACKE + THM_EPSILON=1e-10) + else() + if(OpenMP_FOUND) + target_link_libraries(phphcalc_lib recgrid_lib BLAS::BLAS + LAPACK::LAPACK OpenMP::OpenMP_C) else() - target_compile_definitions(phphcalc_lib PRIVATE THM_EPSILON=1e-10) + target_link_libraries(phphcalc_lib recgrid_lib BLAS::BLAS + LAPACK::LAPACK) + endif() + + if(BLAS_LIBRARIES MATCHES "libmkl") + if(PHONO3PY_USE_MTBLAS) + target_compile_definitions( + phphcalc_lib PRIVATE MKL_BLAS MULTITHREADED_BLAS + THM_EPSILON=1e-10) + else() + target_compile_definitions(phphcalc_lib + PRIVATE MKL_BLAS THM_EPSILON=1e-10) + endif() + endif() + + if(BLAS_LIBRARIES MATCHES "libopenblas") + if(PHONO3PY_USE_MTBLAS) + target_compile_definitions( + phphcalc_lib PRIVATE MULTITHREADED_BLAS THM_EPSILON=1e-10) + else() + target_compile_definitions(phphcalc_lib + PRIVATE THM_EPSILON=1e-10) + endif() endif() endif() + + target_include_directories(phphcalc_lib PRIVATE ${MY_INCLUDES}) endif() if(NOT BUILD_NANOBIND_MODULE) @@ -285,45 +314,67 @@ if(BUILD_PHONONCALC_LIB OR BUILD_NANOBIND_MODULE) # Shared library add_library(phononcalc_lib SHARED ${SOURCES_PHONONCALC}) - if(OpenMP_FOUND) - target_link_libraries(phononcalc_lib BLAS::BLAS LAPACK::LAPACK - OpenMP::OpenMP_C) + if(BUILD_WITHOUT_LAPACKE) + if(OpenMP_FOUND) + target_link_libraries(phononcalc_lib OpenMP::OpenMP_C) + else() + target_link_libraries(phononcalc_lib) + endif() + + target_compile_definitions(phononcalc_lib PRIVATE NO_INCLUDE_LAPACKE) else() - target_link_libraries(phononcalc_lib BLAS::BLAS LAPACK::LAPACK) - endif() + if(OpenMP_FOUND) + target_link_libraries(phononcalc_lib BLAS::BLAS LAPACK::LAPACK + OpenMP::OpenMP_C) + else() + target_link_libraries(phononcalc_lib BLAS::BLAS LAPACK::LAPACK) + endif() - target_include_directories(phononcalc_lib PRIVATE ${MY_INCLUDES}) + if(BLAS_LIBRARIES MATCHES "libmkl") + target_compile_definitions(phononcalc_lib PRIVATE MKL_BLAS + MULTITHREADED_BLAS) + endif() - if(BLAS_LIBRARIES MATCHES "libmkl") - target_compile_definitions(phononcalc_lib PRIVATE MKL_LAPACKE - MULTITHREADED_BLAS) + if(BLAS_LIBRARIES MATCHES "libopenblas") + target_compile_definitions(phononcalc_lib PRIVATE MULTITHREADED_BLAS) + endif() endif() - if(BLAS_LIBRARIES MATCHES "libopenblas") - target_compile_definitions(phononcalc_lib PRIVATE MULTITHREADED_BLAS) - endif() + target_include_directories(phononcalc_lib PRIVATE ${MY_INCLUDES}) else() # Static link library add_library(phononcalc_lib STATIC ${SOURCES_PHONONCALC}) - if(OpenMP_FOUND) - target_link_libraries(phononcalc_lib PRIVATE BLAS::BLAS LAPACK::LAPACK - OpenMP::OpenMP_C) + if(BUILD_WITHOUT_LAPACKE) + if(OpenMP_FOUND) + target_link_libraries(phononcalc_lib PRIVATE OpenMP::OpenMP_C) + else() + target_link_libraries(phononcalc_lib PRIVATE) + endif() + + target_compile_definitions(phononcalc_lib PRIVATE NO_INCLUDE_LAPACKE) else() - target_link_libraries(phononcalc_lib PRIVATE BLAS::BLAS LAPACK::LAPACK) - endif() + if(OpenMP_FOUND) + target_link_libraries( + phononcalc_lib PRIVATE BLAS::BLAS LAPACK::LAPACK + OpenMP::OpenMP_C) + else() + target_link_libraries(phononcalc_lib PRIVATE BLAS::BLAS + LAPACK::LAPACK) + endif() - target_include_directories(phononcalc_lib PRIVATE ${MY_INCLUDES}) + if(BLAS_LIBRARIES MATCHES "libmkl") + target_compile_definitions(phononcalc_lib PRIVATE MKL_BLAS + MULTITHREADED_BLAS) + endif() - if(BLAS_LIBRARIES MATCHES "libmkl") - target_compile_definitions(phononcalc_lib PRIVATE MKL_LAPACKE - MULTITHREADED_BLAS) + if(BLAS_LIBRARIES MATCHES "libopenblas") + target_compile_definitions(phononcalc_lib PRIVATE MULTITHREADED_BLAS) + endif() endif() - if(BLAS_LIBRARIES MATCHES "libopenblas") - target_compile_definitions(phononcalc_lib PRIVATE MULTITHREADED_BLAS) - endif() + target_include_directories(phononcalc_lib PRIVATE ${MY_INCLUDES}) endif() if(NOT BUILD_NANOBIND_MODULE) @@ -432,7 +483,17 @@ if(BUILD_NANOBIND_MODULE) target_link_libraries(_phononcalc PRIVATE phononcalc_lib) target_link_libraries(_recgrid PRIVATE recgrid_lib) - target_compile_definitions(_phono3py PRIVATE THM_EPSILON=1e-10) + if(BUILD_WITHOUT_LAPACKE) + target_compile_definitions(_phono3py PRIVATE NO_INCLUDE_LAPACKE + THM_EPSILON=1e-10) + else() + if(BLAS_LIBRARIES MATCHES "libmkl") + target_compile_definitions(_phono3py PRIVATE MKL_BLAS THM_EPSILON=1e-10) + else() + target_compile_definitions(_phono3py PRIVATE THM_EPSILON=1e-10) + endif() + endif() + install(TARGETS _phono3py LIBRARY DESTINATION ${SKBUILD_PROJECT_NAME}) install(TARGETS _phononcalc LIBRARY DESTINATION ${SKBUILD_PROJECT_NAME}) install(TARGETS _recgrid LIBRARY DESTINATION ${SKBUILD_PROJECT_NAME}) diff --git a/c/_phono3py.c b/c/_phono3py.c deleted file mode 100644 index 27388cb5..00000000 --- a/c/_phono3py.c +++ /dev/null @@ -1,1877 +0,0 @@ -/* Copyright (C) 2015 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 -#include -#include -#include -#include -#include -#include - -// #include "lapack_wrapper.h" -#include "phono3py.h" -#include "phonoc_array.h" - -static PyObject *py_get_interaction(PyObject *self, PyObject *args); -static PyObject *py_get_pp_collision(PyObject *self, PyObject *args); -static PyObject *py_get_pp_collision_with_sigma(PyObject *self, PyObject *args); -static PyObject *py_get_imag_self_energy_with_g(PyObject *self, PyObject *args); -static PyObject *py_get_detailed_imag_self_energy_with_g(PyObject *self, - PyObject *args); -static PyObject *py_get_real_self_energy_at_bands(PyObject *self, - PyObject *args); -static PyObject *py_get_real_self_energy_at_frequency_point(PyObject *self, - PyObject *args); -static PyObject *py_get_collision_matrix(PyObject *self, PyObject *args); -static PyObject *py_get_reducible_collision_matrix(PyObject *self, - PyObject *args); -static PyObject *py_symmetrize_collision_matrix(PyObject *self, PyObject *args); -static PyObject *py_expand_collision_matrix(PyObject *self, PyObject *args); -static PyObject *py_distribute_fc3(PyObject *self, PyObject *args); -static PyObject *py_rotate_delta_fc2s(PyObject *self, PyObject *args); -static PyObject *py_get_isotope_strength(PyObject *self, PyObject *args); -static PyObject *py_get_thm_isotope_strength(PyObject *self, PyObject *args); -static PyObject *py_get_permutation_symmetry_fc3(PyObject *self, - PyObject *args); -static PyObject *py_get_permutation_symmetry_compact_fc3(PyObject *self, - PyObject *args); -static PyObject *py_transpose_compact_fc3(PyObject *self, PyObject *args); -static PyObject *py_get_thm_relative_grid_address(PyObject *self, - PyObject *args); -static PyObject *py_get_neighboring_grid_points(PyObject *self, PyObject *args); -static PyObject *py_get_thm_integration_weights_at_grid_points(PyObject *self, - PyObject *args); -static PyObject *py_tpl_get_triplets_reciprocal_mesh_at_q(PyObject *self, - PyObject *args); -static PyObject *py_tpl_get_BZ_triplets_at_q(PyObject *self, PyObject *args); -static PyObject *py_get_triplets_integration_weights(PyObject *self, - PyObject *args); -static PyObject *py_get_triplets_integration_weights_with_sigma(PyObject *self, - PyObject *args); -static PyObject *py_get_grid_index_from_address(PyObject *self, PyObject *args); -static PyObject *py_get_gr_grid_addresses(PyObject *self, PyObject *args); -static PyObject *py_get_reciprocal_rotations(PyObject *self, PyObject *args); -static PyObject *py_transform_rotations(PyObject *self, PyObject *args); -static PyObject *py_get_snf3x3(PyObject *self, PyObject *args); -static PyObject *py_get_ir_grid_map(PyObject *self, PyObject *args); -static PyObject *py_get_bz_grid_addresses(PyObject *self, PyObject *args); -static PyObject *py_rotate_bz_grid_addresses(PyObject *self, PyObject *args); -static PyObject *py_diagonalize_collision_matrix(PyObject *self, - PyObject *args); -static PyObject *py_pinv_from_eigensolution(PyObject *self, PyObject *args); -static PyObject *py_get_default_colmat_solver(PyObject *self, PyObject *args); -static PyObject *py_lapacke_pinv(PyObject *self, PyObject *args); -static PyObject *py_get_omp_max_threads(PyObject *self, PyObject *args); - -static void show_colmat_info(const PyArrayObject *collision_matrix_py, - const long i_sigma, const long i_temp, - const long adrs_shift); -static Larray *convert_to_larray(PyArrayObject *npyary); -static Darray *convert_to_darray(PyArrayObject *npyary); - -struct module_state { - PyObject *error; -}; - -#if PY_MAJOR_VERSION >= 3 -#define GETSTATE(m) ((struct module_state *)PyModule_GetState(m)) -#else -#define GETSTATE(m) (&_state) -static struct module_state _state; -#endif - -static PyObject *error_out(PyObject *m) { - struct module_state *st = GETSTATE(m); - PyErr_SetString(st->error, "something bad happened"); - return NULL; -} - -static PyMethodDef _phono3py_methods[] = { - {"error_out", (PyCFunction)error_out, METH_NOARGS, NULL}, - {"interaction", (PyCFunction)py_get_interaction, METH_VARARGS, - "Interaction of triplets"}, - {"pp_collision", (PyCFunction)py_get_pp_collision, METH_VARARGS, - "Collision and ph-ph calculation"}, - {"pp_collision_with_sigma", (PyCFunction)py_get_pp_collision_with_sigma, - METH_VARARGS, "Collision and ph-ph calculation for smearing method"}, - {"imag_self_energy_with_g", (PyCFunction)py_get_imag_self_energy_with_g, - METH_VARARGS, "Imaginary part of self energy at frequency points with g"}, - {"detailed_imag_self_energy_with_g", - (PyCFunction)py_get_detailed_imag_self_energy_with_g, METH_VARARGS, - "Detailed contribution to imaginary part of self energy at frequency " - "points with g"}, - {"real_self_energy_at_bands", (PyCFunction)py_get_real_self_energy_at_bands, - METH_VARARGS, "Real part of self energy from third order force constants"}, - {"real_self_energy_at_frequency_point", - (PyCFunction)py_get_real_self_energy_at_frequency_point, METH_VARARGS, - "Real part of self energy from third order force constants at a frequency " - "point"}, - {"collision_matrix", (PyCFunction)py_get_collision_matrix, METH_VARARGS, - "Collision matrix with g"}, - {"reducible_collision_matrix", - (PyCFunction)py_get_reducible_collision_matrix, METH_VARARGS, - "Collision matrix with g for reducible grid points"}, - {"symmetrize_collision_matrix", (PyCFunction)py_symmetrize_collision_matrix, - METH_VARARGS, "Symmetrize collision matrix"}, - {"expand_collision_matrix", (PyCFunction)py_expand_collision_matrix, - METH_VARARGS, "Expand collision matrix"}, - {"distribute_fc3", (PyCFunction)py_distribute_fc3, METH_VARARGS, - "Distribute least fc3 to full fc3"}, - {"rotate_delta_fc2s", (PyCFunction)py_rotate_delta_fc2s, METH_VARARGS, - "Rotate delta fc2s"}, - {"isotope_strength", (PyCFunction)py_get_isotope_strength, METH_VARARGS, - "Isotope scattering strength"}, - {"thm_isotope_strength", (PyCFunction)py_get_thm_isotope_strength, - METH_VARARGS, "Isotope scattering strength for tetrahedron_method"}, - {"permutation_symmetry_fc3", (PyCFunction)py_get_permutation_symmetry_fc3, - METH_VARARGS, "Set permutation symmetry for fc3"}, - {"permutation_symmetry_compact_fc3", - (PyCFunction)py_get_permutation_symmetry_compact_fc3, METH_VARARGS, - "Set permutation symmetry for compact-fc3"}, - {"transpose_compact_fc3", (PyCFunction)py_transpose_compact_fc3, - METH_VARARGS, "Transpose compact fc3"}, - {"tetrahedra_relative_grid_address", - (PyCFunction)py_get_thm_relative_grid_address, METH_VARARGS, - "Relative grid addresses of vertices of 24 tetrahedra"}, - {"neighboring_grid_points", (PyCFunction)py_get_neighboring_grid_points, - METH_VARARGS, "Neighboring grid points by relative grid addresses"}, - {"integration_weights_at_grid_points", - (PyCFunction)py_get_thm_integration_weights_at_grid_points, METH_VARARGS, - "Integration weights of tetrahedron method at grid points"}, - {"triplets_reciprocal_mesh_at_q", - (PyCFunction)py_tpl_get_triplets_reciprocal_mesh_at_q, METH_VARARGS, - "Triplets on reciprocal mesh points at a specific q-point"}, - {"BZ_triplets_at_q", (PyCFunction)py_tpl_get_BZ_triplets_at_q, METH_VARARGS, - "Triplets in reciprocal primitive lattice are transformed to those in " - "BZ."}, - {"triplets_integration_weights", - (PyCFunction)py_get_triplets_integration_weights, METH_VARARGS, - "Integration weights of tetrahedron method for triplets"}, - {"triplets_integration_weights_with_sigma", - (PyCFunction)py_get_triplets_integration_weights_with_sigma, METH_VARARGS, - "Integration weights of smearing method for triplets"}, - {"grid_index_from_address", (PyCFunction)py_get_grid_index_from_address, - METH_VARARGS, "Grid index from grid address"}, - {"ir_grid_map", (PyCFunction)py_get_ir_grid_map, METH_VARARGS, - "Reciprocal mesh points with ir grid mapping table"}, - {"gr_grid_addresses", (PyCFunction)py_get_gr_grid_addresses, METH_VARARGS, - "Get generalized regular grid addresses"}, - {"reciprocal_rotations", (PyCFunction)py_get_reciprocal_rotations, - METH_VARARGS, "Return rotation matrices in reciprocal space"}, - {"transform_rotations", (PyCFunction)py_transform_rotations, METH_VARARGS, - "Transform rotations to those in generalized regular grid"}, - {"snf3x3", (PyCFunction)py_get_snf3x3, METH_VARARGS, - "Get Smith formal form for 3x3 integer matrix"}, - {"bz_grid_addresses", (PyCFunction)py_get_bz_grid_addresses, METH_VARARGS, - "Get grid addresses including Brillouin zone surface"}, - {"rotate_bz_grid_index", (PyCFunction)py_rotate_bz_grid_addresses, - METH_VARARGS, "Rotate grid point considering Brillouin zone surface"}, - {"diagonalize_collision_matrix", - (PyCFunction)py_diagonalize_collision_matrix, METH_VARARGS, - "Diagonalize and optionally pseudo-inverse using Lapack dsyev(d)"}, - {"pinv_from_eigensolution", (PyCFunction)py_pinv_from_eigensolution, - METH_VARARGS, "Pseudo-inverse from eigensolution"}, - {"default_colmat_solver", (PyCFunction)py_get_default_colmat_solver, - METH_VARARGS, "Return default collison matrix solver by integer value"}, - {"lapacke_pinv", (PyCFunction)py_lapacke_pinv, METH_VARARGS, - "Pseudo inversion using lapacke."}, - {"omp_max_threads", py_get_omp_max_threads, METH_VARARGS, - "Return openmp max number of threads. Return 0 unless openmp is " - "activated. "}, - {NULL, NULL, 0, NULL}}; - -#if PY_MAJOR_VERSION >= 3 - -static int _phono3py_traverse(PyObject *m, visitproc visit, void *arg) { - Py_VISIT(GETSTATE(m)->error); - return 0; -} - -static int _phono3py_clear(PyObject *m) { - Py_CLEAR(GETSTATE(m)->error); - return 0; -} - -static struct PyModuleDef moduledef = { - PyModuleDef_HEAD_INIT, "_phono3py", NULL, - sizeof(struct module_state), _phono3py_methods, NULL, - _phono3py_traverse, _phono3py_clear, NULL}; - -#define INITERROR return NULL - -PyObject *PyInit__phono3py(void) -#else -#define INITERROR return - -void init_phono3py(void) -#endif -{ -#if PY_MAJOR_VERSION >= 3 - PyObject *module = PyModule_Create(&moduledef); -#else - PyObject *module = Py_InitModule("_phono3py", _phono3py_methods); -#endif - struct module_state *st; - - if (module == NULL) INITERROR; - st = GETSTATE(module); - - st->error = PyErr_NewException("_phono3py.Error", NULL, NULL); - if (st->error == NULL) { - Py_DECREF(module); - INITERROR; - } - -#if PY_MAJOR_VERSION >= 3 - return module; -#endif -} - -static PyObject *py_get_interaction(PyObject *self, PyObject *args) { - PyArrayObject *py_fc3_normal_squared; - PyArrayObject *py_g_zero; - PyArrayObject *py_frequencies; - PyArrayObject *py_eigenvectors; - PyArrayObject *py_triplets; - PyArrayObject *py_bz_grid_addresses; - PyArrayObject *py_D_diag; - PyArrayObject *py_Q; - PyArrayObject *py_svecs; - PyArrayObject *py_multi; - PyArrayObject *py_fc3; - PyArrayObject *py_masses; - PyArrayObject *py_p2s_map; - PyArrayObject *py_s2p_map; - PyArrayObject *py_band_indices; - PyArrayObject *py_all_shortest; - double cutoff_frequency; - long symmetrize_fc3_q; - long make_r0_average; - long openmp_per_triplets; - - Darray *fc3_normal_squared; - Darray *freqs; - _lapack_complex_double *eigvecs; - long(*triplets)[3]; - long num_triplets; - char *g_zero; - long(*bz_grid_addresses)[3]; - long *D_diag; - long(*Q)[3]; - double *fc3; - double(*svecs)[3]; - long(*multi)[2]; - double *masses; - char *all_shortest; - long *p2s; - long *s2p; - long *band_indices; - long multi_dims[2]; - long i; - long is_compact_fc3; - - if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOOOllOdl", &py_fc3_normal_squared, - &py_g_zero, &py_frequencies, &py_eigenvectors, - &py_triplets, &py_bz_grid_addresses, &py_D_diag, - &py_Q, &py_fc3, &py_svecs, &py_multi, &py_masses, - &py_p2s_map, &py_s2p_map, &py_band_indices, - &symmetrize_fc3_q, &make_r0_average, &py_all_shortest, - &cutoff_frequency, &openmp_per_triplets)) { - return NULL; - } - - fc3_normal_squared = convert_to_darray(py_fc3_normal_squared); - freqs = convert_to_darray(py_frequencies); - /* npy_cdouble and lapack_complex_double may not be compatible. */ - /* So eigenvectors should not be used in Python side */ - eigvecs = (_lapack_complex_double *)PyArray_DATA(py_eigenvectors); - triplets = (long(*)[3])PyArray_DATA(py_triplets); - num_triplets = (long)PyArray_DIMS(py_triplets)[0]; - g_zero = (char *)PyArray_DATA(py_g_zero); - bz_grid_addresses = (long(*)[3])PyArray_DATA(py_bz_grid_addresses); - D_diag = (long *)PyArray_DATA(py_D_diag); - Q = (long(*)[3])PyArray_DATA(py_Q); - fc3 = (double *)PyArray_DATA(py_fc3); - if (PyArray_DIMS(py_fc3)[0] == PyArray_DIMS(py_fc3)[1]) { - is_compact_fc3 = 0; - } else { - is_compact_fc3 = 1; - } - svecs = (double(*)[3])PyArray_DATA(py_svecs); - for (i = 0; i < 2; i++) { - multi_dims[i] = PyArray_DIMS(py_multi)[i]; - } - multi = (long(*)[2])PyArray_DATA(py_multi); - masses = (double *)PyArray_DATA(py_masses); - p2s = (long *)PyArray_DATA(py_p2s_map); - s2p = (long *)PyArray_DATA(py_s2p_map); - band_indices = (long *)PyArray_DATA(py_band_indices); - all_shortest = (char *)PyArray_DATA(py_all_shortest); - - ph3py_get_interaction(fc3_normal_squared, g_zero, freqs, eigvecs, triplets, - num_triplets, bz_grid_addresses, D_diag, Q, fc3, - is_compact_fc3, svecs, multi_dims, multi, masses, p2s, - s2p, band_indices, symmetrize_fc3_q, make_r0_average, - all_shortest, cutoff_frequency, openmp_per_triplets); - - free(fc3_normal_squared); - fc3_normal_squared = NULL; - free(freqs); - freqs = NULL; - - Py_RETURN_NONE; -} - -static PyObject *py_get_pp_collision(PyObject *self, PyObject *args) { - PyArrayObject *py_gamma; - PyArrayObject *py_relative_grid_address; - PyArrayObject *py_frequencies; - PyArrayObject *py_eigenvectors; - PyArrayObject *py_triplets; - PyArrayObject *py_triplet_weights; - PyArrayObject *py_bz_grid_addresses; - PyArrayObject *py_bz_map; - PyArrayObject *py_D_diag; - PyArrayObject *py_Q; - PyArrayObject *py_fc3; - PyArrayObject *py_svecs; - PyArrayObject *py_multi; - PyArrayObject *py_masses; - PyArrayObject *py_p2s_map; - PyArrayObject *py_s2p_map; - PyArrayObject *py_band_indices; - PyArrayObject *py_temperatures; - PyArrayObject *py_all_shortest; - double cutoff_frequency; - long is_NU; - long symmetrize_fc3_q; - long make_r0_average; - long bz_grid_type; - long openmp_per_triplets; - - double *gamma; - long(*relative_grid_address)[4][3]; - double *frequencies; - _lapack_complex_double *eigenvectors; - long(*triplets)[3]; - long num_triplets; - long *triplet_weights; - long(*bz_grid_addresses)[3]; - long *bz_map; - long *D_diag; - long(*Q)[3]; - double *fc3; - double(*svecs)[3]; - long(*multi)[2]; - double *masses; - long *p2s; - long *s2p; - Larray *band_indices; - Darray *temperatures; - char *all_shortest; - long multi_dims[2]; - long i; - long is_compact_fc3; - - if (!PyArg_ParseTuple( - args, "OOOOOOOOlOOOOOOOOOOlllOdl", &py_gamma, - &py_relative_grid_address, &py_frequencies, &py_eigenvectors, - &py_triplets, &py_triplet_weights, &py_bz_grid_addresses, - &py_bz_map, &bz_grid_type, &py_D_diag, &py_Q, &py_fc3, &py_svecs, - &py_multi, &py_masses, &py_p2s_map, &py_s2p_map, &py_band_indices, - &py_temperatures, &is_NU, &symmetrize_fc3_q, &make_r0_average, - &py_all_shortest, &cutoff_frequency, &openmp_per_triplets)) { - return NULL; - } - - gamma = (double *)PyArray_DATA(py_gamma); - relative_grid_address = - (long(*)[4][3])PyArray_DATA(py_relative_grid_address); - frequencies = (double *)PyArray_DATA(py_frequencies); - eigenvectors = (_lapack_complex_double *)PyArray_DATA(py_eigenvectors); - triplets = (long(*)[3])PyArray_DATA(py_triplets); - num_triplets = (long)PyArray_DIMS(py_triplets)[0]; - triplet_weights = (long *)PyArray_DATA(py_triplet_weights); - bz_grid_addresses = (long(*)[3])PyArray_DATA(py_bz_grid_addresses); - bz_map = (long *)PyArray_DATA(py_bz_map); - D_diag = (long *)PyArray_DATA(py_D_diag); - Q = (long(*)[3])PyArray_DATA(py_Q); - fc3 = (double *)PyArray_DATA(py_fc3); - if (PyArray_DIMS(py_fc3)[0] == PyArray_DIMS(py_fc3)[1]) { - is_compact_fc3 = 0; - } else { - is_compact_fc3 = 1; - } - svecs = (double(*)[3])PyArray_DATA(py_svecs); - for (i = 0; i < 2; i++) { - multi_dims[i] = PyArray_DIMS(py_multi)[i]; - } - multi = (long(*)[2])PyArray_DATA(py_multi); - masses = (double *)PyArray_DATA(py_masses); - p2s = (long *)PyArray_DATA(py_p2s_map); - s2p = (long *)PyArray_DATA(py_s2p_map); - band_indices = convert_to_larray(py_band_indices); - temperatures = convert_to_darray(py_temperatures); - all_shortest = (char *)PyArray_DATA(py_all_shortest); - - ph3py_get_pp_collision( - gamma, relative_grid_address, frequencies, eigenvectors, triplets, - num_triplets, triplet_weights, bz_grid_addresses, bz_map, bz_grid_type, - D_diag, Q, fc3, is_compact_fc3, svecs, multi_dims, multi, masses, p2s, - s2p, band_indices, temperatures, is_NU, symmetrize_fc3_q, - make_r0_average, all_shortest, cutoff_frequency, openmp_per_triplets); - - free(band_indices); - band_indices = NULL; - free(temperatures); - temperatures = NULL; - - Py_RETURN_NONE; -} - -static PyObject *py_get_pp_collision_with_sigma(PyObject *self, - PyObject *args) { - PyArrayObject *py_gamma; - PyArrayObject *py_frequencies; - PyArrayObject *py_eigenvectors; - PyArrayObject *py_triplets; - PyArrayObject *py_triplet_weights; - PyArrayObject *py_bz_grid_addresses; - PyArrayObject *py_D_diag; - PyArrayObject *py_Q; - PyArrayObject *py_fc3; - PyArrayObject *py_svecs; - PyArrayObject *py_multi; - PyArrayObject *py_masses; - PyArrayObject *py_p2s_map; - PyArrayObject *py_s2p_map; - PyArrayObject *py_band_indices; - PyArrayObject *py_temperatures; - PyArrayObject *py_all_shortest; - long is_NU; - long symmetrize_fc3_q; - double sigma; - double sigma_cutoff; - long make_r0_average; - double cutoff_frequency; - long openmp_per_triplets; - - double *gamma; - double *frequencies; - _lapack_complex_double *eigenvectors; - long(*triplets)[3]; - long num_triplets; - long *triplet_weights; - long(*bz_grid_addresses)[3]; - long *D_diag; - long(*Q)[3]; - double *fc3; - double(*svecs)[3]; - long(*multi)[2]; - double *masses; - long *p2s; - long *s2p; - Larray *band_indices; - Darray *temperatures; - char *all_shortest; - long multi_dims[2]; - long i; - long is_compact_fc3; - - if (!PyArg_ParseTuple( - args, "OddOOOOOOOOOOOOOOOlllOdl", &py_gamma, &sigma, &sigma_cutoff, - &py_frequencies, &py_eigenvectors, &py_triplets, - &py_triplet_weights, &py_bz_grid_addresses, &py_D_diag, &py_Q, - &py_fc3, &py_svecs, &py_multi, &py_masses, &py_p2s_map, &py_s2p_map, - &py_band_indices, &py_temperatures, &is_NU, &symmetrize_fc3_q, - &make_r0_average, &py_all_shortest, &cutoff_frequency, - &openmp_per_triplets)) { - return NULL; - } - - gamma = (double *)PyArray_DATA(py_gamma); - frequencies = (double *)PyArray_DATA(py_frequencies); - eigenvectors = (_lapack_complex_double *)PyArray_DATA(py_eigenvectors); - triplets = (long(*)[3])PyArray_DATA(py_triplets); - num_triplets = (long)PyArray_DIMS(py_triplets)[0]; - triplet_weights = (long *)PyArray_DATA(py_triplet_weights); - bz_grid_addresses = (long(*)[3])PyArray_DATA(py_bz_grid_addresses); - D_diag = (long *)PyArray_DATA(py_D_diag); - Q = (long(*)[3])PyArray_DATA(py_Q); - fc3 = (double *)PyArray_DATA(py_fc3); - if (PyArray_DIMS(py_fc3)[0] == PyArray_DIMS(py_fc3)[1]) { - is_compact_fc3 = 0; - } else { - is_compact_fc3 = 1; - } - svecs = (double(*)[3])PyArray_DATA(py_svecs); - for (i = 0; i < 2; i++) { - multi_dims[i] = PyArray_DIMS(py_multi)[i]; - } - multi = (long(*)[2])PyArray_DATA(py_multi); - masses = (double *)PyArray_DATA(py_masses); - p2s = (long *)PyArray_DATA(py_p2s_map); - s2p = (long *)PyArray_DATA(py_s2p_map); - band_indices = convert_to_larray(py_band_indices); - temperatures = convert_to_darray(py_temperatures); - all_shortest = (char *)PyArray_DATA(py_all_shortest); - - ph3py_get_pp_collision_with_sigma( - gamma, sigma, sigma_cutoff, frequencies, eigenvectors, triplets, - num_triplets, triplet_weights, bz_grid_addresses, D_diag, Q, fc3, - is_compact_fc3, svecs, multi_dims, multi, masses, p2s, s2p, - band_indices, temperatures, is_NU, symmetrize_fc3_q, make_r0_average, - all_shortest, cutoff_frequency, openmp_per_triplets); - - free(band_indices); - band_indices = NULL; - free(temperatures); - temperatures = NULL; - - Py_RETURN_NONE; -} - -static PyObject *py_get_imag_self_energy_with_g(PyObject *self, - PyObject *args) { - PyArrayObject *py_gamma; - PyArrayObject *py_fc3_normal_squared; - PyArrayObject *py_frequencies; - PyArrayObject *py_triplets; - PyArrayObject *py_triplet_weights; - PyArrayObject *py_g; - PyArrayObject *py_g_zero; - double cutoff_frequency, temperature; - long frequency_point_index; - - Darray *fc3_normal_squared; - double *gamma; - double *g; - char *g_zero; - double *frequencies; - long(*triplets)[3]; - long *triplet_weights; - long num_frequency_points; - - if (!PyArg_ParseTuple(args, "OOOOOdOOdl", &py_gamma, &py_fc3_normal_squared, - &py_triplets, &py_triplet_weights, &py_frequencies, - &temperature, &py_g, &py_g_zero, &cutoff_frequency, - &frequency_point_index)) { - return NULL; - } - - fc3_normal_squared = convert_to_darray(py_fc3_normal_squared); - gamma = (double *)PyArray_DATA(py_gamma); - g = (double *)PyArray_DATA(py_g); - g_zero = (char *)PyArray_DATA(py_g_zero); - frequencies = (double *)PyArray_DATA(py_frequencies); - triplets = (long(*)[3])PyArray_DATA(py_triplets); - triplet_weights = (long *)PyArray_DATA(py_triplet_weights); - num_frequency_points = (long)PyArray_DIMS(py_g)[2]; - - ph3py_get_imag_self_energy_at_bands_with_g( - gamma, fc3_normal_squared, frequencies, triplets, triplet_weights, g, - g_zero, temperature, cutoff_frequency, num_frequency_points, - frequency_point_index); - - free(fc3_normal_squared); - fc3_normal_squared = NULL; - - Py_RETURN_NONE; -} - -static PyObject *py_get_detailed_imag_self_energy_with_g(PyObject *self, - PyObject *args) { - PyArrayObject *py_gamma_detail; - PyArrayObject *py_gamma_N; - PyArrayObject *py_gamma_U; - PyArrayObject *py_fc3_normal_squared; - PyArrayObject *py_frequencies; - PyArrayObject *py_triplets; - PyArrayObject *py_triplet_weights; - PyArrayObject *py_bz_grid_addresses; - PyArrayObject *py_g; - PyArrayObject *py_g_zero; - double cutoff_frequency, temperature; - - Darray *fc3_normal_squared; - double *gamma_detail; - double *gamma_N; - double *gamma_U; - double *g; - char *g_zero; - double *frequencies; - long(*triplets)[3]; - long *triplet_weights; - long(*bz_grid_addresses)[3]; - - if (!PyArg_ParseTuple(args, "OOOOOOOOdOOd", &py_gamma_detail, &py_gamma_N, - &py_gamma_U, &py_fc3_normal_squared, &py_triplets, - &py_triplet_weights, &py_bz_grid_addresses, - &py_frequencies, &temperature, &py_g, &py_g_zero, - &cutoff_frequency)) { - return NULL; - } - - fc3_normal_squared = convert_to_darray(py_fc3_normal_squared); - gamma_detail = (double *)PyArray_DATA(py_gamma_detail); - gamma_N = (double *)PyArray_DATA(py_gamma_N); - gamma_U = (double *)PyArray_DATA(py_gamma_U); - g = (double *)PyArray_DATA(py_g); - g_zero = (char *)PyArray_DATA(py_g_zero); - frequencies = (double *)PyArray_DATA(py_frequencies); - triplets = (long(*)[3])PyArray_DATA(py_triplets); - triplet_weights = (long *)PyArray_DATA(py_triplet_weights); - bz_grid_addresses = (long(*)[3])PyArray_DATA(py_bz_grid_addresses); - - ph3py_get_detailed_imag_self_energy_at_bands_with_g( - gamma_detail, gamma_N, gamma_U, fc3_normal_squared, frequencies, - triplets, triplet_weights, bz_grid_addresses, g, g_zero, temperature, - cutoff_frequency); - - free(fc3_normal_squared); - fc3_normal_squared = NULL; - - Py_RETURN_NONE; -} - -static PyObject *py_get_real_self_energy_at_bands(PyObject *self, - PyObject *args) { - PyArrayObject *py_shift; - PyArrayObject *py_fc3_normal_squared; - PyArrayObject *py_frequencies; - PyArrayObject *py_triplets; - PyArrayObject *py_triplet_weights; - PyArrayObject *py_band_indices; - double epsilon, unit_conversion_factor, cutoff_frequency, temperature; - - Darray *fc3_normal_squared; - double *shift; - double *frequencies; - long *band_indices; - long(*triplets)[3]; - long *triplet_weights; - - if (!PyArg_ParseTuple(args, "OOOOOOdddd", &py_shift, &py_fc3_normal_squared, - &py_triplets, &py_triplet_weights, &py_frequencies, - &py_band_indices, &temperature, &epsilon, - &unit_conversion_factor, &cutoff_frequency)) { - return NULL; - } - - fc3_normal_squared = convert_to_darray(py_fc3_normal_squared); - shift = (double *)PyArray_DATA(py_shift); - frequencies = (double *)PyArray_DATA(py_frequencies); - band_indices = (long *)PyArray_DATA(py_band_indices); - triplets = (long(*)[3])PyArray_DATA(py_triplets); - triplet_weights = (long *)PyArray_DATA(py_triplet_weights); - - ph3py_get_real_self_energy_at_bands( - shift, fc3_normal_squared, band_indices, frequencies, triplets, - triplet_weights, epsilon, temperature, unit_conversion_factor, - cutoff_frequency); - - free(fc3_normal_squared); - fc3_normal_squared = NULL; - - Py_RETURN_NONE; -} - -static PyObject *py_get_real_self_energy_at_frequency_point(PyObject *self, - PyObject *args) { - PyArrayObject *py_shift; - PyArrayObject *py_fc3_normal_squared; - PyArrayObject *py_frequencies; - PyArrayObject *py_triplets; - PyArrayObject *py_triplet_weights; - PyArrayObject *py_band_indices; - double frequency_point, epsilon, unit_conversion_factor, cutoff_frequency; - double temperature; - - Darray *fc3_normal_squared; - double *shift; - double *frequencies; - long *band_indices; - long(*triplets)[3]; - long *triplet_weights; - - if (!PyArg_ParseTuple(args, "OdOOOOOdddd", &py_shift, &frequency_point, - &py_fc3_normal_squared, &py_triplets, - &py_triplet_weights, &py_frequencies, - &py_band_indices, &temperature, &epsilon, - &unit_conversion_factor, &cutoff_frequency)) { - return NULL; - } - - fc3_normal_squared = convert_to_darray(py_fc3_normal_squared); - shift = (double *)PyArray_DATA(py_shift); - frequencies = (double *)PyArray_DATA(py_frequencies); - band_indices = (long *)PyArray_DATA(py_band_indices); - triplets = (long(*)[3])PyArray_DATA(py_triplets); - triplet_weights = (long *)PyArray_DATA(py_triplet_weights); - - ph3py_get_real_self_energy_at_frequency_point( - shift, frequency_point, fc3_normal_squared, band_indices, frequencies, - triplets, triplet_weights, epsilon, temperature, unit_conversion_factor, - cutoff_frequency); - - free(fc3_normal_squared); - fc3_normal_squared = NULL; - - Py_RETURN_NONE; -} - -static PyObject *py_get_collision_matrix(PyObject *self, PyObject *args) { - PyArrayObject *py_collision_matrix; - PyArrayObject *py_fc3_normal_squared; - PyArrayObject *py_frequencies; - PyArrayObject *py_triplets; - PyArrayObject *py_triplets_map; - PyArrayObject *py_map_q; - PyArrayObject *py_g; - PyArrayObject *py_rotated_grid_points; - PyArrayObject *py_rotations_cartesian; - double temperature, unit_conversion_factor, cutoff_frequency; - - Darray *fc3_normal_squared; - double *collision_matrix; - double *g; - double *frequencies; - long(*triplets)[3]; - long *triplets_map; - long *map_q; - long *rotated_grid_points; - long num_gp, num_ir_gp, num_rot; - double *rotations_cartesian; - - if (!PyArg_ParseTuple( - args, "OOOOOOOOOddd", &py_collision_matrix, &py_fc3_normal_squared, - &py_frequencies, &py_g, &py_triplets, &py_triplets_map, &py_map_q, - &py_rotated_grid_points, &py_rotations_cartesian, &temperature, - &unit_conversion_factor, &cutoff_frequency)) { - return NULL; - } - - fc3_normal_squared = convert_to_darray(py_fc3_normal_squared); - collision_matrix = (double *)PyArray_DATA(py_collision_matrix); - g = (double *)PyArray_DATA(py_g); - frequencies = (double *)PyArray_DATA(py_frequencies); - triplets = (long(*)[3])PyArray_DATA(py_triplets); - triplets_map = (long *)PyArray_DATA(py_triplets_map); - num_gp = (long)PyArray_DIMS(py_triplets_map)[0]; - map_q = (long *)PyArray_DATA(py_map_q); - rotated_grid_points = (long *)PyArray_DATA(py_rotated_grid_points); - num_ir_gp = (long)PyArray_DIMS(py_rotated_grid_points)[0]; - num_rot = (long)PyArray_DIMS(py_rotated_grid_points)[1]; - rotations_cartesian = (double *)PyArray_DATA(py_rotations_cartesian); - - assert(num_rot == PyArray_DIMS(py_rotations_cartesian)[0]); - assert(num_gp == PyArray_DIMS(py_frequencies)[0]); - - ph3py_get_collision_matrix(collision_matrix, fc3_normal_squared, - frequencies, triplets, triplets_map, map_q, - rotated_grid_points, rotations_cartesian, g, - num_ir_gp, num_gp, num_rot, temperature, - unit_conversion_factor, cutoff_frequency); - - free(fc3_normal_squared); - fc3_normal_squared = NULL; - - Py_RETURN_NONE; -} - -static PyObject *py_get_reducible_collision_matrix(PyObject *self, - PyObject *args) { - PyArrayObject *py_collision_matrix; - PyArrayObject *py_fc3_normal_squared; - PyArrayObject *py_frequencies; - PyArrayObject *py_triplets; - PyArrayObject *py_triplets_map; - PyArrayObject *py_map_q; - PyArrayObject *py_g; - double temperature, unit_conversion_factor, cutoff_frequency; - - Darray *fc3_normal_squared; - double *collision_matrix; - double *g; - double *frequencies; - long(*triplets)[3]; - long *triplets_map; - long num_gp; - long *map_q; - - if (!PyArg_ParseTuple( - args, "OOOOOOOddd", &py_collision_matrix, &py_fc3_normal_squared, - &py_frequencies, &py_g, &py_triplets, &py_triplets_map, &py_map_q, - &temperature, &unit_conversion_factor, &cutoff_frequency)) { - return NULL; - } - - fc3_normal_squared = convert_to_darray(py_fc3_normal_squared); - collision_matrix = (double *)PyArray_DATA(py_collision_matrix); - g = (double *)PyArray_DATA(py_g); - frequencies = (double *)PyArray_DATA(py_frequencies); - triplets = (long(*)[3])PyArray_DATA(py_triplets); - triplets_map = (long *)PyArray_DATA(py_triplets_map); - num_gp = (long)PyArray_DIMS(py_triplets_map)[0]; - map_q = (long *)PyArray_DATA(py_map_q); - - ph3py_get_reducible_collision_matrix( - collision_matrix, fc3_normal_squared, frequencies, triplets, - triplets_map, map_q, g, num_gp, temperature, unit_conversion_factor, - cutoff_frequency); - - free(fc3_normal_squared); - fc3_normal_squared = NULL; - - Py_RETURN_NONE; -} - -static PyObject *py_symmetrize_collision_matrix(PyObject *self, - PyObject *args) { - PyArrayObject *py_collision_matrix; - - double *collision_matrix; - long num_band, num_grid_points, num_temp, num_sigma; - long num_column; - - if (!PyArg_ParseTuple(args, "O", &py_collision_matrix)) { - return NULL; - } - - collision_matrix = (double *)PyArray_DATA(py_collision_matrix); - num_sigma = (long)PyArray_DIMS(py_collision_matrix)[0]; - num_temp = (long)PyArray_DIMS(py_collision_matrix)[1]; - num_grid_points = (long)PyArray_DIMS(py_collision_matrix)[2]; - num_band = (long)PyArray_DIMS(py_collision_matrix)[3]; - - if (PyArray_NDIM(py_collision_matrix) == 8) { - num_column = num_grid_points * num_band * 3; - } else { - num_column = num_grid_points * num_band; - } - - ph3py_symmetrize_collision_matrix(collision_matrix, num_column, num_temp, - num_sigma); - - Py_RETURN_NONE; -} - -static PyObject *py_expand_collision_matrix(PyObject *self, PyObject *args) { - PyArrayObject *py_collision_matrix; - PyArrayObject *py_ir_grid_points; - PyArrayObject *py_rot_grid_points; - - double *collision_matrix; - long *rot_grid_points; - long *ir_grid_points; - long num_band, num_grid_points, num_temp, num_sigma, num_rot, num_ir_gp; - - if (!PyArg_ParseTuple(args, "OOO", &py_collision_matrix, &py_ir_grid_points, - &py_rot_grid_points)) { - return NULL; - } - - collision_matrix = (double *)PyArray_DATA(py_collision_matrix); - rot_grid_points = (long *)PyArray_DATA(py_rot_grid_points); - ir_grid_points = (long *)PyArray_DATA(py_ir_grid_points); - num_sigma = (long)PyArray_DIMS(py_collision_matrix)[0]; - num_temp = (long)PyArray_DIMS(py_collision_matrix)[1]; - num_grid_points = (long)PyArray_DIMS(py_collision_matrix)[2]; - num_band = (long)PyArray_DIMS(py_collision_matrix)[3]; - num_rot = (long)PyArray_DIMS(py_rot_grid_points)[0]; - num_ir_gp = (long)PyArray_DIMS(py_ir_grid_points)[0]; - - ph3py_expand_collision_matrix(collision_matrix, rot_grid_points, - ir_grid_points, num_ir_gp, num_grid_points, - num_rot, num_sigma, num_temp, num_band); - - Py_RETURN_NONE; -} - -static PyObject *py_get_isotope_strength(PyObject *self, PyObject *args) { - PyArrayObject *py_gamma; - PyArrayObject *py_frequencies; - PyArrayObject *py_eigenvectors; - PyArrayObject *py_band_indices; - PyArrayObject *py_mass_variances; - PyArrayObject *py_ir_grid_points; - PyArrayObject *py_weights; - - long grid_point; - double cutoff_frequency; - double sigma; - - double *gamma; - double *frequencies; - long *ir_grid_points; - double *weights; - _lapack_complex_double *eigenvectors; - long *band_indices; - double *mass_variances; - long num_band, num_band0, num_ir_grid_points; - - if (!PyArg_ParseTuple(args, "OlOOOOOOdd", &py_gamma, &grid_point, - &py_ir_grid_points, &py_weights, &py_mass_variances, - &py_frequencies, &py_eigenvectors, &py_band_indices, - &sigma, &cutoff_frequency)) { - return NULL; - } - - gamma = (double *)PyArray_DATA(py_gamma); - frequencies = (double *)PyArray_DATA(py_frequencies); - eigenvectors = (_lapack_complex_double *)PyArray_DATA(py_eigenvectors); - ir_grid_points = (long *)PyArray_DATA(py_ir_grid_points); - weights = (double *)PyArray_DATA(py_weights); - band_indices = (long *)PyArray_DATA(py_band_indices); - mass_variances = (double *)PyArray_DATA(py_mass_variances); - num_band = (long)PyArray_DIMS(py_frequencies)[1]; - num_band0 = (long)PyArray_DIMS(py_band_indices)[0]; - num_ir_grid_points = (long)PyArray_DIMS(py_ir_grid_points)[0]; - - ph3py_get_isotope_scattering_strength( - gamma, grid_point, ir_grid_points, weights, mass_variances, frequencies, - eigenvectors, num_ir_grid_points, band_indices, num_band, num_band0, - sigma, cutoff_frequency); - - Py_RETURN_NONE; -} - -static PyObject *py_get_thm_isotope_strength(PyObject *self, PyObject *args) { - PyArrayObject *py_gamma; - PyArrayObject *py_frequencies; - PyArrayObject *py_eigenvectors; - PyArrayObject *py_band_indices; - PyArrayObject *py_mass_variances; - PyArrayObject *py_ir_grid_points; - PyArrayObject *py_weights; - PyArrayObject *py_integration_weights; - long grid_point; - double cutoff_frequency; - - double *gamma; - double *frequencies; - long *ir_grid_points; - double *weights; - _lapack_complex_double *eigenvectors; - long *band_indices; - double *mass_variances; - long num_band, num_band0, num_ir_grid_points; - double *integration_weights; - - if (!PyArg_ParseTuple(args, "OlOOOOOOOd", &py_gamma, &grid_point, - &py_ir_grid_points, &py_weights, &py_mass_variances, - &py_frequencies, &py_eigenvectors, &py_band_indices, - &py_integration_weights, &cutoff_frequency)) { - return NULL; - } - - gamma = (double *)PyArray_DATA(py_gamma); - frequencies = (double *)PyArray_DATA(py_frequencies); - ir_grid_points = (long *)PyArray_DATA(py_ir_grid_points); - weights = (double *)PyArray_DATA(py_weights); - eigenvectors = (_lapack_complex_double *)PyArray_DATA(py_eigenvectors); - band_indices = (long *)PyArray_DATA(py_band_indices); - mass_variances = (double *)PyArray_DATA(py_mass_variances); - num_band = (long)PyArray_DIMS(py_frequencies)[1]; - num_band0 = (long)PyArray_DIMS(py_band_indices)[0]; - integration_weights = (double *)PyArray_DATA(py_integration_weights); - num_ir_grid_points = (long)PyArray_DIMS(py_ir_grid_points)[0]; - - ph3py_get_thm_isotope_scattering_strength( - gamma, grid_point, ir_grid_points, weights, mass_variances, frequencies, - eigenvectors, num_ir_grid_points, band_indices, num_band, num_band0, - integration_weights, cutoff_frequency); - - Py_RETURN_NONE; -} - -static PyObject *py_distribute_fc3(PyObject *self, PyObject *args) { - PyArrayObject *force_constants_third; - long target; - long source; - PyArrayObject *rotation_cart_inv; - PyArrayObject *atom_mapping_py; - - double *fc3; - double *rot_cart_inv; - long *atom_mapping; - long num_atom; - - if (!PyArg_ParseTuple(args, "OllOO", &force_constants_third, &target, - &source, &atom_mapping_py, &rotation_cart_inv)) { - return NULL; - } - - fc3 = (double *)PyArray_DATA(force_constants_third); - rot_cart_inv = (double *)PyArray_DATA(rotation_cart_inv); - atom_mapping = (long *)PyArray_DATA(atom_mapping_py); - num_atom = (long)PyArray_DIMS(atom_mapping_py)[0]; - - ph3py_distribute_fc3(fc3, target, source, atom_mapping, num_atom, - rot_cart_inv); - - Py_RETURN_NONE; -} - -static PyObject *py_rotate_delta_fc2s(PyObject *self, PyObject *args) { - PyArrayObject *py_fc3; - PyArrayObject *py_delta_fc2s; - PyArrayObject *py_inv_U; - PyArrayObject *py_site_sym_cart; - PyArrayObject *py_rot_map_syms; - - double(*fc3)[3][3][3]; - double(*delta_fc2s)[3][3]; - double *inv_U; - double(*site_sym_cart)[3][3]; - long *rot_map_syms; - long num_atom, num_disp, num_site_sym; - - if (!PyArg_ParseTuple(args, "OOOOO", &py_fc3, &py_delta_fc2s, &py_inv_U, - &py_site_sym_cart, &py_rot_map_syms)) { - return NULL; - } - - /* (num_atom, num_atom, 3, 3, 3) */ - fc3 = (double(*)[3][3][3])PyArray_DATA(py_fc3); - /* (n_u1, num_atom, num_atom, 3, 3) */ - delta_fc2s = (double(*)[3][3])PyArray_DATA(py_delta_fc2s); - /* (3, n_u1 * n_sym) */ - inv_U = (double *)PyArray_DATA(py_inv_U); - /* (n_sym, 3, 3) */ - site_sym_cart = (double(*)[3][3])PyArray_DATA(py_site_sym_cart); - /* (n_sym, natom) */ - rot_map_syms = (long *)PyArray_DATA(py_rot_map_syms); - - num_atom = (long)PyArray_DIMS(py_fc3)[0]; - num_disp = (long)PyArray_DIMS(py_delta_fc2s)[0]; - num_site_sym = (long)PyArray_DIMS(py_site_sym_cart)[0]; - - ph3py_rotate_delta_fc2(fc3, delta_fc2s, inv_U, site_sym_cart, rot_map_syms, - num_atom, num_site_sym, num_disp); - - Py_RETURN_NONE; -} - -static PyObject *py_get_permutation_symmetry_fc3(PyObject *self, - PyObject *args) { - PyArrayObject *py_fc3; - - double *fc3; - long num_atom; - - if (!PyArg_ParseTuple(args, "O", &py_fc3)) { - return NULL; - } - - fc3 = (double *)PyArray_DATA(py_fc3); - num_atom = (long)PyArray_DIMS(py_fc3)[0]; - - ph3py_get_permutation_symmetry_fc3(fc3, num_atom); - - Py_RETURN_NONE; -} - -static PyObject *py_get_permutation_symmetry_compact_fc3(PyObject *self, - PyObject *args) { - PyArrayObject *py_fc3; - PyArrayObject *py_permutations; - PyArrayObject *py_s2pp_map; - PyArrayObject *py_p2s_map; - PyArrayObject *py_nsym_list; - - double *fc3; - long *s2pp; - long *p2s; - long *nsym_list; - long *perms; - long n_patom, n_satom; - - if (!PyArg_ParseTuple(args, "OOOOO", &py_fc3, &py_permutations, - &py_s2pp_map, &py_p2s_map, &py_nsym_list)) { - return NULL; - } - - fc3 = (double *)PyArray_DATA(py_fc3); - perms = (long *)PyArray_DATA(py_permutations); - s2pp = (long *)PyArray_DATA(py_s2pp_map); - p2s = (long *)PyArray_DATA(py_p2s_map); - nsym_list = (long *)PyArray_DATA(py_nsym_list); - n_patom = (long)PyArray_DIMS(py_fc3)[0]; - n_satom = (long)PyArray_DIMS(py_fc3)[1]; - - ph3py_get_permutation_symmetry_compact_fc3(fc3, p2s, s2pp, nsym_list, perms, - n_satom, n_patom); - - Py_RETURN_NONE; -} - -static PyObject *py_transpose_compact_fc3(PyObject *self, PyObject *args) { - PyArrayObject *py_fc3; - PyArrayObject *py_permutations; - PyArrayObject *py_s2pp_map; - PyArrayObject *py_p2s_map; - PyArrayObject *py_nsym_list; - long t_type; - - double *fc3; - long *s2pp; - long *p2s; - long *nsym_list; - long *perms; - long n_patom, n_satom; - - if (!PyArg_ParseTuple(args, "OOOOOl", &py_fc3, &py_permutations, - &py_s2pp_map, &py_p2s_map, &py_nsym_list, &t_type)) { - return NULL; - } - - fc3 = (double *)PyArray_DATA(py_fc3); - perms = (long *)PyArray_DATA(py_permutations); - s2pp = (long *)PyArray_DATA(py_s2pp_map); - p2s = (long *)PyArray_DATA(py_p2s_map); - nsym_list = (long *)PyArray_DATA(py_nsym_list); - n_patom = (long)PyArray_DIMS(py_fc3)[0]; - n_satom = (long)PyArray_DIMS(py_fc3)[1]; - - ph3py_transpose_compact_fc3(fc3, p2s, s2pp, nsym_list, perms, n_satom, - n_patom, t_type); - - Py_RETURN_NONE; -} - -static PyObject *py_get_thm_relative_grid_address(PyObject *self, - PyObject *args) { - PyArrayObject *py_relative_grid_address; - PyArrayObject *py_reciprocal_lattice_py; - - long(*relative_grid_address)[4][3]; - double(*reciprocal_lattice)[3]; - - if (!PyArg_ParseTuple(args, "OO", &py_relative_grid_address, - &py_reciprocal_lattice_py)) { - return NULL; - } - - relative_grid_address = - (long(*)[4][3])PyArray_DATA(py_relative_grid_address); - reciprocal_lattice = (double(*)[3])PyArray_DATA(py_reciprocal_lattice_py); - - ph3py_get_relative_grid_address(relative_grid_address, reciprocal_lattice); - - Py_RETURN_NONE; -} - -static PyObject *py_get_neighboring_grid_points(PyObject *self, - PyObject *args) { - PyArrayObject *py_relative_grid_points; - PyArrayObject *py_grid_points; - PyArrayObject *py_relative_grid_address; - PyArrayObject *py_D_diag; - PyArrayObject *py_bz_grid_address; - PyArrayObject *py_bz_map; - long bz_grid_type; - - long *relative_grid_points; - long *grid_points; - long num_grid_points, num_relative_grid_address; - long(*relative_grid_address)[3]; - long *D_diag; - long(*bz_grid_address)[3]; - long *bz_map; - - if (!PyArg_ParseTuple(args, "OOOOOOl", &py_relative_grid_points, - &py_grid_points, &py_relative_grid_address, - &py_D_diag, &py_bz_grid_address, &py_bz_map, - &bz_grid_type)) { - return NULL; - } - - relative_grid_points = (long *)PyArray_DATA(py_relative_grid_points); - grid_points = (long *)PyArray_DATA(py_grid_points); - num_grid_points = (long)PyArray_DIMS(py_grid_points)[0]; - relative_grid_address = (long(*)[3])PyArray_DATA(py_relative_grid_address); - num_relative_grid_address = (long)PyArray_DIMS(py_relative_grid_address)[0]; - D_diag = (long *)PyArray_DATA(py_D_diag); - bz_grid_address = (long(*)[3])PyArray_DATA(py_bz_grid_address); - bz_map = (long *)PyArray_DATA(py_bz_map); - - ph3py_get_neighboring_gird_points( - relative_grid_points, grid_points, relative_grid_address, D_diag, - bz_grid_address, bz_map, bz_grid_type, num_grid_points, - num_relative_grid_address); - - Py_RETURN_NONE; -} - -static PyObject *py_get_thm_integration_weights_at_grid_points(PyObject *self, - PyObject *args) { - PyArrayObject *py_iw; - PyArrayObject *py_frequency_points; - PyArrayObject *py_relative_grid_address; - PyArrayObject *py_D_diag; - PyArrayObject *py_grid_points; - PyArrayObject *py_frequencies; - PyArrayObject *py_bz_grid_address; - PyArrayObject *py_gp2irgp_map; - PyArrayObject *py_bz_map; - long bz_grid_type; - char *function; - - double *iw; - double *frequency_points; - long num_frequency_points, num_band, num_gp; - long(*relative_grid_address)[4][3]; - long *D_diag; - long *grid_points; - long(*bz_grid_address)[3]; - long *bz_map; - long *gp2irgp_map; - double *frequencies; - - if (!PyArg_ParseTuple(args, "OOOOOOOOOls", &py_iw, &py_frequency_points, - &py_relative_grid_address, &py_D_diag, - &py_grid_points, &py_frequencies, &py_bz_grid_address, - &py_bz_map, &py_gp2irgp_map, &bz_grid_type, - &function)) { - return NULL; - } - - iw = (double *)PyArray_DATA(py_iw); - frequency_points = (double *)PyArray_DATA(py_frequency_points); - num_frequency_points = (long)PyArray_DIMS(py_frequency_points)[0]; - relative_grid_address = - (long(*)[4][3])PyArray_DATA(py_relative_grid_address); - D_diag = (long *)PyArray_DATA(py_D_diag); - grid_points = (long *)PyArray_DATA(py_grid_points); - num_gp = (long)PyArray_DIMS(py_grid_points)[0]; - bz_grid_address = (long(*)[3])PyArray_DATA(py_bz_grid_address); - bz_map = (long *)PyArray_DATA(py_bz_map); - gp2irgp_map = (long *)PyArray_DATA(py_gp2irgp_map); - frequencies = (double *)PyArray_DATA(py_frequencies); - num_band = (long)PyArray_DIMS(py_frequencies)[1]; - - ph3py_get_thm_integration_weights_at_grid_points( - iw, frequency_points, num_frequency_points, num_band, num_gp, - relative_grid_address, D_diag, grid_points, bz_grid_address, bz_map, - bz_grid_type, frequencies, gp2irgp_map, function[0]); - - Py_RETURN_NONE; -} - -static PyObject *py_tpl_get_triplets_reciprocal_mesh_at_q(PyObject *self, - PyObject *args) { - PyArrayObject *py_map_triplets; - PyArrayObject *py_map_q; - PyArrayObject *py_D_diag; - PyArrayObject *py_rotations; - long fixed_grid_number; - long is_time_reversal; - long swappable; - - long *map_triplets; - long *map_q; - long *D_diag; - long(*rot)[3][3]; - long num_rot; - long num_ir; - - if (!PyArg_ParseTuple(args, "OOlOlOl", &py_map_triplets, &py_map_q, - &fixed_grid_number, &py_D_diag, &is_time_reversal, - &py_rotations, &swappable)) { - return NULL; - } - - map_triplets = (long *)PyArray_DATA(py_map_triplets); - map_q = (long *)PyArray_DATA(py_map_q); - D_diag = (long *)PyArray_DATA(py_D_diag); - rot = (long(*)[3][3])PyArray_DATA(py_rotations); - num_rot = (long)PyArray_DIMS(py_rotations)[0]; - num_ir = ph3py_get_triplets_reciprocal_mesh_at_q( - map_triplets, map_q, fixed_grid_number, D_diag, is_time_reversal, - num_rot, rot, swappable); - - return PyLong_FromLong(num_ir); -} - -static PyObject *py_tpl_get_BZ_triplets_at_q(PyObject *self, PyObject *args) { - PyArrayObject *py_triplets; - PyArrayObject *py_bz_grid_address; - PyArrayObject *py_bz_map; - PyArrayObject *py_map_triplets; - PyArrayObject *py_D_diag; - PyArrayObject *py_Q; - long grid_point; - long bz_grid_type; - - long(*triplets)[3]; - long(*bz_grid_address)[3]; - long *bz_map; - long *map_triplets; - long num_map_triplets; - long *D_diag; - long(*Q)[3]; - long num_ir; - - if (!PyArg_ParseTuple(args, "OlOOOOOl", &py_triplets, &grid_point, - &py_bz_grid_address, &py_bz_map, &py_map_triplets, - &py_D_diag, &py_Q, &bz_grid_type)) { - return NULL; - } - - triplets = (long(*)[3])PyArray_DATA(py_triplets); - bz_grid_address = (long(*)[3])PyArray_DATA(py_bz_grid_address); - bz_map = (long *)PyArray_DATA(py_bz_map); - map_triplets = (long *)PyArray_DATA(py_map_triplets); - num_map_triplets = (long)PyArray_DIMS(py_map_triplets)[0]; - D_diag = (long *)PyArray_DATA(py_D_diag); - Q = (long(*)[3])PyArray_DATA(py_Q); - - num_ir = ph3py_get_BZ_triplets_at_q(triplets, grid_point, bz_grid_address, - bz_map, map_triplets, num_map_triplets, - D_diag, Q, bz_grid_type); - - return PyLong_FromLong(num_ir); -} - -static PyObject *py_get_triplets_integration_weights(PyObject *self, - PyObject *args) { - PyArrayObject *py_iw; - PyArrayObject *py_iw_zero; - PyArrayObject *py_frequency_points; - PyArrayObject *py_relative_grid_address; - PyArrayObject *py_D_diag; - PyArrayObject *py_triplets; - PyArrayObject *py_frequencies1; - PyArrayObject *py_frequencies2; - PyArrayObject *py_bz_grid_addresses; - PyArrayObject *py_bz_map; - long bz_grid_type; - long tp_type; - - double *iw; - char *iw_zero; - double *frequency_points; - long(*relative_grid_address)[4][3]; - long *D_diag; - long(*triplets)[3]; - long(*bz_grid_addresses)[3]; - long *bz_map; - double *frequencies1, *frequencies2; - long num_band0, num_band1, num_band2, num_triplets; - - if (!PyArg_ParseTuple(args, "OOOOOOOOOOll", &py_iw, &py_iw_zero, - &py_frequency_points, &py_relative_grid_address, - &py_D_diag, &py_triplets, &py_frequencies1, - &py_frequencies2, &py_bz_grid_addresses, &py_bz_map, - &bz_grid_type, &tp_type)) { - return NULL; - } - - iw = (double *)PyArray_DATA(py_iw); - iw_zero = (char *)PyArray_DATA(py_iw_zero); - frequency_points = (double *)PyArray_DATA(py_frequency_points); - num_band0 = (long)PyArray_DIMS(py_frequency_points)[0]; - relative_grid_address = - (long(*)[4][3])PyArray_DATA(py_relative_grid_address); - D_diag = (long *)PyArray_DATA(py_D_diag); - triplets = (long(*)[3])PyArray_DATA(py_triplets); - num_triplets = (long)PyArray_DIMS(py_triplets)[0]; - bz_grid_addresses = (long(*)[3])PyArray_DATA(py_bz_grid_addresses); - bz_map = (long *)PyArray_DATA(py_bz_map); - frequencies1 = (double *)PyArray_DATA(py_frequencies1); - frequencies2 = (double *)PyArray_DATA(py_frequencies2); - num_band1 = (long)PyArray_DIMS(py_frequencies1)[1]; - num_band2 = (long)PyArray_DIMS(py_frequencies2)[1]; - - ph3py_get_integration_weight( - iw, iw_zero, frequency_points, num_band0, relative_grid_address, D_diag, - triplets, num_triplets, bz_grid_addresses, bz_map, bz_grid_type, - frequencies1, num_band1, frequencies2, num_band2, tp_type, 1); - - Py_RETURN_NONE; -} - -static PyObject *py_get_triplets_integration_weights_with_sigma( - PyObject *self, PyObject *args) { - PyArrayObject *py_iw; - PyArrayObject *py_iw_zero; - PyArrayObject *py_frequency_points; - PyArrayObject *py_triplets; - PyArrayObject *py_frequencies; - double sigma, sigma_cutoff; - - double *iw; - char *iw_zero; - double *frequency_points; - long(*triplets)[3]; - double *frequencies; - long num_band0, num_band, num_iw, num_triplets; - - if (!PyArg_ParseTuple(args, "OOOOOdd", &py_iw, &py_iw_zero, - &py_frequency_points, &py_triplets, &py_frequencies, - &sigma, &sigma_cutoff)) { - return NULL; - } - - iw = (double *)PyArray_DATA(py_iw); - iw_zero = (char *)PyArray_DATA(py_iw_zero); - frequency_points = (double *)PyArray_DATA(py_frequency_points); - num_band0 = (long)PyArray_DIMS(py_frequency_points)[0]; - triplets = (long(*)[3])PyArray_DATA(py_triplets); - num_triplets = (long)PyArray_DIMS(py_triplets)[0]; - frequencies = (double *)PyArray_DATA(py_frequencies); - num_band = (long)PyArray_DIMS(py_frequencies)[1]; - num_iw = (long)PyArray_DIMS(py_iw)[0]; - - ph3py_get_integration_weight_with_sigma( - iw, iw_zero, sigma, sigma_cutoff, frequency_points, num_band0, triplets, - num_triplets, frequencies, num_band, num_iw); - - Py_RETURN_NONE; -} - -static PyObject *py_get_grid_index_from_address(PyObject *self, - PyObject *args) { - PyArrayObject *py_address; - PyArrayObject *py_D_diag; - - long *address; - long *D_diag; - long gp; - - if (!PyArg_ParseTuple(args, "OO", &py_address, &py_D_diag)) { - return NULL; - } - - address = (long *)PyArray_DATA(py_address); - D_diag = (long *)PyArray_DATA(py_D_diag); - - gp = ph3py_get_grid_index_from_address(address, D_diag); - - return PyLong_FromLong(gp); -} - -static PyObject *py_get_gr_grid_addresses(PyObject *self, PyObject *args) { - PyArrayObject *py_gr_grid_addresses; - PyArrayObject *py_D_diag; - - long(*gr_grid_addresses)[3]; - long *D_diag; - - if (!PyArg_ParseTuple(args, "OO", &py_gr_grid_addresses, &py_D_diag)) { - return NULL; - } - - gr_grid_addresses = (long(*)[3])PyArray_DATA(py_gr_grid_addresses); - D_diag = (long *)PyArray_DATA(py_D_diag); - - ph3py_get_gr_grid_addresses(gr_grid_addresses, D_diag); - - Py_RETURN_NONE; -} - -static PyObject *py_get_reciprocal_rotations(PyObject *self, PyObject *args) { - PyArrayObject *py_rec_rotations; - PyArrayObject *py_rotations; - long is_time_reversal; - - long(*rec_rotations)[3][3]; - long(*rotations)[3][3]; - long num_rot, num_rec_rot; - - if (!PyArg_ParseTuple(args, "OOl", &py_rec_rotations, &py_rotations, - &is_time_reversal)) { - return NULL; - } - - rec_rotations = (long(*)[3][3])PyArray_DATA(py_rec_rotations); - rotations = (long(*)[3][3])PyArray_DATA(py_rotations); - num_rot = (long)PyArray_DIMS(py_rotations)[0]; - - num_rec_rot = ph3py_get_reciprocal_rotations(rec_rotations, rotations, - num_rot, is_time_reversal); - - return PyLong_FromLong(num_rec_rot); -} - -static PyObject *py_transform_rotations(PyObject *self, PyObject *args) { - PyArrayObject *py_transformed_rotations; - PyArrayObject *py_rotations; - PyArrayObject *py_D_diag; - PyArrayObject *py_Q; - - long(*transformed_rotations)[3][3]; - long(*rotations)[3][3]; - long *D_diag; - long(*Q)[3]; - long num_rot, succeeded; - - if (!PyArg_ParseTuple(args, "OOOO", &py_transformed_rotations, - &py_rotations, &py_D_diag, &py_Q)) { - return NULL; - } - - transformed_rotations = - (long(*)[3][3])PyArray_DATA(py_transformed_rotations); - rotations = (long(*)[3][3])PyArray_DATA(py_rotations); - D_diag = (long *)PyArray_DATA(py_D_diag); - Q = (long(*)[3])PyArray_DATA(py_Q); - num_rot = (long)PyArray_DIMS(py_transformed_rotations)[0]; - - succeeded = ph3py_transform_rotations(transformed_rotations, rotations, - num_rot, D_diag, Q); - if (succeeded) { - Py_RETURN_TRUE; - } else { - Py_RETURN_FALSE; - } -} - -static PyObject *py_get_snf3x3(PyObject *self, PyObject *args) { - PyArrayObject *py_D_diag; - PyArrayObject *py_P; - PyArrayObject *py_Q; - PyArrayObject *py_A; - - long *D_diag; - long(*P)[3]; - long(*Q)[3]; - long(*A)[3]; - long succeeded; - - if (!PyArg_ParseTuple(args, "OOOO", &py_D_diag, &py_P, &py_Q, &py_A)) { - return NULL; - } - - D_diag = (long *)PyArray_DATA(py_D_diag); - P = (long(*)[3])PyArray_DATA(py_P); - Q = (long(*)[3])PyArray_DATA(py_Q); - A = (long(*)[3])PyArray_DATA(py_A); - - succeeded = ph3py_get_snf3x3(D_diag, P, Q, A); - if (succeeded) { - Py_RETURN_TRUE; - } else { - Py_RETURN_FALSE; - } -} - -static PyObject *py_get_ir_grid_map(PyObject *self, PyObject *args) { - PyArrayObject *py_grid_mapping_table; - PyArrayObject *py_D_diag; - PyArrayObject *py_is_shift; - PyArrayObject *py_rotations; - - long *D_diag; - long *is_shift; - long(*rot)[3][3]; - long num_rot; - - long *grid_mapping_table; - long num_ir; - - if (!PyArg_ParseTuple(args, "OOOO", &py_grid_mapping_table, &py_D_diag, - &py_is_shift, &py_rotations)) { - return NULL; - } - - D_diag = (long *)PyArray_DATA(py_D_diag); - is_shift = (long *)PyArray_DATA(py_is_shift); - rot = (long(*)[3][3])PyArray_DATA(py_rotations); - num_rot = (long)PyArray_DIMS(py_rotations)[0]; - grid_mapping_table = (long *)PyArray_DATA(py_grid_mapping_table); - - num_ir = ph3py_get_ir_grid_map(grid_mapping_table, D_diag, is_shift, rot, - num_rot); - return PyLong_FromLong(num_ir); -} - -static PyObject *py_get_bz_grid_addresses(PyObject *self, PyObject *args) { - PyArrayObject *py_bz_grid_addresses; - PyArrayObject *py_bz_map; - PyArrayObject *py_bzg2grg; - PyArrayObject *py_D_diag; - PyArrayObject *py_Q; - PyArrayObject *py_PS; - PyArrayObject *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; - - if (!PyArg_ParseTuple(args, "OOOOOOOl", &py_bz_grid_addresses, &py_bz_map, - &py_bzg2grg, &py_D_diag, &py_Q, &py_PS, - &py_reciprocal_lattice, &type)) { - return NULL; - } - - bz_grid_addresses = (long(*)[3])PyArray_DATA(py_bz_grid_addresses); - bz_map = (long *)PyArray_DATA(py_bz_map); - bzg2grg = (long *)PyArray_DATA(py_bzg2grg); - D_diag = (long *)PyArray_DATA(py_D_diag); - Q = (long(*)[3])PyArray_DATA(py_Q); - PS = (long *)PyArray_DATA(py_PS); - reciprocal_lattice = (double(*)[3])PyArray_DATA(py_reciprocal_lattice); - - num_total_gp = - ph3py_get_bz_grid_addresses(bz_grid_addresses, bz_map, bzg2grg, D_diag, - Q, PS, reciprocal_lattice, type); - - return PyLong_FromLong(num_total_gp); -} - -static PyObject *py_rotate_bz_grid_addresses(PyObject *self, PyObject *args) { - PyArrayObject *py_bz_grid_addresses; - PyArrayObject *py_rotation; - PyArrayObject *py_bz_map; - PyArrayObject *py_D_diag; - PyArrayObject *py_PS; - long bz_grid_index; - long type; - - long(*bz_grid_addresses)[3]; - long(*rotation)[3]; - long *bz_map; - long *D_diag; - long *PS; - long ret_bz_gp; - - if (!PyArg_ParseTuple(args, "lOOOOOl", &bz_grid_index, &py_rotation, - &py_bz_grid_addresses, &py_bz_map, &py_D_diag, &py_PS, - &type)) { - return NULL; - } - - bz_grid_addresses = (long(*)[3])PyArray_DATA(py_bz_grid_addresses); - rotation = (long(*)[3])PyArray_DATA(py_rotation); - bz_map = (long *)PyArray_DATA(py_bz_map); - D_diag = (long *)PyArray_DATA(py_D_diag); - PS = (long *)PyArray_DATA(py_PS); - - ret_bz_gp = ph3py_rotate_bz_grid_index( - bz_grid_index, rotation, bz_grid_addresses, bz_map, D_diag, PS, type); - - return PyLong_FromLong(ret_bz_gp); -} - -static PyObject *py_diagonalize_collision_matrix(PyObject *self, - PyObject *args) { - PyArrayObject *py_collision_matrix; - PyArrayObject *py_eigenvalues; - double cutoff; - long i_sigma, i_temp, is_pinv, solver; - - double *collision_matrix; - double *eigvals; - long num_temp, num_grid_point, num_band; - long num_column, adrs_shift; - long info; - - if (!PyArg_ParseTuple(args, "OOlldll", &py_collision_matrix, - &py_eigenvalues, &i_sigma, &i_temp, &cutoff, &solver, - &is_pinv)) { - return NULL; - } - - collision_matrix = (double *)PyArray_DATA(py_collision_matrix); - eigvals = (double *)PyArray_DATA(py_eigenvalues); - - if (PyArray_NDIM(py_collision_matrix) == 2) { - num_temp = 1; - num_column = (long)PyArray_DIM(py_collision_matrix, 1); - } else { - num_temp = (long)PyArray_DIM(py_collision_matrix, 1); - num_grid_point = (long)PyArray_DIM(py_collision_matrix, 2); - num_band = (long)PyArray_DIM(py_collision_matrix, 3); - if (PyArray_NDIM(py_collision_matrix) == 8) { - num_column = num_grid_point * num_band * 3; - } else { - num_column = num_grid_point * num_band; - } - } - adrs_shift = (i_sigma * num_column * num_column * num_temp + - i_temp * num_column * num_column); - - /* show_colmat_info(py_collision_matrix, i_sigma, i_temp, adrs_shift); */ - - info = ph3py_phonopy_dsyev(collision_matrix + adrs_shift, eigvals, - num_column, solver); - if (is_pinv) { - ph3py_pinv_from_eigensolution(collision_matrix + adrs_shift, eigvals, - num_column, cutoff, 0); - } - - return PyLong_FromLong(info); -} - -static PyObject *py_pinv_from_eigensolution(PyObject *self, PyObject *args) { - PyArrayObject *py_collision_matrix; - PyArrayObject *py_eigenvalues; - double cutoff; - long i_sigma, i_temp, pinv_method; - - double *collision_matrix; - double *eigvals; - long num_temp, num_grid_point, num_band; - long num_column, adrs_shift; - - if (!PyArg_ParseTuple(args, "OOlldl", &py_collision_matrix, &py_eigenvalues, - &i_sigma, &i_temp, &cutoff, &pinv_method)) { - return NULL; - } - - collision_matrix = (double *)PyArray_DATA(py_collision_matrix); - eigvals = (double *)PyArray_DATA(py_eigenvalues); - num_temp = (long)PyArray_DIMS(py_collision_matrix)[1]; - num_grid_point = (long)PyArray_DIMS(py_collision_matrix)[2]; - num_band = (long)PyArray_DIMS(py_collision_matrix)[3]; - - if (PyArray_NDIM(py_collision_matrix) == 8) { - num_column = num_grid_point * num_band * 3; - } else { - num_column = num_grid_point * num_band; - } - adrs_shift = (i_sigma * num_column * num_column * num_temp + - i_temp * num_column * num_column); - - /* show_colmat_info(py_collision_matrix, i_sigma, i_temp, adrs_shift); */ - - ph3py_pinv_from_eigensolution(collision_matrix + adrs_shift, eigvals, - num_column, cutoff, pinv_method); - - Py_RETURN_NONE; -} - -static PyObject *py_get_default_colmat_solver(PyObject *self, PyObject *args) { - if (!PyArg_ParseTuple(args, "")) { - return NULL; - } - -#if defined(MKL_LAPACKE) || defined(SCIPY_MKL_H) - return PyLong_FromLong((long)1); -#else - return PyLong_FromLong((long)4); -#endif -} - -static PyObject *py_lapacke_pinv(PyObject *self, PyObject *args) { - PyArrayObject *data_in_py; - PyArrayObject *data_out_py; - double cutoff; - - int m; - int n; - double *data_in; - double *data_out; - int info; - - if (!PyArg_ParseTuple(args, "OOd", &data_out_py, &data_in_py, &cutoff)) { - return NULL; - } - - m = (long)PyArray_DIMS(data_in_py)[0]; - n = (long)PyArray_DIMS(data_in_py)[1]; - data_in = (double *)PyArray_DATA(data_in_py); - data_out = (double *)PyArray_DATA(data_out_py); - - info = ph3py_phonopy_pinv(data_out, data_in, m, n, cutoff); - - return PyLong_FromLong((long)info); -} - -static PyObject *py_get_omp_max_threads(PyObject *self, PyObject *args) { - return PyLong_FromLong(ph3py_get_max_threads()); -} - -static void show_colmat_info(const PyArrayObject *py_collision_matrix, - const long i_sigma, const long i_temp, - const long adrs_shift) { - long i; - - printf(" Array_shape:("); - for (i = 0; i < PyArray_NDIM(py_collision_matrix); i++) { - printf("%d", (int)PyArray_DIM(py_collision_matrix, i)); - if (i < PyArray_NDIM(py_collision_matrix) - 1) { - printf(","); - } else { - printf("), "); - } - } - printf("Data shift:%lu [%lu, %lu]\n", adrs_shift, i_sigma, i_temp); -} - -/** - * @brief Convert numpy "int_" array to phono3py long array structure. - * - * @param npyary - * @return Larray* - */ -static Larray *convert_to_larray(PyArrayObject *npyary) { - long i; - Larray *ary; - - ary = (Larray *)malloc(sizeof(Larray)); - for (i = 0; i < PyArray_NDIM(npyary); i++) { - ary->dims[i] = PyArray_DIMS(npyary)[i]; - } - ary->data = (long *)PyArray_DATA(npyary); - return ary; -} - -/** - * @brief Convert numpy "double" array to phono3py double array structure. - * - * @param npyary - * @return Darray* - * @note PyArray_NDIM receives non-const (PyArrayObject *). - */ -static Darray *convert_to_darray(PyArrayObject *npyary) { - int i; - Darray *ary; - - ary = (Darray *)malloc(sizeof(Darray)); - for (i = 0; i < PyArray_NDIM(npyary); i++) { - ary->dims[i] = PyArray_DIMS(npyary)[i]; - } - ary->data = (double *)PyArray_DATA(npyary); - return ary; -} diff --git a/c/_phono3py.cpp b/c/_phono3py.cpp index 061631a5..a7017d66 100644 --- a/c/_phono3py.cpp +++ b/c/_phono3py.cpp @@ -1049,6 +1049,20 @@ long py_rotate_bz_grid_addresses(long bz_grid_index, nb::ndarray<> py_rotation, return ret_bz_gp; } +long py_get_default_colmat_solver() { +#if defined(MKL_BLAS) || defined(SCIPY_MKL_H) || !defined(NO_INCLUDE_LAPACKE) + return (long)1; +#else + return (long)4; +#endif +} + +long py_get_omp_max_threads() { return ph3py_get_max_threads(); } + +#ifdef NO_INCLUDE_LAPACKE +long py_include_lapacke() { return 0; } +#else +long py_include_lapacke() { return 1; } long py_diagonalize_collision_matrix(nb::ndarray<> py_collision_matrix, nb::ndarray<> py_eigenvalues, long i_sigma, long i_temp, double cutoff, long solver, @@ -1118,14 +1132,6 @@ void py_pinv_from_eigensolution(nb::ndarray<> py_collision_matrix, num_column, cutoff, pinv_method); } -long py_get_default_colmat_solver() { -#if defined(MKL_LAPACKE) || defined(SCIPY_MKL_H) - return (long)1; -#else - return (long)4; -#endif -} - long py_lapacke_pinv(nb::ndarray<> data_out_py, nb::ndarray<> data_in_py, double cutoff) { long m; @@ -1143,8 +1149,7 @@ long py_lapacke_pinv(nb::ndarray<> data_out_py, nb::ndarray<> data_in_py, return info; } - -long py_get_omp_max_threads() { return ph3py_get_max_threads(); } +#endif NB_MODULE(_phono3py, m) { m.def("interaction", &py_get_interaction); @@ -1179,9 +1184,12 @@ 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("default_colmat_solver", &py_get_default_colmat_solver); + m.def("omp_max_threads", &py_get_omp_max_threads); + m.def("include_lapacke", &py_include_lapacke); +#ifndef NO_INCLUDE_LAPACKE 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); m.def("lapacke_pinv", &py_lapacke_pinv); - m.def("omp_max_threads", &py_get_omp_max_threads); +#endif } diff --git a/c/_phononcalc.c b/c/_phononcalc.c deleted file mode 100644 index 941b70da..00000000 --- a/c/_phononcalc.c +++ /dev/null @@ -1,232 +0,0 @@ -/* 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 -#include - -// #include "lapack_wrapper.h" -#include "phononcalc.h" - -static PyObject *py_get_phonons_at_gridpoints(PyObject *self, PyObject *args); - -struct module_state { - PyObject *error; -}; - -#if PY_MAJOR_VERSION >= 3 -#define GETSTATE(m) ((struct module_state *)PyModule_GetState(m)) -#else -#define GETSTATE(m) (&_state) -static struct module_state _state; -#endif - -static PyObject *error_out(PyObject *m) { - struct module_state *st = GETSTATE(m); - PyErr_SetString(st->error, "something bad happened"); - return NULL; -} - -static PyMethodDef _phononcalc_methods[] = { - {"error_out", (PyCFunction)error_out, METH_NOARGS, NULL}, - {"phonons_at_gridpoints", py_get_phonons_at_gridpoints, METH_VARARGS, - "Set phonons at grid points"}, - {NULL, NULL, 0, NULL}}; - -#if PY_MAJOR_VERSION >= 3 - -static int _phononcalc_traverse(PyObject *m, visitproc visit, void *arg) { - Py_VISIT(GETSTATE(m)->error); - return 0; -} - -static int _phononcalc_clear(PyObject *m) { - Py_CLEAR(GETSTATE(m)->error); - return 0; -} - -static struct PyModuleDef moduledef = { - PyModuleDef_HEAD_INIT, "_phononcalc", NULL, - sizeof(struct module_state), _phononcalc_methods, NULL, - _phononcalc_traverse, _phononcalc_clear, NULL}; - -#define INITERROR return NULL - -PyObject *PyInit__phononcalc(void) -#else -#define INITERROR return - -void init_phononcalc(void) -#endif -{ -#if PY_MAJOR_VERSION >= 3 - PyObject *module = PyModule_Create(&moduledef); -#else - PyObject *module = Py_InitModule("_phononcalc", _phononcalc_methods); -#endif - struct module_state *st; - - if (module == NULL) INITERROR; - st = GETSTATE(module); - - st->error = PyErr_NewException("_phononcalc.Error", NULL, NULL); - if (st->error == NULL) { - Py_DECREF(module); - INITERROR; - } - -#if PY_MAJOR_VERSION >= 3 - return module; -#endif -} - -static PyObject *py_get_phonons_at_gridpoints(PyObject *self, PyObject *args) { - PyArrayObject *py_frequencies; - PyArrayObject *py_eigenvectors; - PyArrayObject *py_phonon_done; - PyArrayObject *py_grid_points; - PyArrayObject *py_grid_address; - PyArrayObject *py_QDinv; - PyArrayObject *py_shortest_vectors_fc2; - PyArrayObject *py_multiplicity_fc2; - PyArrayObject *py_positions_fc2; - PyArrayObject *py_fc2; - PyArrayObject *py_masses_fc2; - PyArrayObject *py_p2s_map_fc2; - PyArrayObject *py_s2p_map_fc2; - PyArrayObject *py_reciprocal_lattice; - PyArrayObject *py_born_effective_charge; - PyArrayObject *py_q_direction; - PyArrayObject *py_dielectric_constant; - PyArrayObject *py_dd_q0; - PyArrayObject *py_G_list; - double nac_factor; - double unit_conversion_factor; - double lambda; - char *uplo; - - double(*born)[3][3]; - double(*dielectric)[3]; - double *q_dir; - double *freqs; - _lapack_complex_double *eigvecs; - char *phonon_done; - long *grid_points; - long(*grid_address)[3]; - double(*QDinv)[3]; - double *fc2; - double(*svecs_fc2)[3]; - long(*multi_fc2)[2]; - double(*positions_fc2)[3]; - double *masses_fc2; - long *p2s_fc2; - long *s2p_fc2; - double(*rec_lat)[3]; - double(*dd_q0)[2]; - double(*G_list)[3]; - long num_patom, num_satom, num_phonons, num_grid_points, num_G_points; - - if (!PyArg_ParseTuple( - args, "OOOOOOOOOOOOOdOOOOdOOds", &py_frequencies, &py_eigenvectors, - &py_phonon_done, &py_grid_points, &py_grid_address, &py_QDinv, - &py_fc2, &py_shortest_vectors_fc2, &py_multiplicity_fc2, - &py_positions_fc2, &py_masses_fc2, &py_p2s_map_fc2, &py_s2p_map_fc2, - &unit_conversion_factor, &py_born_effective_charge, - &py_dielectric_constant, &py_reciprocal_lattice, &py_q_direction, - &nac_factor, &py_dd_q0, &py_G_list, &lambda, &uplo)) { - return NULL; - } - - freqs = (double *)PyArray_DATA(py_frequencies); - eigvecs = (_lapack_complex_double *)PyArray_DATA(py_eigenvectors); - phonon_done = (char *)PyArray_DATA(py_phonon_done); - grid_points = (long *)PyArray_DATA(py_grid_points); - grid_address = (long(*)[3])PyArray_DATA(py_grid_address); - QDinv = (double(*)[3])PyArray_DATA(py_QDinv); - fc2 = (double *)PyArray_DATA(py_fc2); - svecs_fc2 = (double(*)[3])PyArray_DATA(py_shortest_vectors_fc2); - multi_fc2 = (long(*)[2])PyArray_DATA(py_multiplicity_fc2); - masses_fc2 = (double *)PyArray_DATA(py_masses_fc2); - p2s_fc2 = (long *)PyArray_DATA(py_p2s_map_fc2); - s2p_fc2 = (long *)PyArray_DATA(py_s2p_map_fc2); - rec_lat = (double(*)[3])PyArray_DATA(py_reciprocal_lattice); - num_patom = (long)PyArray_DIMS(py_multiplicity_fc2)[1]; - num_satom = (long)PyArray_DIMS(py_multiplicity_fc2)[0]; - num_phonons = (long)PyArray_DIMS(py_frequencies)[0]; - num_grid_points = (long)PyArray_DIMS(py_grid_points)[0]; - if ((PyObject *)py_born_effective_charge == Py_None) { - born = NULL; - } else { - born = (double(*)[3][3])PyArray_DATA(py_born_effective_charge); - } - if ((PyObject *)py_dielectric_constant == Py_None) { - dielectric = NULL; - } else { - dielectric = (double(*)[3])PyArray_DATA(py_dielectric_constant); - } - if ((PyObject *)py_q_direction == Py_None) { - q_dir = NULL; - } else { - q_dir = (double *)PyArray_DATA(py_q_direction); - if (fabs(q_dir[0]) < 1e-10 && fabs(q_dir[1]) < 1e-10 && - fabs(q_dir[2]) < 1e-10) { - q_dir = NULL; - } - } - if ((PyObject *)py_dd_q0 == Py_None) { - dd_q0 = NULL; - } else { - dd_q0 = (double(*)[2])PyArray_DATA(py_dd_q0); - } - if ((PyObject *)py_G_list == Py_None) { - G_list = NULL; - num_G_points = 0; - } else { - G_list = (double(*)[3])PyArray_DATA(py_G_list); - num_G_points = (long)PyArray_DIMS(py_G_list)[0]; - } - if ((PyObject *)py_positions_fc2 == Py_None) { - positions_fc2 = NULL; - } else { - positions_fc2 = (double(*)[3])PyArray_DATA(py_positions_fc2); - } - - phcalc_get_phonons_at_gridpoints( - freqs, eigvecs, phonon_done, num_phonons, grid_points, num_grid_points, - grid_address, QDinv, fc2, svecs_fc2, multi_fc2, positions_fc2, - num_patom, num_satom, masses_fc2, p2s_fc2, s2p_fc2, - unit_conversion_factor, born, dielectric, rec_lat, q_dir, nac_factor, - dd_q0, G_list, num_G_points, lambda, uplo[0]); - - Py_RETURN_NONE; -} diff --git a/c/lapack_wrapper.c b/c/lapack_wrapper.c index 0b6704b9..166e2000 100644 --- a/c/lapack_wrapper.c +++ b/c/lapack_wrapper.c @@ -38,13 +38,16 @@ #define min(a, b) ((a) > (b) ? (b) : (a)) -#if defined(MKL_LAPACKE) || defined(SCIPY_MKL_H) -MKL_Complex16 lapack_make_complex_double(double re, double im) { - MKL_Complex16 z; +#if (defined(MKL_BLAS) || defined(SCIPY_MKL_H)) && !defined(NO_INCLUDE_LAPACKE) +lapack_complex_double lapack_make_complex_double(double re, double im) { + lapack_complex_double z; z.real = re; z.imag = im; return z; } +#endif + +#if defined(MKL_BLAS) || defined(SCIPY_MKL_H) || defined(NO_INCLUDE_LAPACKE) #ifndef LAPACKE_malloc #define LAPACKE_malloc(size) malloc(size) #endif @@ -53,6 +56,18 @@ MKL_Complex16 lapack_make_complex_double(double re, double im) { #endif #endif +lapack_complex_double phonoc_complex_prod(const lapack_complex_double a, + const lapack_complex_double b) { + lapack_complex_double c; + c = lapack_make_complex_double( + lapack_complex_double_real(a) * lapack_complex_double_real(b) - + lapack_complex_double_imag(a) * lapack_complex_double_imag(b), + lapack_complex_double_imag(a) * lapack_complex_double_real(b) + + lapack_complex_double_real(a) * lapack_complex_double_imag(b)); + return c; +} + +#ifndef NO_INCLUDE_LAPACKE int phonopy_zheev(double *w, lapack_complex_double *a, const int n, const char uplo) { lapack_int info; @@ -186,17 +201,6 @@ int phonopy_dsyev(double *data, double *eigvals, const int size, return (int)info; } -lapack_complex_double phonoc_complex_prod(const lapack_complex_double a, - const lapack_complex_double b) { - lapack_complex_double c; - c = lapack_make_complex_double( - lapack_complex_double_real(a) * lapack_complex_double_real(b) - - lapack_complex_double_imag(a) * lapack_complex_double_imag(b), - lapack_complex_double_imag(a) * lapack_complex_double_real(b) + - lapack_complex_double_real(a) * lapack_complex_double_imag(b)); - return c; -} - void pinv_from_eigensolution(double *data, const double *eigvals, const long size, const double cutoff, const long pinv_method) { @@ -284,3 +288,4 @@ void pinv_from_eigensolution(double *data, const double *eigvals, free(tmp_data); tmp_data = NULL; } +#endif diff --git a/c/lapack_wrapper.h b/c/lapack_wrapper.h index 13d398c7..b18db54a 100644 --- a/c/lapack_wrapper.h +++ b/c/lapack_wrapper.h @@ -35,16 +35,34 @@ #ifndef __lapack_wrapper_H__ #define __lapack_wrapper_H__ -#if defined(MKL_LAPACKE) || defined(SCIPY_MKL_H) +#ifdef NO_INCLUDE_LAPACKE +#include +#define lapack_complex_double double _Complex +#ifdef CMPLX +#define lapack_make_complex_double(re, im) CMPLX(re, im) +#else +#define lapack_make_complex_double(re, im) ((double _Complex)((re) + (im) * I)) +#endif +#define lapack_complex_double_real(z) (creal(z)) +#define lapack_complex_double_imag(z) (cimag(z)) +#endif + +#if defined(MKL_BLAS) || defined(SCIPY_MKL_H) #include #define lapack_complex_double MKL_Complex16 +MKL_Complex16 lapack_make_complex_double(double re, double im); #define lapack_complex_double_real(z) ((z).real) #define lapack_complex_double_imag(z) ((z).imag) -MKL_Complex16 lapack_make_complex_double(double re, double im); -#else +#endif + +#if !defined(MKL_BLAS) && !defined(SCIPY_MKL_H) && !defined(NO_INCLUDE_LAPACKE) #include #endif +lapack_complex_double phonoc_complex_prod(const lapack_complex_double a, + const lapack_complex_double b); + +#ifndef NO_INCLUDE_LAPACKE int phonopy_zheev(double *w, lapack_complex_double *a, const int n, const char uplo); int phonopy_pinv(double *data_out, const double *data_in, const int m, @@ -56,10 +74,9 @@ void phonopy_pinv_mt(double *data_out, int *info_out, const double *data_in, int phonopy_dsyev(double *data, double *eigvals, const int size, const int algorithm); -lapack_complex_double phonoc_complex_prod(const lapack_complex_double a, - const lapack_complex_double b); - void pinv_from_eigensolution(double *data, const double *eigvals, const long size, const double cutoff, const long pinv_method); #endif + +#endif diff --git a/c/phono3py.c b/c/phono3py.c index 4f7f625b..39e7050e 100644 --- a/c/phono3py.c +++ b/c/phono3py.c @@ -697,6 +697,15 @@ long ph3py_get_thm_integration_weights_at_grid_points( return 1; } +long ph3py_get_max_threads(void) { +#ifdef _OPENMP + return omp_get_max_threads(); +#else + return 0; +#endif +} + +#ifndef NO_INCLUDE_LAPACKE long ph3py_phonopy_dsyev(double *data, double *eigvals, const long size, const long algorithm) { return (long)phonopy_dsyev(data, eigvals, (int)size, (int)algorithm); @@ -712,11 +721,4 @@ void ph3py_pinv_from_eigensolution(double *data, const double *eigvals, const long pinv_method) { pinv_from_eigensolution(data, eigvals, size, cutoff, pinv_method); } - -long ph3py_get_max_threads(void) { -#ifdef _OPENMP - return omp_get_max_threads(); -#else - return 0; #endif -} diff --git a/c/phono3py.h b/c/phono3py.h index d6a2972f..54247e24 100644 --- a/c/phono3py.h +++ b/c/phono3py.h @@ -233,6 +233,9 @@ long ph3py_get_thm_integration_weights_at_grid_points( const long *grid_points, const long (*bz_grid_addresses)[3], const long *bz_map, const long bz_grid_type, const double *frequencies, const long *gp2irgp_map, const char function); +long ph3py_get_max_threads(void); + +#ifndef NO_INCLUDE_LAPACKE long ph3py_phonopy_dsyev(double *data, double *eigvals, const long size, const long algorithm); long ph3py_phonopy_pinv(double *data_out, const double *data_in, const long m, @@ -240,7 +243,7 @@ long ph3py_phonopy_pinv(double *data_out, const double *data_in, const long m, void ph3py_pinv_from_eigensolution(double *data, const double *eigvals, const long size, const double cutoff, const long pinv_method); -long ph3py_get_max_threads(void); +#endif #ifdef __cplusplus } diff --git a/c/phonon.c b/c/phonon.c index 9282093b..f4e8a36d 100644 --- a/c/phonon.c +++ b/c/phonon.c @@ -36,6 +36,7 @@ #include #include +#include #include #include "dynmat.h" @@ -66,16 +67,6 @@ static void get_gonze_undone_phonons( 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 char uplo); -static void get_phonons(lapack_complex_double *eigvecs, const double q[3], - const double *fc2, const double *masses, - const long *p2s, const long *s2p, - const long (*multi)[2], const long num_patom, - const long num_satom, const double (*svecs)[3], - const long is_nac, const double (*born)[3][3], - const double dielectric[3][3], - const double reciprocal_lattice[3][3], - const double *q_direction, const double nac_factor, - const double unit_conversion_factor); static void get_gonze_phonons( lapack_complex_double *eigvecs, const double q[3], const double *fc2, const double *masses, const long *p2s, const long *s2p, @@ -208,12 +199,13 @@ static void get_undone_phonons( } is_nac = needs_nac(born, grid_address, gp, q_direction); - get_phonons(eigenvectors + num_band * num_band * gp, q, fc2, masses_fc2, - p2s_fc2, s2p_fc2, multi_fc2, num_patom, num_satom, - svecs_fc2, is_nac, born, dielectric, reciprocal_lattice, - q_direction, nac_factor, unit_conversion_factor); + get_dynamical_matrix(eigenvectors + num_band * num_band * gp, q, fc2, + masses_fc2, p2s_fc2, s2p_fc2, multi_fc2, num_patom, + num_satom, svecs_fc2, is_nac, born, dielectric, + reciprocal_lattice, q_direction, nac_factor); } +#ifndef NO_INCLUDE_LAPACKE /* To avoid multithreaded BLAS in OpenMP loop */ #ifndef MULTITHREADED_BLAS #ifdef _OPENMP @@ -235,6 +227,7 @@ static void get_undone_phonons( unit_conversion_factor; } } +#endif } static void get_gonze_undone_phonons( @@ -274,6 +267,7 @@ static void get_gonze_undone_phonons( nac_factor, dd_q0, G_list, num_G_points, lambda); } +#ifndef NO_INCLUDE_LAPACKE /* To avoid multithreaded BLAS in OpenMP loop */ #ifndef MULTITHREADED_BLAS #ifdef _OPENMP @@ -295,22 +289,7 @@ static void get_gonze_undone_phonons( unit_conversion_factor; } } -} - -static void get_phonons(lapack_complex_double *eigvecs, const double q[3], - const double *fc2, const double *masses, - const long *p2s, const long *s2p, - const long (*multi)[2], const long num_patom, - const long num_satom, const double (*svecs)[3], - const long is_nac, const double (*born)[3][3], - const double dielectric[3][3], - const double reciprocal_lattice[3][3], - const double *q_direction, const double nac_factor, - const double unit_conversion_factor) { - /* Store dynamical matrix in eigvecs array. */ - get_dynamical_matrix(eigvecs, q, fc2, masses, p2s, s2p, multi, num_patom, - num_satom, svecs, is_nac, born, dielectric, - reciprocal_lattice, q_direction, nac_factor); +#endif } static void get_gonze_phonons( diff --git a/c/reciprocal_to_normal.c b/c/reciprocal_to_normal.c index d061999c..02d6ffbd 100644 --- a/c/reciprocal_to_normal.c +++ b/c/reciprocal_to_normal.c @@ -35,7 +35,7 @@ #include "reciprocal_to_normal.h" #ifdef MULTITHREADED_BLAS -#if defined(MKL_LAPACKE) || defined(SCIPY_MKL_H) +#if defined(MKL_BLAS) || defined(SCIPY_MKL_H) #include #else #include diff --git a/doc/changelog.md b/doc/changelog.md index fea78c71..dd7516ef 100644 --- a/doc/changelog.md +++ b/doc/changelog.md @@ -2,6 +2,12 @@ # Change Log +## Nov-25-2024: Version 3.8.0 + +- Follow the change due to phonopy's refactoring of MLP interface. +- Experimental option (`BUILD_WITHOUT_LAPACKE=ON`) to compile phono3py without + LAPACKE in C + ## Nov-13-2024: Version 3.7.0 - Update to follow the change of phonopy's internal functions diff --git a/doc/conf.py b/doc/conf.py index e42ff2fa..8a8a5292 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -58,9 +58,9 @@ # built documents. # # The short X.Y version. -version = "3.7" +version = "3.8" # The full version, including alpha/beta/rc tags. -release = "3.7.0" +release = "3.8.0" # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. diff --git a/doc/direct-solution.md b/doc/direct-solution.md index a13c3fe1..2d716feb 100644 --- a/doc/direct-solution.md +++ b/doc/direct-solution.md @@ -307,7 +307,7 @@ is used through them. Multithreaded OpenBLAS is installed by conda and can be used via LAPACKE in phono3py. MKL LAPACK and BLAS are also able to be used via LAPACKE in phono3py with appropriate setting in `setup.py`. -Using `--pinv-solver=[number]`, one of the following solver is chosen: +Using `--pinv-solver NUMBER`, one of the following solver is chosen: 1. Lapacke `dsyev`: Smaller memory consumption than `dsyevd`, but slower. This is the default solver when MKL LAPACKE is integrated or scipy is not @@ -323,3 +323,8 @@ Using `--pinv-solver=[number]`, one of the following solver is chosen: LAPACKE is not integrated. 5. Scipy's `dsyevd`: Similar to solver (2), this solver should be used carefully. + +```{note} +`pinv-solver` =3, 4, and 5 are only available when phono3py is compiled +without LAPACKE. +``` diff --git a/doc/install.md b/doc/install.md index 04a10c1f..13c965f3 100644 --- a/doc/install.md +++ b/doc/install.md @@ -49,7 +49,7 @@ conda environment. Automatic search of required libraries and flags that are already on the system is performed by cmake. (install_with_cmake)= -### Build with automatic search of library configurations by cmake +### Building with automatic search of library configurations by cmake With installed cmake and required libraries on the system, cmake tries to find libraries to be linked and the compiler flags that are required during the build. @@ -61,6 +61,21 @@ See an example at {ref}`install_an_example`. In the standard output, flags and libraries found by cmake are shown. Please carefully check if those configurations are expected ones or not. +### Building without linking LAPACKE + +**Experimental** + +To compile phono3py without linking LAPACKE in C, use the following command: + +``` +% BUILD_WITHOUT_LAPACKE=ON pip install -e . -vvv +``` + +When this option is enabled, linking to BLAS and LAPACKE libraries is not +required, so installing these libraries for C compilation may not be necessary. +However, since numpy and scipy rely on BLAS and LAPACK libraries, their runtime +versions are still required. + (install_an_example)= ## Installation instruction of latest development version of phono3py diff --git a/phono3py/api_phono3py.py b/phono3py/api_phono3py.py index c8cf0600..8f112e10 100644 --- a/phono3py/api_phono3py.py +++ b/phono3py/api_phono3py.py @@ -54,14 +54,9 @@ symmetrize_force_constants, ) from phonopy.interface.fc_calculator import get_fc2 +from phonopy.interface.mlp import PhonopyMLP from phonopy.interface.pypolymlp import ( - PypolymlpData, PypolymlpParams, - develop_mlp_by_pypolymlp, - develop_polymlp, - evalulate_polymlp, - load_polymlp, - parse_mlp_params, ) from phonopy.structure.atoms import PhonopyAtoms from phonopy.structure.cells import ( @@ -721,6 +716,15 @@ def phonon_mlp_dataset(self, mlp_dataset: dict): self._check_mlp_dataset(mlp_dataset) self._phonon_mlp_dataset = mlp_dataset + @property + def mlp(self) -> Optional[PhonopyMLP]: + """Setter and getter of PhonopyMLP dataclass.""" + return self._mlp + + @mlp.setter + def mlp(self, mlp) -> Optional[PhonopyMLP]: + self._mlp = mlp + @property def band_indices(self) -> list[np.ndarray]: """Setter and getter of band indices. @@ -2176,7 +2180,7 @@ def develop_mlp( params: Optional[Union[PypolymlpParams, dict, str]] = None, test_size: float = 0.1, ): - """Develop MLP of pypolymlp. + """Develop machine learning potential. Parameters ---------- @@ -2192,20 +2196,28 @@ def develop_mlp( if self._mlp_dataset is None: raise RuntimeError("MLP dataset is not set.") - self._mlp = develop_mlp_by_pypolymlp( + self._mlp = PhonopyMLP(log_level=self._log_level) + self._mlp.develop( self._mlp_dataset, self._supercell, params=params, test_size=test_size, - log_level=self._log_level, ) - def load_mlp(self, filename: str = "phono3py.pmlp"): - """Load machine learning potential of pypolymlp.""" - self._mlp = load_polymlp(filename=filename) + def save_mlp(self, filename: Optional[str] = None): + """Save machine learning potential.""" + if self._mlp is None: + raise RuntimeError("MLP is not developed yet.") + + self._mlp.save(filename=filename) + + def load_mlp(self, filename: Optional[str] = None): + """Load machine learning potential.""" + self._mlp = PhonopyMLP(log_level=self._log_level) + self._mlp.load(filename=filename) def evaluate_mlp(self): - """Evaluate machine learning potential of pypolymlp. + """Evaluate machine learning potential. This method calculates the supercell energies and forces from the MLP for the displacements in self._dataset of type 2. The results are stored @@ -2223,9 +2235,7 @@ def evaluate_mlp(self): if self.supercells_with_displacements is None: raise RuntimeError("Displacements are not set. Run generate_displacements.") - energies, forces, _ = evalulate_polymlp( - self._mlp, self.supercells_with_displacements - ) + energies, forces, _ = self._mlp.evaluate(self.supercells_with_displacements) self.supercell_energies = energies self.forces = forces @@ -2234,7 +2244,7 @@ def develop_phonon_mlp( params: Optional[Union[PypolymlpParams, dict, str]] = None, test_size: float = 0.1, ): - """Develop MLP of pypolymlp for fc2. + """Develop MLP for fc2. Parameters ---------- @@ -2250,31 +2260,28 @@ def develop_phonon_mlp( if self._phonon_mlp_dataset is None: raise RuntimeError("MLP dataset is not set.") - if params is not None: - _params = parse_mlp_params(params) - else: - _params = params - - disps = self._phonon_mlp_dataset["displacements"] - forces = self._phonon_mlp_dataset["forces"] - energies = self._phonon_mlp_dataset["supercell_energies"] - n = int(len(disps) * (1 - test_size)) - train_data = PypolymlpData( - displacements=disps[:n], forces=forces[:n], supercell_energies=energies[:n] - ) - test_data = PypolymlpData( - displacements=disps[n:], forces=forces[n:], supercell_energies=energies[n:] - ) - self._phonon_mlp = develop_polymlp( + self._phonon_mlp = PhonopyMLP(log_level=self._log_level) + self._phonon_mlp.develop( + self._phonon_mlp_dataset, self._phonon_supercell, - train_data, - test_data, - params=_params, - verbose=self._log_level - 1 > 0, + params=params, + test_size=test_size, ) + def save_phonon_mlp(self, filename: Optional[str] = None): + """Save machine learning potential.""" + if self._mlp is None: + raise RuntimeError("MLP is not developed yet.") + + self._phonon_mlp.save(filename=filename) + + def load_phonon_mlp(self, filename: Optional[str] = None): + """Load machine learning potential.""" + self._phonon_mlp = PhonopyMLP(log_level=self._log_level) + self._phonon_mlp.load(filename=filename) + def evaluate_phonon_mlp(self): - """Evaluate the machine learning potential of pypolymlp. + """Evaluate the machine learning potential. This method calculates the supercell energies and forces from the MLP for the displacements in self._dataset of type 2. The results are stored @@ -2298,9 +2305,7 @@ def evaluate_phonon_mlp(self): mlp = self._mlp else: mlp = self._phonon_mlp - energies, forces, _ = evalulate_polymlp( - mlp, self.phonon_supercells_with_displacements - ) + energies, forces, _ = mlp.evaluate(self.phonon_supercells_with_displacements) self.phonon_supercell_energies = energies self.phonon_forces = forces diff --git a/phono3py/conductivity/utils.py b/phono3py/conductivity/utils.py index ea5797bc..ebdf67c6 100644 --- a/phono3py/conductivity/utils.py +++ b/phono3py/conductivity/utils.py @@ -686,6 +686,13 @@ def select_colmat_solver(pinv_solver): print("Phono3py C-routine is not compiled correctly.") default_solver = 4 + if not phono3c.include_lapacke(): + if pinv_solver in (1, 2, 6): + raise RuntimeError( + "Use pinv-solver 3, 4, or 5 because " + "phono3py is not compiled with LAPACKE." + ) + solver_numbers = (1, 2, 3, 4, 5, 6) solver = pinv_solver diff --git a/phono3py/cui/create_force_constants.py b/phono3py/cui/create_force_constants.py index 94c82b46..c1a97efa 100644 --- a/phono3py/cui/create_force_constants.py +++ b/phono3py/cui/create_force_constants.py @@ -533,7 +533,7 @@ def run_pypolymlp_to_compute_forces( if log_level: print("Developing MLPs by pypolymlp...", flush=True) ph3py.develop_mlp(params=mlp_params) - ph3py.mlp.save_mlp(filename=mlp_filename) + ph3py.save_mlp(filename=mlp_filename) if log_level: print(f'MLPs were written into "{mlp_filename}"', flush=True) else: diff --git a/phono3py/cui/kaccum_script.py b/phono3py/cui/kaccum_script.py index 149bf244..4eee31e5 100644 --- a/phono3py/cui/kaccum_script.py +++ b/phono3py/cui/kaccum_script.py @@ -85,7 +85,6 @@ def _get_ir_grid_info(bz_grid: BZGrid, weights, qpoints=None, ir_grid_points=Non """ (ir_grid_points_ref, weights_ref, ir_grid_map) = get_ir_grid_points(bz_grid) - ir_grid_points_ref = np.array(bz_grid.grg2bzg[ir_grid_points_ref], dtype="int_") _assert_grid_in_hdf5( weights, qpoints, ir_grid_points, weights_ref, ir_grid_points_ref, bz_grid @@ -300,13 +299,13 @@ def main(): if "weight" in f_kappa: frequencies = f_kappa["frequency"][:] ir_weights = f_kappa["weight"][:] - ir_grid_points = None + ir_grid_points_BZ = None qpoints = f_kappa["qpoint"][:] elif "ir_grid_weights" in f_kappa and not args.no_gridsym: ir_weights = f_kappa["ir_grid_weights"][:] - ir_grid_points = f_kappa["ir_grid_points"][:] + ir_grid_points_BZ = f_kappa["ir_grid_points"][:] qpoints = None - frequencies = np.array(f_kappa["frequency"][ir_grid_points], dtype="double") + frequencies = np.array(f_kappa["frequency"][ir_grid_points_BZ], dtype="double") else: frequencies = f_kappa["frequency"][:] ir_weights = np.ones(len(frequencies), dtype="int_") @@ -318,14 +317,14 @@ def main(): store_dense_gp_map=True, ) if args.no_gridsym or (ir_weights == 1).all(): - ir_grid_points = bz_grid.grg2bzg + ir_grid_points = None ir_grid_map = None else: ir_grid_points, ir_grid_map = _get_ir_grid_info( bz_grid, ir_weights, qpoints=qpoints, - ir_grid_points=ir_grid_points, + ir_grid_points=bz_grid.bzg2grg[ir_grid_points_BZ], ) conditions = frequencies > 0 if np.logical_not(conditions).sum() > 3: diff --git a/phono3py/cui/phono3py_script.py b/phono3py/cui/phono3py_script.py index 2f699db2..804aac86 100644 --- a/phono3py/cui/phono3py_script.py +++ b/phono3py/cui/phono3py_script.py @@ -243,6 +243,8 @@ def start_phono3py(**argparse_control) -> tuple[argparse.Namespace, int]: max_threads = phono3c.omp_max_threads() if max_threads > 0: print(f"Compiled with OpenMP support (max {max_threads} threads).") + if phono3c.include_lapacke(): + print("Compiled with LAPACKE.") if argparse_control.get("load_phono3py_yaml", False): print("Running in phono3py.load mode.") diff --git a/phono3py/other/kaccum.py b/phono3py/other/kaccum.py index edc81f70..86811cc2 100644 --- a/phono3py/other/kaccum.py +++ b/phono3py/other/kaccum.py @@ -17,17 +17,35 @@ class KappaDOSTHM: """Class to calculate DOS like spectram with tetrahedron method. - To compute usual DOS, + To compute usual DOS on all GR grid points: + ``` + freqs, _, _ = ph3.get_phonon_data() + freqs_grg = freqs[bzgrid.grg2bzg] kappados = KappaDOSTHM( - np.ones(freqs.shape)[None, :, :, None], - freqs, + np.ones(freqs_grg.shape, dtype=float)[None, :, :, None], + freqs_grg, bzgrid, - ir_grid_points, - ir_grid_weights=ir_weights, - ir_grid_map=ir_grid_map, num_sampling_points=201 ) + ``` + + To compute DOS on ir-grid points: + + ``` + freqs, _, _ = ph3.get_phonon_data() + ir_grid_points, ir_grid_weights, ir_grid_map = get_ir_grid_points(bzgrid) + freqs_ir = freqs[bzgrid.grg2bzg[ir_grid_points]] + kappados = KappaDOSTHM( + np.ones(freqs_ir.shape, dtype=float)[None, :, :, None], + freqs_ir, + bzgrid, + ir_grid_points=ir_grid_points, + ir_grid_weights=ir_grid_weights, + ir_grid_map=ir_grid_map, + num_sampling_points=201, + ) + ``` """ @@ -36,7 +54,7 @@ def __init__( mode_kappa: np.ndarray, frequencies: np.ndarray, bz_grid: BZGrid, - ir_grid_points: np.ndarray, + ir_grid_points: Optional[np.ndarray] = None, ir_grid_weights: Optional[np.ndarray] = None, ir_grid_map: Optional[np.ndarray] = None, frequency_points: Optional[Union[np.ndarray, Sequence]] = None, @@ -54,20 +72,24 @@ def __init__( Frequencies at ir-grid points. shape=(ir_grid_points, num_band), dtype='double' bz_grid : BZGrid + BZGrid instance. ir_grid_points : ndarray - Irreducible grid point indices in GR-grid. + Irreducible grid point indices in GR-grid (as obtained by + get_ir_grid_points). shape=(num_ir_grid_points, ), dtype='int_' ir_grid_weights : ndarray - Weights of irreducible grid points. Its sum is the number of - grid points in GR-grid (prod(D_diag)). + Weights of irreducible grid points. Its sum is the number of grid + points in GR-grid (prod(D_diag)) (as obtained by + get_ir_grid_points). shape=(num_ir_grid_points, ), dtype='int_' ir_grid_map : ndarray Index mapping table to irreducible grid points from all grid points - such as, [0, 0, 2, 3, 3, ...]. + in GR-grid such as, [0, 0, 2, 3, 3, ...]. (as obtained by + get_ir_grid_points). shape=(prod(D_diag), ), dtype='int_' frequency_points : array_like, optional, default=None - This is used as the frequency points. When None, - frequency points are created from `num_sampling_points`. + This is used as the frequency points. When None, frequency points + are created from `num_sampling_points`. num_sampling_points : int, optional, default=100 Number of uniform sampling points. @@ -85,13 +107,20 @@ def __init__( self._kdos = np.zeros( (n_temp, len(self._frequency_points), 2, n_elem), dtype="double" ) - + if ir_grid_points is None: + _ir_grid_points = np.arange(len(frequencies), dtype=int) + else: + _ir_grid_points = ir_grid_points + grid_points = bz_grid.grg2bzg[_ir_grid_points] if ir_grid_map is None: - bzgp2irgp_map = None + _ir_grid_map = np.arange(len(frequencies), dtype=int) else: - bzgp2irgp_map = self._get_bzgp2irgp_map(bz_grid, ir_grid_map) + _ir_grid_map = ir_grid_map + bzgp2irgp_map = self._get_bzgp2irgp_map( + bz_grid.bzg2grg, _ir_grid_map, _ir_grid_points + ) if ir_grid_weights is None: - grid_weights = np.ones(mode_kappa.shape[1]) + grid_weights = np.ones(mode_kappa.shape[1], dtype="int_") else: grid_weights = ir_grid_weights for j, function in enumerate(("J", "I")): @@ -99,7 +128,7 @@ def __init__( self._frequency_points, frequencies, bz_grid, - grid_points=bz_grid.grg2bzg[ir_grid_points], + grid_points=grid_points, bzgp2irgp_map=bzgp2irgp_map, function=function, ) @@ -124,11 +153,38 @@ def get_kdos(self): """ return self._frequency_points, self._kdos - def _get_bzgp2irgp_map(self, bz_grid, ir_grid_map): + def _get_bzgp2irgp_map(self, bzg2grg, ir_grid_map, ir_grid_points): + """Return mapping table from BZ-grid indices to ir-grid point indices. + + More precisely, return mapping table from grid points in BZ-grid to + indices of ir-grid points. The length of the set of the indices is the + number of the ir-grid points. + + Parameters + ---------- + bzg2grg : ndarray + Mapping table from BZ-grid to GR-grid. + shape=(len(all-BZ-grid-points), ), dtype='int_' + ir_grid_map : ndarray + Mapping table from all grid points to ir-grid points in GR-grid. + shape=(np.prod(D_diag), ), dtype='int_' + ir_grid_points : ndarray + Irreducible grid points in GR-grid. shape=(num_ir_grid_points, ), + dtype='int_' + + Returns + ------- + np.ndarray + Mapping table from BZ-grid to indices of ir-grid points. + shape=(len(ir-grid-points), ), dtype=' + + """ unique_gps = np.unique(ir_grid_map) + assert np.array_equal(unique_gps, ir_grid_points) + # ir-grid points in GR-grid to the index of unique grid points. gp_map = {j: i for i, j in enumerate(unique_gps)} bzgp2irgp_map = np.array( - [gp_map[ir_grid_map[grgp]] for grgp in bz_grid.bzg2grg], dtype="int_" + [gp_map[ir_grid_map[grgp]] for grgp in bzg2grg], dtype="int_" ) return bzgp2irgp_map @@ -236,7 +292,7 @@ def run_prop_dos( mode_prop, frequencies, bz_grid, - bz_grid.bzg2grg[ir_grid_points], + ir_grid_points=ir_grid_points, ir_grid_map=ir_grid_map, num_sampling_points=num_sampling_points, ) @@ -266,7 +322,7 @@ def run_mfp_dos( mode_prop[i : i + 1, :, :], mean_freepath[i], bz_grid, - bz_grid.bzg2grg[ir_grid_points], + ir_grid_points=ir_grid_points, ir_grid_map=ir_grid_map, num_sampling_points=num_sampling_points, ) diff --git a/phono3py/phonon/solver.py b/phono3py/phonon/solver.py index a694c5bd..6c216a49 100644 --- a/phono3py/phonon/solver.py +++ b/phono3py/phonon/solver.py @@ -53,6 +53,14 @@ def run_phonon_solver_c( ): """Bulid and solve dynamical matrices on grid in C-API. + Note + ---- + When LAPACKE is linked in C, `phononcalc.phonons_at_gridpoints` constucts + and solves dynamical matrices on grid points. Otherwise, it only constructs + dynamical matrices and solves them in python. + + Parameters + ---------- dm : DynamicalMatrix DynamicalMatrix instance. frequencies, eigenvectors, phonon_done : @@ -65,14 +73,15 @@ def run_phonon_solver_c( QDinv : ndarray See BZGrid.QDinv. frequency_conversion_factor : float, optional - Frequency convertion factor that is multiplied with - sqrt or eigenvalue of dynamical matrix. Default is VaspToTHz. + Frequency convertion factor that is multiplied with sqrt or eigenvalue + of dynamical matrix. Default is VaspToTHz. nac_q_direction : array_like, optional See Interaction.nac_q_direction. Default is None. lapack_zheev_uplo : str, optional 'U' or 'L' for lapack zheev solver. Default is 'L'. """ + import phono3py._phono3py as phono3c import phono3py._phononcalc as phononcalc ( @@ -124,6 +133,12 @@ def run_phonon_solver_c( assert QDinv.flags.c_contiguous assert lapack_zheev_uplo in ("L", "U") + if not phono3c.include_lapacke(): + # phonon_done is set even with phono3c.include_lapacke() == 0 for which + # dynamical matrices are not diagonalized in + # phononcalc.phonons_at_gridpoints. + phonon_undone = np.where(phonon_done == 0)[0] + fc_p2s, fc_s2p = _get_fc_elements_mapping(dm, fc) phononcalc.phonons_at_gridpoints( frequencies, @@ -154,6 +169,19 @@ def run_phonon_solver_c( lapack_zheev_uplo, ) + if not phono3c.include_lapacke(): + # The variable `eigenvectors` contains dynamical matrices. + # They are diagonalized in python as follows. + for gp in phonon_undone: + frequencies[gp], eigenvectors[gp] = np.linalg.eigh( + eigenvectors[gp], UPLO=lapack_zheev_uplo + ) + frequencies[phonon_undone] = ( + np.sign(frequencies[phonon_undone]) + * np.sqrt(np.abs(frequencies[phonon_undone])) + * frequency_conversion_factor + ) + def run_phonon_solver_py( grid_point, diff --git a/phono3py/version.py b/phono3py/version.py index 6dc2d7e2..23d47634 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.7.0" +__version__ = "3.8.0" diff --git a/pyproject.toml b/pyproject.toml index 094ad8eb..5836dc48 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -16,7 +16,7 @@ dependencies = [ "matplotlib>=2.2.2", "h5py>=3.0", "spglib>=2.3", - "phonopy>=2.30,<2.31", + "phonopy>=2.31,<2.32", ] license = { file = "LICENSE" } @@ -52,6 +52,7 @@ sdist.include = [ PHONO3PY_USE_MTBLAS = {env="PHONO3PY_USE_MTBLAS", default="ON"} USE_CONDA_PATH = {env="USE_CONDA_PATH", default="ON"} PHONO3PY_USE_OMP = {env="PHONO3PY_USE_OMP", default="ON"} +BUILD_WITHOUT_LAPACKE = {env="BUILD_WITHOUT_LAPACKE", default="OFF"} [tool.setuptools_scm] write_to = "phono3py/_version.py" diff --git a/test/conductivity/test_kappa_LBTE.py b/test/conductivity/test_kappa_LBTE.py index f494b0e0..10326417 100644 --- a/test/conductivity/test_kappa_LBTE.py +++ b/test/conductivity/test_kappa_LBTE.py @@ -1,12 +1,28 @@ """Tests for direct solution of LBTE.""" import numpy as np +import pytest +import phono3py._phono3py as phono3c from phono3py.api_phono3py import Phono3py -def test_kappa_LBTE(si_pbesol: Phono3py): +@pytest.mark.skipif( + not phono3c.include_lapacke(), reason="test for phono3py compliled with lapacke" +) +@pytest.mark.parametrize("pinv_solver", [1, 2]) +def test_kappa_LBTE_12(si_pbesol: Phono3py, pinv_solver: int): """Test for symmetry reduced collision matrix.""" + _test_kappa_LBTE(si_pbesol, pinv_solver) + + +@pytest.mark.parametrize("pinv_solver", [3, 4, 5]) +def test_kappa_LBTE_345(si_pbesol: Phono3py, pinv_solver: int): + """Test for symmetry reduced collision matrix.""" + _test_kappa_LBTE(si_pbesol, pinv_solver) + + +def _test_kappa_LBTE(si_pbesol: Phono3py, pinv_solver: int): if si_pbesol._make_r0_average: ref_kappa = [110.896, 110.896, 110.896, 0, 0, 0] else: @@ -18,11 +34,30 @@ def test_kappa_LBTE(si_pbesol: Phono3py): temperatures=[ 300, ], + pinv_solver=pinv_solver, ) kappa = si_pbesol.thermal_conductivity.kappa.ravel() np.testing.assert_allclose(ref_kappa, kappa, atol=0.3) +@pytest.mark.skipif( + phono3c.include_lapacke(), reason="test for phono3py compliled without lapacke" +) +@pytest.mark.parametrize("pinv_solver", [1, 2]) +def test_kappa_LBTE_witout_lapacke(si_pbesol: Phono3py, pinv_solver: int): + """Test for symmetry reduced collision matrix.""" + si_pbesol.mesh_numbers = [9, 9, 9] + si_pbesol.init_phph_interaction() + with pytest.raises(RuntimeError): + si_pbesol.run_thermal_conductivity( + is_LBTE=True, + temperatures=[ + 300, + ], + pinv_solver=pinv_solver, + ) + + def test_kappa_LBTE_full_colmat(si_pbesol: Phono3py): """Test for full collision matrix.""" if si_pbesol._make_r0_average: @@ -40,7 +75,7 @@ def test_kappa_LBTE_full_colmat(si_pbesol: Phono3py): is_reducible_collision_matrix=True, ) kappa = si_pbesol.thermal_conductivity.kappa.ravel() - np.testing.assert_allclose(ref_kappa, kappa, atol=0.3) + np.testing.assert_allclose(ref_kappa, kappa, atol=0.5) def test_kappa_LBTE_aln(aln_lda: Phono3py): diff --git a/test/other/test_kaccum.py b/test/other/test_kaccum.py index becd2b13..67cee18a 100644 --- a/test/other/test_kaccum.py +++ b/test/other/test_kaccum.py @@ -323,51 +323,51 @@ def test_run_prop_dos(si_pbesol: Phono3py): ref_kdos = [ 0.00000, 0.00000, - 2.19162, - 5.16697, - 28.22125, - 18.97280, - 58.56343, - 12.19206, - 69.05896, - 3.47035, - 73.17626, - 1.48915, - 74.74544, - 0.43485, - 75.87064, - 1.74135, - 79.08179, - 2.30428, - 81.21678, + 2.19027, + 5.16378, + 25.73771, + 15.50548, + 56.63845, + 19.48634, + 68.64061, + 3.16016, + 72.58422, + 1.62370, + 74.53304, + 0.78480, + 76.94719, + 5.37456, + 80.63569, + 0.61517, + 81.15021, 0.00000, ] # print(",".join([f"{v:10.5f}" for v in mfp[0, :, :, 0].ravel()])) ref_mfp = [ 0.00000, 0.00000, - 29.19150, - 0.02604, - 42.80717, - 0.01202, - 52.09457, - 0.01158, - 61.79908, - 0.01140, - 69.49177, - 0.00784, - 74.57499, - 0.00501, - 77.99145, - 0.00364, - 80.33477, - 0.00210, - 81.21678, + 33.66617, + 0.02216, + 44.81836, + 0.01032, + 53.09252, + 0.01062, + 62.18077, + 0.01085, + 69.59018, + 0.00765, + 74.58994, + 0.00496, + 77.95139, + 0.00359, + 80.27019, + 0.00209, + 81.15021, 0.00000, ] # print(",".join([f"{v:10.5f}" for v in sampling_points[0]])) ref_sampling_points = [ - -0.00000, + 0.00000, 1.69664, 3.39328, 5.08992, @@ -381,15 +381,15 @@ def test_run_prop_dos(si_pbesol: Phono3py): # print(",".join([f"{v:10.5f}" for v in sampling_points_mfp[0]])) ref_sampling_points_mfp = [ 0.00000, - 803.91710, - 1607.83420, - 2411.75130, - 3215.66841, - 4019.58551, - 4823.50261, - 5627.41971, - 6431.33681, - 7235.25391, + 803.28312, + 1606.56625, + 2409.84937, + 3213.13249, + 4016.41562, + 4819.69874, + 5622.98186, + 6426.26499, + 7229.54811, ] np.testing.assert_allclose(ref_kdos, kdos[0, :, :, 0].ravel(), atol=0.5) np.testing.assert_allclose(ref_mfp, mfp[0, :, :, 0].ravel(), atol=0.5) diff --git a/test/phonon3/test_fc3.py b/test/phonon3/test_fc3.py index 5f358277..bfe9e675 100644 --- a/test/phonon3/test_fc3.py +++ b/test/phonon3/test_fc3.py @@ -3,6 +3,7 @@ import numpy as np import pytest +import phono3py._phono3py as phono3c from phono3py import Phono3py from phono3py.phonon3.fc3 import ( cutoff_fc3_by_zero, @@ -257,40 +258,43 @@ def test_phonon_smat_alm_cutoff_fc3(si_pbesol_111_222_alm_cutoff_fc3: Phono3py): np.testing.assert_allclose(ph.fc2[0, 33], fc2_ref, atol=1e-6, rtol=0) -@pytest.mark.parametrize("pinv_solver", ["numpy", "lapacke"]) -def test_fc3_lapacke_solver(si_pbesol_111: Phono3py, pinv_solver: str): +@pytest.mark.skipif( + not phono3c.include_lapacke(), reason="requires to complile with lapacke" +) +def test_fc3_lapacke_solver(si_pbesol_111: Phono3py): """Test fc3 with Si PBEsol 1x1x1 using lapacke solver.""" - ph = si_pbesol_111 - _, fc3 = get_fc3( - ph.supercell, - ph.primitive, - ph.dataset, - ph.symmetry, - pinv_solver=pinv_solver, - verbose=True, - ) - set_translational_invariance_fc3(fc3) - set_permutation_symmetry_fc3(fc3) - - fc3_ref = [ - [ - [1.07250822e-01, 1.86302073e-17, -4.26452855e-18], - [8.96414569e-03, -1.43046911e-01, -1.38498937e-01], - [-8.96414569e-03, -1.38498937e-01, -1.43046911e-01], - ], - [ - [-8.96414569e-03, -1.43046911e-01, -1.38498937e-01], - [-3.39457157e-02, -4.63315728e-17, -4.17779237e-17], - [-3.31746167e-01, -2.60025724e-02, -2.60025724e-02], - ], - [ - [8.96414569e-03, -1.38498937e-01, -1.43046911e-01], - [-3.31746167e-01, 2.60025724e-02, 2.60025724e-02], - [-3.39457157e-02, 3.69351540e-17, 5.94504191e-18], - ], - ] - - np.testing.assert_allclose(fc3[0, 1, 7], fc3_ref, atol=1e-8, rtol=0) + for pinv_solver in ["lapacke", "numpy"]: + ph = si_pbesol_111 + _, fc3 = get_fc3( + ph.supercell, + ph.primitive, + ph.dataset, + ph.symmetry, + pinv_solver=pinv_solver, + verbose=True, + ) + set_translational_invariance_fc3(fc3) + set_permutation_symmetry_fc3(fc3) + + fc3_ref = [ + [ + [1.07250822e-01, 1.86302073e-17, -4.26452855e-18], + [8.96414569e-03, -1.43046911e-01, -1.38498937e-01], + [-8.96414569e-03, -1.38498937e-01, -1.43046911e-01], + ], + [ + [-8.96414569e-03, -1.43046911e-01, -1.38498937e-01], + [-3.39457157e-02, -4.63315728e-17, -4.17779237e-17], + [-3.31746167e-01, -2.60025724e-02, -2.60025724e-02], + ], + [ + [8.96414569e-03, -1.38498937e-01, -1.43046911e-01], + [-3.31746167e-01, 2.60025724e-02, 2.60025724e-02], + [-3.39457157e-02, 3.69351540e-17, 5.94504191e-18], + ], + ] + + np.testing.assert_allclose(fc3[0, 1, 7], fc3_ref, atol=1e-8, rtol=0) def test_fc3_symfc(si_pbesol_111_symfc: Phono3py):