Use nanobind for C binding

This commit is contained in:
Atsushi Togo 2024-07-04 17:55:36 +09:00
parent a0873ec6cd
commit 2731d8c3f4
26 changed files with 940 additions and 118 deletions

View File

@ -1,73 +1,128 @@
cmake_minimum_required(VERSION 3.0)
cmake_policy(SET CMP0007 NEW)
cmake_minimum_required(VERSION 3.20)
project(phonopy C)
if(NOT SKBUILD_PROJECT_NAME)
set(SKBUILD_PROJECT_NAME phonopy)
endif()
if(PHONOPY_LIBS)
project(phonopy LANGUAGES C)
set(CMAKE_MACOSX_RPATH 1)
set(CMAKE_C_FLAGS_RELEASE "-Wall -O2")
set(CMAKE_C_FLAGS_DEBUG "-g -DTHMWARNING")
if(NOT CMAKE_BUILD_TYPE)
set(CMAKE_BUILD_TYPE Release)
endif(NOT CMAKE_BUILD_TYPE)
message(STATUS "Build type: ${CMAKE_BUILD_TYPE}")
message(STATUS "CMAKE_SYSTEM_PREFIX_PATH: ${CMAKE_SYSTEM_PREFIX_PATH}")
include(GNUInstallDirs)
# set(CMAKE_POSITION_INDEPENDENT_CODE ON)
# Version numbers
file(READ ${PROJECT_SOURCE_DIR}/phonopy/version.py version_file)
string(REGEX MATCH "__version__.*([0-9]+)[.]([0-9]+)[.]([0-9]+)" phono3py_version ${version_file})
string(REGEX MATCH "__version__.*([0-9]+)[.]([0-9]+)[.]([0-9]+)"
phono3py_version ${version_file})
set(phonopy_major_version ${CMAKE_MATCH_1})
set(phonopy_minor_version ${CMAKE_MATCH_2})
set(phonopy_micro_version ${CMAKE_MATCH_3})
set(serial "${phonopy_major_version}.${phonopy_minor_version}.${phonopy_micro_version}")
set(serial
"${phonopy_major_version}.${phonopy_minor_version}.${phonopy_micro_version}"
)
set(soserial "1")
else()
project(${SKBUILD_PROJECT_NAME})
set(DEV_MODULE Development.Module)
find_package(
Python 3.8
REQUIRED COMPONENTS Interpreter ${DEV_MODULE}
OPTIONAL_COMPONENTS Development.SABIModule
)
endif()
cmake_policy(SET CMP0007 NEW)
if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
set(CMAKE_BUILD_TYPE
Release
CACHE STRING "Choose the type of build." FORCE)
set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS "Debug" "Release"
"MinSizeRel" "RelWithDebInfo")
endif()
message(STATUS "Build type: ${CMAKE_BUILD_TYPE}")
message(STATUS "CMAKE_SYSTEM_PREFIX_PATH: ${CMAKE_SYSTEM_PREFIX_PATH}")
find_package(OpenMP)
if(OpenMP_FOUND)
message(STATUS "OpenMP libs: ${OpenMP_C_LIBRARIES}")
message(STATUS "OpenMP flags: ${OpenMP_C_FLAGS}")
endif()
if (PHONOPY)
set(SOURCES_PHONOPY
${PROJECT_SOURCE_DIR}/c/phonopy.c
${PROJECT_SOURCE_DIR}/c/dynmat.c
${PROJECT_SOURCE_DIR}/c/derivative_dynmat.c
${PROJECT_SOURCE_DIR}/c/rgrid.c
${PROJECT_SOURCE_DIR}/c/phonopy.c ${PROJECT_SOURCE_DIR}/c/dynmat.c
${PROJECT_SOURCE_DIR}/c/derivative_dynmat.c ${PROJECT_SOURCE_DIR}/c/rgrid.c
${PROJECT_SOURCE_DIR}/c/tetrahedron_method.c)
if (BUILD_SHARED_LIBRARIES)
if(BUILD_SHARED_LIBS)
# Shared library
add_library(phpy SHARED ${SOURCES_PHONOPY})
add_library(phonopy_libs SHARED ${SOURCES_PHONOPY})
if(OpenMP_FOUND)
target_link_libraries(phpy m ${OpenMP_C_LIBRARIES})
target_compile_options(phpy PRIVATE ${OpenMP_C_FLAGS})
target_link_libraries(
phonopy_libs
PRIVATE m
PRIVATE OpenMP::OpenMP_C)
else()
target_link_libraries(phpy m)
endif()
target_include_directories(phpy PRIVATE ${MY_INCLUDES})
target_compile_definitions(phpy PRIVATE THM_EPSILON=1e-10)
set_property(TARGET phpy PROPERTY VERSION ${serial})
set_property(TARGET phpy PROPERTY SOVERSION ${soserial})
install(TARGETS phpy LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR})
target_link_libraries(phonopy_libs PRIVATE m)
endif()
target_compile_definitions(phonopy_libs PRIVATE THM_EPSILON=1e-10)
else()
# Static link library
add_library(phpy_static STATIC ${SOURCES_PHONOPY})
if (OpenMP_FOUND)
target_link_libraries(phpy_static m ${OpenMP_C_LIBRARIES})
target_compile_options(phpy_static PRIVATE ${OpenMP_C_FLAGS})
else()
target_link_libraries(phpy_static m)
endif()
target_include_directories(phpy_static PRIVATE ${MY_INCLUDES})
target_compile_definitions(phpy_static PRIVATE THM_EPSILON=1e-10)
set_property(TARGET phpy_static PROPERTY VERSION ${serial})
set_property(TARGET phpy_static PROPERTY SOVERSION ${soserial})
set_property(TARGET phpy_static PROPERTY OUTPUT_NAME phpy)
install(TARGETS phpy_static ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR})
add_library(phonopy_libs STATIC ${SOURCES_PHONOPY})
# Header file
install(FILES ${PROJECT_SOURCE_DIR}/c/phonopy.h DESTINATION ${CMAKE_INSTALL_INCLUDEDIR})
if(OpenMP_FOUND)
target_link_libraries(
phonopy_libs
PRIVATE m
PRIVATE OpenMP::OpenMP_C)
else()
target_link_libraries(phonopy_libs PRIVATE m)
endif()
target_compile_definitions(phonopy_libs PRIVATE THM_EPSILON=1e-10)
endif()
if(PHONOPY_LIBS)
install(FILES ${PROJECT_SOURCE_DIR}/c/phonopy.h
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR})
if(BUILD_SHARED_LIBS)
install(TARGETS phonopy_libs LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR})
set_property(TARGET phonopy_libs PROPERTY VERSION ${serial})
set_property(TARGET phonopy_libs PROPERTY SOVERSION ${soserial})
else()
install(TARGETS phonopy_libs ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR})
set_property(TARGET phonopy_libs PROPERTY VERSION ${serial})
set_property(TARGET phonopy_libs PROPERTY SOVERSION ${soserial})
set_property(TARGET phonopy_libs PROPERTY OUTPUT_NAME phonopy_libs)
endif()
else()
set_target_properties(phonopy_libs PROPERTIES POSITION_INDEPENDENT_CODE ON)
execute_process(
COMMAND "${Python_EXECUTABLE}" -m nanobind --cmake_dir
OUTPUT_STRIP_TRAILING_WHITESPACE
OUTPUT_VARIABLE NB_DIR)
list(APPEND CMAKE_PREFIX_PATH "${NB_DIR}")
find_package(nanobind CONFIG REQUIRED)
nanobind_add_module(
_phonopy
STABLE_ABI
${PROJECT_SOURCE_DIR}/c/phonopy.h
${PROJECT_SOURCE_DIR}/c/_phonopy.cpp
)
if(OpenMP_FOUND)
target_link_libraries(_phonopy PRIVATE phonopy_libs PRIVATE OpenMP::OpenMP_C)
else()
target_link_libraries(_phonopy PRIVATE phonopy_libs)
endif()
target_compile_definitions(_phonopy PRIVATE THM_EPSILON=1e-10)
install(TARGETS _phonopy LIBRARY DESTINATION ${SKBUILD_PROJECT_NAME})
endif()

681
c/_phonopy.cpp Normal file
View File

@ -0,0 +1,681 @@
#include <nanobind/nanobind.h>
#include <nanobind/ndarray.h>
#include "phonopy.h"
namespace nb = nanobind;
void py_transform_dynmat_to_fc(nb::ndarray<> py_force_constants,
nb::ndarray<> py_dynamical_matrices,
nb::ndarray<> py_commensurate_points,
nb::ndarray<> py_svecs, nb::ndarray<> py_multi,
nb::ndarray<> py_masses,
nb::ndarray<> py_s2pp_map,
nb::ndarray<> py_fc_index_map, long use_openmp) {
double *fc;
double(*dm)[2];
double(*comm_points)[3];
double(*svecs)[3];
double *masses;
long(*multi)[2];
long *s2pp_map;
long *fc_index_map;
long num_patom;
long num_satom;
fc = (double *)py_force_constants.data();
dm = (double(*)[2])py_dynamical_matrices.data();
comm_points = (double(*)[3])py_commensurate_points.data();
svecs = (double(*)[3])py_svecs.data();
masses = (double *)py_masses.data();
multi = (long(*)[2])py_multi.data();
s2pp_map = (long *)py_s2pp_map.data();
fc_index_map = (long *)py_fc_index_map.data();
num_patom = py_multi.shape(1);
num_satom = py_multi.shape(0);
phpy_transform_dynmat_to_fc(fc, dm, comm_points, svecs, multi, masses,
s2pp_map, fc_index_map, num_patom, num_satom,
use_openmp);
};
void py_perm_trans_symmetrize_fc(nb::ndarray<> py_force_constants, long level) {
double *fc;
int n_satom;
fc = (double *)py_force_constants.data();
n_satom = py_force_constants.shape(0);
phpy_perm_trans_symmetrize_fc(fc, n_satom, level);
}
void py_perm_trans_symmetrize_compact_fc(nb::ndarray<> py_force_constants,
nb::ndarray<> py_permutations,
nb::ndarray<> py_s2pp_map,
nb::ndarray<> py_p2s_map,
nb::ndarray<> py_nsym_list,
long level) {
double *fc;
int *perms;
int *s2pp;
int *p2s;
int *nsym_list;
int n_patom, n_satom;
fc = (double *)py_force_constants.data();
perms = (int *)py_permutations.data();
s2pp = (int *)py_s2pp_map.data();
p2s = (int *)py_p2s_map.data();
nsym_list = (int *)py_nsym_list.data();
n_patom = py_force_constants.shape(0);
n_satom = py_force_constants.shape(1);
phpy_perm_trans_symmetrize_compact_fc(fc, p2s, s2pp, nsym_list, perms,
n_satom, n_patom, level);
}
void py_transpose_compact_fc(nb::ndarray<> py_force_constants,
nb::ndarray<> py_permutations,
nb::ndarray<> py_s2pp_map,
nb::ndarray<> py_p2s_map,
nb::ndarray<> py_nsym_list) {
double *fc;
int *s2pp;
int *p2s;
int *nsym_list;
int *perms;
int n_patom, n_satom;
fc = (double *)py_force_constants.data();
perms = (int *)py_permutations.data();
s2pp = (int *)py_s2pp_map.data();
p2s = (int *)py_p2s_map.data();
nsym_list = (int *)py_nsym_list.data();
n_patom = py_force_constants.shape(0);
n_satom = py_force_constants.shape(1);
phpy_set_index_permutation_symmetry_compact_fc(fc, p2s, s2pp, nsym_list,
perms, n_satom, n_patom, 1);
}
void py_get_dynamical_matrices_with_dd_openmp_over_qpoints(
nb::ndarray<> py_dynamical_matrix, nb::ndarray<> py_qpoints,
nb::ndarray<> py_force_constants, nb::ndarray<> py_svecs,
nb::ndarray<> py_multi, nb::ndarray<> py_positions, nb::ndarray<> py_masses,
nb::ndarray<> py_s2p_map, nb::ndarray<> py_p2s_map,
nb::ndarray<> py_q_direction, nb::ndarray<> py_born,
nb::ndarray<> py_dielectric, nb::ndarray<> py_reciprocal_lattice,
double nac_factor, nb::ndarray<> py_dd_q0, nb::ndarray<> py_G_list,
double lambda, long use_Wang_NAC) {
double(*dm)[2];
double *fc;
double *q_direction;
double(*qpoints)[3];
double(*svecs)[3];
long(*multi)[2];
double(*positions)[3];
double *masses;
double(*born)[3][3];
double(*dielectric)[3];
double(*reciprocal_lattice)[3];
double(*dd_q0)[2];
double(*G_list)[3];
long *s2p_map;
long *p2s_map;
long num_patom;
long num_satom;
long n_qpoints;
long n_Gpoints;
dm = (double(*)[2])py_dynamical_matrix.data();
qpoints = (double(*)[3])py_qpoints.data();
n_qpoints = py_qpoints.shape(0);
fc = (double *)py_force_constants.data();
svecs = (double(*)[3])py_svecs.data();
multi = (long(*)[2])py_multi.data();
masses = (double *)py_masses.data();
s2p_map = (long *)py_s2p_map.data();
p2s_map = (long *)py_p2s_map.data();
q_direction = (double *)py_q_direction.data();
born = (double(*)[3][3])py_born.data();
dielectric = (double(*)[3])py_dielectric.data();
reciprocal_lattice = (double(*)[3])py_reciprocal_lattice.data();
if (use_Wang_NAC) {
positions = NULL;
dd_q0 = NULL;
G_list = NULL;
n_Gpoints = 0;
} else {
positions = (double(*)[3])py_positions.data();
dd_q0 = (double(*)[2])py_dd_q0.data();
G_list = (double(*)[3])py_G_list.data();
n_Gpoints = py_G_list.shape(0);
}
num_patom = py_p2s_map.shape(0);
num_satom = py_s2p_map.shape(0);
phpy_dynamical_matrices_with_dd_openmp_over_qpoints(
dm, qpoints, n_qpoints, fc, svecs, multi, positions, num_patom,
num_satom, masses, p2s_map, s2p_map, born, dielectric,
reciprocal_lattice, q_direction, nac_factor, dd_q0, G_list, n_Gpoints,
lambda, use_Wang_NAC);
}
void py_get_dynamical_matrix(nb::ndarray<> py_dynamical_matrix,
nb::ndarray<> py_force_constants,
nb::ndarray<> py_q, nb::ndarray<> py_svecs,
nb::ndarray<> py_multi, nb::ndarray<> py_masses,
nb::ndarray<> py_s2p_map, nb::ndarray<> py_p2s_map,
long use_openmp) {
double(*dm)[2];
double *fc;
double *q;
double(*svecs)[3];
double *m;
long(*multi)[2];
long *s2p_map;
long *p2s_map;
long num_patom;
long num_satom;
dm = (double(*)[2])py_dynamical_matrix.data();
fc = (double *)py_force_constants.data();
q = (double *)py_q.data();
svecs = (double(*)[3])py_svecs.data();
m = (double *)py_masses.data();
multi = (long(*)[2])py_multi.data();
s2p_map = (long *)py_s2p_map.data();
p2s_map = (long *)py_p2s_map.data();
num_patom = py_p2s_map.shape(0);
num_satom = py_s2p_map.shape(0);
if (py_q.ndim() == 1) {
phpy_get_dynamical_matrix_at_q(dm, num_patom, num_satom, fc, q, svecs,
multi, m, s2p_map, p2s_map, NULL,
use_openmp);
} else {
phpy_get_dynamical_matrices_openmp_over_qpoints(
dm, num_patom, num_satom, fc, (double(*)[3])q, py_q.shape(0), svecs,
multi, m, s2p_map, p2s_map, NULL);
}
}
void py_get_nac_dynamical_matrix(
nb::ndarray<> py_dynamical_matrix, nb::ndarray<> py_force_constants,
nb::ndarray<> py_q, nb::ndarray<> py_svecs, nb::ndarray<> py_multi,
nb::ndarray<> py_masses, nb::ndarray<> py_s2p_map, nb::ndarray<> py_p2s_map,
nb::ndarray<> py_q_cart, nb::ndarray<> py_born, double factor,
long use_openmp) {
double(*dm)[2];
double *fc;
double *q_cart;
double *q;
double(*svecs)[3];
double *m;
double(*born)[3][3];
long(*multi)[2];
long *s2p_map;
long *p2s_map;
long num_patom;
long num_satom;
long n;
double(*charge_sum)[3][3];
dm = (double(*)[2])py_dynamical_matrix.data();
fc = (double *)py_force_constants.data();
q_cart = (double *)py_q_cart.data();
q = (double *)py_q.data();
svecs = (double(*)[3])py_svecs.data();
m = (double *)py_masses.data();
born = (double(*)[3][3])py_born.data();
multi = (long(*)[2])py_multi.data();
s2p_map = (long *)py_s2p_map.data();
p2s_map = (long *)py_p2s_map.data();
num_patom = py_p2s_map.shape(0);
num_satom = py_s2p_map.shape(0);
charge_sum =
(double(*)[3][3])malloc(sizeof(double[3][3]) * num_patom * num_patom);
n = num_satom / num_patom;
phpy_get_charge_sum(charge_sum, num_patom, factor / n, q_cart, born);
if (py_q.ndim() == 1) {
phpy_get_dynamical_matrix_at_q(dm, num_patom, num_satom, fc, q, svecs,
multi, m, s2p_map, p2s_map, charge_sum,
use_openmp);
} else {
phpy_get_dynamical_matrices_openmp_over_qpoints(
dm, num_patom, num_satom, fc, (double(*)[3])q, py_q.shape(0), svecs,
multi, m, s2p_map, p2s_map, charge_sum);
}
free(charge_sum);
charge_sum = NULL;
}
void py_get_recip_dipole_dipole(
nb::ndarray<> py_dd, nb::ndarray<> py_dd_q0, nb::ndarray<> py_G_list,
nb::ndarray<> py_q_cart, nb::ndarray<> py_q_direction,
nb::ndarray<> py_born, nb::ndarray<> py_dielectric,
nb::ndarray<> py_positions, long is_q_zero, double factor, double lambda,
double tolerance, long use_openmp) {
double(*dd)[2];
double(*dd_q0)[2];
double(*G_list)[3];
double *q_vector;
double *q_direction;
double(*born)[3][3];
double(*dielectric)[3];
double(*pos)[3];
long num_patom, num_G;
dd = (double(*)[2])py_dd.data();
dd_q0 = (double(*)[2])py_dd_q0.data();
G_list = (double(*)[3])py_G_list.data();
if (is_q_zero) {
q_direction = NULL;
} else {
q_direction = (double *)py_q_direction.data();
}
q_vector = (double *)py_q_cart.data();
born = (double(*)[3][3])py_born.data();
dielectric = (double(*)[3])py_dielectric.data();
pos = (double(*)[3])py_positions.data();
num_G = py_G_list.shape(0);
num_patom = py_positions.shape(0);
phpy_get_recip_dipole_dipole(dd, /* [natom, 3, natom, 3, (real, imag)] */
dd_q0, /* [natom, 3, 3, (real, imag)] */
G_list, /* [num_kvec, 3] */
num_G, num_patom, q_vector, q_direction, born,
dielectric, pos, /* [natom, 3] */
factor, /* 4pi/V*unit-conv */
lambda, /* 4 * Lambda^2 */
tolerance, use_openmp);
}
void py_get_recip_dipole_dipole_q0(nb::ndarray<> py_dd_q0,
nb::ndarray<> py_G_list,
nb::ndarray<> py_born,
nb::ndarray<> py_dielectric,
nb::ndarray<> py_positions, double lambda,
double tolerance, long use_openmp) {
double(*dd_q0)[2];
double(*G_list)[3];
double(*born)[3][3];
double(*dielectric)[3];
double(*pos)[3];
long num_patom, num_G;
dd_q0 = (double(*)[2])py_dd_q0.data();
G_list = (double(*)[3])py_G_list.data();
born = (double(*)[3][3])py_born.data();
dielectric = (double(*)[3])py_dielectric.data();
pos = (double(*)[3])py_positions.data();
num_G = py_G_list.shape(0);
num_patom = py_positions.shape(0);
phpy_get_recip_dipole_dipole_q0(dd_q0, /* [natom, 3, 3, (real, imag)] */
G_list, /* [num_kvec, 3] */
num_G, num_patom, born, dielectric,
pos, /* [natom, 3] */
lambda, /* 4 * Lambda^2 */
tolerance, use_openmp);
}
void py_get_derivative_dynmat(
nb::ndarray<> py_derivative_dynmat, nb::ndarray<> py_force_constants,
nb::ndarray<> py_q_vector, nb::ndarray<> py_lattice, nb::ndarray<> py_svecs,
nb::ndarray<> py_multi, nb::ndarray<> py_masses, nb::ndarray<> py_s2p_map,
nb::ndarray<> py_p2s_map, double nac_factor, nb::ndarray<> py_born,
nb::ndarray<> py_dielectric, nb::ndarray<> py_q_direction, long is_nac,
long is_q_zero, long use_openmp) {
double(*ddm)[2];
double *fc;
double *q_vector;
double *lat;
double(*svecs)[3];
double *masses;
long(*multi)[2];
long *s2p_map;
long *p2s_map;
long num_patom;
long num_satom;
double *born;
double *epsilon;
double *q_dir;
ddm = (double(*)[2])py_derivative_dynmat.data();
fc = (double *)py_force_constants.data();
q_vector = (double *)py_q_vector.data();
lat = (double *)py_lattice.data();
svecs = (double(*)[3])py_svecs.data();
masses = (double *)py_masses.data();
multi = (long(*)[2])py_multi.data();
s2p_map = (long *)py_s2p_map.data();
p2s_map = (long *)py_p2s_map.data();
num_patom = py_p2s_map.shape(0);
num_satom = py_s2p_map.shape(0);
epsilon = (double *)py_dielectric.data();
born = (double *)py_born.data();
if (is_q_zero) {
q_dir = NULL;
} else {
q_dir = (double *)py_q_direction.data();
}
phpy_get_derivative_dynmat_at_q(
ddm, num_patom, num_satom, fc, q_vector, lat, svecs, multi, masses,
s2p_map, p2s_map, nac_factor, born, epsilon, q_dir, is_nac, use_openmp);
}
void py_get_thermal_properties(nb::ndarray<> py_thermal_props,
nb::ndarray<> py_temperatures,
nb::ndarray<> py_frequencies,
nb::ndarray<> py_weights,
double cutoff_frequency) {
double *temperatures;
double *freqs;
double *thermal_props;
long *weights;
long num_qpoints;
long num_bands;
long num_temp;
thermal_props = (double *)py_thermal_props.data();
temperatures = (double *)py_temperatures.data();
num_temp = (long)py_temperatures.shape(0);
freqs = (double *)py_frequencies.data();
num_qpoints = (long)py_frequencies.shape(0);
weights = (long *)py_weights.data();
num_bands = (long)py_frequencies.shape(1);
phpy_get_thermal_properties(thermal_props, temperatures, freqs, weights,
num_temp, num_qpoints, num_bands,
cutoff_frequency);
}
void py_distribute_fc2(nb::ndarray<> py_force_constants,
nb::ndarray<> py_atom_list,
nb::ndarray<> py_fc_indices_of_atom_list,
nb::ndarray<> py_rotations_cart,
nb::ndarray<> py_permutations,
nb::ndarray<> py_map_atoms, nb::ndarray<> py_map_syms) {
double(*r_carts)[3][3];
double(*fc2)[3][3];
int *permutations;
int *map_atoms;
int *map_syms;
int *atom_list;
int *fc_indices_of_atom_list;
long num_pos, num_rot, len_atom_list;
fc2 = (double(*)[3][3])py_force_constants.data();
atom_list = (int *)py_atom_list.data();
len_atom_list = py_atom_list.shape(0);
fc_indices_of_atom_list = (int *)py_fc_indices_of_atom_list.data();
permutations = (int *)py_permutations.data();
map_atoms = (int *)py_map_atoms.data();
map_syms = (int *)py_map_syms.data();
r_carts = (double(*)[3][3])py_rotations_cart.data();
num_rot = py_permutations.shape(0);
num_pos = py_permutations.shape(1);
phpy_distribute_fc2(fc2, atom_list, len_atom_list, fc_indices_of_atom_list,
r_carts, permutations, map_atoms, map_syms, num_rot,
num_pos);
}
bool py_compute_permutation(nb::ndarray<> permutation, nb::ndarray<> lattice,
nb::ndarray<> positions,
nb::ndarray<> permuted_positions, double symprec) {
int *rot_atoms;
double(*lat)[3];
double(*pos)[3];
double(*rot_pos)[3];
int num_pos;
int is_found;
rot_atoms = (int *)permutation.data();
lat = (double(*)[3])lattice.data();
pos = (double(*)[3])positions.data();
rot_pos = (double(*)[3])permuted_positions.data();
num_pos = positions.shape(0);
is_found = phpy_compute_permutation(rot_atoms, lat, pos, rot_pos, num_pos,
symprec);
if (is_found) {
return true;
} else {
return false;
}
}
void py_gsv_set_smallest_vectors_sparse(
nb::ndarray<> py_smallest_vectors, nb::ndarray<> py_multiplicity,
nb::ndarray<> py_pos_to, nb::ndarray<> py_pos_from,
nb::ndarray<> py_lattice_points, nb::ndarray<> py_reduced_basis,
nb::ndarray<> py_trans_mat, double symprec) {
double(*smallest_vectors)[27][3];
int *multiplicity;
double(*pos_to)[3];
double(*pos_from)[3];
int(*lattice_points)[3];
double(*reduced_basis)[3];
int(*trans_mat)[3];
int num_pos_to, num_pos_from, num_lattice_points;
smallest_vectors = (double(*)[27][3])py_smallest_vectors.data();
multiplicity = (int *)py_multiplicity.data();
pos_to = (double(*)[3])py_pos_to.data();
pos_from = (double(*)[3])py_pos_from.data();
num_pos_to = py_pos_to.shape(0);
num_pos_from = py_pos_from.shape(0);
lattice_points = (int(*)[3])py_lattice_points.data();
num_lattice_points = py_lattice_points.shape(0);
reduced_basis = (double(*)[3])py_reduced_basis.data();
trans_mat = (int(*)[3])py_trans_mat.data();
phpy_set_smallest_vectors_sparse(smallest_vectors, multiplicity, pos_to,
num_pos_to, pos_from, num_pos_from,
lattice_points, num_lattice_points,
reduced_basis, trans_mat, symprec);
}
void py_gsv_set_smallest_vectors_dense(
nb::ndarray<> py_smallest_vectors, nb::ndarray<> py_multiplicity,
nb::ndarray<> py_pos_to, nb::ndarray<> py_pos_from,
nb::ndarray<> py_lattice_points, nb::ndarray<> py_reduced_basis,
nb::ndarray<> py_trans_mat, long initialize, double symprec) {
double(*smallest_vectors)[3];
long(*multiplicity)[2];
double(*pos_to)[3];
double(*pos_from)[3];
long(*lattice_points)[3];
double(*reduced_basis)[3];
long(*trans_mat)[3];
long num_pos_to, num_pos_from, num_lattice_points;
smallest_vectors = (double(*)[3])py_smallest_vectors.data();
multiplicity = (long(*)[2])py_multiplicity.data();
pos_to = (double(*)[3])py_pos_to.data();
pos_from = (double(*)[3])py_pos_from.data();
num_pos_to = py_pos_to.shape(0);
num_pos_from = py_pos_from.shape(0);
lattice_points = (long(*)[3])py_lattice_points.data();
num_lattice_points = py_lattice_points.shape(0);
reduced_basis = (double(*)[3])py_reduced_basis.data();
trans_mat = (long(*)[3])py_trans_mat.data();
phpy_set_smallest_vectors_dense(
smallest_vectors, multiplicity, pos_to, num_pos_to, pos_from,
num_pos_from, lattice_points, num_lattice_points, reduced_basis,
trans_mat, initialize, symprec);
}
void py_thm_relative_grid_address(nb::ndarray<> py_relative_grid_address,
nb::ndarray<> py_reciprocal_lattice_py) {
long(*relative_grid_address)[4][3];
double(*reciprocal_lattice)[3];
relative_grid_address = (long(*)[4][3])py_relative_grid_address.data();
reciprocal_lattice = (double(*)[3])py_reciprocal_lattice_py.data();
phpy_get_relative_grid_address(relative_grid_address, reciprocal_lattice);
}
void py_thm_all_relative_grid_address(nb::ndarray<> py_relative_grid_address) {
long(*relative_grid_address)[24][4][3];
relative_grid_address = (long(*)[24][4][3])py_relative_grid_address.data();
phpy_get_all_relative_grid_address(relative_grid_address);
}
double py_thm_integration_weight(double omega,
nb::ndarray<> py_tetrahedra_omegas,
const char *function) {
double(*tetrahedra_omegas)[4];
double iw;
tetrahedra_omegas = (double(*)[4])py_tetrahedra_omegas.data();
iw = phpy_get_integration_weight(omega, tetrahedra_omegas, function[0]);
return iw;
}
void py_thm_integration_weight_at_omegas(nb::ndarray<> py_integration_weights,
nb::ndarray<> py_omegas,
nb::ndarray<> py_tetrahedra_omegas,
const char *function) {
double *omegas;
double *iw;
long num_omegas;
double(*tetrahedra_omegas)[4];
long i;
omegas = (double *)py_omegas.data();
iw = (double *)py_integration_weights.data();
num_omegas = (long)py_omegas.shape(0);
tetrahedra_omegas = (double(*)[4])py_tetrahedra_omegas.data();
#ifdef _OPENMP
#pragma omp parallel for
#endif
for (i = 0; i < num_omegas; i++) {
iw[i] = phpy_get_integration_weight(omegas[i], tetrahedra_omegas,
function[0]);
}
}
void py_get_tetrahedra_frequenies(nb::ndarray<> py_freq_tetras,
nb::ndarray<> py_grid_points,
nb::ndarray<> py_mesh,
nb::ndarray<> py_grid_address,
nb::ndarray<> py_gp_ir_index,
nb::ndarray<> py_relative_grid_address,
nb::ndarray<> py_frequencies) {
double *freq_tetras;
long *grid_points;
long *mesh;
long(*grid_address)[3];
long *gp_ir_index;
long(*relative_grid_address)[3];
double *frequencies;
long num_gp_in, num_band;
freq_tetras = (double *)py_freq_tetras.data();
grid_points = (long *)py_grid_points.data();
num_gp_in = py_grid_points.shape(0);
mesh = (long *)py_mesh.data();
grid_address = (long(*)[3])py_grid_address.data();
gp_ir_index = (long *)py_gp_ir_index.data();
relative_grid_address = (long(*)[3])py_relative_grid_address.data();
frequencies = (double *)py_frequencies.data();
num_band = py_frequencies.shape(1);
phpy_get_tetrahedra_frequenies(freq_tetras, mesh, grid_points, grid_address,
relative_grid_address, gp_ir_index,
frequencies, num_band, num_gp_in);
}
void py_tetrahedron_method_dos(nb::ndarray<> py_dos, nb::ndarray<> py_mesh,
nb::ndarray<> py_freq_points,
nb::ndarray<> py_frequencies,
nb::ndarray<> py_coef,
nb::ndarray<> py_grid_address,
nb::ndarray<> py_grid_mapping_table,
nb::ndarray<> py_relative_grid_address) {
double *dos;
long *mesh;
double *freq_points;
double *frequencies;
double *coef;
long(*grid_address)[3];
long num_gp, num_ir_gp, num_band, num_freq_points, num_coef;
long *grid_mapping_table;
long(*relative_grid_address)[4][3];
/* dos[num_ir_gp][num_band][num_freq_points][num_coef] */
dos = (double *)py_dos.data();
mesh = (long *)py_mesh.data();
freq_points = (double *)py_freq_points.data();
num_freq_points = (long)py_freq_points.shape(0);
frequencies = (double *)py_frequencies.data();
num_ir_gp = (long)py_frequencies.shape(0);
num_band = (long)py_frequencies.shape(1);
coef = (double *)py_coef.data();
num_coef = (long)py_coef.shape(1);
grid_address = (long(*)[3])py_grid_address.data();
num_gp = (long)py_grid_address.shape(0);
grid_mapping_table = (long *)py_grid_mapping_table.data();
relative_grid_address = (long(*)[4][3])py_relative_grid_address.data();
phpy_tetrahedron_method_dos(dos, mesh, grid_address, relative_grid_address,
grid_mapping_table, freq_points, frequencies,
coef, num_freq_points, num_ir_gp, num_band,
num_coef, num_gp);
}
NB_MODULE(_phonopy, m) {
m.def("transform_dynmat_to_fc", &py_transform_dynmat_to_fc);
m.def("perm_trans_symmetrize_fc", &py_perm_trans_symmetrize_fc);
m.def("perm_trans_symmetrize_compact_fc",
&py_perm_trans_symmetrize_compact_fc);
m.def("transpose_compact_fc", &py_transpose_compact_fc);
m.def("dynamical_matrices_with_dd_openmp_over_qpoints",
&py_get_dynamical_matrices_with_dd_openmp_over_qpoints);
m.def("dynamical_matrix", &py_get_dynamical_matrix);
m.def("nac_dynamical_matrix", &py_get_nac_dynamical_matrix);
m.def("recip_dipole_dipole", &py_get_recip_dipole_dipole);
m.def("recip_dipole_dipole_q0", &py_get_recip_dipole_dipole_q0);
m.def("derivative_dynmat", &py_get_derivative_dynmat);
m.def("thermal_properties", &py_get_thermal_properties);
m.def("distribute_fc2", &py_distribute_fc2);
m.def("compute_permutation", &py_compute_permutation);
m.def("gsv_set_smallest_vectors_sparse",
&py_gsv_set_smallest_vectors_sparse);
m.def("gsv_set_smallest_vectors_dense", &py_gsv_set_smallest_vectors_dense);
m.def("tetrahedra_relative_grid_address", &py_thm_relative_grid_address);
m.def("all_tetrahedra_relative_grid_address",
&py_thm_all_relative_grid_address);
m.def("tetrahedra_integration_weight", &py_thm_integration_weight);
m.def("tetrahedra_integration_weight_at_omegas",
&py_thm_integration_weight_at_omegas);
m.def("tetrahedra_frequencies", &py_get_tetrahedra_frequenies);
m.def("tetrahedron_method_dos", &py_tetrahedron_method_dos);
m.def("use_openmp", &phpy_use_openmp);
m.def("get_max_threads", &phpy_get_max_threads);
}

View File

@ -45,7 +45,7 @@ static void get_derivative_dynmat_at_q(
const double *q, const double *lattice, /* column vector */
const double (*svecs)[3], const long (*multi)[2], const double *mass,
const long *s2p_map, const long *p2s_map, const double nac_factor,
const double *born, const double *dielectric, const double *q_direction);
const double *born, const double *dielectric);
static void get_derivative_nac(double *ddnac, double *dnac,
const long num_patom, const double *lattice,
const double *mass, const double *q,
@ -66,27 +66,11 @@ void ddm_get_derivative_dynmat_at_q(
const double (*svecs)[3], const long (*multi)[2], const double *mass,
const long *s2p_map, const long *p2s_map, const double nac_factor,
const double *born, const double *dielectric, const double *q_direction,
const long use_openmp) {
long i, j, k, ij, adrs, adrsT, is_nac;
const long is_nac, const long use_openmp) {
long i, j, k, ij, adrs, adrsT;
double factor;
double *ddnac, *dnac;
if (born) {
is_nac = 1;
if (q_direction) {
if (fabs(q_direction[0]) < 1e-5 && fabs(q_direction[1]) < 1e-5 &&
fabs(q_direction[2]) < 1e-5) {
is_nac = 0;
}
} else {
if (fabs(q[0]) < 1e-5 && fabs(q[1]) < 1e-5 && fabs(q[2]) < 1e-5) {
is_nac = 0;
}
}
} else {
is_nac = 0;
}
if (is_nac) {
ddnac = (double *)malloc(sizeof(double) * num_patom * num_patom * 27);
dnac = (double *)malloc(sizeof(double) * num_patom * num_patom * 9);
@ -102,10 +86,10 @@ void ddm_get_derivative_dynmat_at_q(
for (ij = 0; ij < num_patom * num_patom; ij++) {
i = ij / num_patom;
j = ij % num_patom;
get_derivative_dynmat_at_q(
derivative_dynmat, i, j, ddnac, dnac, is_nac, num_patom,
num_satom, fc, q, lattice, svecs, multi, mass, s2p_map, p2s_map,
nac_factor, born, dielectric, q_direction);
get_derivative_dynmat_at_q(derivative_dynmat, i, j, ddnac, dnac,
is_nac, num_patom, num_satom, fc, q,
lattice, svecs, multi, mass, s2p_map,
p2s_map, nac_factor, born, dielectric);
}
} else {
for (i = 0; i < num_patom; i++) {
@ -113,7 +97,7 @@ void ddm_get_derivative_dynmat_at_q(
get_derivative_dynmat_at_q(
derivative_dynmat, i, j, ddnac, dnac, is_nac, num_patom,
num_satom, fc, q, lattice, svecs, multi, mass, s2p_map,
p2s_map, nac_factor, born, dielectric, q_direction);
p2s_map, nac_factor, born, dielectric);
}
}
}
@ -149,7 +133,7 @@ void get_derivative_dynmat_at_q(
const double *q, const double *lattice, /* column vector */
const double (*svecs)[3], const long (*multi)[2], const double *mass,
const long *s2p_map, const long *p2s_map, const double nac_factor,
const double *born, const double *dielectric, const double *q_direction) {
const double *born, const double *dielectric) {
long k, l, m, n, adrs, m_pair, i_pair, svecs_adrs;
double coef[3], real_coef[3], imag_coef[3];
double c, s, phase, mass_sqrt, fc_elem, real_phase, imag_phase;

View File

@ -42,6 +42,6 @@ void ddm_get_derivative_dynmat_at_q(
const double (*svecs)[3], const long (*multi)[2], const double *mass,
const long *s2p_map, const long *p2s_map, const double nac_factor,
const double *born, const double *dielectric, const double *q_direction,
const long use_openmp);
const long is_nac, const long use_openmp);
#endif

View File

@ -161,11 +161,11 @@ void phpy_get_derivative_dynmat_at_q(
const double (*svecs)[3], const long (*multi)[2], const double *mass,
const long *s2p_map, const long *p2s_map, const double nac_factor,
const double *born, const double *dielectric, const double *q_direction,
const long use_openmp) {
const long is_nac, const long use_openmp) {
ddm_get_derivative_dynmat_at_q(derivative_dynmat, num_patom, num_satom, fc,
q, lattice, svecs, multi, mass, s2p_map,
p2s_map, nac_factor, born, dielectric,
q_direction, use_openmp);
q_direction, is_nac, use_openmp);
}
void phpy_get_relative_grid_address(long relative_grid_address[24][4][3],

View File

@ -37,6 +37,10 @@
#ifndef __phonopy_H__
#define __phonopy_H__
#ifdef __cplusplus
extern "C" {
#endif
void phpy_transform_dynmat_to_fc(double *fc, const double (*dm)[2],
const double (*comm_points)[3],
const double (*svecs)[3],
@ -94,7 +98,7 @@ void phpy_get_derivative_dynmat_at_q(
const double (*svecs)[3], const long (*multi)[2], const double *mass,
const long *s2p_map, const long *p2s_map, const double nac_factor,
const double *born, const double *dielectric, const double *q_direction,
const long use_openmp);
const long is_nac, const long use_openmp);
void phpy_get_neighboring_grid_points(
size_t neighboring_grid_points[], const size_t grid_point,
const int relative_grid_address[][3], const int num_relative_grid_address,
@ -162,4 +166,9 @@ void phpy_set_index_permutation_symmetry_compact_fc(
const int is_transpose);
long phpy_use_openmp(void);
long phpy_get_max_threads(void);
#ifdef __cplusplus
}
#endif
#endif

View File

@ -141,18 +141,23 @@ class DerivativeOfDynamicalMatrix:
dtype=("c%d" % (np.dtype("double").itemsize * 2)),
)
if isinstance(self._dynmat, DynamicalMatrixNAC):
is_nac = True
born = self._dynmat.born
dielectric = self._dynmat.dielectric_constant
nac_factor = self._dynmat.nac_factor
if q_direction is None:
q_dir = None
q_dir = np.zerod(3)
is_q_zero = True
else:
q_dir = np.array(q_direction, dtype="double", order="C")
is_q_zero = False
else:
born = None
dielectric = None
nac_factor = 0
q_dir = None
born = np.zeros(0) # dummy value
dielectric = np.zeros(0) # dummy value
nac_factor = 0 # dummy value
q_dir = np.zerod(3) # dummy value
is_nac = False
is_q_zero = True
if fc.shape[0] == fc.shape[1]: # full fc
phonoc.derivative_dynmat(
@ -169,6 +174,8 @@ class DerivativeOfDynamicalMatrix:
born,
dielectric,
q_dir,
is_nac * 1,
is_q_zero * 1,
self._dynmat.use_openmp * 1,
)
else:
@ -186,6 +193,8 @@ class DerivativeOfDynamicalMatrix:
born,
dielectric,
q_dir,
is_nac * 1,
is_q_zero * 1,
self._dynmat.use_openmp * 1,
)

View File

@ -802,7 +802,9 @@ class DynamicalMatrixGL(DynamicalMatrixNAC):
return C_dd
def _get_c_recip_dipole_dipole(self, q_cart, q_dir_cart):
def _get_c_recip_dipole_dipole(
self, q_cart: np.ndarray, q_dir_cart: Optional[np.ndarray]
) -> np.ndarray:
"""Reciprocal part of Eq.(71) on the right hand side.
This is subtracted from supercell force constants to create
@ -828,15 +830,23 @@ class DynamicalMatrixGL(DynamicalMatrixNAC):
else:
dd_q0 = self._dd_q0
if q_dir_cart is None:
is_q_zero = True
_q_dir_cart = np.zeros(3, dtype="double")
else:
is_q_zero = False
_q_dir_cart = q_dir_cart
phonoc.recip_dipole_dipole(
dd.view(dtype="double"),
dd_q0.view(dtype="double"),
self._G_list,
q_cart,
q_dir_cart,
_q_dir_cart,
self._born,
self._dielectric,
np.array(pos, dtype="double", order="C"),
is_q_zero * 1,
self._unit_conversion * 4.0 * np.pi / volume,
self._Lambda,
self.Q_DIRECTION_TOLERANCE,
@ -1224,10 +1234,9 @@ def run_dynamical_matrix_solver_c(
) = gonze_nac_dataset # Convergence parameter
fc = gonze_fc
else:
positions = None
dd_q0 = None
G_list = None
Lambda = 0
dd_q0 = np.zeros(2) # dummy value
G_list = np.zeros(3) # dummy value
Lambda = 0 # dummy value
fc = dm.force_constants
if isinstance(dm, DynamicalMatrixWang):
use_Wang_NAC = True

View File

@ -293,7 +293,7 @@ def distribute_force_constants(
atom_list_done,
lattice, # column vectors
rotations, # scaled (fractional)
permutations,
permutations: np.ndarray,
atom_list=None,
fc_indices_of_atom_list=None,
):
@ -322,6 +322,16 @@ def distribute_force_constants(
_fc_indices_of_atom_list = np.arange(len(targets), dtype="intc")
else:
_fc_indices_of_atom_list = np.array(fc_indices_of_atom_list, dtype="intc")
if map_atoms.ndim != 1 or map_atoms.shape[0] != permutations.shape[1]:
raise ValueError("wrong shape for map_atoms")
if map_syms.ndim != 1 or map_syms.shape[0] != permutations.shape[1]:
raise ValueError("wrong shape for map_syms")
if rots_cartesian.shape[0] != permutations.shape[0]:
raise ValueError("permutations and rotations are different length")
import phonopy._phonopy as phonoc
phonoc.distribute_fc2(

View File

@ -36,7 +36,9 @@
from phonopy.cui.phonopy_script import main
if __name__ == "__main__":
def run():
"""Run phonopy script."""
argparse_control = {
"fc_symmetry": False,
"is_nac": False,

View File

@ -36,7 +36,9 @@
from phonopy.cui.phonopy_script import main
if __name__ == "__main__":
def run():
"""Run phonopy-load script."""
argparse_control = {
"fc_symmetry": True,
"is_nac": True,

View File

@ -38,14 +38,6 @@ import warnings
import numpy as np
try:
import phonopy._phonopy as phonoc
except ImportError:
import sys
print("Phonopy C-extension has to be built properly.")
sys.exit(1)
def get_tetrahedra_relative_grid_address(microzone_lattice, lang="C"):
"""Return relative (differences of) grid addresses from the central.
@ -57,6 +49,14 @@ def get_tetrahedra_relative_grid_address(microzone_lattice, lang="C"):
microzone_lattice = np.linalg.inv(cell.get_cell()) / mesh
"""
try:
import phonopy._phonopy as phonoc
except ImportError:
import sys
print("Phonopy C-extension has to be built properly.")
sys.exit(1)
relative_grid_address = np.zeros((24, 4, 3), dtype="int_", order="C")
if lang == "C":
phonoc.tetrahedra_relative_grid_address(
@ -76,6 +76,14 @@ def get_all_tetrahedra_relative_grid_address(lang="C"):
This exists only for the test.
"""
try:
import phonopy._phonopy as phonoc
except ImportError:
import sys
print("Phonopy C-extension has to be built properly.")
sys.exit(1)
relative_grid_address = np.zeros((4, 24, 4, 3), dtype="int_", order="C")
if lang == "C":
phonoc.all_tetrahedra_relative_grid_address(relative_grid_address)
@ -103,6 +111,14 @@ def get_tetrahedra_integration_weight(omegas, tetrahedra_omegas, function="I"):
'J' is for intetration and 'I' is for its derivative.
"""
try:
import phonopy._phonopy as phonoc
except ImportError:
import sys
print("Phonopy C-extension has to be built properly.")
sys.exit(1)
if isinstance(omegas, float):
return phonoc.tetrahedra_integration_weight(
omegas, np.array(tetrahedra_omegas, dtype="double", order="C"), function

View File

@ -1,5 +1,53 @@
[build-system]
requires = ["setuptools", "wheel", "numpy"]
requires = ["scikit-build-core", "nanobind", "numpy"]
build-backend = "scikit_build_core.build"
[project]
name = "phonopy"
dynamic = ["version"]
readme = { file = "README.md", content-type = "text/markdown" }
description = "This is the phonopy module."
authors = [{ name = "Atsushi Togo", email = "atztogo@gmail.com" }]
requires-python = ">=3.8"
dependencies = [
"numpy>=1.17.0",
"PyYAML>=5.3",
"matplotlib>=2.2.2",
"h5py>=3.0",
"spglib>=2.3",
]
license = { file = "LICENSE" }
[project.optional-dependencies]
cp2k = ["cp2k-input-tools"]
plot = ["seekpath"]
[project.scripts]
phonopy = "phonopy.scripts.phonopy:run"
phonopy-load = "phonopy.scripts.phonopy_load:run"
# phonopy-qha = "phonopy.scripts.phonopy-qha"
# phonopy-bandplot = "phonopy.scripts.phonopy-bandplot"
# phonopy-vasp-born = "phonopy.scripts.phonopy-vasp-born"
# phonopy-vasp-efe = "phonopy.scripts.phonopy-vasp-efe"
# phonopy-crystal-born = "phonopy.scripts.phonopy-crystal-born"
# phonopy-calc-convert = "phonopy.scripts.phonopy-calc-convert"
# phonopy-propplot = "phonopy.scripts.phonopy-propplot"
# phonopy-tdplot = "phonopy.scripts.phonopy-tdplot"
# phonopy-gruneisen = "phonopy.scripts.phonopy-gruneisen"
# phonopy-gruneisenplot = "phonopy.scripts.phonopy-gruneisenplot"
# phonopy-pdosplot = "phonopy.scripts.phonopy-pdosplot"
[tool.sciki-build]
cmake.verbose = true
logging.level = "INFO"
[tool.scikit-build]
metadata.version.provider = "scikit_build_core.metadata.setuptools_scm"
sdist.include = ["phonopy/_version.py"]
[tool.setuptools_scm]
write_to = "phonopy/_version.py"
[tool.ruff]
line-length = 88
@ -10,10 +58,7 @@ lint.select = [
"E", # pycodestyle-errors
"D", # pydocstyle
]
lint.extend-ignore = [
"D417",
"D100",
]
lint.extend-ignore = ["D417", "D100"]
exclude = [
"test/phonon/test_irreps.py",
"test/qha/test_electron.py",