scnc
/
phono3py
mirror of https://github.com/phonopy/phono3py.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

size_t -> long for triplets related integer values

This commit is contained in:
Atsushi Togo 2021-02-22 12:46:31 +09:00
parent 545af58938
commit 38ef33df25
40 changed files with 1296 additions and 1644 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
CMakeLists.txt
View File

@ -58,6 +58,7 @@ set(SOURCES_PHONO3PY
${PROJECT_SOURCE_DIR}/c/real_self_energy.c
${PROJECT_SOURCE_DIR}/c/real_to_reciprocal.c
${PROJECT_SOURCE_DIR}/c/reciprocal_to_normal.c
${PROJECT_SOURCE_DIR}/c/rgrid.c
${PROJECT_SOURCE_DIR}/c/tetrahedron_method.c
${PROJECT_SOURCE_DIR}/c/triplet.c
${PROJECT_SOURCE_DIR}/c/triplet_iw.c

318
c/_phono3py.c
View File

@ -89,13 +89,13 @@ static PyObject * py_get_default_colmat_solver(PyObject *self, PyObject *args);
static void pinv_from_eigensolution(double *data,
const double *eigvals,
const size_t size,
const long size,
const double cutoff,
const int pinv_method);
static void show_colmat_info(const PyArrayObject *collision_matrix_py,
const size_t i_sigma,
const size_t i_temp,
const size_t adrs_shift);
const long i_sigma,
const long i_temp,
const long adrs_shift);
static Iarray* convert_to_iarray(const PyArrayObject* npyary);
static Darray* convert_to_darray(const PyArrayObject* npyary);
@ -120,7 +120,7 @@ error_out(PyObject *m) {
static PyMethodDef _phono3py_methods[] = {
{"error_out", (PyCFunction)error_out, METH_NOARGS, NULL},
{"interaction",
{"interaction",
(PyCFunction)py_get_interaction,
METH_VARARGS,
"Interaction of triplets"},
@ -311,8 +311,8 @@ static PyObject * py_get_interaction(PyObject *self, PyObject *args)
Darray *fc3_normal_squared;
Darray *freqs;
lapack_complex_double *eigvecs;
size_t (*triplets)[3];
npy_intp num_triplets;
long (*triplets)[3];
long num_triplets;
char* g_zero;
int *grid_address;
int *mesh;
@ -353,8 +353,8 @@ static PyObject * py_get_interaction(PyObject *self, PyObject *args)
/* 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 = (size_t(*)[3])PyArray_DATA(py_triplets);
num_triplets = PyArray_DIMS(py_triplets)[0];
triplets = (long(*)[3])PyArray_DATA(py_triplets);
num_triplets = (long)PyArray_DIMS(py_triplets)[0];
g_zero = (char*)PyArray_DATA(py_g_zero);
grid_address = (int*)PyArray_DATA(py_grid_address);
mesh = (int*)PyArray_DATA(py_mesh);
@ -429,11 +429,11 @@ static PyObject * py_get_pp_collision(PyObject *self, PyObject *args)
int (*relative_grid_address)[4][3];
double *frequencies;
lapack_complex_double *eigenvectors;
size_t (*triplets)[3];
npy_intp num_triplets;
int *triplet_weights;
long (*triplets)[3];
long num_triplets;
long *triplet_weights;
int *grid_address;
size_t *bz_map;
long *bz_map;
int *mesh;
double *fc3;
double *svecs;
@ -475,11 +475,11 @@ static PyObject * py_get_pp_collision(PyObject *self, PyObject *args)
relative_grid_address = (int(*)[4][3])PyArray_DATA(py_relative_grid_address);
frequencies = (double*)PyArray_DATA(py_frequencies);
eigenvectors = (lapack_complex_double*)PyArray_DATA(py_eigenvectors);
triplets = (size_t(*)[3])PyArray_DATA(py_triplets);
num_triplets = PyArray_DIMS(py_triplets)[0];
triplet_weights = (int*)PyArray_DATA(py_triplet_weights);
triplets = (long(*)[3])PyArray_DATA(py_triplets);
num_triplets = (long)PyArray_DIMS(py_triplets)[0];
triplet_weights = (long*)PyArray_DATA(py_triplet_weights);
grid_address = (int*)PyArray_DATA(py_grid_address);
bz_map = (size_t*)PyArray_DATA(py_bz_map);
bz_map = (long*)PyArray_DATA(py_bz_map);
mesh = (int*)PyArray_DATA(py_mesh);
fc3 = (double*)PyArray_DATA(py_fc3);
if (PyArray_DIMS(py_fc3)[0] == PyArray_DIMS(py_fc3)[1]) {
@ -556,9 +556,9 @@ static PyObject * py_get_pp_collision_with_sigma(PyObject *self, PyObject *args)
double *gamma;
double *frequencies;
lapack_complex_double *eigenvectors;
size_t (*triplets)[3];
npy_intp num_triplets;
int *triplet_weights;
long (*triplets)[3];
long num_triplets;
long *triplet_weights;
int *grid_address;
int *mesh;
double *fc3;
@ -600,9 +600,9 @@ static PyObject * py_get_pp_collision_with_sigma(PyObject *self, PyObject *args)
gamma = (double*)PyArray_DATA(py_gamma);
frequencies = (double*)PyArray_DATA(py_frequencies);
eigenvectors = (lapack_complex_double*)PyArray_DATA(py_eigenvectors);
triplets = (size_t(*)[3])PyArray_DATA(py_triplets);
num_triplets = PyArray_DIMS(py_triplets)[0];
triplet_weights = (int*)PyArray_DATA(py_triplet_weights);
triplets = (long(*)[3])PyArray_DATA(py_triplets);
num_triplets = (long)PyArray_DIMS(py_triplets)[0];
triplet_weights = (long*)PyArray_DATA(py_triplet_weights);
grid_address = (int*)PyArray_DATA(py_grid_address);
mesh = (int*)PyArray_DATA(py_mesh);
fc3 = (double*)PyArray_DATA(py_fc3);
@ -671,9 +671,9 @@ static PyObject * py_get_imag_self_energy_with_g(PyObject *self, PyObject *args)
double *g;
char* g_zero;
double *frequencies;
size_t (*triplets)[3];
int *triplet_weights;
int num_frequency_points;
long (*triplets)[3];
long *triplet_weights;
long num_frequency_points;
if (!PyArg_ParseTuple(args, "OOOOOdOOdi",
&py_gamma,
@ -694,9 +694,9 @@ static PyObject * py_get_imag_self_energy_with_g(PyObject *self, PyObject *args)
g = (double*)PyArray_DATA(py_g);
g_zero = (char*)PyArray_DATA(py_g_zero);
frequencies = (double*)PyArray_DATA(py_frequencies);
triplets = (size_t(*)[3])PyArray_DATA(py_triplets);
triplet_weights = (int*)PyArray_DATA(py_triplet_weights);
num_frequency_points = PyArray_DIMS(py_g)[2];
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,
@ -738,8 +738,8 @@ py_get_detailed_imag_self_energy_with_g(PyObject *self, PyObject *args)
double *g;
char* g_zero;
double *frequencies;
size_t (*triplets)[3];
int *triplet_weights;
long (*triplets)[3];
long *triplet_weights;
int *grid_address;
if (!PyArg_ParseTuple(args, "OOOOOOOOdOOd",
@ -765,8 +765,8 @@ py_get_detailed_imag_self_energy_with_g(PyObject *self, PyObject *args)
g = (double*)PyArray_DATA(py_g);
g_zero = (char*)PyArray_DATA(py_g_zero);
frequencies = (double*)PyArray_DATA(py_frequencies);
triplets = (size_t(*)[3])PyArray_DATA(py_triplets);
triplet_weights = (int*)PyArray_DATA(py_triplet_weights);
triplets = (long(*)[3])PyArray_DATA(py_triplets);
triplet_weights = (long*)PyArray_DATA(py_triplet_weights);
grid_address = (int*)PyArray_DATA(py_grid_address);
ph3py_get_detailed_imag_self_energy_at_bands_with_g(gamma_detail,
@ -803,8 +803,8 @@ static PyObject * py_get_real_self_energy_at_bands(PyObject *self,
double *shift;
double *frequencies;
int *band_indices;
size_t (*triplets)[3];
int *triplet_weights;
long (*triplets)[3];
long *triplet_weights;
if (!PyArg_ParseTuple(args, "OOOOOOdddd",
&py_shift,
@ -825,8 +825,8 @@ static PyObject * py_get_real_self_energy_at_bands(PyObject *self,
shift = (double*)PyArray_DATA(py_shift);
frequencies = (double*)PyArray_DATA(py_frequencies);
band_indices = (int*)PyArray_DATA(py_band_indices);
triplets = (size_t(*)[3])PyArray_DATA(py_triplets);
triplet_weights = (int*)PyArray_DATA(py_triplet_weights);
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,
@ -861,8 +861,8 @@ static PyObject * py_get_real_self_energy_at_frequency_point(PyObject *self,
double *shift;
double *frequencies;
int *band_indices;
size_t (*triplets)[3];
int *triplet_weights;
long (*triplets)[3];
long *triplet_weights;
if (!PyArg_ParseTuple(args, "OdOOOOOdddd",
&py_shift,
@ -884,8 +884,8 @@ static PyObject * py_get_real_self_energy_at_frequency_point(PyObject *self,
shift = (double*)PyArray_DATA(py_shift);
frequencies = (double*)PyArray_DATA(py_frequencies);
band_indices = (int*)PyArray_DATA(py_band_indices);
triplets = (size_t(*)[3])PyArray_DATA(py_triplets);
triplet_weights = (int*)PyArray_DATA(py_triplet_weights);
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,
@ -922,11 +922,11 @@ static PyObject * py_get_collision_matrix(PyObject *self, PyObject *args)
double *collision_matrix;
double *g;
double *frequencies;
size_t (*triplets)[3];
size_t *triplets_map;
size_t *map_q;
size_t *rotated_grid_points;
npy_intp num_gp, num_ir_gp, num_rot;
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",
@ -949,13 +949,13 @@ static PyObject * py_get_collision_matrix(PyObject *self, PyObject *args)
collision_matrix = (double*)PyArray_DATA(py_collision_matrix);
g = (double*)PyArray_DATA(py_g);
frequencies = (double*)PyArray_DATA(py_frequencies);
triplets = (size_t(*)[3])PyArray_DATA(py_triplets);
triplets_map = (size_t*)PyArray_DATA(py_triplets_map);
num_gp = PyArray_DIMS(py_triplets_map)[0];
map_q = (size_t*)PyArray_DATA(py_map_q);
rotated_grid_points = (size_t*)PyArray_DATA(py_rotated_grid_points);
num_ir_gp = PyArray_DIMS(py_rotated_grid_points)[0];
num_rot = PyArray_DIMS(py_rotated_grid_points)[1];
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]);
@ -998,10 +998,10 @@ static PyObject * py_get_reducible_collision_matrix(PyObject *self, PyObject *ar
double *collision_matrix;
double *g;
double *frequencies;
size_t (*triplets)[3];
size_t *triplets_map;
npy_intp num_gp;
size_t *map_q;
long (*triplets)[3];
long *triplets_map;
long num_gp;
long *map_q;
if (!PyArg_ParseTuple(args, "OOOOOOOddd",
&py_collision_matrix,
@ -1021,10 +1021,10 @@ static PyObject * py_get_reducible_collision_matrix(PyObject *self, PyObject *ar
collision_matrix = (double*)PyArray_DATA(py_collision_matrix);
g = (double*)PyArray_DATA(py_g);
frequencies = (double*)PyArray_DATA(py_frequencies);
triplets = (size_t(*)[3])PyArray_DATA(py_triplets);
triplets_map = (size_t*)PyArray_DATA(py_triplets_map);
num_gp = PyArray_DIMS(py_triplets_map)[0];
map_q = (size_t*)PyArray_DATA(py_map_q);
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,
@ -1084,8 +1084,8 @@ static PyObject * py_expand_collision_matrix(PyObject *self, PyObject *args)
PyArrayObject *py_rot_grid_points;
double *collision_matrix;
size_t *rot_grid_points;
size_t *ir_grid_points;
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",
@ -1096,8 +1096,8 @@ static PyObject * py_expand_collision_matrix(PyObject *self, PyObject *args)
}
collision_matrix = (double*)PyArray_DATA(py_collision_matrix);
rot_grid_points = (size_t*)PyArray_DATA(py_rot_grid_points);
ir_grid_points = (size_t*)PyArray_DATA(py_ir_grid_points);
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];
@ -1135,7 +1135,7 @@ static PyObject * py_get_isotope_strength(PyObject *self, PyObject *args)
lapack_complex_double *eigenvectors;
int *band_indices;
double *mass_variances;
npy_intp num_band, num_band0;
long num_band, num_band0;
if (!PyArg_ParseTuple(args, "OlOOOOldd",
&py_gamma,
@ -1156,8 +1156,8 @@ static PyObject * py_get_isotope_strength(PyObject *self, PyObject *args)
eigenvectors = (lapack_complex_double*)PyArray_DATA(py_eigenvectors);
band_indices = (int*)PyArray_DATA(py_band_indices);
mass_variances = (double*)PyArray_DATA(py_mass_variances);
num_band = PyArray_DIMS(py_frequencies)[1];
num_band0 = PyArray_DIMS(py_band_indices)[0];
num_band = (long)PyArray_DIMS(py_frequencies)[1];
num_band0 = (long)PyArray_DIMS(py_band_indices)[0];
ph3py_get_isotope_scattering_strength(gamma,
grid_point,
@ -1189,12 +1189,12 @@ static PyObject * py_get_thm_isotope_strength(PyObject *self, PyObject *args)
double *gamma;
double *frequencies;
size_t *ir_grid_points;
int *weights;
long *ir_grid_points;
long *weights;
lapack_complex_double *eigenvectors;
int *band_indices;
double *mass_variances;
npy_intp num_band, num_band0, num_ir_grid_points;
long num_band, num_band0, num_ir_grid_points;
double *integration_weights;
if (!PyArg_ParseTuple(args, "OlOOOOOOOd",
@ -1214,15 +1214,15 @@ static PyObject * py_get_thm_isotope_strength(PyObject *self, PyObject *args)
gamma = (double*)PyArray_DATA(py_gamma);
frequencies = (double*)PyArray_DATA(py_frequencies);
ir_grid_points = (size_t*)PyArray_DATA(py_ir_grid_points);
weights = (int*)PyArray_DATA(py_weights);
ir_grid_points = (long*)PyArray_DATA(py_ir_grid_points);
weights = (long*)PyArray_DATA(py_weights);
eigenvectors = (lapack_complex_double*)PyArray_DATA(py_eigenvectors);
band_indices = (int*)PyArray_DATA(py_band_indices);
mass_variances = (double*)PyArray_DATA(py_mass_variances);
num_band = PyArray_DIMS(py_frequencies)[1];
num_band0 = PyArray_DIMS(py_band_indices)[0];
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 = PyArray_DIMS(py_ir_grid_points)[0];
num_ir_grid_points = (long)PyArray_DIMS(py_ir_grid_points)[0];
ph3py_get_thm_isotope_scattering_strength(gamma,
grid_point,
@ -1252,7 +1252,7 @@ static PyObject * py_distribute_fc3(PyObject *self, PyObject *args)
double *fc3;
double *rot_cart_inv;
int *atom_mapping;
npy_intp num_atom;
long num_atom;
if (!PyArg_ParseTuple(args, "OiiOO",
&force_constants_third,
@ -1266,7 +1266,7 @@ static PyObject * py_distribute_fc3(PyObject *self, PyObject *args)
fc3 = (double*)PyArray_DATA(force_constants_third);
rot_cart_inv = (double*)PyArray_DATA(rotation_cart_inv);
atom_mapping = (int*)PyArray_DATA(atom_mapping_py);
num_atom = PyArray_DIMS(atom_mapping_py)[0];
num_atom = (long)PyArray_DIMS(atom_mapping_py)[0];
ph3py_distribute_fc3(fc3,
target,
@ -1291,7 +1291,7 @@ static PyObject * py_rotate_delta_fc2s(PyObject *self, PyObject *args)
double *inv_U;
double (*site_sym_cart)[3][3];
int *rot_map_syms;
npy_intp num_atom, num_disp, num_site_sym;
long num_atom, num_disp, num_site_sym;
if (!PyArg_ParseTuple(args, "OOOOO",
&py_fc3,
@ -1313,9 +1313,9 @@ static PyObject * py_rotate_delta_fc2s(PyObject *self, PyObject *args)
/* (n_sym, natom) */
rot_map_syms = (int*)PyArray_DATA(py_rot_map_syms);
num_atom = PyArray_DIMS(py_fc3)[0];
num_disp = PyArray_DIMS(py_delta_fc2s)[0];
num_site_sym = PyArray_DIMS(py_site_sym_cart)[0];
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,
@ -1335,14 +1335,14 @@ py_set_permutation_symmetry_fc3(PyObject *self, PyObject *args)
PyArrayObject *py_fc3;
double *fc3;
npy_intp num_atom;
long num_atom;
if (!PyArg_ParseTuple(args, "O", &py_fc3)) {
return NULL;
}
fc3 = (double*)PyArray_DATA(py_fc3);
num_atom = PyArray_DIMS(py_fc3)[0];
num_atom = (long)PyArray_DIMS(py_fc3)[0];
ph3py_set_permutation_symmetry_fc3(fc3, num_atom);
@ -1363,7 +1363,7 @@ py_set_permutation_symmetry_compact_fc3(PyObject *self, PyObject *args)
int *p2s;
int *nsym_list;
int *perms;
npy_intp n_patom, n_satom;
long n_patom, n_satom;
if (!PyArg_ParseTuple(args, "OOOOO",
&py_fc3,
@ -1379,8 +1379,8 @@ py_set_permutation_symmetry_compact_fc3(PyObject *self, PyObject *args)
s2pp = (int*)PyArray_DATA(py_s2pp_map);
p2s = (int*)PyArray_DATA(py_p2s_map);
nsym_list = (int*)PyArray_DATA(py_nsym_list);
n_patom = PyArray_DIMS(py_fc3)[0];
n_satom = PyArray_DIMS(py_fc3)[1];
n_patom = (long)PyArray_DIMS(py_fc3)[0];
n_satom = (long)PyArray_DIMS(py_fc3)[1];
ph3py_set_permutation_symmetry_compact_fc3(fc3,
p2s,
@ -1407,7 +1407,7 @@ static PyObject * py_transpose_compact_fc3(PyObject *self, PyObject *args)
int *p2s;
int *nsym_list;
int *perms;
npy_intp n_patom, n_satom;
long n_patom, n_satom;
if (!PyArg_ParseTuple(args, "OOOOOi",
&py_fc3,
@ -1424,8 +1424,8 @@ static PyObject * py_transpose_compact_fc3(PyObject *self, PyObject *args)
s2pp = (int*)PyArray_DATA(py_s2pp_map);
p2s = (int*)PyArray_DATA(py_p2s_map);
nsym_list = (int*)PyArray_DATA(py_nsym_list);
n_patom = PyArray_DIMS(py_fc3)[0];
n_satom = PyArray_DIMS(py_fc3)[1];
n_patom = (long)PyArray_DIMS(py_fc3)[0];
n_satom = (long)PyArray_DIMS(py_fc3)[1];
ph3py_transpose_compact_fc3(fc3,
p2s,
@ -1448,13 +1448,13 @@ static PyObject * py_get_neighboring_gird_points(PyObject *self, PyObject *args)
PyArrayObject *py_bz_grid_address;
PyArrayObject *py_bz_map;
size_t *relative_grid_points;
size_t *grid_points;
npy_intp num_grid_points, num_relative_grid_address;
long *relative_grid_points;
long *grid_points;
long num_grid_points, num_relative_grid_address;
int (*relative_grid_address)[3];
int *mesh;
int (*bz_grid_address)[3];
size_t *bz_map;
long *bz_map;
if (!PyArg_ParseTuple(args, "OOOOOO",
&py_relative_grid_points,
@ -1466,14 +1466,14 @@ static PyObject * py_get_neighboring_gird_points(PyObject *self, PyObject *args)
return NULL;
}
relative_grid_points = (size_t*)PyArray_DATA(py_relative_grid_points);
grid_points = (size_t*)PyArray_DATA(py_grid_points);
num_grid_points = PyArray_DIMS(py_grid_points)[0];
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 = (int(*)[3])PyArray_DATA(py_relative_grid_address);
num_relative_grid_address = PyArray_DIMS(py_relative_grid_address)[0];
num_relative_grid_address = (long)PyArray_DIMS(py_relative_grid_address)[0];
mesh = (int*)PyArray_DATA(py_mesh);
bz_grid_address = (int(*)[3])PyArray_DATA(py_bz_grid_address);
bz_map = (size_t*)PyArray_DATA(py_bz_map);
bz_map = (long*)PyArray_DATA(py_bz_map);
ph3py_get_neighboring_gird_points(relative_grid_points,
grid_points,
@ -1500,12 +1500,12 @@ static PyObject * py_set_integration_weights(PyObject *self, PyObject *args)
double *iw;
double *frequency_points;
npy_intp num_band0, num_band, num_gp;
long num_band0, num_band, num_gp;
int (*relative_grid_address)[4][3];
int *mesh;
size_t *grid_points;
long *grid_points;
int (*bz_grid_address)[3];
size_t *bz_map;
long *bz_map;
double *frequencies;
if (!PyArg_ParseTuple(args, "OOOOOOOO",
@ -1522,15 +1522,15 @@ static PyObject * py_set_integration_weights(PyObject *self, PyObject *args)
iw = (double*)PyArray_DATA(py_iw);
frequency_points = (double*)PyArray_DATA(py_frequency_points);
num_band0 = PyArray_DIMS(py_frequency_points)[0];
num_band0 = (long)PyArray_DIMS(py_frequency_points)[0];
relative_grid_address = (int(*)[4][3])PyArray_DATA(py_relative_grid_address);
mesh = (int*)PyArray_DATA(py_mesh);
grid_points = (size_t*)PyArray_DATA(py_grid_points);
num_gp = PyArray_DIMS(py_grid_points)[0];
grid_points = (long*)PyArray_DATA(py_grid_points);
num_gp = (long)PyArray_DIMS(py_grid_points)[0];
bz_grid_address = (int(*)[3])PyArray_DATA(py_bz_grid_address);
bz_map = (size_t*)PyArray_DATA(py_bz_map);
bz_map = (long*)PyArray_DATA(py_bz_map);
frequencies = (double*)PyArray_DATA(py_frequencies);
num_band = PyArray_DIMS(py_frequencies)[1];
num_band = (long)PyArray_DIMS(py_frequencies)[1];
ph3py_set_integration_weights(iw,
frequency_points,
@ -1560,12 +1560,12 @@ py_tpl_get_triplets_reciprocal_mesh_at_q(PyObject *self, PyObject *args)
int swappable;
int (*grid_address)[3];
size_t *map_triplets;
size_t *map_q;
long *map_triplets;
long *map_q;
int *mesh;
int (*rot)[3][3];
npy_intp num_rot;
size_t num_ir;
long num_rot;
long num_ir;
if (!PyArg_ParseTuple(args, "OOOlOiOi",
&py_map_triplets,
@ -1580,11 +1580,11 @@ py_tpl_get_triplets_reciprocal_mesh_at_q(PyObject *self, PyObject *args)
}
grid_address = (int(*)[3])PyArray_DATA(py_grid_address);
map_triplets = (size_t*)PyArray_DATA(py_map_triplets);
map_q = (size_t*)PyArray_DATA(py_map_q);
map_triplets = (long*)PyArray_DATA(py_map_triplets);
map_q = (long*)PyArray_DATA(py_map_q);
mesh = (int*)PyArray_DATA(py_mesh);
rot = (int(*)[3][3])PyArray_DATA(py_rotations);
num_rot = PyArray_DIMS(py_rotations)[0];
num_rot = (long)PyArray_DIMS(py_rotations)[0];
num_ir = ph3py_get_triplets_reciprocal_mesh_at_q(map_triplets,
map_q,
grid_address,
@ -1595,7 +1595,7 @@ py_tpl_get_triplets_reciprocal_mesh_at_q(PyObject *self, PyObject *args)
rot,
swappable);
return PyLong_FromSize_t(num_ir);
return PyLong_FromLong(num_ir);
}
static PyObject * py_tpl_get_BZ_triplets_at_q(PyObject *self, PyObject *args)
@ -1607,13 +1607,13 @@ static PyObject * py_tpl_get_BZ_triplets_at_q(PyObject *self, PyObject *args)
PyArrayObject *py_mesh;
long grid_point;
size_t (*triplets)[3];
long (*triplets)[3];
int (*bz_grid_address)[3];
size_t *bz_map;
size_t *map_triplets;
npy_intp num_map_triplets;
long *bz_map;
long *map_triplets;
long num_map_triplets;
int *mesh;
size_t num_ir;
long num_ir;
if (!PyArg_ParseTuple(args, "OlOOOO",
&py_triplets,
@ -1625,11 +1625,11 @@ static PyObject * py_tpl_get_BZ_triplets_at_q(PyObject *self, PyObject *args)
return NULL;
}
triplets = (size_t(*)[3])PyArray_DATA(py_triplets);
triplets = (long(*)[3])PyArray_DATA(py_triplets);
bz_grid_address = (int(*)[3])PyArray_DATA(py_bz_grid_address);
bz_map = (size_t*)PyArray_DATA(py_bz_map);
map_triplets = (size_t*)PyArray_DATA(py_map_triplets);
num_map_triplets = PyArray_DIMS(py_map_triplets)[0];
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];
mesh = (int*)PyArray_DATA(py_mesh);
num_ir = ph3py_get_BZ_triplets_at_q(triplets,
@ -1640,7 +1640,7 @@ static PyObject * py_tpl_get_BZ_triplets_at_q(PyObject *self, PyObject *args)
num_map_triplets,
mesh);
return PyLong_FromSize_t(num_ir);
return PyLong_FromLong(num_ir);
}
static PyObject *
@ -1663,11 +1663,11 @@ py_set_triplets_integration_weights(PyObject *self, PyObject *args)
double *frequency_points;
int (*relative_grid_address)[4][3];
int *mesh;
size_t (*triplets)[3];
long (*triplets)[3];
int (*bz_grid_address)[3];
size_t *bz_map;
long *bz_map;
double *frequencies1, *frequencies2;
npy_intp num_band0, num_band1, num_band2, num_triplets;
long num_band0, num_band1, num_band2, num_triplets;
if (!PyArg_ParseTuple(args, "OOOOOOOOOOi",
&py_iw,
@ -1687,17 +1687,17 @@ py_set_triplets_integration_weights(PyObject *self, PyObject *args)
iw = (double*)PyArray_DATA(py_iw);
iw_zero = (char*)PyArray_DATA(py_iw_zero);
frequency_points = (double*)PyArray_DATA(py_frequency_points);
num_band0 = PyArray_DIMS(py_frequency_points)[0];
num_band0 = (long)PyArray_DIMS(py_frequency_points)[0];
relative_grid_address = (int(*)[4][3])PyArray_DATA(py_relative_grid_address);
mesh = (int*)PyArray_DATA(py_mesh);
triplets = (size_t(*)[3])PyArray_DATA(py_triplets);
num_triplets = PyArray_DIMS(py_triplets)[0];
triplets = (long(*)[3])PyArray_DATA(py_triplets);
num_triplets = (long)PyArray_DIMS(py_triplets)[0];
bz_grid_address = (int(*)[3])PyArray_DATA(py_bz_grid_address);
bz_map = (size_t*)PyArray_DATA(py_bz_map);
bz_map = (long*)PyArray_DATA(py_bz_map);
frequencies1 = (double*)PyArray_DATA(py_frequencies1);
frequencies2 = (double*)PyArray_DATA(py_frequencies2);
num_band1 = PyArray_DIMS(py_frequencies1)[1];
num_band2 = PyArray_DIMS(py_frequencies2)[1];
num_band1 = (long)PyArray_DIMS(py_frequencies1)[1];
num_band2 = (long)PyArray_DIMS(py_frequencies2)[1];
ph3py_get_integration_weight(iw,
iw_zero,
@ -1733,9 +1733,9 @@ py_set_triplets_integration_weights_with_sigma(PyObject *self, PyObject *args)
double *iw;
char *iw_zero;
double *frequency_points;
size_t (*triplets)[3];
long (*triplets)[3];
double *frequencies;
npy_intp num_band0, num_band, num_iw, num_triplets;
long num_band0, num_band, num_iw, num_triplets;
if (!PyArg_ParseTuple(args, "OOOOOdd",
&py_iw,
@ -1751,12 +1751,12 @@ py_set_triplets_integration_weights_with_sigma(PyObject *self, PyObject *args)
iw = (double*)PyArray_DATA(py_iw);
iw_zero = (char*)PyArray_DATA(py_iw_zero);
frequency_points = (double*)PyArray_DATA(py_frequency_points);
num_band0 = PyArray_DIMS(py_frequency_points)[0];
triplets = (size_t(*)[3])PyArray_DATA(py_triplets);
num_triplets = PyArray_DIMS(py_triplets)[0];
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 = PyArray_DIMS(py_frequencies)[1];
num_iw = PyArray_DIMS(py_iw)[0];
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,
@ -1783,8 +1783,8 @@ py_diagonalize_collision_matrix(PyObject *self, PyObject *args)
double *collision_matrix;
double *eigvals;
npy_intp num_temp, num_grid_point, num_band;
size_t num_column, adrs_shift;
long num_temp, num_grid_point, num_band;
long num_column, adrs_shift;
int info;
if (!PyArg_ParseTuple(args, "OOiidii",
@ -1803,11 +1803,11 @@ py_diagonalize_collision_matrix(PyObject *self, PyObject *args)
if (PyArray_NDIM(py_collision_matrix) == 2) {
num_temp = 1;
num_column = PyArray_DIM(py_collision_matrix, 1);
num_column = (long)PyArray_DIM(py_collision_matrix, 1);
} else {
num_temp = PyArray_DIM(py_collision_matrix, 1);
num_grid_point = PyArray_DIM(py_collision_matrix, 2);
num_band = PyArray_DIM(py_collision_matrix, 3);
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 {
@ -1838,8 +1838,8 @@ static PyObject * py_pinv_from_eigensolution(PyObject *self, PyObject *args)
double *collision_matrix;
double *eigvals;
npy_intp num_temp, num_grid_point, num_band;
size_t num_column, adrs_shift;
long num_temp, num_grid_point, num_band;
long num_column, adrs_shift;
if (!PyArg_ParseTuple(args, "OOiidi",
&py_collision_matrix,
@ -1853,9 +1853,9 @@ static PyObject * py_pinv_from_eigensolution(PyObject *self, PyObject *args)
collision_matrix = (double*)PyArray_DATA(py_collision_matrix);
eigvals = (double*)PyArray_DATA(py_eigenvalues);
num_temp = PyArray_DIMS(py_collision_matrix)[1];
num_grid_point = PyArray_DIMS(py_collision_matrix)[2];
num_band = PyArray_DIMS(py_collision_matrix)[3];
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;
@ -1889,14 +1889,14 @@ static PyObject * py_get_default_colmat_solver(PyObject *self, PyObject *args)
static void pinv_from_eigensolution(double *data,
const double *eigvals,
const size_t size,
const long size,
const double cutoff,
const int pinv_method)
{
size_t i, ib, j, k, max_l, i_s, j_s;
long i, ib, j, k, max_l, i_s, j_s;
double *tmp_data;
double e, sum;
size_t *l;
long *l;
l = NULL;
tmp_data = NULL;
@ -1908,7 +1908,7 @@ static void pinv_from_eigensolution(double *data,
tmp_data[i] = data[i];
}
l = (size_t*)malloc(sizeof(size_t) * size);
l = (long*)malloc(sizeof(long) * size);
max_l = 0;
for (i = 0; i < size; i++) {
if (pinv_method == 0) {
@ -1972,9 +1972,9 @@ static void pinv_from_eigensolution(double *data,
}
static void show_colmat_info(const PyArrayObject *py_collision_matrix,
const size_t i_sigma,
const size_t i_temp,
const size_t adrs_shift)
const long i_sigma,
const long i_temp,
const long adrs_shift)
{
int i;

132
c/collision_matrix.c
View File

@ -41,16 +41,16 @@
static void get_collision_matrix(double *collision_matrix,
const double *fc3_normal_squared,
const size_t num_band0,
const size_t num_band,
const long num_band0,
const long num_band,
const double *frequencies,
const size_t (*triplets)[3],
const size_t *triplets_map,
const size_t num_gp,
const size_t *map_q,
const size_t *rot_grid_points,
const size_t num_ir_gp,
const size_t num_rot,
const long (*triplets)[3],
const long *triplets_map,
const long num_gp,
const long *map_q,
const long *rot_grid_points,
const long num_ir_gp,
const long num_rot,
const double *rotations_cartesian,
const double *g,
const double temperature,
@ -59,46 +59,46 @@ static void get_collision_matrix(double *collision_matrix,
static void
get_reducible_collision_matrix(double *collision_matrix,
const double *fc3_normal_squared,
const size_t num_band0,
const size_t num_band,
const long num_band0,
const long num_band,
const double *frequencies,
const size_t (*triplets)[3],
const size_t *triplets_map,
const size_t num_gp,
const size_t *map_q,
const long (*triplets)[3],
const long *triplets_map,
const long num_gp,
const long *map_q,
const double *g,
const double temperature,
const double unit_conversion_factor,
const double cutoff_frequency);
static void get_inv_sinh(double *inv_sinh,
const size_t gp,
const long gp,
const double temperature,
const double *frequencies,
const size_t triplet[3],
const size_t *triplets_map,
const size_t *map_q,
const size_t num_band,
const long triplet[3],
const long *triplets_map,
const long *map_q,
const long num_band,
const double cutoff_frequency);
static size_t *create_gp2tp_map(const size_t *triplets_map,
const size_t num_gp);
static long *create_gp2tp_map(const long *triplets_map,
const long num_gp);
void col_get_collision_matrix(double *collision_matrix,
const Darray *fc3_normal_squared,
const double *frequencies,
const size_t (*triplets)[3],
const size_t *triplets_map,
const size_t *map_q,
const size_t *rot_grid_points,
const long (*triplets)[3],
const long *triplets_map,
const long *map_q,
const long *rot_grid_points,
const double *rotations_cartesian,
const double *g,
const size_t num_ir_gp,
const size_t num_gp,
const size_t num_rot,
const long num_ir_gp,
const long num_gp,
const long num_rot,
const double temperature,
const double unit_conversion_factor,
const double cutoff_frequency)
{
size_t num_triplets, num_band0, num_band;
long num_triplets, num_band0, num_band;
num_triplets = fc3_normal_squared->dims[0];
num_band0 = fc3_normal_squared->dims[1];
@ -127,16 +127,16 @@ void col_get_collision_matrix(double *collision_matrix,
void col_get_reducible_collision_matrix(double *collision_matrix,
const Darray *fc3_normal_squared,
const double *frequencies,
const size_t (*triplets)[3],
const size_t *triplets_map,
const size_t *map_q,
const long (*triplets)[3],
const long *triplets_map,
const long *map_q,
const double *g,
const size_t num_gp,
const long num_gp,
const double temperature,
const double unit_conversion_factor,
const double cutoff_frequency)
{
size_t num_triplets, num_band, num_band0;
long num_triplets, num_band, num_band0;
num_triplets = fc3_normal_squared->dims[0];
num_band0 = fc3_normal_squared->dims[1];
@ -160,24 +160,24 @@ void col_get_reducible_collision_matrix(double *collision_matrix,
static void get_collision_matrix(double *collision_matrix,
const double *fc3_normal_squared,
const size_t num_band0,
const size_t num_band,
const long num_band0,
const long num_band,
const double *frequencies,
const size_t (*triplets)[3],
const size_t *triplets_map,
const size_t num_gp,
const size_t *map_q,
const size_t *rot_grid_points,
const size_t num_ir_gp,
const size_t num_rot,
const long (*triplets)[3],
const long *triplets_map,
const long num_gp,
const long *map_q,
const long *rot_grid_points,
const long num_ir_gp,
const long num_rot,
const double *rotations_cartesian,
const double *g,
const double temperature,
const double unit_conversion_factor,
const double cutoff_frequency)
{
size_t i, j, k, l, m, n, ti, r_gp;
size_t *gp2tp_map;
long i, j, k, l, m, n, ti, r_gp;
long *gp2tp_map;
double collision;
double *inv_sinh;
@ -234,20 +234,20 @@ static void get_collision_matrix(double *collision_matrix,
static void
get_reducible_collision_matrix(double *collision_matrix,
const double *fc3_normal_squared,
const size_t num_band0,
const size_t num_band,
const long num_band0,
const long num_band,
const double *frequencies,
const size_t (*triplets)[3],
const size_t *triplets_map,
const size_t num_gp,
const size_t *map_q,
const long (*triplets)[3],
const long *triplets_map,
const long num_gp,
const long *map_q,
const double *g,
const double temperature,
const double unit_conversion_factor,
const double cutoff_frequency)
{
size_t i, j, k, l, ti;
size_t *gp2tp_map;
long i, j, k, l, ti;
long *gp2tp_map;
double collision;
double *inv_sinh;
@ -293,16 +293,16 @@ get_reducible_collision_matrix(double *collision_matrix,
}
static void get_inv_sinh(double *inv_sinh,
const size_t gp,
const long gp,
const double temperature,
const double *frequencies,
const size_t triplet[3],
const size_t *triplets_map,
const size_t *map_q,
const size_t num_band,
const long triplet[3],
const long *triplets_map,
const long *map_q,
const long num_band,
const double cutoff_frequency)
{
size_t i, gp2;
long i, gp2;
double f;
if (triplets_map[gp] == map_q[gp]) {
@ -320,11 +320,11 @@ static void get_inv_sinh(double *inv_sinh,
}
}
static size_t *create_gp2tp_map(const size_t *triplets_map,
const size_t num_gp)
static long *create_gp2tp_map(const long *triplets_map,
const long num_gp)
{
size_t i, max_i, count;
size_t *gp2tp_map;
long i, max_i, count;
long *gp2tp_map;
max_i = 0;
for (i = 0; i < num_gp; i++) {
@ -333,7 +333,7 @@ static size_t *create_gp2tp_map(const size_t *triplets_map,
}
}
gp2tp_map = (size_t*)malloc(sizeof(size_t) * (max_i + 1));
gp2tp_map = (long*)malloc(sizeof(long) * (max_i + 1));
for (i = 0; i < max_i + 1; i++) {
gp2tp_map[i] = 0;
}
@ -344,7 +344,7 @@ static size_t *create_gp2tp_map(const size_t *triplets_map,
gp2tp_map[i] = count;
count++;
}
}
}
return gp2tp_map;
}

22
c/collision_matrix.h
View File

@ -40,26 +40,26 @@
void col_get_collision_matrix(double *collision_matrix,
const Darray *fc3_normal_squared,
const double *frequencies,
const size_t (*triplets)[3],
const size_t *triplets_map,
const size_t *map_q,
const size_t *rot_grid_points,
const long (*triplets)[3],
const long *triplets_map,
const long *map_q,
const long *rot_grid_points,
const double *rotations_cartesian,
const double *g,
const size_t num_ir_gp,
const size_t num_gp,
const size_t num_rot,
const long num_ir_gp,
const long num_gp,
const long num_rot,
const double temperature,
const double unit_conversion_factor,
const double cutoff_frequency);
void col_get_reducible_collision_matrix(double *collision_matrix,
const Darray *fc3_normal_squared,
const double *frequencies,
const size_t (*triplets)[3],
const size_t *triplets_map,
const size_t *map_q,
const long (*triplets)[3],
const long *triplets_map,
const long *map_q,
const double *g,
const size_t num_gp,
const long num_gp,
const double temperature,
const double unit_conversion_factor,
const double cutoff_frequency);

102
c/fc3.c
View File

@ -37,14 +37,14 @@
#include "fc3.h"
static void rotate_delta_fc2s(double (*rot_delta_fc2s)[3][3],
const size_t i_atom,
const size_t j_atom,
const long i_atom,
const long j_atom,
PHPYCONST double (*delta_fc2s)[3][3],
PHPYCONST double (*site_sym_cart)[3][3],
const int *rot_map_sym,
const size_t num_atom,
const size_t num_site_sym,
const size_t num_disp);
const long num_atom,
const long num_site_sym,
const long num_disp);
static void tensor2_rotation(double rot_tensor[3][3],
PHPYCONST double tensor[3][3],
PHPYCONST double r[3][3]);
@ -59,44 +59,44 @@ static void copy_permutation_symmetry_fc3_elem(double *fc3,
const int a,
const int b,
const int c,
const size_t num_atom);
const long num_atom);
static void set_permutation_symmetry_fc3_elem(double *fc3_elem,
const double *fc3,
const int a,
const int b,
const int c,
const size_t num_atom);
const long num_atom);
static void set_permutation_symmetry_compact_fc3(double * fc3,
const int p2s[],
const int s2pp[],
const int nsym_list[],
const int perms[],
const size_t n_satom,
const size_t n_patom);
const long n_satom,
const long n_patom);
static void transpose_compact_fc3_type01(double * fc3,
const int p2s[],
const int s2pp[],
const int nsym_list[],
const int perms[],
const size_t n_satom,
const size_t n_patom,
const long n_satom,
const long n_patom,
const int t_type);
static void transpose_compact_fc3_type2(double * fc3,
const int p2s[],
const int s2pp[],
const int nsym_list[],
const int perms[],
const size_t n_satom,
const size_t n_patom);
const long n_satom,
const long n_patom);
void fc3_distribute_fc3(double *fc3,
const int target,
const int source,
const int *atom_mapping,
const size_t num_atom,
const long num_atom,
const double *rot_cart)
{
size_t i, j, adrs_out, adrs_in;
long i, j, adrs_out, adrs_in;
for (i = 0; i < num_atom; i++) {
for (j = 0; j < num_atom; j++) {
@ -114,11 +114,11 @@ void fc3_rotate_delta_fc2(double (*fc3)[3][3][3],
const double *inv_U,
PHPYCONST double (*site_sym_cart)[3][3],
const int *rot_map_syms,
const size_t num_atom,
const size_t num_site_sym,
const size_t num_disp)
const long num_atom,
const long num_site_sym,
const long num_disp)
{
size_t i_atoms, i, j, k, l, m, n;
long i_atoms, i, j, k, l, m, n;
double (*rot_delta_fc2s)[3][3];
rot_delta_fc2s = (double(*)[3][3]) malloc(sizeof(double[3][3])
@ -153,10 +153,10 @@ void fc3_rotate_delta_fc2(double (*fc3)[3][3][3],
rot_delta_fc2s = NULL;
}
void fc3_set_permutation_symmetry_fc3(double *fc3, const size_t num_atom)
void fc3_set_permutation_symmetry_fc3(double *fc3, const long num_atom)
{
double fc3_elem[27];
size_t i, j, k;
long i, j, k;
#pragma omp parallel for private(j, k, fc3_elem)
for (i = 0; i < num_atom; i++) {
@ -175,8 +175,8 @@ void fc3_set_permutation_symmetry_compact_fc3(double * fc3,
const int s2pp[],
const int nsym_list[],
const int perms[],
const size_t n_satom,
const size_t n_patom)
const long n_satom,
const long n_patom)
{
set_permutation_symmetry_compact_fc3(fc3,
p2s,
@ -192,8 +192,8 @@ void fc3_transpose_compact_fc3(double * fc3,
const int s2pp[],
const int nsym_list[],
const int perms[],
const size_t n_satom,
const size_t n_patom,
const long n_satom,
const long n_patom,
const int t_type)
{
/* Three types of index permutations */
@ -223,16 +223,16 @@ void fc3_transpose_compact_fc3(double * fc3,
}
static void rotate_delta_fc2s(double (*rot_delta_fc2s)[3][3],
const size_t i_atom,
const size_t j_atom,
const long i_atom,
const long j_atom,
PHPYCONST double (*delta_fc2s)[3][3],
PHPYCONST double (*site_sym_cart)[3][3],
const int *rot_map_sym,
const size_t num_atom,
const size_t num_site_sym,
const size_t num_disp)
const long num_atom,
const long num_site_sym,
const long num_disp)
{
size_t i, j;
long i, j;
for (i = 0; i < num_disp; i++) {
for (j = 0; j < num_site_sym; j++) {
@ -249,7 +249,7 @@ static void tensor2_rotation(double rot_tensor[3][3],
PHPYCONST double tensor[3][3],
PHPYCONST double r[3][3])
{
size_t i, j, k, l;
long i, j, k, l;
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
@ -272,7 +272,7 @@ static void tensor3_rotation(double *rot_tensor,
const double *tensor,
const double *rot_cartesian)
{
size_t l;
long l;
for (l = 0; l < 27; l++) {
rot_tensor[l] = tensor3_rotation_elem(tensor, rot_cartesian, l);
@ -283,7 +283,7 @@ static double tensor3_rotation_elem(const double *tensor,
const double *r,
const int pos)
{
size_t i, j, k, l, m, n;
long i, j, k, l, m, n;
double sum;
l = pos / 9;
@ -307,9 +307,9 @@ static void copy_permutation_symmetry_fc3_elem(double *fc3,
const int a,
const int b,
const int c,
const size_t num_atom)
const long num_atom)
{
size_t i, j, k;
long i, j, k;
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
@ -348,9 +348,9 @@ static void set_permutation_symmetry_fc3_elem(double *fc3_elem,
const int a,
const int b,
const int c,
const size_t num_atom)
const long num_atom)
{
size_t i, j, k;
long i, j, k;
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
@ -384,17 +384,17 @@ static void set_permutation_symmetry_compact_fc3(double * fc3,
const int s2pp[],
const int nsym_list[],
const int perms[],
const size_t n_satom,
const size_t n_patom)
const long n_satom,
const long n_patom)
{
/* fc3 shape=(n_patom, n_satom, n_satom, 3, 3, 3) */
/* 1D indexing: */
/* i * n_satom * n_satom * 27 + j * n_satom * 27 + */
/* k * 27 + l * 9 + m * 3 + n */
size_t i, j, k, l, m, n, i_p, j_p, k_p;
long i, j, k, l, m, n, i_p, j_p, k_p;
int done_any;
size_t i_trans_j, k_trans_j, i_trans_k, j_trans_k;
size_t adrs[6];
long i_trans_j, k_trans_j, i_trans_k, j_trans_k;
long adrs[6];
double fc3_elem[3][3][3];
char *done;
@ -479,16 +479,16 @@ static void transpose_compact_fc3_type01(double * fc3,
const int s2pp[],
const int nsym_list[],
const int perms[],
const size_t n_satom,
const size_t n_patom,
const long n_satom,
const long n_patom,
const int t_type)
{
/* Three types of index permutations */
/* t_type=0: dim[0] <-> dim[1] */
/* t_type=1: dim[0] <-> dim[2] */
/* t_type=2: dim[1] <-> dim[2] */
size_t i, j, k, l, m, n, i_p, j_p, i_trans, k_trans;
size_t adrs, adrs_t;
long i, j, k, l, m, n, i_p, j_p, i_trans, k_trans;
long adrs, adrs_t;
double fc3_elem[3][3][3];
char *done;
@ -585,11 +585,11 @@ static void transpose_compact_fc3_type2(double * fc3,
const int s2pp[],
const int nsym_list[],
const int perms[],
const size_t n_satom,
const size_t n_patom)
const long n_satom,
const long n_patom)
{
size_t j, k, l, m, n, i_p;
size_t adrs, adrs_t;
long j, k, l, m, n, i_p;
long adrs, adrs_t;
double fc3_elem[3][3][3];
for (i_p = 0; i_p < n_patom; i_p++) {

18
c/fc3.h
View File

@ -39,31 +39,31 @@ void fc3_distribute_fc3(double *fc3,
const int target,
const int source,
const int *atom_mapping,
const size_t num_atom,
const long num_atom,
const double *rot_cart);
void fc3_rotate_delta_fc2(double (*fc3)[3][3][3],
PHPYCONST double (*delta_fc2s)[3][3],
const double *inv_U,
PHPYCONST double (*site_sym_cart)[3][3],
const int *rot_map_syms,
const size_t num_atom,
const size_t num_site_sym,
const size_t num_disp);
void fc3_set_permutation_symmetry_fc3(double *fc3, const size_t num_atom);
const long num_atom,
const long num_site_sym,
const long num_disp);
void fc3_set_permutation_symmetry_fc3(double *fc3, const long num_atom);
void fc3_set_permutation_symmetry_compact_fc3(double * fc3,
const int p2s[],
const int s2pp[],
const int nsym_list[],
const int perms[],
const size_t n_satom,
const size_t n_patom);
const long n_satom,
const long n_patom);
void fc3_transpose_compact_fc3(double * fc3,
const int p2s[],
const int s2pp[],
const int nsym_list[],
const int perms[],
const size_t n_satom,
const size_t n_patom,
const long n_satom,
const long n_patom,
const int t_type);
#endif

88
c/imag_self_energy_with_g.c
View File

@ -42,26 +42,26 @@
#include "triplet.h"
static int ise_set_g_pos_frequency_point(int (*g_pos)[4],
const size_t num_band0,
const size_t num_band,
const long num_band0,
const long num_band,
const char *g_zero);
static void
detailed_imag_self_energy_at_triplet(double *detailed_imag_self_energy,
double *imag_self_energy,
const size_t num_band0,
const size_t num_band,
const long num_band0,
const long num_band,
const double *fc3_normal_squared,
const double *frequencies,
const size_t triplet[3],
const long triplet[3],
const double *g1,
const double *g2_3,
const char *g_zero,
const double *temperatures,
const size_t num_temps,
const long num_temps,
const double cutoff_frequency);
static double
collect_detailed_imag_self_energy(double *imag_self_energy,
const size_t num_band,
const long num_band,
const double *fc3_normal_squared,
const double *n1,
const double *n2,
@ -70,7 +70,7 @@ collect_detailed_imag_self_energy(double *imag_self_energy,
const char *g_zero);
static double
collect_detailed_imag_self_energy_0K(double *imag_self_energy,
const size_t num_band,
const long num_band,
const double *fc3_normal_squared,
const double *n1,
const double *n2,
@ -78,17 +78,17 @@ collect_detailed_imag_self_energy_0K(double *imag_self_energy,
const char *g_zero);
static void set_occupations(double *n1,
double *n2,
const size_t num_band,
const long num_band,
const double temperature,
const size_t triplet[3],
const long triplet[3],
const double *frequencies,
const double cutoff_frequency);
void ise_get_imag_self_energy_at_bands_with_g(double *imag_self_energy,
const Darray *fc3_normal_squared,
const double *frequencies,
const size_t (*triplets)[3],
const int *weights,
const long (*triplets)[3],
const long *triplet_weights,
const double *g,
const char *g_zero,
const double temperature,
@ -96,8 +96,8 @@ void ise_get_imag_self_energy_at_bands_with_g(double *imag_self_energy,
const int num_frequency_points,
const int frequency_point_index)
{
size_t i, j, num_triplets, num_band0, num_band, num_band_prod;
size_t num_g_pos, g_index_dims, g_index_shift;
long i, j, num_triplets, num_band0, num_band, num_band_prod;
long num_g_pos, g_index_dims, g_index_shift;
int (*g_pos)[4];
double *ise;
int at_a_frequency_point;
@ -150,7 +150,7 @@ void ise_get_imag_self_energy_at_bands_with_g(double *imag_self_energy,
fc3_normal_squared->data + i * num_band_prod,
frequencies,
triplets[i],
weights[i],
triplet_weights[i],
g + i * g_index_dims + g_index_shift,
g + (i + num_triplets) * g_index_dims + g_index_shift,
g_pos,
@ -185,8 +185,8 @@ void ise_get_detailed_imag_self_energy_at_bands_with_g
double *imag_self_energy_U,
const Darray *fc3_normal_squared,
const double *frequencies,
const size_t (*triplets)[3],
const int *weights,
const long (*triplets)[3],
const long *triplet_weights,
const int *grid_address,
const double *g,
const char *g_zero,
@ -194,7 +194,7 @@ void ise_get_detailed_imag_self_energy_at_bands_with_g
const double cutoff_frequency)
{
double *ise;
size_t i, j, num_triplets, num_band0, num_band, num_band_prod;
long i, j, num_triplets, num_band0, num_band, num_band_prod;
int *is_N;
double ise_tmp, N, U;
@ -237,7 +237,7 @@ void ise_get_detailed_imag_self_energy_at_bands_with_g
U = 0;
/* #pragma omp parallel for private(ise_tmp) reduction(+:N,U) */
for (j = 0; j < num_triplets; j++) {
ise_tmp = ise[j * num_band0 + i] * weights[j];
ise_tmp = ise[j * num_band0 + i] * triplet_weights[j];
if (is_N[j]) {
N += ise_tmp;
} else {
@ -255,23 +255,23 @@ void ise_get_detailed_imag_self_energy_at_bands_with_g
}
void ise_imag_self_energy_at_triplet(double *imag_self_energy,
const size_t num_band0,
const size_t num_band,
const long num_band0,
const long num_band,
const double *fc3_normal_squared,
const double *frequencies,
const size_t triplet[3],
const int triplet_weight,
const long triplet[3],
const long triplet_weight,
const double *g1,
const double *g2_3,
PHPYCONST int (*g_pos)[4],
const size_t num_g_pos,
const long num_g_pos,
const double *temperatures,
const size_t num_temps,
const long num_temps,
const double cutoff_frequency,
const int openmp_at_bands,
const int at_a_frequency_point)
{
size_t i, j;
long i, j;
double *n1, *n2;
int g_pos_3;
@ -329,11 +329,11 @@ void ise_imag_self_energy_at_triplet(double *imag_self_energy,
}
int ise_set_g_pos(int (*g_pos)[4],
const size_t num_band0,
const size_t num_band,
const long num_band0,
const long num_band,
const char *g_zero)
{
size_t num_g_pos, j, k, l, jkl;
long num_g_pos, j, k, l, jkl;
num_g_pos = 0;
jkl = 0;
@ -355,11 +355,11 @@ int ise_set_g_pos(int (*g_pos)[4],
}
static int ise_set_g_pos_frequency_point(int (*g_pos)[4],
const size_t num_band0,
const size_t num_band,
const long num_band0,
const long num_band,
const char *g_zero)
{
size_t num_g_pos, j, k, l, kl, jkl;
long num_g_pos, j, k, l, kl, jkl;
num_g_pos = 0;
jkl = 0;
@ -385,19 +385,19 @@ static int ise_set_g_pos_frequency_point(int (*g_pos)[4],
static void
detailed_imag_self_energy_at_triplet(double *detailed_imag_self_energy,
double *imag_self_energy,
const size_t num_band0,
const size_t num_band,
const long num_band0,
const long num_band,
const double *fc3_normal_squared,
const double *frequencies,
const size_t triplet[3],
const long triplet[3],
const double *g1,
const double *g2_3,
const char *g_zero,
const double *temperatures,
const size_t num_temps,
const long num_temps,
const double cutoff_frequency)
{
size_t i, j, adrs_shift;
long i, j, adrs_shift;
double *n1, *n2;
n1 = NULL;
@ -450,7 +450,7 @@ detailed_imag_self_energy_at_triplet(double *detailed_imag_self_energy,
static double
collect_detailed_imag_self_energy(double *imag_self_energy,
const size_t num_band,
const long num_band,
const double *fc3_normal_squared,
const double *n1,
const double *n2,
@ -458,7 +458,7 @@ collect_detailed_imag_self_energy(double *imag_self_energy,
const double *g2_3,
const char *g_zero)
{
size_t ij, i, j;
long ij, i, j;
double sum_g;
sum_g = 0;
@ -479,14 +479,14 @@ collect_detailed_imag_self_energy(double *imag_self_energy,
static double
collect_detailed_imag_self_energy_0K(double *imag_self_energy,
const size_t num_band,
const long num_band,
const double *fc3_normal_squared,
const double *n1,
const double *n2,
const double *g1,
const char *g_zero)
{
size_t ij, i, j;
long ij, i, j;
double sum_g;
sum_g = 0;
@ -505,13 +505,13 @@ collect_detailed_imag_self_energy_0K(double *imag_self_energy,
static void set_occupations(double *n1,
double *n2,
const size_t num_band,
const long num_band,
const double temperature,
const size_t triplet[3],
const long triplet[3],
const double *frequencies,
const double cutoff_frequency)
{
size_t j;
long j;
double f1, f2;
for (j = 0; j < num_band; j++) {

24
c/imag_self_energy_with_g.h
View File

@ -41,8 +41,8 @@
void ise_get_imag_self_energy_at_bands_with_g(double *imag_self_energy,
const Darray *fc3_normal_squared,
const double *frequencies,
const size_t (*triplets)[3],
const int *weights,
const long (*triplets)[3],
const long *triplet_weights,
const double *g,
const char *g_zero,
const double temperature,
@ -55,32 +55,32 @@ void ise_get_detailed_imag_self_energy_at_bands_with_g
double *imag_self_energy_U,
const Darray *fc3_normal_squared,
const double *frequencies,
const size_t (*triplets)[3],
const int *weights,
const long (*triplets)[3],
const long *triplet_weights,
const int *grid_address,
const double *g,
const char *g_zero,
const double temperature,
const double cutoff_frequency);
void ise_imag_self_energy_at_triplet(double *imag_self_energy,
const size_t num_band0,
const size_t num_band,
const long num_band0,
const long num_band,
const double *fc3_normal_squared,
const double *frequencies,
const size_t triplet[3],
const int triplet_weight,
const long triplet[3],
const long triplet_weight,
const double *g1,
const double *g2_3,
PHPYCONST int (*g_pos)[4],
const size_t num_g_pos,
const long num_g_pos,
const double *temperatures,
const size_t num_temps,
const long num_temps,
const double cutoff_frequency,
const int openmp_at_bands,
const int at_a_frequency_point);
int ise_set_g_pos(int (*g_pos)[4],
const size_t num_band0,
const size_t num_band,
const long num_band0,
const long num_band,
const char *g_zero);
#endif

68
c/interaction.c
View File

@ -50,7 +50,7 @@ static const int index_exchange[6][3] = {{0, 1, 2},
{1, 0, 2}};
static void real_to_normal(double *fc3_normal_squared,
PHPYCONST int (*g_pos)[4],
const size_t num_g_pos,
const long num_g_pos,
const double *freqs0,
const double *freqs1,
const double *freqs2,
@ -67,15 +67,15 @@ static void real_to_normal(double *fc3_normal_squared,
const int *p2s_map,
const int *s2p_map,
const int *band_indices,
const size_t num_band0,
const size_t num_band,
const long num_band0,
const long num_band,
const double cutoff_frequency,
const size_t triplet_index,
const size_t num_triplets,
const long triplet_index,
const long num_triplets,
const int openmp_at_bands);
static void real_to_normal_sym_q(double *fc3_normal_squared,
PHPYCONST int (*g_pos)[4],
const size_t num_g_pos,
const long num_g_pos,
PHPYCONST double *freqs[3],
PHPYCONST lapack_complex_double *eigvecs[3],
const double *fc3,
@ -88,11 +88,11 @@ static void real_to_normal_sym_q(double *fc3_normal_squared,
const int *p2s_map,
const int *s2p_map,
const int *band_indices,
const size_t num_band0,
const size_t num_band,
const long num_band0,
const long num_band,
const double cutoff_frequency,
const size_t triplet_index,
const size_t num_triplets,
const long triplet_index,
const long num_triplets,
const int openmp_at_bands);
/* fc3_normal_squared[num_triplets, num_band0, num_band, num_band] */
@ -100,8 +100,8 @@ void itr_get_interaction(Darray *fc3_normal_squared,
const char *g_zero,
const Darray *frequencies,
const lapack_complex_double *eigenvectors,
const size_t (*triplets)[3],
const size_t num_triplets,
const long (*triplets)[3],
const long num_triplets,
const int *grid_address,
const int *mesh,
const double *fc3,
@ -118,8 +118,8 @@ void itr_get_interaction(Darray *fc3_normal_squared,
{
int openmp_per_triplets;
int (*g_pos)[4];
size_t i;
size_t num_band, num_band0, num_band_prod, num_g_pos;
long i;
long num_band, num_band0, num_band_prod, num_g_pos;
g_pos = NULL;
@ -173,13 +173,13 @@ void itr_get_interaction(Darray *fc3_normal_squared,
}
void itr_get_interaction_at_triplet(double *fc3_normal_squared,
const size_t num_band0,
const size_t num_band,
const long num_band0,
const long num_band,
PHPYCONST int (*g_pos)[4],
const size_t num_g_pos,
const long num_g_pos,
const double *frequencies,
const lapack_complex_double *eigenvectors,
const size_t triplet[3],
const long triplet[3],
const int *grid_address,
const int *mesh,
const double *fc3,
@ -193,11 +193,11 @@ void itr_get_interaction_at_triplet(double *fc3_normal_squared,
const int *band_indices,
const int symmetrize_fc3_q,
const double cutoff_frequency,
const size_t triplet_index, /* only for print */
const size_t num_triplets, /* only for print */
const long triplet_index, /* only for print */
const long num_triplets, /* only for print */
const int openmp_at_bands)
{
size_t j, k;
long j, k;
double *freqs[3];
lapack_complex_double *eigvecs[3];
double q[9];
@ -278,7 +278,7 @@ void itr_get_interaction_at_triplet(double *fc3_normal_squared,
static void real_to_normal(double *fc3_normal_squared,
PHPYCONST int (*g_pos)[4],
const size_t num_g_pos,
const long num_g_pos,
const double *freqs0,
const double *freqs1,
const double *freqs2,
@ -295,14 +295,14 @@ static void real_to_normal(double *fc3_normal_squared,
const int *p2s_map,
const int *s2p_map,
const int *band_indices,
const size_t num_band0,
const size_t num_band,
const long num_band0,
const long num_band,
const double cutoff_frequency,
const size_t triplet_index,
const size_t num_triplets,
const long triplet_index,
const long num_triplets,
const int openmp_at_bands)
{
size_t num_patom;
long num_patom;
lapack_complex_double *fc3_reciprocal;
num_patom = num_band / 3;
@ -350,7 +350,7 @@ static void real_to_normal(double *fc3_normal_squared,
static void real_to_normal_sym_q(double *fc3_normal_squared,
PHPYCONST int (*g_pos)[4],
const size_t num_g_pos,
const long num_g_pos,
PHPYCONST double *freqs[3],
PHPYCONST lapack_complex_double *eigvecs[3],
const double *fc3,
@ -363,15 +363,15 @@ static void real_to_normal_sym_q(double *fc3_normal_squared,
const int *p2s_map,
const int *s2p_map,
const int *band_indices,
const size_t num_band0,
const size_t num_band,
const long num_band0,
const long num_band,
const double cutoff_frequency,
const size_t triplet_index,
const size_t num_triplets,
const long triplet_index,
const long num_triplets,
const int openmp_at_bands)
{
size_t i, j, k, l;
size_t band_ex[3];
long i, j, k, l;
long band_ex[3];
double q_ex[9];
double *fc3_normal_squared_ex;

16
c/interaction.h
View File

@ -42,8 +42,8 @@ void itr_get_interaction(Darray *fc3_normal_squared,
const char *g_zero,
const Darray *frequencies,
const lapack_complex_double *eigenvectors,
const size_t (*triplets)[3],
const size_t num_triplets,
const long (*triplets)[3],
const long num_triplets,
const int *grid_address,
const int *mesh,
const double *fc3,
@ -58,13 +58,13 @@ void itr_get_interaction(Darray *fc3_normal_squared,
const int symmetrize_fc3_q,
const double cutoff_frequency);
void itr_get_interaction_at_triplet(double *fc3_normal_squared,
const size_t num_band0,
const size_t num_band,
const long num_band0,
const long num_band,
PHPYCONST int (*g_pos)[4],
const size_t num_g_pos,
const long num_g_pos,
const double *frequencies,
const lapack_complex_double *eigenvectors,
const size_t triplet[3],
const long triplet[3],
const int *grid_address,
const int *mesh,
const double *fc3,
@ -78,8 +78,8 @@ void itr_get_interaction_at_triplet(double *fc3_normal_squared,
const int *band_indices,
const int symmetrize_fc3_q,
const double cutoff_frequency,
const size_t triplet_index, /* only for print */
const size_t num_triplets, /* only for print */
const long triplet_index, /* only for print */
const long num_triplets, /* only for print */
const int openmp_at_bands);
#endif

24
c/isotope.c
View File

@ -40,18 +40,18 @@
void
iso_get_isotope_scattering_strength(double *gamma,
const size_t grid_point,
const long grid_point,
const double *mass_variances,
const double *frequencies,
const lapack_complex_double *eigenvectors,
const size_t num_grid_points,
const long num_grid_points,
const int *band_indices,
const size_t num_band,
const size_t num_band0,
const long num_band,
const long num_band0,
const double sigma,
const double cutoff_frequency)
{
size_t i, j, k, l, m;
long i, j, k, l, m;
double *e0_r, *e0_i, e1_r, e1_i, a, b, f, *f0, dist, sum_g, sum_g_k;
e0_r = (double*)malloc(sizeof(double) * num_band * num_band0);
@ -128,20 +128,20 @@ iso_get_isotope_scattering_strength(double *gamma,
void iso_get_thm_isotope_scattering_strength
(double *gamma,
const size_t grid_point,
const size_t *ir_grid_points,
const int *weights,
const long grid_point,
const long *ir_grid_points,
const long *weights,
const double *mass_variances,
const double *frequencies,
const lapack_complex_double *eigenvectors,
const size_t num_grid_points,
const long num_grid_points,
const int *band_indices,
const size_t num_band,
const size_t num_band0,
const long num_band,
const long num_band0,
const double *integration_weights,
const double cutoff_frequency)
{
size_t i, j, k, l, m, gp;
long i, j, k, l, m, gp;
double *e0_r, *e0_i, *f0, *gamma_ij;
double e1_r, e1_i, a, b, f, dist, sum_g_k;

20
c/isotope.h
View File

@ -39,28 +39,28 @@
void
iso_get_isotope_scattering_strength(double *gamma,
const size_t grid_point,
const long grid_point,
const double *mass_variances,
const double *frequencies,
const lapack_complex_double *eigenvectors,
const size_t num_grid_points,
const long num_grid_points,
const int *band_indices,
const size_t num_band,
const size_t num_band0,
const long num_band,
const long num_band0,
const double sigma,
const double cutoff_frequency);
void iso_get_thm_isotope_scattering_strength
(double *gamma,
const size_t grid_point,
const size_t *ir_grid_points,
const int *weights,
const long grid_point,
const long *ir_grid_points,
const long *weights,
const double *mass_variances,
const double *frequencies,
const lapack_complex_double *eigenvectors,
const size_t num_grid_points,
const long num_grid_points,
const int *band_indices,
const size_t num_band,
const size_t num_band0,
const long num_band,
const long num_band0,
const double *integration_weights,
const double cutoff_frequency);
#endif

546
c/kpoint.c
View File

@ -46,221 +46,41 @@
#define warning_print(...)
#endif
#define KPT_NUM_BZ_SEARCH_SPACE 125
static int bz_search_space[KPT_NUM_BZ_SEARCH_SPACE][3] = {
{ 0, 0, 0},
{ 0, 0, 1},
{ 0, 0, 2},
{ 0, 0, -2},
{ 0, 0, -1},
{ 0, 1, 0},
{ 0, 1, 1},
{ 0, 1, 2},
{ 0, 1, -2},
{ 0, 1, -1},
{ 0, 2, 0},
{ 0, 2, 1},
{ 0, 2, 2},
{ 0, 2, -2},
{ 0, 2, -1},
{ 0, -2, 0},
{ 0, -2, 1},
{ 0, -2, 2},
{ 0, -2, -2},
{ 0, -2, -1},
{ 0, -1, 0},
{ 0, -1, 1},
{ 0, -1, 2},
{ 0, -1, -2},
{ 0, -1, -1},
{ 1, 0, 0},
{ 1, 0, 1},
{ 1, 0, 2},
{ 1, 0, -2},
{ 1, 0, -1},
{ 1, 1, 0},
{ 1, 1, 1},
{ 1, 1, 2},
{ 1, 1, -2},
{ 1, 1, -1},
{ 1, 2, 0},
{ 1, 2, 1},
{ 1, 2, 2},
{ 1, 2, -2},
{ 1, 2, -1},
{ 1, -2, 0},
{ 1, -2, 1},
{ 1, -2, 2},
{ 1, -2, -2},
{ 1, -2, -1},
{ 1, -1, 0},
{ 1, -1, 1},
{ 1, -1, 2},
{ 1, -1, -2},
{ 1, -1, -1},
{ 2, 0, 0},
{ 2, 0, 1},
{ 2, 0, 2},
{ 2, 0, -2},
{ 2, 0, -1},
{ 2, 1, 0},
{ 2, 1, 1},
{ 2, 1, 2},
{ 2, 1, -2},
{ 2, 1, -1},
{ 2, 2, 0},
{ 2, 2, 1},
{ 2, 2, 2},
{ 2, 2, -2},
{ 2, 2, -1},
{ 2, -2, 0},
{ 2, -2, 1},
{ 2, -2, 2},
{ 2, -2, -2},
{ 2, -2, -1},
{ 2, -1, 0},
{ 2, -1, 1},
{ 2, -1, 2},
{ 2, -1, -2},
{ 2, -1, -1},
{-2, 0, 0},
{-2, 0, 1},
{-2, 0, 2},
{-2, 0, -2},
{-2, 0, -1},
{-2, 1, 0},
{-2, 1, 1},
{-2, 1, 2},
{-2, 1, -2},
{-2, 1, -1},
{-2, 2, 0},
{-2, 2, 1},
{-2, 2, 2},
{-2, 2, -2},
{-2, 2, -1},
{-2, -2, 0},
{-2, -2, 1},
{-2, -2, 2},
{-2, -2, -2},
{-2, -2, -1},
{-2, -1, 0},
{-2, -1, 1},
{-2, -1, 2},
{-2, -1, -2},
{-2, -1, -1},
{-1, 0, 0},
{-1, 0, 1},
{-1, 0, 2},
{-1, 0, -2},
{-1, 0, -1},
{-1, 1, 0},
{-1, 1, 1},
{-1, 1, 2},
{-1, 1, -2},
{-1, 1, -1},
{-1, 2, 0},
{-1, 2, 1},
{-1, 2, 2},
{-1, 2, -2},
{-1, 2, -1},
{-1, -2, 0},
{-1, -2, 1},
{-1, -2, 2},
{-1, -2, -2},
{-1, -2, -1},
{-1, -1, 0},
{-1, -1, 1},
{-1, -1, 2},
{-1, -1, -2},
{-1, -1, -1}
};
static MatINT *get_point_group_reciprocal(const MatINT * rotations,
const int is_time_reversal);
static MatINT *get_point_group_reciprocal_with_q(const MatINT * rot_reciprocal,
const double symprec,
const size_t num_q,
const long num_q,
SPGCONST double qpoints[][3]);
static size_t get_dense_ir_reciprocal_mesh(int grid_address[][3],
size_t ir_mapping_table[],
const int mesh[3],
const int is_shift[3],
const MatINT *rot_reciprocal);
static size_t get_dense_ir_reciprocal_mesh_normal(int grid_address[][3],
size_t ir_mapping_table[],
const int mesh[3],
const int is_shift[3],
const MatINT *rot_reciprocal);
static size_t get_dense_ir_reciprocal_mesh_distortion(int grid_address[][3],
size_t ir_mapping_table[],
const int mesh[3],
const int is_shift[3],
const MatINT *rot_reciprocal);
static size_t get_dense_num_ir(size_t ir_mapping_table[], const int mesh[3]);
static size_t relocate_dense_BZ_grid_address(int bz_grid_address[][3],
size_t bz_map[],
SPGCONST int grid_address[][3],
const int mesh[3],
SPGCONST double rec_lattice[3][3],
const int is_shift[3]);
static double get_tolerance_for_BZ_reduction(SPGCONST double rec_lattice[3][3],
const int mesh[3]);
static long get_dense_ir_reciprocal_mesh(int grid_address[][3],
long ir_mapping_table[],
const int mesh[3],
const int is_shift[3],
const MatINT *rot_reciprocal);
static long get_dense_ir_reciprocal_mesh_normal(int grid_address[][3],
long ir_mapping_table[],
const int mesh[3],
const int is_shift[3],
const MatINT *rot_reciprocal);
static long get_dense_ir_reciprocal_mesh_distortion(int grid_address[][3],
long ir_mapping_table[],
const int mesh[3],
const int is_shift[3],
const MatINT *rot_reciprocal);
static long get_dense_num_ir(long ir_mapping_table[], const int mesh[3]);
static int check_mesh_symmetry(const int mesh[3],
const int is_shift[3],
const MatINT *rot_reciprocal);
/* grid_address (e.g. 4x4x4 mesh, unless GRID_ORDER_XYZ is defined) */
/* [[ 0 0 0] */
/* [ 1 0 0] */
/* [ 2 0 0] */
/* [-1 0 0] */
/* [ 0 1 0] */
/* [ 1 1 0] */
/* [ 2 1 0] */
/* [-1 1 0] */
/* .... ] */
/* */
/* Each value of 'map' correspnds to the index of grid_point. */
int kpt_get_irreducible_reciprocal_mesh(int grid_address[][3],
int ir_mapping_table[],
const int mesh[3],
const int is_shift[3],
const MatINT *rot_reciprocal)
long kpt_get_dense_irreducible_reciprocal_mesh(int grid_address[][3],
long ir_mapping_table[],
const int mesh[3],
const int is_shift[3],
const MatINT *rot_reciprocal)
{
int num_ir;
size_t i;
size_t *dense_ir_mapping_table;
if ((dense_ir_mapping_table =
(size_t*)malloc(sizeof(size_t) * mesh[0] * mesh[1] * mesh[2])) == NULL) {
warning_print("spglib: Memory of unique_rot could not be allocated.");
return 0;
}
num_ir = kpt_get_dense_irreducible_reciprocal_mesh(grid_address,
dense_ir_mapping_table,
mesh,
is_shift,
rot_reciprocal);
for (i = 0; i < mesh[0] * mesh[1] * mesh[2]; i++) {
ir_mapping_table[i] = dense_ir_mapping_table[i];
}
free(dense_ir_mapping_table);
dense_ir_mapping_table = NULL;
return num_ir;
}
size_t kpt_get_dense_irreducible_reciprocal_mesh(int grid_address[][3],
size_t ir_mapping_table[],
const int mesh[3],
const int is_shift[3],
const MatINT *rot_reciprocal)
{
size_t num_ir;
long num_ir;
num_ir = get_dense_ir_reciprocal_mesh(grid_address,
ir_mapping_table,
@ -271,181 +91,6 @@ size_t kpt_get_dense_irreducible_reciprocal_mesh(int grid_address[][3],
return num_ir;
}
int kpt_get_stabilized_reciprocal_mesh(int grid_address[][3],
int ir_mapping_table[],
const int mesh[3],
const int is_shift[3],
const int is_time_reversal,
const MatINT * rotations,
const size_t num_q,
SPGCONST double qpoints[][3])
{
int num_ir;
size_t i;
size_t *dense_ir_mapping_table;
if ((dense_ir_mapping_table =
(size_t*)malloc(sizeof(size_t) * mesh[0] * mesh[1] * mesh[2])) == NULL) {
warning_print("spglib: Memory of unique_rot could not be allocated.");
return 0;
}
num_ir = kpt_get_dense_stabilized_reciprocal_mesh(grid_address,
dense_ir_mapping_table,
mesh,
is_shift,
is_time_reversal,
rotations,
num_q,
qpoints);
for (i = 0; i < mesh[0] * mesh[1] * mesh[2]; i++) {
ir_mapping_table[i] = dense_ir_mapping_table[i];
}
free(dense_ir_mapping_table);
dense_ir_mapping_table = NULL;
return num_ir;
}
size_t kpt_get_dense_stabilized_reciprocal_mesh(int grid_address[][3],
size_t ir_mapping_table[],
const int mesh[3],
const int is_shift[3],
const int is_time_reversal,
const MatINT * rotations,
const size_t num_q,
SPGCONST double qpoints[][3])
{
size_t num_ir;
MatINT *rot_reciprocal, *rot_reciprocal_q;
double tolerance;
rot_reciprocal = NULL;
rot_reciprocal_q = NULL;
rot_reciprocal = get_point_group_reciprocal(rotations, is_time_reversal);
tolerance = 0.01 / (mesh[0] + mesh[1] + mesh[2]);
rot_reciprocal_q = get_point_group_reciprocal_with_q(rot_reciprocal,
tolerance,
num_q,
qpoints);
num_ir = get_dense_ir_reciprocal_mesh(grid_address,
ir_mapping_table,
mesh,
is_shift,
rot_reciprocal_q);
mat_free_MatINT(rot_reciprocal_q);
rot_reciprocal_q = NULL;
mat_free_MatINT(rot_reciprocal);
rot_reciprocal = NULL;
return num_ir;
}
void
kpt_get_dense_grid_points_by_rotations(size_t rot_grid_points[],
const int address_orig[3],
SPGCONST int (*rot_reciprocal)[3][3],
const int num_rot,
const int mesh[3],
const int is_shift[3])
{
int i;
int address_double_orig[3], address_double[3];
for (i = 0; i < 3; i++) {
address_double_orig[i] = address_orig[i] * 2 + is_shift[i];
}
for (i = 0; i < num_rot; i++) {
mat_multiply_matrix_vector_i3(address_double,
rot_reciprocal[i],
address_double_orig);
rot_grid_points[i] = kgd_get_dense_grid_point_double_mesh(address_double, mesh);
}
}
void
kpt_get_dense_BZ_grid_points_by_rotations(size_t rot_grid_points[],
const int address_orig[3],
SPGCONST int (*rot_reciprocal)[3][3],
const int num_rot,
const int mesh[3],
const int is_shift[3],
const size_t bz_map[])
{
int i;
int address_double_orig[3], address_double[3], bzmesh[3];
for (i = 0; i < 3; i++) {
bzmesh[i] = mesh[i] * 2;
address_double_orig[i] = address_orig[i] * 2 + is_shift[i];
}
for (i = 0; i < num_rot; i++) {
mat_multiply_matrix_vector_i3(address_double,
rot_reciprocal[i],
address_double_orig);
rot_grid_points[i] =
bz_map[kgd_get_dense_grid_point_double_mesh(address_double, bzmesh)];
}
}
int kpt_relocate_BZ_grid_address(int bz_grid_address[][3],
int bz_map[],
SPGCONST int grid_address[][3],
const int mesh[3],
SPGCONST double rec_lattice[3][3],
const int is_shift[3])
{
int i, num_bz_map, num_bzgp;
size_t *dense_bz_map;
num_bz_map = mesh[0] * mesh[1] * mesh[2] * 8;
if ((dense_bz_map =
(size_t*)malloc(sizeof(size_t) * num_bz_map)) == NULL) {
warning_print("spglib: Memory of unique_rot could not be allocated.");
return 0;
}
num_bzgp = kpt_relocate_dense_BZ_grid_address(bz_grid_address,
dense_bz_map,
grid_address,
mesh,
rec_lattice,
is_shift);
for (i = 0; i < num_bz_map; i++) {
if (dense_bz_map[i] == num_bz_map) {
bz_map[i] = -1;
} else {
bz_map[i] = dense_bz_map[i];
}
}
free(dense_bz_map);
dense_bz_map = NULL;
return num_bzgp;
}
size_t kpt_relocate_dense_BZ_grid_address(int bz_grid_address[][3],
size_t bz_map[],
SPGCONST int grid_address[][3],
const int mesh[3],
SPGCONST double rec_lattice[3][3],
const int is_shift[3])
{
return relocate_dense_BZ_grid_address(bz_grid_address,
bz_map,
grid_address,
mesh,
rec_lattice,
is_shift);
}
MatINT *kpt_get_point_group_reciprocal(const MatINT * rotations,
const int is_time_reversal)
{
@ -454,7 +99,7 @@ MatINT *kpt_get_point_group_reciprocal(const MatINT * rotations,
MatINT *kpt_get_point_group_reciprocal_with_q(const MatINT * rot_reciprocal,
const double symprec,
const size_t num_q,
const long num_q,
SPGCONST double qpoints[][3])
{
return get_point_group_reciprocal_with_q(rot_reciprocal,
@ -542,7 +187,7 @@ static MatINT *get_point_group_reciprocal(const MatINT * rotations,
/* Return NULL if failed */
static MatINT *get_point_group_reciprocal_with_q(const MatINT * rot_reciprocal,
const double symprec,
const size_t num_q,
const long num_q,
SPGCONST double qpoints[][3])
{
int i, j, k, l, is_all_ok, num_rot;
@ -608,11 +253,11 @@ static MatINT *get_point_group_reciprocal_with_q(const MatINT * rot_reciprocal,
return rot_reciprocal_q;
}
static size_t get_dense_ir_reciprocal_mesh(int grid_address[][3],
size_t ir_mapping_table[],
const int mesh[3],
const int is_shift[3],
const MatINT *rot_reciprocal)
static long get_dense_ir_reciprocal_mesh(int grid_address[][3],
long ir_mapping_table[],
const int mesh[3],
const int is_shift[3],
const MatINT *rot_reciprocal)
{
if (check_mesh_symmetry(mesh, is_shift, rot_reciprocal)) {
return get_dense_ir_reciprocal_mesh_normal(grid_address,
@ -629,11 +274,11 @@ static size_t get_dense_ir_reciprocal_mesh(int grid_address[][3],
}
}
static size_t get_dense_ir_reciprocal_mesh_normal(int grid_address[][3],
size_t ir_mapping_table[],
const int mesh[3],
const int is_shift[3],
const MatINT *rot_reciprocal)
static long get_dense_ir_reciprocal_mesh_normal(int grid_address[][3],
long ir_mapping_table[],
const int mesh[3],
const int is_shift[3],
const MatINT *rot_reciprocal)
{
/* In the following loop, mesh is doubled. */
/* Even and odd mesh numbers correspond to */
@ -642,14 +287,14 @@ static size_t get_dense_ir_reciprocal_mesh_normal(int grid_address[][3],
/* grid: reducible grid points */
/* ir_mapping_table: the mapping from each point to ir-point. */
size_t i, grid_point_rot;
long i, grid_point_rot;
int j;
int address_double[3], address_double_rot[3];
kgd_get_all_grid_addresses(grid_address, mesh);
#pragma omp parallel for private(j, grid_point_rot, address_double, address_double_rot)
for (i = 0; i < mesh[0] * mesh[1] * (size_t)(mesh[2]); i++) {
for (i = 0; i < mesh[0] * mesh[1] * (long)(mesh[2]); i++) {
kgd_get_grid_address_double_mesh(address_double,
grid_address[i],
mesh,
@ -674,14 +319,14 @@ static size_t get_dense_ir_reciprocal_mesh_normal(int grid_address[][3],
return get_dense_num_ir(ir_mapping_table, mesh);
}
static size_t
static long
get_dense_ir_reciprocal_mesh_distortion(int grid_address[][3],
size_t ir_mapping_table[],
long ir_mapping_table[],
const int mesh[3],
const int is_shift[3],
const MatINT *rot_reciprocal)
{
size_t i, grid_point_rot;
long i, grid_point_rot;
int j, k, indivisible;
int address_double[3], address_double_rot[3];
long long_address_double[3], long_address_double_rot[3], divisor[3];
@ -696,7 +341,7 @@ get_dense_ir_reciprocal_mesh_distortion(int grid_address[][3],
}
#pragma omp parallel for private(j, k, grid_point_rot, address_double, address_double_rot, long_address_double, long_address_double_rot)
for (i = 0; i < mesh[0] * mesh[1] * (size_t)(mesh[2]); i++) {
for (i = 0; i < mesh[0] * mesh[1] * (long)(mesh[2]); i++) {
kgd_get_grid_address_double_mesh(address_double,
grid_address[i],
mesh,
@ -742,21 +387,21 @@ get_dense_ir_reciprocal_mesh_distortion(int grid_address[][3],
return get_dense_num_ir(ir_mapping_table, mesh);
}
static size_t get_dense_num_ir(size_t ir_mapping_table[], const int mesh[3])
static long get_dense_num_ir(long ir_mapping_table[], const int mesh[3])
{
size_t i, num_ir;
long i, num_ir;
num_ir = 0;
#pragma omp parallel for reduction(+:num_ir)
for (i = 0; i < mesh[0] * mesh[1] * (size_t)(mesh[2]); i++) {
for (i = 0; i < mesh[0] * mesh[1] * (long)(mesh[2]); i++) {
if (ir_mapping_table[i] == i) {
num_ir++;
}
}
#ifdef _OPENMP
for (i = 0; i < mesh[0] * mesh[1] * (size_t)(mesh[2]); i++) {
for (i = 0; i < mesh[0] * mesh[1] * (long)(mesh[2]); i++) {
ir_mapping_table[i] = ir_mapping_table[ir_mapping_table[i]];
}
#endif
@ -764,103 +409,6 @@ static size_t get_dense_num_ir(size_t ir_mapping_table[], const int mesh[3])
return num_ir;
}
static size_t relocate_dense_BZ_grid_address(int bz_grid_address[][3],
size_t bz_map[],
SPGCONST int grid_address[][3],
const int mesh[3],
SPGCONST double rec_lattice[3][3],
const int is_shift[3])
{
double tolerance, min_distance;
double q_vector[3], distance[KPT_NUM_BZ_SEARCH_SPACE];
int bzmesh[3], bz_address_double[3];
size_t i, boundary_num_gp, total_num_gp, bzgp, gp, num_bzmesh;
int j, k, min_index;
tolerance = get_tolerance_for_BZ_reduction(rec_lattice, mesh);
for (j = 0; j < 3; j++) {
bzmesh[j] = mesh[j] * 2;
}
num_bzmesh = bzmesh[0] * bzmesh[1] * (size_t)(bzmesh[2]);
for (i = 0; i < num_bzmesh; i++) {
bz_map[i] = num_bzmesh;
}
boundary_num_gp = 0;
total_num_gp = mesh[0] * mesh[1] * (size_t)(mesh[2]);
/* Multithreading doesn't work for this loop since gp calculated */
/* with boundary_num_gp is unstable to store bz_grid_address. */
for (i = 0; i < total_num_gp; i++) {
for (j = 0; j < KPT_NUM_BZ_SEARCH_SPACE; j++) {
for (k = 0; k < 3; k++) {
q_vector[k] =
((grid_address[i][k] + bz_search_space[j][k] * mesh[k]) * 2 +
is_shift[k]) / ((double)mesh[k]) / 2;
}
mat_multiply_matrix_vector_d3(q_vector, rec_lattice, q_vector);
distance[j] = mat_norm_squared_d3(q_vector);
}
min_distance = distance[0];
min_index = 0;
for (j = 1; j < KPT_NUM_BZ_SEARCH_SPACE; j++) {
if (distance[j] < min_distance) {
min_distance = distance[j];
min_index = j;
}
}
for (j = 0; j < KPT_NUM_BZ_SEARCH_SPACE; j++) {
if (distance[j] < min_distance + tolerance) {
if (j == min_index) {
gp = i;
} else {
gp = boundary_num_gp + total_num_gp;
}
for (k = 0; k < 3; k++) {
bz_grid_address[gp][k] =
grid_address[i][k] + bz_search_space[j][k] * mesh[k];
bz_address_double[k] = bz_grid_address[gp][k] * 2 + is_shift[k];
}
bzgp = kgd_get_dense_grid_point_double_mesh(bz_address_double, bzmesh);
bz_map[bzgp] = gp;
if (j != min_index) {
boundary_num_gp++;
}
}
}
}
return boundary_num_gp + total_num_gp;
}
static double get_tolerance_for_BZ_reduction(SPGCONST double rec_lattice[3][3],
const int mesh[3])
{
int i, j;
double tolerance;
double length[3];
for (i = 0; i < 3; i++) {
length[i] = 0;
for (j = 0; j < 3; j++) {
length[i] += rec_lattice[j][i] * rec_lattice[j][i];
}
length[i] /= mesh[i] * mesh[i];
}
tolerance = length[0];
for (i = 1; i < 3; i++) {
if (tolerance < length[i]) {
tolerance = length[i];
}
}
tolerance *= 0.01;
return tolerance;
}
static int check_mesh_symmetry(const int mesh[3],
const int is_shift[3],
const MatINT *rot_reciprocal)

60
c/kpoint.h
View File

@ -38,64 +38,16 @@
#include <stddef.h>
#include "mathfunc.h"
int kpt_get_irreducible_reciprocal_mesh(int grid_address[][3],
int ir_mapping_table[],
const int mesh[3],
const int is_shift[3],
const MatINT *rot_reciprocal);
size_t kpt_get_dense_irreducible_reciprocal_mesh(int grid_address[][3],
size_t ir_mapping_table[],
const int mesh[3],
const int is_shift[3],
const MatINT *rot_reciprocal);
int kpt_get_stabilized_reciprocal_mesh(int grid_address[][3],
int ir_mapping_table[],
const int mesh[3],
const int is_shift[3],
const int is_time_reversal,
const MatINT * rotations,
const size_t num_q,
SPGCONST double qpoints[][3]);
size_t kpt_get_dense_stabilized_reciprocal_mesh(int grid_address[][3],
size_t ir_mapping_table[],
const int mesh[3],
const int is_shift[3],
const int is_time_reversal,
const MatINT * rotations,
const size_t num_q,
SPGCONST double qpoints[][3]);
void
kpt_get_dense_grid_points_by_rotations(size_t rot_grid_points[],
const int address_orig[3],
SPGCONST int (*rot_reciprocal)[3][3],
const int num_rot,
const int mesh[3],
const int is_shift[3]);
void
kpt_get_dense_BZ_grid_points_by_rotations(size_t rot_grid_points[],
const int address_orig[3],
SPGCONST int (*rot_reciprocal)[3][3],
const int num_rot,
const int mesh[3],
const int is_shift[3],
const size_t bz_map[]);
int kpt_relocate_BZ_grid_address(int bz_grid_address[][3],
int bz_map[],
SPGCONST int grid_address[][3],
const int mesh[3],
SPGCONST double rec_lattice[3][3],
const int is_shift[3]);
size_t kpt_relocate_dense_BZ_grid_address(int bz_grid_address[][3],
size_t bz_map[],
SPGCONST int grid_address[][3],
const int mesh[3],
SPGCONST double rec_lattice[3][3],
const int is_shift[3]);
long kpt_get_dense_irreducible_reciprocal_mesh(int grid_address[][3],
long ir_mapping_table[],
const int mesh[3],
const int is_shift[3],
const MatINT *rot_reciprocal);
MatINT *kpt_get_point_group_reciprocal(const MatINT * rotations,
const int is_time_reversal);
MatINT *kpt_get_point_group_reciprocal_with_q(const MatINT * rot_reciprocal,
const double symprec,
const size_t num_q,
const long num_q,
SPGCONST double qpoints[][3]);
#endif

174
c/phono3py.c
View File

@ -45,15 +45,19 @@
#include "fc3.h"
#include "tetrahedron_method.h"
#include "triplet.h"
#include "triplet_iw.h"
#include "kgrid.h"
#include <stdio.h>
void ph3py_get_interaction(Darray *fc3_normal_squared,
const char *g_zero,
const Darray *frequencies,
const lapack_complex_double *eigenvectors,
const size_t (*triplets)[3],
const size_t num_triplets,
const long (*triplets)[3],
const long num_triplets,
const int *grid_address,
const int *mesh,
const double *fc3,
@ -94,11 +98,11 @@ void ph3py_get_pp_collision(double *imag_self_energy,
PHPYCONST int relative_grid_address[24][4][3], /* thm */
const double *frequencies,
const lapack_complex_double *eigenvectors,
const size_t (*triplets)[3],
const size_t num_triplets,
const int *triplet_weights,
const long (*triplets)[3],
const long num_triplets,
const long *triplet_weights,
const int *grid_address, /* thm */
const size_t *bz_map, /* thm */
const long *bz_map, /* thm */
const int *mesh, /* thm */
const double *fc3,
const int is_compact_fc3,
@ -146,9 +150,9 @@ void ph3py_get_pp_collision_with_sigma(
const double sigma_cutoff,
const double *frequencies,
const lapack_complex_double *eigenvectors,
const size_t (*triplets)[3],
const size_t num_triplets,
const int *triplet_weights,
const long (*triplets)[3],
const long num_triplets,
const long *triplet_weights,
const int *grid_address,
const int *mesh,
const double *fc3,
@ -195,8 +199,8 @@ void ph3py_get_imag_self_energy_at_bands_with_g(
double *imag_self_energy,
const Darray *fc3_normal_squared,
const double *frequencies,
const size_t (*triplets)[3],
const int *triplet_weights,
const long (*triplets)[3],
const long *triplet_weights,
const double *g,
const char *g_zero,
const double temperature,
@ -224,8 +228,8 @@ void ph3py_get_detailed_imag_self_energy_at_bands_with_g(
double *imag_self_energy_U,
const Darray *fc3_normal_squared,
const double *frequencies,
const size_t (*triplets)[3],
const int *triplet_weights,
const long (*triplets)[3],
const long *triplet_weights,
const int *grid_address,
const double *g,
const char *g_zero,
@ -251,8 +255,8 @@ void ph3py_get_real_self_energy_at_bands(double *real_self_energy,
const Darray *fc3_normal_squared,
const int *band_indices,
const double *frequencies,
const size_t (*triplets)[3],
const int *triplet_weights,
const long (*triplets)[3],
const long *triplet_weights,
const double epsilon,
const double temperature,
const double unit_conversion_factor,
@ -277,8 +281,8 @@ void ph3py_get_real_self_energy_at_frequency_point(
const Darray *fc3_normal_squared,
const int *band_indices,
const double *frequencies,
const size_t (*triplets)[3],
const int *triplet_weights,
const long (*triplets)[3],
const long *triplet_weights,
const double epsilon,
const double temperature,
const double unit_conversion_factor,
@ -301,15 +305,15 @@ void ph3py_get_real_self_energy_at_frequency_point(
void ph3py_get_collision_matrix(double *collision_matrix,
const Darray *fc3_normal_squared,
const double *frequencies,
const size_t (*triplets)[3],
const size_t *triplets_map,
const size_t *map_q,
const size_t *rotated_grid_points,
const long (*triplets)[3],
const long *triplets_map,
const long *map_q,
const long *rotated_grid_points,
const double *rotations_cartesian,
const double *g,
const size_t num_ir_gp,
const size_t num_gp,
const size_t num_rot,
const long num_ir_gp,
const long num_gp,
const long num_rot,
const double temperature,
const double unit_conversion_factor,
const double cutoff_frequency)
@ -335,11 +339,11 @@ void ph3py_get_collision_matrix(double *collision_matrix,
void ph3py_get_reducible_collision_matrix(double *collision_matrix,
const Darray *fc3_normal_squared,
const double *frequencies,
const size_t (*triplets)[3],
const size_t *triplets_map,
const size_t *map_q,
const long (*triplets)[3],
const long *triplets_map,
const long *map_q,
const double *g,
const size_t num_gp,
const long num_gp,
const double temperature,
const double unit_conversion_factor,
const double cutoff_frequency)
@ -360,14 +364,14 @@ void ph3py_get_reducible_collision_matrix(double *collision_matrix,
void ph3py_get_isotope_scattering_strength(
double *gamma,
const size_t grid_point,
const long grid_point,
const double *mass_variances,
const double *frequencies,
const lapack_complex_double *eigenvectors,
const size_t num_grid_points,
const long num_grid_points,
const int *band_indices,
const size_t num_band,
const size_t num_band0,
const long num_band,
const long num_band0,
const double sigma,
const double cutoff_frequency)
{
@ -387,16 +391,16 @@ void ph3py_get_isotope_scattering_strength(
void ph3py_get_thm_isotope_scattering_strength
(double *gamma,
const size_t grid_point,
const size_t *ir_grid_points,
const int *weights,
const long grid_point,
const long *ir_grid_points,
const long *weights,
const double *mass_variances,
const double *frequencies,
const lapack_complex_double *eigenvectors,
const size_t num_ir_grid_points,
const long num_ir_grid_points,
const int *band_indices,
const size_t num_band,
const size_t num_band0,
const long num_band,
const long num_band0,
const double *integration_weights,
const double cutoff_frequency)
{
@ -419,7 +423,7 @@ void ph3py_distribute_fc3(double *fc3,
const int target,
const int source,
const int *atom_mapping,
const size_t num_atom,
const long num_atom,
const double *rot_cart)
{
fc3_distribute_fc3(fc3,
@ -436,9 +440,9 @@ void ph3py_rotate_delta_fc2(double (*fc3)[3][3][3],
const double *inv_U,
PHPYCONST double (*site_sym_cart)[3][3],
const int *rot_map_syms,
const size_t num_atom,
const size_t num_site_sym,
const size_t num_disp)
const long num_atom,
const long num_site_sym,
const long num_disp)
{
fc3_rotate_delta_fc2(fc3,
delta_fc2s,
@ -451,7 +455,7 @@ void ph3py_rotate_delta_fc2(double (*fc3)[3][3][3],
}
void ph3py_set_permutation_symmetry_fc3(double *fc3, const size_t num_atom)
void ph3py_set_permutation_symmetry_fc3(double *fc3, const long num_atom)
{
fc3_set_permutation_symmetry_fc3(fc3, num_atom);
}
@ -462,8 +466,8 @@ void ph3py_set_permutation_symmetry_compact_fc3(double * fc3,
const int s2pp[],
const int nsym_list[],
const int perms[],
const size_t n_satom,
const size_t n_patom)
const long n_satom,
const long n_patom)
{
fc3_set_permutation_symmetry_compact_fc3(fc3,
p2s,
@ -479,8 +483,8 @@ void ph3py_transpose_compact_fc3(double * fc3,
const int s2pp[],
const int nsym_list[],
const int perms[],
const size_t n_satom,
const size_t n_patom,
const long n_satom,
const long n_patom,
const int t_type)
{
fc3_transpose_compact_fc3(fc3,
@ -494,15 +498,15 @@ void ph3py_transpose_compact_fc3(double * fc3,
}
size_t ph3py_get_triplets_reciprocal_mesh_at_q(size_t *map_triplets,
size_t *map_q,
int (*grid_address)[3],
const size_t grid_point,
const int mesh[3],
const int is_time_reversal,
const int num_rot,
PHPYCONST int (*rotations)[3][3],
const int swappable)
long ph3py_get_triplets_reciprocal_mesh_at_q(long *map_triplets,
long *map_q,
int (*grid_address)[3],
const long grid_point,
const int mesh[3],
const int is_time_reversal,
const long num_rot,
PHPYCONST int (*rotations)[3][3],
const int swappable)
{
return tpl_get_triplets_reciprocal_mesh_at_q(map_triplets,
map_q,
@ -516,13 +520,13 @@ size_t ph3py_get_triplets_reciprocal_mesh_at_q(size_t *map_triplets,
}
size_t ph3py_get_BZ_triplets_at_q(size_t (*triplets)[3],
const size_t grid_point,
PHPYCONST int (*bz_grid_address)[3],
const size_t *bz_map,
const size_t *map_triplets,
const size_t num_map_triplets,
const int mesh[3])
long ph3py_get_BZ_triplets_at_q(long (*triplets)[3],
const long grid_point,
PHPYCONST int (*bz_grid_address)[3],
const long *bz_map,
const long *map_triplets,
const long num_map_triplets,
const int mesh[3])
{
return tpl_get_BZ_triplets_at_q(triplets,
grid_point,
@ -537,18 +541,18 @@ size_t ph3py_get_BZ_triplets_at_q(size_t (*triplets)[3],
void ph3py_get_integration_weight(double *iw,
char *iw_zero,
const double *frequency_points,
const size_t num_band0,
const long num_band0,
PHPYCONST int relative_grid_address[24][4][3],
const int mesh[3],
PHPYCONST size_t (*triplets)[3],
const size_t num_triplets,
PHPYCONST long (*triplets)[3],
const long num_triplets,
PHPYCONST int (*bz_grid_address)[3],
const size_t *bz_map,
const long *bz_map,
const double *frequencies1,
const size_t num_band1,
const long num_band1,
const double *frequencies2,
const size_t num_band2,
const size_t tp_type,
const long num_band2,
const long tp_type,
const int openmp_per_triplets,
const int openmp_per_bands)
{
@ -577,12 +581,12 @@ void ph3py_get_integration_weight_with_sigma(double *iw,
const double sigma,
const double sigma_cutoff,
const double *frequency_points,
const size_t num_band0,
PHPYCONST size_t (*triplets)[3],
const size_t num_triplets,
const long num_band0,
PHPYCONST long (*triplets)[3],
const long num_triplets,
const double *frequencies,
const size_t num_band,
const size_t tp_type)
const long num_band,
const long tp_type)
{
tpl_get_integration_weight_with_sigma(iw,
iw_zero,
@ -625,8 +629,8 @@ void ph3py_symmetrize_collision_matrix(double *collision_matrix,
void ph3py_expand_collision_matrix(double *collision_matrix,
const size_t *rot_grid_points,
const size_t *ir_grid_points,
const long *rot_grid_points,
const long *ir_grid_points,
const long num_ir_gp,
const long num_grid_points,
const long num_rot,
@ -694,12 +698,12 @@ void ph3py_expand_collision_matrix(double *collision_matrix,
}
void ph3py_get_neighboring_gird_points(size_t *relative_grid_points,
const size_t *grid_points,
void ph3py_get_neighboring_gird_points(long *relative_grid_points,
const long *grid_points,
PHPYCONST int (*relative_grid_address)[3],
const int mesh[3],
PHPYCONST int (*bz_grid_address)[3],
const size_t *bz_map,
const long *bz_map,
const long num_grid_points,
const long num_relative_grid_address)
{
@ -707,7 +711,7 @@ void ph3py_get_neighboring_gird_points(size_t *relative_grid_points,
#pragma omp parallel for
for (i = 0; i < num_grid_points; i++) {
thm_get_dense_neighboring_grid_points
tpi_get_dense_neighboring_grid_points
(relative_grid_points + i * num_relative_grid_address,
grid_points[i],
relative_grid_address,
@ -726,19 +730,19 @@ void ph3py_set_integration_weights(double *iw,
const long num_gp,
PHPYCONST int (*relative_grid_address)[4][3],
const int mesh[3],
const size_t *grid_points,
const long *grid_points,
PHPYCONST int (*bz_grid_address)[3],
const size_t *bz_map,
const long *bz_map,
const double *frequencies)
{
long i, j, k, bi;
size_t vertices[24][4];
long vertices[24][4];
double freq_vertices[24][4];
#pragma omp parallel for private(j, k, bi, vertices, freq_vertices)
for (i = 0; i < num_gp; i++) {
for (j = 0; j < 24; j++) {
thm_get_dense_neighboring_grid_points(vertices[j],
tpi_get_dense_neighboring_grid_points(vertices[j],
grid_points[i],
relative_grid_address[j],
4,

164
c/phono3py.h
View File

@ -46,8 +46,8 @@ void ph3py_get_interaction(Darray *fc3_normal_squared,
const char *g_zero,
const Darray *frequencies,
const lapack_complex_double *eigenvectors,
const size_t (*triplets)[3],
const size_t num_triplets,
const long (*triplets)[3],
const long num_triplets,
const int *grid_address,
const int *mesh,
const double *fc3,
@ -65,11 +65,11 @@ void ph3py_get_pp_collision(double *imag_self_energy,
PHPYCONST int relative_grid_address[24][4][3], /* thm */
const double *frequencies,
const lapack_complex_double *eigenvectors,
const size_t (*triplets)[3],
const size_t num_triplets,
const int *triplet_weights,
const long (*triplets)[3],
const long num_triplets,
const long *triplet_weights,
const int *grid_address, /* thm */
const size_t *bz_map, /* thm */
const long *bz_map, /* thm */
const int *mesh, /* thm */
const double *fc3,
const int is_compact_fc3,
@ -90,9 +90,9 @@ void ph3py_get_pp_collision_with_sigma(
const double sigma_cutoff,
const double *frequencies,
const lapack_complex_double *eigenvectors,
const size_t (*triplets)[3],
const size_t num_triplets,
const int *triplet_weights,
const long (*triplets)[3],
const long num_triplets,
const long *triplet_weights,
const int *grid_address,
const int *mesh,
const double *fc3,
@ -112,8 +112,8 @@ void ph3py_get_imag_self_energy_at_bands_with_g(
double *imag_self_energy,
const Darray *fc3_normal_squared,
const double *frequencies,
const size_t (*triplets)[3],
const int *triplet_weights,
const long (*triplets)[3],
const long *triplet_weights,
const double *g,
const char *g_zero,
const double temperature,
@ -126,8 +126,8 @@ void ph3py_get_detailed_imag_self_energy_at_bands_with_g(
double *imag_self_energy_U,
const Darray *fc3_normal_squared,
const double *frequencies,
const size_t (*triplets)[3],
const int *triplet_weights,
const long (*triplets)[3],
const long *triplet_weights,
const int *grid_address,
const double *g,
const char *g_zero,
@ -137,8 +137,8 @@ void ph3py_get_real_self_energy_at_bands(double *real_self_energy,
const Darray *fc3_normal_squared,
const int *band_indices,
const double *frequencies,
const size_t (*triplets)[3],
const int *triplet_weights,
const long (*triplets)[3],
const long *triplet_weights,
const double epsilon,
const double temperature,
const double unit_conversion_factor,
@ -149,8 +149,8 @@ void ph3py_get_real_self_energy_at_frequency_point(
const Darray *fc3_normal_squared,
const int *band_indices,
const double *frequencies,
const size_t (*triplets)[3],
const int *triplet_weights,
const long (*triplets)[3],
const long *triplet_weights,
const double epsilon,
const double temperature,
const double unit_conversion_factor,
@ -158,116 +158,116 @@ void ph3py_get_real_self_energy_at_frequency_point(
void ph3py_get_collision_matrix(double *collision_matrix,
const Darray *fc3_normal_squared,
const double *frequencies,
const size_t (*triplets)[3],
const size_t *triplets_map,
const size_t *map_q,
const size_t *rotated_grid_points,
const long (*triplets)[3],
const long *triplets_map,
const long *map_q,
const long *rotated_grid_points,
const double *rotations_cartesian,
const double *g,
const size_t num_ir_gp,
const size_t num_gp,
const size_t num_rot,
const long num_ir_gp,
const long num_gp,
const long num_rot,
const double temperature,
const double unit_conversion_factor,
const double cutoff_frequency);
void ph3py_get_reducible_collision_matrix(double *collision_matrix,
const Darray *fc3_normal_squared,
const double *frequencies,
const size_t (*triplets)[3],
const size_t *triplets_map,
const size_t *map_q,
const long (*triplets)[3],
const long *triplets_map,
const long *map_q,
const double *g,
const size_t num_gp,
const long num_gp,
const double temperature,
const double unit_conversion_factor,
const double cutoff_frequency);
void ph3py_get_isotope_scattering_strength(
double *gamma,
const size_t grid_point,
const long grid_point,
const double *mass_variances,
const double *frequencies,
const lapack_complex_double *eigenvectors,
const size_t num_grid_points,
const long num_grid_points,
const int *band_indices,
const size_t num_band,
const size_t num_band0,
const long num_band,
const long num_band0,
const double sigma,
const double cutoff_frequency);
void ph3py_get_thm_isotope_scattering_strength(
double *gamma,
const size_t grid_point,
const size_t *ir_grid_points,
const int *weights,
const long grid_point,
const long *ir_grid_points,
const long *weights,
const double *mass_variances,
const double *frequencies,
const lapack_complex_double *eigenvectors,
const size_t num_ir_grid_points,
const long num_ir_grid_points,
const int *band_indices,
const size_t num_band,
const size_t num_band0,
const long num_band,
const long num_band0,
const double *integration_weights,
const double cutoff_frequency);
void ph3py_distribute_fc3(double *fc3,
const int target,
const int source,
const int *atom_mapping,
const size_t num_atom,
const long num_atom,
const double *rot_cart);
void ph3py_rotate_delta_fc2(double (*fc3)[3][3][3],
PHPYCONST double (*delta_fc2s)[3][3],
const double *inv_U,
PHPYCONST double (*site_sym_cart)[3][3],
const int *rot_map_syms,
const size_t num_atom,
const size_t num_site_sym,
const size_t num_disp);
void ph3py_set_permutation_symmetry_fc3(double *fc3, const size_t num_atom);
const long num_atom,
const long num_site_sym,
const long num_disp);
void ph3py_set_permutation_symmetry_fc3(double *fc3, const long num_atom);
void ph3py_set_permutation_symmetry_compact_fc3(double * fc3,
const int p2s[],
const int s2pp[],
const int nsym_list[],
const int perms[],
const size_t n_satom,
const size_t n_patom);
const long n_satom,
const long n_patom);
void ph3py_transpose_compact_fc3(double * fc3,
const int p2s[],
const int s2pp[],
const int nsym_list[],
const int perms[],
const size_t n_satom,
const size_t n_patom,
const long n_satom,
const long n_patom,
const int t_type);
size_t ph3py_get_triplets_reciprocal_mesh_at_q(size_t *map_triplets,
size_t *map_q,
int (*grid_address)[3],
const size_t grid_point,
const int mesh[3],
const int is_time_reversal,
const int num_rot,
PHPYCONST int (*rotations)[3][3],
const int swappable);
size_t ph3py_get_BZ_triplets_at_q(size_t (*triplets)[3],
const size_t grid_point,
PHPYCONST int (*bz_grid_address)[3],
const size_t *bz_map,
const size_t *map_triplets,
const size_t num_map_triplets,
const int mesh[3]);
long ph3py_get_triplets_reciprocal_mesh_at_q(long *map_triplets,
long *map_q,
int (*grid_address)[3],
const long grid_point,
const int mesh[3],
const int is_time_reversal,
const long num_rot,
PHPYCONST int (*rotations)[3][3],
const int swappable);
long ph3py_get_BZ_triplets_at_q(long (*triplets)[3],
const long grid_point,
PHPYCONST int (*bz_grid_address)[3],
const long *bz_map,
const long *map_triplets,
const long num_map_triplets,
const int mesh[3]);
void ph3py_get_integration_weight(double *iw,
char *iw_zero,
const double *frequency_points,
const size_t num_band0,
const long num_band0,
PHPYCONST int relative_grid_address[24][4][3],
const int mesh[3],
PHPYCONST size_t (*triplets)[3],
const size_t num_triplets,
PHPYCONST long (*triplets)[3],
const long num_triplets,
PHPYCONST int (*bz_grid_address)[3],
const size_t *bz_map,
const long *bz_map,
const double *frequencies1,
const size_t num_band1,
const long num_band1,
const double *frequencies2,
const size_t num_band2,
const size_t tp_type,
const long num_band2,
const long tp_type,
const int openmp_per_triplets,
const int openmp_per_bands);
void ph3py_get_integration_weight_with_sigma(double *iw,
@ -275,12 +275,12 @@ void ph3py_get_integration_weight_with_sigma(double *iw,
const double sigma,
const double sigma_cutoff,
const double *frequency_points,
const size_t num_band0,
PHPYCONST size_t (*triplets)[3],
const size_t num_triplets,
const long num_band0,
PHPYCONST long (*triplets)[3],
const long num_triplets,
const double *frequencies,
const size_t num_band,
const size_t tp_type);
const long num_band,
const long tp_type);
void ph3py_symmetrize_collision_matrix(double *collision_matrix,
@ -288,20 +288,20 @@ void ph3py_symmetrize_collision_matrix(double *collision_matrix,
const long num_temp,
const long num_sigma);
void ph3py_expand_collision_matrix(double *collision_matrix,
const size_t *rot_grid_points,
const size_t *ir_grid_points,
const long *rot_grid_points,
const long *ir_grid_points,
const long num_ir_gp,
const long num_grid_points,
const long num_rot,
const long num_sigma,
const long num_temp,
const long num_band);
void ph3py_get_neighboring_gird_points(size_t *relative_grid_points,
const size_t *grid_points,
void ph3py_get_neighboring_gird_points(long *relative_grid_points,
const long *grid_points,
PHPYCONST int (*relative_grid_address)[3],
const int mesh[3],
PHPYCONST int (*bz_grid_address)[3],
const size_t *bz_map,
const long *bz_map,
const long num_grid_points,
const long num_relative_grid_address);
void ph3py_set_integration_weights(double *iw,
@ -311,9 +311,9 @@ void ph3py_set_integration_weights(double *iw,
const long num_gp,
PHPYCONST int (*relative_grid_address)[4][3],
const int mesh[3],
const size_t *grid_points,
const long *grid_points,
PHPYCONST int (*bz_grid_address)[3],
const size_t *bz_map,
const long *bz_map,
const double *frequencies);
#endif

70
c/pp_collision.c
View File

@ -45,16 +45,16 @@
#include "lapack_wrapper.h"
static void get_collision(double *ise,
const size_t num_band0,
const size_t num_band,
const size_t num_temps,
const long num_band0,
const long num_band,
const long num_temps,
const double *temperatures,
const double *g,
const char *g_zero,
const double *frequencies,
const lapack_complex_double *eigenvectors,
const size_t triplet[3],
const int weight,
const long triplet[3],
const long triplet_weight,
const int *grid_address,
const int *mesh,
const double *fc3,
@ -72,21 +72,21 @@ static void get_collision(double *ise,
static void finalize_ise(double *imag_self_energy,
const double *ise,
const int *grid_address,
const size_t (*triplets)[3],
const size_t num_triplets,
const size_t num_temps,
const size_t num_band0,
const long (*triplets)[3],
const long num_triplets,
const long num_temps,
const long num_band0,
const int is_NU);
void ppc_get_pp_collision(double *imag_self_energy,
PHPYCONST int relative_grid_address[24][4][3], /* thm */
const double *frequencies,
const lapack_complex_double *eigenvectors,
const size_t (*triplets)[3],
const size_t num_triplets,
const int *weights,
const long (*triplets)[3],
const long num_triplets,
const long *triplet_weights,
const int *grid_address, /* thm */
const size_t *bz_map, /* thm */
const long *bz_map, /* thm */
const int *mesh, /* thm */
const double *fc3,
const int is_compact_fc3,
@ -102,8 +102,8 @@ void ppc_get_pp_collision(double *imag_self_energy,
const int symmetrize_fc3_q,
const double cutoff_frequency)
{
size_t i;
size_t num_band, num_band0, num_band_prod, num_temps;
long i;
long num_band, num_band0, num_band_prod, num_temps;
int openmp_per_triplets;
double *ise, *freqs_at_gp, *g;
char *g_zero;
@ -166,7 +166,7 @@ void ppc_get_pp_collision(double *imag_self_energy,
frequencies,
eigenvectors,
triplets[i],
weights[i],
triplet_weights[i],
grid_address,
mesh,
fc3,
@ -209,9 +209,9 @@ void ppc_get_pp_collision_with_sigma(
const double sigma_cutoff,
const double *frequencies,
const lapack_complex_double *eigenvectors,
const size_t (*triplets)[3],
const size_t num_triplets,
const int *weights,
const long (*triplets)[3],
const long num_triplets,
const long *triplet_weights,
const int *grid_address,
const int *mesh,
const double *fc3,
@ -228,8 +228,8 @@ void ppc_get_pp_collision_with_sigma(
const int symmetrize_fc3_q,
const double cutoff_frequency)
{
size_t i;
size_t num_band, num_band0, num_band_prod, num_temps;
long i;
long num_band, num_band0, num_band_prod, num_temps;
int openmp_per_triplets, const_adrs_shift;
double cutoff;
double *ise, *freqs_at_gp, *g;
@ -288,7 +288,7 @@ void ppc_get_pp_collision_with_sigma(
frequencies,
eigenvectors,
triplets[i],
weights[i],
triplet_weights[i],
grid_address,
mesh,
fc3,
@ -326,16 +326,16 @@ void ppc_get_pp_collision_with_sigma(
}
static void get_collision(double *ise,
const size_t num_band0,
const size_t num_band,
const size_t num_temps,
const long num_band0,
const long num_band,
const long num_temps,
const double *temperatures,
const double *g,
const char *g_zero,
const double *frequencies,
const lapack_complex_double *eigenvectors,
const size_t triplet[3],
const int weight,
const long triplet[3],
const long triplet_weight,
const int *grid_address,
const int *mesh,
const double *fc3,
@ -351,8 +351,8 @@ static void get_collision(double *ise,
const double cutoff_frequency,
const int openmp_per_triplets)
{
size_t i;
size_t num_band_prod, num_g_pos;
long i;
long num_band_prod, num_g_pos;
double *fc3_normal_squared;
int (*g_pos)[4];
@ -405,7 +405,7 @@ static void get_collision(double *ise,
fc3_normal_squared,
frequencies,
triplet,
weight,
triplet_weight,
g,
g + num_band_prod,
g_pos,
@ -425,13 +425,13 @@ static void get_collision(double *ise,
static void finalize_ise(double *imag_self_energy,
const double *ise,
const int *grid_address,
const size_t (*triplets)[3],
const size_t num_triplets,
const size_t num_temps,
const size_t num_band0,
const long (*triplets)[3],
const long num_triplets,
const long num_temps,
const long num_band0,
const int is_NU)
{
size_t i, j, k;
long i, j, k;
int is_N;
if (is_NU) {

14
c/pp_collision.h
View File

@ -43,11 +43,11 @@ void ppc_get_pp_collision(double *imag_self_energy,
PHPYCONST int relative_grid_address[24][4][3],
const double *frequencies,
const lapack_complex_double *eigenvectors,
const size_t (*triplets)[3],
const size_t num_triplets,
const int *weights,
const long (*triplets)[3],
const long num_triplets,
const long *triplet_weights,
const int *grid_address,
const size_t *bz_map,
const long *bz_map,
const int *mesh,
const double *fc3,
const int is_compact_fc3,
@ -69,9 +69,9 @@ void ppc_get_pp_collision_with_sigma(
const double sigma_cutoff,
const double *frequencies,
const lapack_complex_double *eigenvectors,
const size_t (*triplets)[3],
const size_t num_triplets,
const int *weights,
const long (*triplets)[3],
const long num_triplets,
const long *triplet_weights,
const int *grid_address,
const int *mesh,
const double *fc3,

16
c/real_self_energy.c
View File

@ -43,8 +43,8 @@ static double get_real_self_energy_at_band(const int band_index,
const Darray *fc3_normal_squared,
const double fpoint,
const double *frequencies,
const size_t (*triplets)[3],
const int *triplet_weights,
const long (*triplets)[3],
const long *triplet_weights,
const double epsilon,
const double temperature,
const double unit_conversion_factor,
@ -69,8 +69,8 @@ void rse_get_real_self_energy_at_bands(double *real_self_energy,
const Darray *fc3_normal_squared,
const int *band_indices,
const double *frequencies,
const size_t (*triplets)[3],
const int *triplet_weights,
const long (*triplets)[3],
const long *triplet_weights,
const double epsilon,
const double temperature,
const double unit_conversion_factor,
@ -110,8 +110,8 @@ void rse_get_real_self_energy_at_frequency_point(
const Darray *fc3_normal_squared,
const int *band_indices,
const double *frequencies,
const size_t (*triplets)[3],
const int *triplet_weights,
const long (*triplets)[3],
const long *triplet_weights,
const double epsilon,
const double temperature,
const double unit_conversion_factor,
@ -145,8 +145,8 @@ static double get_real_self_energy_at_band(const int band_index,
const Darray *fc3_normal_squared,
const double fpoint,
const double *frequencies,
const size_t (*triplets)[3],
const int *triplet_weights,
const long (*triplets)[3],
const long *triplet_weights,
const double epsilon,
const double temperature,
const double unit_conversion_factor,

8
c/real_self_energy.h
View File

@ -43,8 +43,8 @@ void rse_get_real_self_energy_at_bands(double *real_self_energy,
const Darray *fc3_normal_squared,
const int *band_indices,
const double *frequencies,
const size_t (*triplets)[3],
const int *triplet_weights,
const long (*triplets)[3],
const long *triplet_weights,
const double epsilon,
const double temperature,
const double unit_conversion_factor,
@ -55,8 +55,8 @@ void rse_get_real_self_energy_at_frequency_point(
const Darray *fc3_normal_squared,
const int *band_indices,
const double *frequencies,
const size_t (*triplets)[3],
const int *triplet_weights,
const long (*triplets)[3],
const long *triplet_weights,
const double epsilon,
const double temperature,
const double unit_conversion_factor,

28
c/real_to_reciprocal.c
View File

@ -71,15 +71,15 @@ static void real_to_reciprocal_elements(lapack_complex_double *fc3_rec_elem,
const int *multiplicity,
const int *p2s,
const int *s2p,
const size_t pi0,
const size_t pi1,
const size_t pi2);
const long pi0,
const long pi1,
const long pi2);
static lapack_complex_double get_phase_factor(const double q[],
const int qi,
const double *shortest_vectors,
const int multi);
static lapack_complex_double
get_pre_phase_factor(const size_t i,
get_pre_phase_factor(const long i,
const double q[9],
const double *shortest_vectors,
const int svecs_dims[3],
@ -133,8 +133,8 @@ real_to_reciprocal_single_thread(lapack_complex_double *fc3_reciprocal,
const int *p2s_map,
const int *s2p_map)
{
size_t i, j, k;
size_t num_patom, adrs_shift;
long i, j, k;
long num_patom, adrs_shift;
lapack_complex_double pre_phase_factor;
num_patom = svecs_dims[1];
@ -179,8 +179,8 @@ real_to_reciprocal_openmp(lapack_complex_double *fc3_reciprocal,
const int *p2s_map,
const int *s2p_map)
{
size_t i, j, k, jk;
size_t num_patom, adrs_shift;
long i, j, k, jk;
long num_patom, adrs_shift;
lapack_complex_double pre_phase_factor;
num_patom = svecs_dims[1];
@ -225,12 +225,12 @@ static void real_to_reciprocal_elements(lapack_complex_double *fc3_rec_elem,
const int *multiplicity,
const int *p2s,
const int *s2p,
const size_t pi0,
const size_t pi1,
const size_t pi2)
const long pi0,
const long pi1,
const long pi2)
{
size_t i, j, k, l;
size_t num_satom, adrs_shift, adrs_vec1, adrs_vec2 ;
long i, j, k, l;
long num_satom, adrs_shift, adrs_vec1, adrs_vec2 ;
lapack_complex_double phase_factor, phase_factor1, phase_factor2;
double fc3_rec_real[27], fc3_rec_imag[27];
@ -286,7 +286,7 @@ static void real_to_reciprocal_elements(lapack_complex_double *fc3_rec_elem,
}
static lapack_complex_double
get_pre_phase_factor(const size_t i,
get_pre_phase_factor(const long i,
const double q[9],
const double *shortest_vectors,
const int svecs_dims[3],

42
c/reciprocal_to_normal.c
View File

@ -46,20 +46,20 @@
#endif
static lapack_complex_double fc3_sum_in_reciprocal_to_normal
(const size_t bi0,
const size_t bi1,
const size_t bi2,
(const long bi0,
const long bi1,
const long bi2,
const lapack_complex_double *eigvecs0,
const lapack_complex_double *eigvecs1,
const lapack_complex_double *eigvecs2,
const lapack_complex_double *fc3_reciprocal,
const double *masses,
const size_t num_atom);
const long num_atom);
static double get_fc3_sum
(const size_t j,
const size_t k,
const size_t bi,
(const long j,
const long k,
const long bi,
const double *freqs0,
const double *freqs1,
const double *freqs2,
@ -68,13 +68,13 @@ static double get_fc3_sum
const lapack_complex_double *eigvecs2,
const lapack_complex_double *fc3_reciprocal,
const double *masses,
const size_t num_atom,
const long num_atom,
const double cutoff_frequency);
void reciprocal_to_normal_squared
(double *fc3_normal_squared,
PHPYCONST int (*g_pos)[4],
const size_t num_g_pos,
const long num_g_pos,
const lapack_complex_double *fc3_reciprocal,
const double *freqs0,
const double *freqs1,
@ -84,12 +84,12 @@ void reciprocal_to_normal_squared
const lapack_complex_double *eigvecs2,
const double *masses,
const int *band_indices,
const size_t num_band0,
const size_t num_band,
const long num_band0,
const long num_band,
const double cutoff_frequency,
const int openmp_at_bands)
{
size_t i, num_atom;
long i, num_atom;
#ifdef MEASURE_R2N
double loopTotalCPUTime, loopTotalWallTime;
@ -132,9 +132,9 @@ void reciprocal_to_normal_squared
}
static double get_fc3_sum
(const size_t j,
const size_t k,
const size_t bi,
(const long j,
const long k,
const long bi,
const double *freqs0,
const double *freqs1,
const double *freqs2,
@ -143,7 +143,7 @@ static double get_fc3_sum
const lapack_complex_double *eigvecs2,
const lapack_complex_double *fc3_reciprocal,
const double *masses,
const size_t num_atom,
const long num_atom,
const double cutoff_frequency)
{
double fff, sum_real, sum_imag;
@ -166,17 +166,17 @@ static double get_fc3_sum
}
static lapack_complex_double fc3_sum_in_reciprocal_to_normal
(const size_t bi0,
const size_t bi1,
const size_t bi2,
(const long bi0,
const long bi1,
const long bi2,
const lapack_complex_double *eigvecs0,
const lapack_complex_double *eigvecs1,
const lapack_complex_double *eigvecs2,
const lapack_complex_double *fc3_reciprocal,
const double *masses,
const size_t num_atom)
const long num_atom)
{
size_t baseIndex, index_l, index_lm, i, j, k, l, m, n;
long baseIndex, index_l, index_lm, i, j, k, l, m, n;
double sum_real, sum_imag, mmm, mass_l, mass_lm;
lapack_complex_double eig_prod, eig_prod1;

6
c/reciprocal_to_normal.h
View File

@ -42,7 +42,7 @@
void reciprocal_to_normal_squared
(double *fc3_normal_squared,
PHPYCONST int (*g_pos)[4],
const size_t num_g_pos,
const long num_g_pos,
const lapack_complex_double *fc3_reciprocal,
const double *freqs0,
const double *freqs1,
@ -52,8 +52,8 @@ void reciprocal_to_normal_squared
const lapack_complex_double *eigvecs2,
const double *masses,
const int *band_indices,
const size_t num_band0,
const size_t num_band,
const long num_band0,
const long num_band,
const double cutoff_frequency,
const int openmp_at_bands);

163
c/rgrid.c Normal file
View File

@ -0,0 +1,163 @@
/* Copyright (C) 2015 Atsushi Togo */
/* All rights reserved. */
/* This file was originally part of spglib and is part of kspclib. */
/* 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 kspclib 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 <stdio.h>
#include <stddef.h>
#include <assert.h>
#include "rgrid.h"
static void get_all_grid_addresses(long grid_address[][3], const long mesh[3]);
static long get_double_grid_index(const long address_double[3],
const long mesh[3]);
static long get_grid_index_single_mesh(const long address[3],
const long mesh[3]);
static void reduce_grid_address(long address[3], const long mesh[3]);
static void reduce_double_grid_address(long address[3], const long mesh[3]);
static long mat_modulo_l(const long a, const long b);
void rgd_get_all_grid_addresses(long grid_address[][3], const long mesh[3])
{
get_all_grid_addresses(grid_address, mesh);
}
long rgd_get_double_grid_index(const long address_double[3],
const long mesh[3])
{
return get_double_grid_index(address_double, mesh);
}
void rgd_get_double_grid_address(long address_double[3],
const long address[3],
const long mesh[3],
const long is_shift[3])
{
long i;
for (i = 0; i < 3; i++) {
address_double[i] = address[i] * 2 + (is_shift[i] != 0);
}
reduce_double_grid_address(address_double, mesh);
}
static void get_all_grid_addresses(long grid_address[][3], const long mesh[3])
{
long i, j, k;
long grid_index;
long address[3];
for (i = 0; i < mesh[0]; i++) {
address[0] = i;
for (j = 0; j < mesh[1]; j++) {
address[1] = j;
for (k = 0; k < mesh[2]; k++) {
address[2] = k;
grid_index = get_grid_index_single_mesh(address, mesh);
assert(mesh[0] * mesh[1] * mesh[2] > grid_index);
grid_address[grid_index][0] = address[0];
grid_address[grid_index][1] = address[1];
grid_address[grid_index][2] = address[2];
reduce_grid_address(grid_address[grid_index], mesh);
}
}
}
}
static long get_double_grid_index(const long address_double[3],
const long mesh[3])
{
long i;
long address[3];
for (i = 0; i < 3; i++) {
if (address_double[i] % 2 == 0) {
address[i] = address_double[i] / 2;
} else {
address[i] = (address_double[i] - 1) / 2;
}
address[i] = mat_modulo_l(address[i], mesh[i]);
}
return get_grid_index_single_mesh(address, mesh);
}
static long get_grid_index_single_mesh(const long address[3],
const long mesh[3])
{
#ifndef GRID_ORDER_XYZ
return (address[2] * mesh[0] * (long)(mesh[1])
+ address[1] * mesh[0] + address[0]);
#else
return (address[0] * mesh[1] * (long)(mesh[2])
+ address[1] * mesh[2] + address[2]);
#endif
}
static void reduce_grid_address(long address[3], const long mesh[3])
{
long i;
for (i = 0; i < 3; i++) {
#ifndef GRID_BOUNDARY_AS_NEGATIVE
address[i] -= mesh[i] * (address[i] > mesh[i] / 2);
#else
address[i] -= mesh[i] * (address[i] > (mesh[i] - 1) / 2);
#endif
}
}
static void reduce_double_grid_address(long address[3], const long mesh[3])
{
long i;
for (i = 0; i < 3; i++) {
#ifndef GRID_BOUNDARY_AS_NEGATIVE
address[i] -= 2 * mesh[i] * (address[i] > mesh[i]);
#else
address[i] -= 2 * mesh[i] * (address[i] > mesh[i] - 1);
#endif
}
}
static long mat_modulo_l(const long a, const long b)
{
long c;
c = a % b;
if (c < 0) {
c += b;
}
return c;
}

81
c/rgrid.h Normal file
View File

@ -0,0 +1,81 @@
/* Copyright (C) 2015 Atsushi Togo */
/* All rights reserved. */
/* This file was originally part of spglib and is part of kspclib. */
/* 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 kspclib project nor the names of its */
/* contributors may be used to endorse or promote products derived */
/* from this software without specific prior written permission. */
/* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS */
/* "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT */
/* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS */
/* FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE */
/* COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, */
/* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, */
/* BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; */
/* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER */
/* CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT */
/* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN */
/* ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE */
/* POSSIBILITY OF SUCH DAMAGE. */
#ifndef __rgrid_H__
#define __rgrid_H__
/* #define GRID_ORDER_XYZ */
/* This changes behaviour of index order of address. */
/* Without GRID_ORDER_XYZ, left most element of address runs first. */
/* grid_address (e.g. 4x4x4 mesh, unless GRID_ORDER_XYZ is defined) */
/* [[ 0 0 0] */
/* [ 1 0 0] */
/* [ 2 0 0] */
/* [-1 0 0] */
/* [ 0 1 0] */
/* [ 1 1 0] */
/* [ 2 1 0] */
/* [-1 1 0] */
/* .... ] */
/* */
/* With GRID_ORDER_XYZ, right most element of address runs first. */
/* grid_address (e.g. 4x4x4 mesh, if GRID_ORDER_XYZ is defined) */
/* [[ 0 0 0] */
/* [ 0 0 1] */
/* [ 0 0 2] */
/* [ 0 0 -1] */
/* [ 0 1 0] */
/* [ 0 1 1] */
/* [ 0 1 2] */
/* [ 0 1 -1] */
/* .... ] */
/* #define GRID_BOUNDARY_AS_NEGATIVE */
/* This changes the behaviour of address elements on the surface of */
/* parallelepiped. */
/* For odd mesh number, this affects nothing, e.g., [-2, -1, 0, 1, 2]. */
/* regardless of with and without GRID_BOUNDARY_AS_NEGATIVE. */
/* For even mesh number, this affects as follows: */
/* without GRID_BOUNDARY_AS_NEGATIVE, e.g., [-2, -1, 0, 1, 2, 3]. */
/* with GRID_BOUNDARY_AS_NEGATIVE, e.g., [-3, -2, -1, 0, 1, 2]. */
void rgd_get_all_grid_addresses(long grid_address[][3], const long mesh[3]);
long rgd_get_double_grid_index(const long address_double[3],
const long mesh[3]);
void rgd_get_double_grid_address(long address_double[3],
const long address[3],
const long mesh[3],
const long is_shift[3]);
#endif

213
c/tetrahedron_method.c
View File

@ -15,7 +15,7 @@
/* the documentation and/or other materials provided with the */
/* distribution. */
/* * Neither the name of the phonopy project nor the names of its */
/* * Neither the name of the kspclib project nor the names of its */
/* contributors may be used to endorse or promote products derived */
/* from this software without specific prior written permission. */
@ -35,7 +35,6 @@
/* Copyright (C) 2014 Atsushi Togo */
#include <stddef.h>
#include "kgrid.h"
#include "tetrahedron_method.h"
#ifdef THMWARNING
@ -65,12 +64,12 @@
/* 7: c + a + b 3, 5, 6 */
static int main_diagonals[4][3] = {{ 1, 1, 1}, /* 0-7 */
{-1, 1, 1}, /* 1-6 */
{ 1,-1, 1}, /* 2-5 */
{ 1, 1,-1}}; /* 3-4 */
static long main_diagonals[4][3] = {{ 1, 1, 1}, /* 0-7 */
{-1, 1, 1}, /* 1-6 */
{ 1,-1, 1}, /* 2-5 */
{ 1, 1,-1}}; /* 3-4 */
static int db_relative_grid_address[4][24][4][3] = {
static long db_relative_grid_address[4][24][4][3] = {
{
{ { 0, 0, 0}, { 1, 0, 0}, { 1, 1, 0}, { 1, 1, 1}, },
{ { 0, 0, 0}, { 1, 0, 0}, { 1, 0, 1}, { 1, 1, 1}, },
@ -177,50 +176,38 @@ static int db_relative_grid_address[4][24][4][3] = {
},
};
static void
get_integration_weight_at_omegas(double *integration_weights,
const int num_omegas,
const double *omegas,
THMCONST double tetrahedra_omegas[24][4],
double (*gn)(const int,
const double,
const double[4]),
double (*IJ)(const int,
const int,
const double,
const double[4]));
static double
get_integration_weight(const double omega,
THMCONST double tetrahedra_omegas[24][4],
double (*gn)(const int,
double (*gn)(const long,
const double,
const double[4]),
double (*IJ)(const int,
const int,
double (*IJ)(const long,
const long,
const double,
const double[4]));
static int get_main_diagonal(THMCONST double rec_lattice[3][3]);
static int sort_omegas(double v[4]);
static long get_main_diagonal(THMCONST double rec_lattice[3][3]);
static long sort_omegas(double v[4]);
static double norm_squared_d3(const double a[3]);
static void multiply_matrix_vector_di3(double v[3],
THMCONST double a[3][3],
const int b[3]);
static double _f(const int n,
const int m,
const long b[3]);
static double _f(const long n,
const long m,
const double omega,
const double vertices_omegas[4]);
static double _J(const int i,
const int ci,
static double _J(const long i,
const long ci,
const double omega,
const double vertices_omegas[4]);
static double _I(const int i,
const int ci,
static double _I(const long i,
const long ci,
const double omega,
const double vertices_omegas[4]);
static double _n(const int i,
static double _n(const long i,
const double omega,
const double vertices_omegas[4]);
static double _g(const int i,
static double _g(const long i,
const double omega,
const double vertices_omegas[4]);
static double _n_0(void);
@ -293,10 +280,10 @@ static double _I_33(const double omega,
static double _I_4(void);
void thm_get_relative_grid_address(int relative_grid_address[24][4][3],
void thm_get_relative_grid_address(long relative_grid_address[24][4][3],
THMCONST double rec_lattice[3][3])
{
int i, j, k, main_diag_index;
long i, j, k, main_diag_index;
main_diag_index = get_main_diagonal(rec_lattice);
@ -310,9 +297,9 @@ void thm_get_relative_grid_address(int relative_grid_address[24][4][3],
}
}
void thm_get_all_relative_grid_address(int relative_grid_address[4][24][4][3])
void thm_get_all_relative_grid_address(long relative_grid_address[4][24][4][3])
{
int i, j, k, main_diag_index;
long i, j, k, main_diag_index;
for (main_diag_index = 0; main_diag_index < 4; main_diag_index++) {
for (i = 0; i < 24; i++) {
@ -341,126 +328,18 @@ double thm_get_integration_weight(const double omega,
}
}
void
thm_get_integration_weight_at_omegas(double *integration_weights,
const int num_omegas,
const double *omegas,
THMCONST double tetrahedra_omegas[24][4],
const char function)
{
if (function == 'I') {
get_integration_weight_at_omegas(integration_weights,
num_omegas,
omegas,
tetrahedra_omegas,
_g, _I);
} else {
get_integration_weight_at_omegas(integration_weights,
num_omegas,
omegas,
tetrahedra_omegas,
_n, _J);
}
}
void thm_get_neighboring_grid_points(int neighboring_grid_points[],
const int grid_point,
THMCONST int relative_grid_address[][3],
const int num_relative_grid_address,
const int mesh[3],
THMCONST int bz_grid_address[][3],
const int bz_map[])
{
int bzmesh[3], address_double[3], bz_address_double[3];
int i, j, bz_gp;
for (i = 0; i < 3; i++) {
bzmesh[i] = mesh[i] * 2;
}
for (i = 0; i < num_relative_grid_address; i++) {
for (j = 0; j < 3; j++) {
address_double[j] = (bz_grid_address[grid_point][j] +
relative_grid_address[i][j]) * 2;
bz_address_double[j] = address_double[j];
}
bz_gp = bz_map[kgd_get_grid_point_double_mesh(bz_address_double, bzmesh)];
if (bz_gp == -1) {
neighboring_grid_points[i] =
kgd_get_grid_point_double_mesh(address_double, mesh);
} else {
neighboring_grid_points[i] = bz_gp;
}
}
}
void
thm_get_dense_neighboring_grid_points(size_t neighboring_grid_points[],
const size_t grid_point,
THMCONST int relative_grid_address[][3],
const int num_relative_grid_address,
const int mesh[3],
THMCONST int bz_grid_address[][3],
const size_t bz_map[])
{
int bzmesh[3], address_double[3], bz_address_double[3];
size_t i, j, bz_gp, prod_bz_mesh;
prod_bz_mesh = 1;
for (i = 0; i < 3; i++) {
bzmesh[i] = mesh[i] * 2;
prod_bz_mesh *= bzmesh[i];
}
for (i = 0; i < num_relative_grid_address; i++) {
for (j = 0; j < 3; j++) {
address_double[j] = (bz_grid_address[grid_point][j] +
relative_grid_address[i][j]) * 2;
bz_address_double[j] = address_double[j];
}
bz_gp = bz_map[kgd_get_dense_grid_point_double_mesh(bz_address_double, bzmesh)];
if (bz_gp == prod_bz_mesh) {
neighboring_grid_points[i] =
kgd_get_dense_grid_point_double_mesh(address_double, mesh);
} else {
neighboring_grid_points[i] = bz_gp;
}
}
}
static void
get_integration_weight_at_omegas(double *integration_weights,
const int num_omegas,
const double *omegas,
THMCONST double tetrahedra_omegas[24][4],
double (*gn)(const int,
const double,
const double[4]),
double (*IJ)(const int,
const int,
const double,
const double[4]))
{
int i;
#pragma omp parallel for
for (i = 0; i < num_omegas; i++) {
integration_weights[i] = get_integration_weight(omegas[i],
tetrahedra_omegas,
gn, IJ);
}
}
static double
get_integration_weight(const double omega,
THMCONST double tetrahedra_omegas[24][4],
double (*gn)(const int,
double (*gn)(const long,
const double,
const double[4]),
double (*IJ)(const int,
const int,
double (*IJ)(const long,
const long,
const double,
const double[4]))
{
int i, j, ci;
long i, j, ci;
double sum;
double v[4];
@ -493,9 +372,9 @@ get_integration_weight(const double omega,
return sum / 6;
}
static int sort_omegas(double v[4])
static long sort_omegas(double v[4])
{
int i;
long i;
double w[4];
i = 0;
@ -563,9 +442,9 @@ static int sort_omegas(double v[4])
return i;
}
static int get_main_diagonal(THMCONST double rec_lattice[3][3])
static long get_main_diagonal(THMCONST double rec_lattice[3][3])
{
int i, shortest;
long i, shortest;
double length, min_length;
double main_diag[3];
@ -590,9 +469,9 @@ static double norm_squared_d3(const double a[3])
static void multiply_matrix_vector_di3(double v[3],
THMCONST double a[3][3],
const int b[3])
const long b[3])
{
int i;
long i;
double c[3];
for (i = 0; i < 3; i++) {
@ -604,8 +483,8 @@ static void multiply_matrix_vector_di3(double v[3],
}
}
static double _f(const int n,
const int m,
static double _f(const long n,
const long m,
const double omega,
const double vertices_omegas[4])
{
@ -613,8 +492,8 @@ static double _f(const int n,
(vertices_omegas[n] - vertices_omegas[m]));
}
static double _J(const int i,
const int ci,
static double _J(const long i,
const long ci,
const double omega,
const double vertices_omegas[4])
{
@ -667,8 +546,8 @@ static double _J(const int i,
}
static double _I(const int i,
const int ci,
static double _I(const long i,
const long ci,
const double omega,
const double vertices_omegas[4])
{
@ -720,7 +599,7 @@ static double _I(const int i,
return 0;
}
static double _n(const int i,
static double _n(const long i,
const double omega,
const double vertices_omegas[4])
{
@ -745,7 +624,7 @@ static double _n(const int i,
return 0;
}
static double _g(const int i,
static double _g(const long i,
const double omega,
const double vertices_omegas[4])
{
@ -847,10 +726,10 @@ static double _g_2(const double omega,
static double _g_3(const double omega,
const double vertices_omegas[4])
{
return (3 *
_f(1, 3, omega, vertices_omegas) *
_f(2, 3, omega, vertices_omegas) /
(vertices_omegas[3] - vertices_omegas[0]));
return (3 *
_f(1, 3, omega, vertices_omegas) *
_f(2, 3, omega, vertices_omegas) /
(vertices_omegas[3] - vertices_omegas[0]));
}
/* omega4 < omega */

27
c/tetrahedron_method.h
View File

@ -15,7 +15,7 @@
/* the documentation and/or other materials provided with the */
/* distribution. */
/* * Neither the name of the phonopy project nor the names of its */
/* * Neither the name of the kspclib project nor the names of its */
/* contributors may be used to endorse or promote products derived */
/* from this software without specific prior written permission. */
@ -41,32 +41,11 @@
#include <stddef.h>
void thm_get_relative_grid_address(int relative_grid_address[24][4][3],
void thm_get_relative_grid_address(long relative_grid_address[24][4][3],
THMCONST double rec_lattice[3][3]);
void thm_get_all_relative_grid_address(int relative_grid_address[4][24][4][3]);
void thm_get_all_relative_grid_address(long relative_grid_address[4][24][4][3]);
double thm_get_integration_weight(const double omega,
THMCONST double tetrahedra_omegas[24][4],
const char function);
void
thm_get_integration_weight_at_omegas(double *integration_weights,
const int num_omegas,
const double *omegas,
THMCONST double tetrahedra_omegas[24][4],
const char function);
void thm_get_neighboring_grid_points(int neighboring_grid_points[],
const int grid_point,
THMCONST int relative_grid_address[][3],
const int num_relative_grid_address,
const int mesh[3],
THMCONST int bz_grid_address[][3],
const int bz_map[]);
void
thm_get_dense_neighboring_grid_points(size_t neighboring_grid_points[],
const size_t grid_point,
THMCONST int relative_grid_address[][3],
const int num_relative_grid_address,
const int mesh[3],
THMCONST int bz_grid_address[][3],
const size_t bz_map[]);
#endif

105
c/triplet.c
View File

@ -40,23 +40,23 @@
#include "triplet_iw.h"
#include "triplet_kpoint.h"
static size_t get_triplets_reciprocal_mesh_at_q(size_t *map_triplets,
size_t *map_q,
int (*grid_address)[3],
const int grid_point,
const int mesh[3],
const int is_time_reversal,
const int num_rot,
TPLCONST int (*rotations)[3][3],
const int swappable);
static long get_triplets_reciprocal_mesh_at_q(long *map_triplets,
long *map_q,
int (*grid_address)[3],
const int grid_point,
const int mesh[3],
const int is_time_reversal,
const long num_rot,
TPLCONST int (*rotations)[3][3],
const int swappable);
size_t tpl_get_BZ_triplets_at_q(size_t (*triplets)[3],
const size_t grid_point,
TPLCONST int (*bz_grid_address)[3],
const size_t *bz_map,
const size_t *map_triplets,
const size_t num_map_triplets,
const int mesh[3])
long tpl_get_BZ_triplets_at_q(long (*triplets)[3],
const long grid_point,
TPLCONST int (*bz_grid_address)[3],
const long *bz_map,
const long *map_triplets,
const long num_map_triplets,
const int mesh[3])
{
return tpk_get_BZ_triplets_at_q(triplets,
grid_point,
@ -67,15 +67,15 @@ size_t tpl_get_BZ_triplets_at_q(size_t (*triplets)[3],
mesh);
}
size_t tpl_get_triplets_reciprocal_mesh_at_q(size_t *map_triplets,
size_t *map_q,
int (*grid_address)[3],
const size_t grid_point,
const int mesh[3],
const int is_time_reversal,
const int num_rot,
TPLCONST int (*rotations)[3][3],
const int swappable)
long tpl_get_triplets_reciprocal_mesh_at_q(long *map_triplets,
long *map_q,
int (*grid_address)[3],
const long grid_point,
const int mesh[3],
const int is_time_reversal,
const long num_rot,
TPLCONST int (*rotations)[3][3],
const int swappable)
{
return get_triplets_reciprocal_mesh_at_q(map_triplets,
map_q,
@ -91,22 +91,22 @@ size_t tpl_get_triplets_reciprocal_mesh_at_q(size_t *map_triplets,
void tpl_get_integration_weight(double *iw,
char *iw_zero,
const double *frequency_points,
const size_t num_band0,
const long num_band0,
TPLCONST int relative_grid_address[24][4][3],
const int mesh[3],
TPLCONST size_t (*triplets)[3],
const size_t num_triplets,
TPLCONST long (*triplets)[3],
const long num_triplets,
TPLCONST int (*bz_grid_address)[3],
const size_t *bz_map,
const long *bz_map,
const double *frequencies1,
const size_t num_band1,
const long num_band1,
const double *frequencies2,
const size_t num_band2,
const size_t tp_type,
const long num_band2,
const long tp_type,
const int openmp_per_triplets,
const int openmp_per_bands)
{
size_t i, num_band_prod;
long i, num_band_prod;
int tp_relative_grid_address[2][24][4][3];
tpl_set_relative_grid_address(tp_relative_grid_address,
@ -141,14 +141,14 @@ void tpl_get_integration_weight_with_sigma(double *iw,
const double sigma,
const double sigma_cutoff,
const double *frequency_points,
const size_t num_band0,
TPLCONST size_t (*triplets)[3],
const size_t num_triplets,
const long num_band0,
TPLCONST long (*triplets)[3],
const long num_triplets,
const double *frequencies,
const size_t num_band,
const size_t tp_type)
const long num_band,
const long tp_type)
{
size_t i, num_band_prod, const_adrs_shift;
long i, num_band_prod, const_adrs_shift;
double cutoff;
cutoff = sigma * sigma_cutoff;
@ -174,7 +174,7 @@ void tpl_get_integration_weight_with_sigma(double *iw,
}
int tpl_is_N(const size_t triplet[3], const int *grid_address)
int tpl_is_N(const long triplet[3], const int *grid_address)
{
int i, j, sum_q, is_N;
@ -195,9 +195,9 @@ int tpl_is_N(const size_t triplet[3], const int *grid_address)
void tpl_set_relative_grid_address(
int tp_relative_grid_address[2][24][4][3],
TPLCONST int relative_grid_address[24][4][3],
const size_t tp_type)
const long tp_type)
{
size_t i, j, k, l;
long i, j, k, l;
int signs[2];
signs[0] = 1;
@ -221,19 +221,18 @@ void tpl_set_relative_grid_address(
}
}
static size_t get_triplets_reciprocal_mesh_at_q(size_t *map_triplets,
size_t *map_q,
int (*grid_address)[3],
const int grid_point,
const int mesh[3],
const int is_time_reversal,
const int num_rot,
TPLCONST int (*rotations)[3][3],
const int swappable)
static long get_triplets_reciprocal_mesh_at_q(long *map_triplets,
long *map_q,
int (*grid_address)[3],
const int grid_point,
const int mesh[3],
const int is_time_reversal,
const long num_rot,
TPLCONST int (*rotations)[3][3],
const int swappable)
{
MatINT *rot_real;
int i;
size_t num_ir;
long i, num_ir;
rot_real = mat_alloc_MatINT(num_rot);
for (i = 0; i < num_rot; i++) {

60
c/triplet.h
View File

@ -49,43 +49,43 @@
/* and map_q[prod(mesh)] are required. rotations are point-group- */
/* operations in real space for which duplicate operations are allowed */
/* in the input. */
size_t tpl_get_triplets_reciprocal_mesh_at_q(size_t *map_triplets,
size_t *map_q,
int (*grid_address)[3],
const size_t grid_point,
const int mesh[3],
const int is_time_reversal,
const int num_rot,
TPLCONST int (*rotations)[3][3],
const int swappable);
long tpl_get_triplets_reciprocal_mesh_at_q(long *map_triplets,
long *map_q,
int (*grid_address)[3],
const long grid_point,
const int mesh[3],
const int is_time_reversal,
const long num_rot,
TPLCONST int (*rotations)[3][3],
const int swappable);
/* Irreducible grid-point-triplets in BZ are stored. */
/* triplets are recovered from grid_point and triplet_weights. */
/* BZ boundary is considered in this recovery. Therefore grid addresses */
/* are given not by grid_address, but by bz_grid_address. */
/* triplets[num_ir_triplets][3] = number of non-zero triplets weights*/
/* Number of ir-triplets is returned. */
size_t tpl_get_BZ_triplets_at_q(size_t (*triplets)[3],
const size_t grid_point,
TPLCONST int (*bz_grid_address)[3],
const size_t *bz_map,
const size_t *map_triplets,
const size_t num_map_triplets,
const int mesh[3]);
long tpl_get_BZ_triplets_at_q(long (*triplets)[3],
const long grid_point,
TPLCONST int (*bz_grid_address)[3],
const long *bz_map,
const long *map_triplets,
const long num_map_triplets,
const int mesh[3]);
void tpl_get_integration_weight(double *iw,
char *iw_zero,
const double *frequency_points,
const size_t num_band0,
const long num_band0,
TPLCONST int relative_grid_address[24][4][3],
const int mesh[3],
TPLCONST size_t (*triplets)[3],
const size_t num_triplets,
TPLCONST long (*triplets)[3],
const long num_triplets,
TPLCONST int (*bz_grid_address)[3],
const size_t *bz_map,
const long *bz_map,
const double *frequencies1,
const size_t num_band1,
const long num_band1,
const double *frequencies2,
const size_t num_band2,
const size_t tp_type,
const long num_band2,
const long tp_type,
const int openmp_per_triplets,
const int openmp_per_bands);
void tpl_get_integration_weight_with_sigma(double *iw,
@ -93,17 +93,17 @@ void tpl_get_integration_weight_with_sigma(double *iw,
const double sigma,
const double sigma_cutoff,
const double *frequency_points,
const size_t num_band0,
TPLCONST size_t (*triplets)[3],
const size_t num_triplets,
const long num_band0,
TPLCONST long (*triplets)[3],
const long num_triplets,
const double *frequencies,
const size_t num_band,
const size_t tp_type);
const long num_band,
const long tp_type);
int tpl_is_N(const size_t triplet[3], const int *grid_address);
int tpl_is_N(const long triplet[3], const int *grid_address);
void tpl_set_relative_grid_address(
int tp_relative_grid_address[2][24][4][3],
TPLCONST int relative_grid_address[24][4][3],
const size_t tp_type);
const long tp_type);
#endif

92
c/triplet_iw.c
View File

@ -34,6 +34,7 @@
#include <stddef.h>
#include <math.h>
#include "kgrid.h"
#include "phonoc_utils.h"
#include "triplet.h"
#include "triplet_iw.h"
@ -42,45 +43,45 @@
static void set_freq_vertices(double freq_vertices[3][24][4],
const double *frequencies1,
const double *frequencies2,
TPLCONST size_t vertices[2][24][4],
TPLCONST long vertices[2][24][4],
const int num_band1,
const int num_band2,
const int b1,
const int b2,
const size_t tp_type);
const long tp_type);
static int set_g(double g[3],
const double f0,
TPLCONST double freq_vertices[3][24][4],
const size_t max_i);
const long max_i);
static int in_tetrahedra(const double f0, TPLCONST double freq_vertices[24][4]);
static void get_triplet_tetrahedra_vertices(
size_t vertices[2][24][4],
long vertices[2][24][4],
TPLCONST int tp_relative_grid_address[2][24][4][3],
const int mesh[3],
const size_t triplet[3],
const long triplet[3],
TPLCONST int (*bz_grid_address)[3],
const size_t *bz_map);
const long *bz_map);
void
tpi_get_integration_weight(double *iw,
char *iw_zero,
const double *frequency_points,
const size_t num_band0,
const long num_band0,
TPLCONST int tp_relative_grid_address[2][24][4][3],
const int mesh[3],
const size_t triplets[3],
const size_t num_triplets,
const long triplets[3],
const long num_triplets,
TPLCONST int (*bz_grid_address)[3],
const size_t *bz_map,
const long *bz_map,
const double *frequencies1,
const size_t num_band1,
const long num_band1,
const double *frequencies2,
const size_t num_band2,
const size_t tp_type,
const long num_band2,
const long tp_type,
const int openmp_per_bands)
{
size_t max_i, j, b1, b2, b12, num_band_prod, adrs_shift;
size_t vertices[2][24][4];
long max_i, j, b1, b2, b12, num_band_prod, adrs_shift;
long vertices[2][24][4];
double g[3];
double freq_vertices[3][24][4];
@ -144,15 +145,15 @@ void tpi_get_integration_weight_with_sigma(double *iw,
const double sigma,
const double cutoff,
const double *frequency_points,
const size_t num_band0,
const size_t triplet[3],
const size_t const_adrs_shift,
const long num_band0,
const long triplet[3],
const long const_adrs_shift,
const double *frequencies,
const size_t num_band,
const size_t tp_type,
const long num_band,
const long tp_type,
const int openmp_per_bands)
{
size_t j, b12, b1, b2, adrs_shift;
long j, b12, b1, b2, adrs_shift;
double f0, f1, f2, g0, g1, g2;
#pragma omp parallel for private(j, b1, b2, f0, f1, f2, g0, g1, g2, adrs_shift) if (openmp_per_bands)
@ -206,15 +207,50 @@ void tpi_get_integration_weight_with_sigma(double *iw,
}
}
void
tpi_get_dense_neighboring_grid_points(long neighboring_grid_points[],
const long grid_point,
TPLCONST int relative_grid_address[][3],
const int num_relative_grid_address,
const int mesh[3],
TPLCONST int bz_grid_address[][3],
const long bz_map[])
{
int bzmesh[3], address_double[3], bz_address_double[3];
long i, j, bz_gp, prod_bz_mesh;
prod_bz_mesh = 1;
for (i = 0; i < 3; i++) {
bzmesh[i] = mesh[i] * 2;
prod_bz_mesh *= bzmesh[i];
}
for (i = 0; i < num_relative_grid_address; i++) {
for (j = 0; j < 3; j++) {
address_double[j] = (bz_grid_address[grid_point][j] +
relative_grid_address[i][j]) * 2;
bz_address_double[j] = address_double[j];
}
bz_gp = bz_map[kgd_get_dense_grid_point_double_mesh(bz_address_double, bzmesh)];
if (bz_gp == prod_bz_mesh) {
neighboring_grid_points[i] =
kgd_get_dense_grid_point_double_mesh(address_double, mesh);
} else {
neighboring_grid_points[i] = bz_gp;
}
}
}
static void set_freq_vertices(double freq_vertices[3][24][4],
const double *frequencies1,
const double *frequencies2,
TPLCONST size_t vertices[2][24][4],
TPLCONST long vertices[2][24][4],
const int num_band1,
const int num_band2,
const int b1,
const int b2,
const size_t tp_type)
const long tp_type)
{
int i, j;
double f1, f2;
@ -248,7 +284,7 @@ static void set_freq_vertices(double freq_vertices[3][24][4],
static int set_g(double g[3],
const double f0,
TPLCONST double freq_vertices[3][24][4],
const size_t max_i)
const long max_i)
{
int i, iw_zero;
@ -293,18 +329,18 @@ static int in_tetrahedra(const double f0, TPLCONST double freq_vertices[24][4])
}
static void get_triplet_tetrahedra_vertices(
size_t vertices[2][24][4],
long vertices[2][24][4],
TPLCONST int tp_relative_grid_address[2][24][4][3],
const int mesh[3],
const size_t triplet[3],
const long triplet[3],
TPLCONST int (*bz_grid_address)[3],
const size_t *bz_map)
const long *bz_map)
{
int i, j;
for (i = 0; i < 2; i++) {
for (j = 0; j < 24; j++) {
thm_get_dense_neighboring_grid_points(vertices[i][j],
tpi_get_dense_neighboring_grid_points(vertices[i][j],
triplet[i + 1],
tp_relative_grid_address[i][j],
4,

32
c/triplet_iw.h
View File

@ -42,30 +42,38 @@ void
tpi_get_integration_weight(double *iw,
char *iw_zero,
const double *frequency_points,
const size_t num_band0,
const long num_band0,
TPLCONST int tp_relative_grid_address[2][24][4][3],
const int mesh[3],
const size_t triplets[3],
const size_t num_triplets,
const long triplets[3],
const long num_triplets,
TPLCONST int (*bz_grid_address)[3],
const size_t *bz_map,
const long *bz_map,
const double *frequencies1,
const size_t num_band1,
const long num_band1,
const double *frequencies2,
const size_t num_band2,
const size_t tp_type,
const long num_band2,
const long tp_type,
const int openmp_per_bands);
void tpi_get_integration_weight_with_sigma(double *iw,
char *iw_zero,
const double sigma,
const double cutoff,
const double *frequency_points,
const size_t num_band0,
const size_t triplet[3],
const size_t const_adrs_shift,
const long num_band0,
const long triplet[3],
const long const_adrs_shift,
const double *frequencies,
const size_t num_band,
const size_t tp_type,
const long num_band,
const long tp_type,
const int openmp_per_bands);
void
tpi_get_dense_neighboring_grid_points(long neighboring_grid_points[],
const long grid_point,
TPLCONST int relative_grid_address[][3],
const int num_relative_grid_address,
const int mesh[3],
TPLCONST int bz_grid_address[][3],
const long bz_map[]);
#endif

118
c/triplet_kpoint.c
View File

@ -172,38 +172,38 @@ static int bz_search_space[KPT_NUM_BZ_SEARCH_SPACE][3] = {
};
static void grid_point_to_address_double(int address_double[3],
const size_t grid_point,
const long grid_point,
const int mesh[3],
const int is_shift[3]);
static size_t get_ir_triplets_at_q(size_t *map_triplets,
size_t *map_q,
int (*grid_address)[3],
const size_t grid_point,
const int mesh[3],
const MatINT * rot_reciprocal,
const int swappable);
static size_t get_BZ_triplets_at_q(size_t (*triplets)[3],
const size_t grid_point,
TPLCONST int (*bz_grid_address)[3],
const size_t *bz_map,
const size_t *map_triplets,
const size_t num_map_triplets,
const int mesh[3]);
static long get_ir_triplets_at_q(long *map_triplets,
long *map_q,
int (*grid_address)[3],
const long grid_point,
const int mesh[3],
const MatINT * rot_reciprocal,
const int swappable);
static long get_BZ_triplets_at_q(long (*triplets)[3],
const long grid_point,
TPLCONST int (*bz_grid_address)[3],
const long *bz_map,
const long *map_triplets,
const long num_map_triplets,
const int mesh[3]);
static int get_third_q_of_triplets_at_q(int bz_address[3][3],
const int q_index,
const size_t *bz_map,
const long *bz_map,
const int mesh[3],
const int bzmesh[3]);
static void modulo_i3(int v[3], const int m[3]);
size_t tpk_get_ir_triplets_at_q(size_t *map_triplets,
size_t *map_q,
int (*grid_address)[3],
const int grid_point,
const int mesh[3],
const int is_time_reversal,
const MatINT * rotations,
const int swappable)
long tpk_get_ir_triplets_at_q(long *map_triplets,
long *map_q,
int (*grid_address)[3],
const int grid_point,
const int mesh[3],
const int is_time_reversal,
const MatINT * rotations,
const int swappable)
{
int num_ir;
MatINT *rot_reciprocal;
@ -220,13 +220,13 @@ size_t tpk_get_ir_triplets_at_q(size_t *map_triplets,
return num_ir;
}
size_t tpk_get_BZ_triplets_at_q(size_t (*triplets)[3],
const size_t grid_point,
TPLCONST int (*bz_grid_address)[3],
const size_t *bz_map,
const size_t *map_triplets,
const size_t num_map_triplets,
const int mesh[3])
long tpk_get_BZ_triplets_at_q(long (*triplets)[3],
const long grid_point,
TPLCONST int (*bz_grid_address)[3],
const long *bz_map,
const long *map_triplets,
const long num_map_triplets,
const int mesh[3])
{
return get_BZ_triplets_at_q(triplets,
grid_point,
@ -237,18 +237,18 @@ size_t tpk_get_BZ_triplets_at_q(size_t (*triplets)[3],
mesh);
}
static size_t get_ir_triplets_at_q(size_t *map_triplets,
size_t *map_q,
int (*grid_address)[3],
const size_t grid_point,
const int mesh[3],
const MatINT * rot_reciprocal,
const int swappable)
static long get_ir_triplets_at_q(long *map_triplets,
long *map_q,
int (*grid_address)[3],
const long grid_point,
const int mesh[3],
const MatINT * rot_reciprocal,
const int swappable)
{
size_t i, j, num_grid, q_2, num_ir_q, num_ir_triplets, ir_grid_point;
long i, j, num_grid, q_2, num_ir_q, num_ir_triplets, ir_grid_point;
int mesh_double[3], is_shift[3];
int address_double0[3], address_double1[3], address_double2[3];
size_t *ir_grid_points, *third_q;
long *ir_grid_points, *third_q;
double tolerance;
double stabilizer_q[1][3];
MatINT *rot_reciprocal_q;
@ -258,7 +258,7 @@ static size_t get_ir_triplets_at_q(size_t *map_triplets,
rot_reciprocal_q = NULL;
tolerance = 0.01 / (mesh[0] + mesh[1] + mesh[2]);
num_grid = mesh[0] * mesh[1] * (size_t)mesh[2];
num_grid = mesh[0] * mesh[1] * (long)mesh[2];
for (i = 0; i < 3; i++) {
/* Only consider the gamma-point */
@ -286,8 +286,8 @@ static size_t get_ir_triplets_at_q(size_t *map_triplets,
mat_free_MatINT(rot_reciprocal_q);
rot_reciprocal_q = NULL;
third_q = (size_t*) malloc(sizeof(size_t) * num_ir_q);
ir_grid_points = (size_t*) malloc(sizeof(size_t) * num_ir_q);
third_q = (long*) malloc(sizeof(long) * num_ir_q);
ir_grid_points = (long*) malloc(sizeof(long) * num_ir_q);
num_ir_q = 0;
for (i = 0; i < num_grid; i++) {
if (map_q[i] == i) {
@ -346,18 +346,18 @@ static size_t get_ir_triplets_at_q(size_t *map_triplets,
return num_ir_triplets;
}
static size_t get_BZ_triplets_at_q(size_t (*triplets)[3],
const size_t grid_point,
TPLCONST int (*bz_grid_address)[3],
const size_t *bz_map,
const size_t *map_triplets,
const size_t num_map_triplets,
const int mesh[3])
static long get_BZ_triplets_at_q(long (*triplets)[3],
const long grid_point,
TPLCONST int (*bz_grid_address)[3],
const long *bz_map,
const long *map_triplets,
const long num_map_triplets,
const int mesh[3])
{
size_t i, num_ir;
long i, num_ir;
int j, k;
int bz_address[3][3], bz_address_double[3], bzmesh[3];
size_t *ir_grid_points;
long *ir_grid_points;
ir_grid_points = NULL;
@ -366,7 +366,7 @@ static size_t get_BZ_triplets_at_q(size_t (*triplets)[3],
}
num_ir = 0;
ir_grid_points = (size_t*) malloc(sizeof(size_t) * num_map_triplets);
ir_grid_points = (long*) malloc(sizeof(long) * num_map_triplets);
for (i = 0; i < num_map_triplets; i++) {
if (map_triplets[i] == i) {
ir_grid_points[num_ir] = i;
@ -407,16 +407,16 @@ static size_t get_BZ_triplets_at_q(size_t (*triplets)[3],
static int get_third_q_of_triplets_at_q(int bz_address[3][3],
const int q_index,
const size_t *bz_map,
const long *bz_map,
const int mesh[3],
const int bzmesh[3])
{
int i, j, smallest_g, smallest_index, sum_g, delta_g[3];
size_t prod_bzmesh;
size_t bzgp[KPT_NUM_BZ_SEARCH_SPACE];
long prod_bzmesh;
long bzgp[KPT_NUM_BZ_SEARCH_SPACE];
int bz_address_double[3];
prod_bzmesh = (size_t)bzmesh[0] * bzmesh[1] * bzmesh[2];
prod_bzmesh = (long)bzmesh[0] * bzmesh[1] * bzmesh[2];
modulo_i3(bz_address[q_index], mesh);
for (i = 0; i < 3; i++) {
@ -442,7 +442,7 @@ static int get_third_q_of_triplets_at_q(int bz_address[3][3],
}
}
escape:
escape:
smallest_g = 4;
smallest_index = 0;
@ -467,7 +467,7 @@ static int get_third_q_of_triplets_at_q(int bz_address[3][3],
}
static void grid_point_to_address_double(int address_double[3],
const size_t grid_point,
const long grid_point,
const int mesh[3],
const int is_shift[3])
{

30
c/triplet_kpoint.h
View File

@ -41,20 +41,20 @@
#include "mathfunc.h"
#include "triplet.h"
size_t tpk_get_ir_triplets_at_q(size_t *map_triplets,
size_t *map_q,
int (*grid_address)[3],
const int grid_point,
const int mesh[3],
const int is_time_reversal,
const MatINT * rotations,
const int swappable);
size_t tpk_get_BZ_triplets_at_q(size_t (*triplets)[3],
const size_t grid_point,
TPLCONST int (*bz_grid_address)[3],
const size_t *bz_map,
const size_t *map_triplets,
const size_t num_map_triplets,
const int mesh[3]);
long tpk_get_ir_triplets_at_q(long *map_triplets,
long *map_q,
int (*grid_address)[3],
const int grid_point,
const int mesh[3],
const int is_time_reversal,
const MatINT * rotations,
const int swappable);
long tpk_get_BZ_triplets_at_q(long (*triplets)[3],
const long grid_point,
TPLCONST int (*bz_grid_address)[3],
const long *bz_map,
const long *map_triplets,
const long num_map_triplets,
const int mesh[3]);
#endif

8
phono3py/other/isotope.py
View File

@ -103,7 +103,7 @@ class Isotope(object):
def set_grid_point(self, grid_point):
self._grid_point = grid_point
self._grid_points = np.arange(np.prod(self._mesh), dtype='uintp')
self._grid_points = np.arange(np.prod(self._mesh), dtype='int_')
if self._grid_address is None:
primitive_lattice = np.linalg.inv(self._primitive.cell)
@ -208,7 +208,7 @@ class Isotope(object):
gamma = np.zeros(len(self._band_indices), dtype='double')
if self._sigma is None:
self._set_integration_weights()
weights = np.ones(len(self._grid_points), dtype='intc')
weights = np.ones(len(self._grid_points), dtype='int_')
phono3c.thm_isotope_strength(gamma,
self._grid_point,
self._grid_points,
@ -246,7 +246,7 @@ class Isotope(object):
import phono3py._phono3py as phono3c
unique_vertices = thm.get_unique_tetrahedra_vertices()
neighboring_grid_points = np.zeros(
len(unique_vertices) * len(self._grid_points), dtype='uintp')
len(unique_vertices) * len(self._grid_points), dtype='int_')
phono3c.neighboring_grid_points(
neighboring_grid_points,
self._grid_points,
@ -255,7 +255,7 @@ class Isotope(object):
self._grid_address,
self._bz_map)
unique_grid_points = np.array(np.unique(neighboring_grid_points),
dtype='uintp')
dtype='int_')
self._run_phonon_solver_c(unique_grid_points)
freq_points = np.array(
self._frequencies[self._grid_point, self._band_indices],

4
phono3py/phonon/solver.py
View File

@ -75,8 +75,8 @@ def run_phonon_solver_c(dm,
Lambda = 0
fc = dm.force_constants
# assert grid_points.dtype == 'int_'
# assert grid_points.flags.c_contiguous
assert grid_points.dtype == 'int_'
assert grid_points.flags.c_contiguous
fc_p2s, fc_s2p = _get_fc_elements_mapping(dm, fc)
phononmod.phonons_at_gridpoints(

2
phono3py/phonon3/conductivity.py
View File

@ -340,7 +340,7 @@ class Conductivity(object):
elif not self._is_kappa_star: # All grid points
coarse_grid_address = get_grid_address(self._coarse_mesh)
coarse_grid_points = np.arange(np.prod(self._coarse_mesh),
dtype='uintp')
dtype='int_')
self._grid_points = from_coarse_to_dense_grid_points(
self._mesh,
self._mesh_divisors,

6
phono3py/phonon3/conductivity_LBTE.py
View File

@ -1023,7 +1023,7 @@ class Conductivity_LBTE(Conductivity):
6), dtype='double')
self._rot_grid_points = np.zeros(
(len(self._ir_grid_points), len(self._point_operations)),
dtype='uintp')
dtype='int_')
for i, ir_gp in enumerate(self._ir_grid_points):
self._rot_grid_points[i] = get_grid_points_by_rotations(
self._grid_address[ir_gp],
@ -1199,7 +1199,7 @@ class Conductivity_LBTE(Conductivity):
num_mesh_points = np.prod(self._mesh)
num_rot = len(self._point_operations)
rot_grid_points = np.zeros((num_rot, num_mesh_points), dtype='uintp')
rot_grid_points = np.zeros((num_rot, num_mesh_points), dtype='int_')
for i in range(num_mesh_points):
rot_grid_points[:, i] = get_grid_points_by_rotations(
@ -1253,7 +1253,7 @@ class Conductivity_LBTE(Conductivity):
"""Returns weights used for collision matrix and |X> and |f>
self._rot_grid_points : ndarray
shape=(ir_grid_points, point_operations), dtype='uintp'
shape=(ir_grid_points, point_operations), dtype='int_'
r_gps : grid points of arms of k-star with duplicates
len(r_gps) == order of crystallographic point group

6
phono3py/phonon3/interaction.py
View File

@ -331,7 +331,7 @@ class Interaction(object):
if self._nac_q_direction is not None:
if (grid_address[grid_point] == 0).all():
self._phonon_done[grid_point] = 0
self.run_phonon_solver(np.array([grid_point], dtype='uintp'))
self.run_phonon_solver(np.array([grid_point], dtype='int_'))
rotations = []
for r in self._symmetry.get_pointgroup_operations():
dq = self._nac_q_direction
@ -397,7 +397,7 @@ class Interaction(object):
if solve_dynamical_matrices:
self.run_phonon_solver(verbose=verbose)
else:
self.run_phonon_solver(np.array([0], dtype='uintp'),
self.run_phonon_solver(np.array([0], dtype='int_'),
verbose=verbose)
if (self._grid_address[0] == 0).all():
@ -437,7 +437,7 @@ class Interaction(object):
def run_phonon_solver(self, grid_points=None, verbose=False):
if grid_points is None:
_grid_points = np.arange(len(self._grid_address), dtype='uintp')
_grid_points = np.arange(len(self._grid_address), dtype='int_')
else:
_grid_points = grid_points
self._run_phonon_solver_c(_grid_points, verbose=verbose)

4
phono3py/phonon3/joint_dos.py
View File

@ -109,7 +109,7 @@ class JointDos(object):
def run(self):
self.run_phonon_solver(
np.arange(len(self._grid_address), dtype='uintp'))
np.arange(len(self._grid_address), dtype='int_'))
try:
import phono3py._phono3py as phono3c
self._run_c()
@ -189,7 +189,7 @@ class JointDos(object):
self._joint_dos = None
self._frequency_points = None
self._phonon_done[0] = 0
self.run_phonon_solver(np.array([grid_point], dtype='uintp'))
self.run_phonon_solver(np.array([grid_point], dtype='int_'))
def get_triplets_at_q(self):
return self._triplets_at_q, self._weights_at_q

28
phono3py/phonon3/triplets.py
View File

@ -71,7 +71,7 @@ def get_triplets_at_q(grid_point,
triplets_at_q : ndarray
Symmetry reduced number of triplets are stored as grid point
integer numbers.
shape=(n_triplets, 3), dtype='uintp'
shape=(n_triplets, 3), dtype='int_'
weights : ndarray
Weights of triplets in Brillouin zone
shape=(n_triplets,), dtype='intc'
@ -84,7 +84,7 @@ def get_triplets_at_q(grid_point,
bz_map : ndarray
Grid point mapping table containing BZ surface. See more
detail in spglib docstring.
shape=(prod(mesh*2),), dtype='uintp'
shape=(prod(mesh*2),), dtype='int_'
map_tripelts : ndarray or None
Returns when stores_triplets_map=True, otherwise None is
returned. Mapping table of all triplets to symmetrically
@ -93,12 +93,12 @@ def get_triplets_at_q(grid_point,
q+q'+q''=G. Considering q' is enough because q is fixed and
q''=G-q-q' where G is automatically determined to choose
smallest |G|.
shape=(prod(mesh),), dtype='uintp'
shape=(prod(mesh),), dtype='int_'
map_q : ndarray or None
Returns when stores_triplets_map=True, otherwise None is
returned. Irreducible q-points stabilized by q-point of
specified grid_point.
shape=(prod(mesh),), dtype='uintp'
shape=(prod(mesh),), dtype='int_'
"""
@ -113,6 +113,7 @@ def get_triplets_at_q(grid_point,
mesh,
reciprocal_lattice,
is_dense=True)
bz_map = np.array(bz_map, dtype='int_')
triplets_at_q, weights = _get_BZ_triplets_at_q(
grid_point,
bz_grid_address,
@ -156,6 +157,7 @@ def get_nosym_triplets_at_q(grid_point,
mesh,
reciprocal_lattice,
is_dense=True)
bz_map = np.array(bz_map, dtype='int_')
map_triplets = np.arange(len(grid_address), dtype=bz_map.dtype)
triplets_at_q, weights = _get_BZ_triplets_at_q(
grid_point,
@ -262,7 +264,7 @@ def reduce_grid_points(mesh_divisors,
coarse_mesh_shifts=None):
divisors = np.array(mesh_divisors, dtype='intc')
if (divisors == 1).all():
coarse_grid_points = np.array(dense_grid_points, dtype='uintp')
coarse_grid_points = np.array(dense_grid_points, dtype='int_')
if dense_grid_weights is not None:
coarse_grid_weights = np.array(dense_grid_weights, dtype='intc')
else:
@ -295,7 +297,7 @@ def from_coarse_to_dense_grid_points(dense_mesh,
dense_address = cga * mesh_divisors + shifts * (mesh_divisors // 2)
dense_grid_points.append(get_grid_point_from_address(dense_address,
dense_mesh))
return np.array(dense_grid_points, dtype='uintp')
return np.array(dense_grid_points, dtype='int_')
def get_coarse_ir_grid_points(primitive,
@ -322,7 +324,7 @@ def get_coarse_ir_grid_points(primitive,
if not is_kappa_star:
coarse_grid_address = get_grid_address(coarse_mesh)
coarse_grid_points = np.arange(np.prod(coarse_mesh), dtype='uintp')
coarse_grid_points = np.arange(np.prod(coarse_mesh), dtype='int_')
else:
(coarse_ir_grid_points,
coarse_ir_grid_weights,
@ -450,7 +452,7 @@ def get_tetrahedra_vertices(relative_address,
grid_order = [1, mesh[0], mesh[0] * mesh[1]]
bz_grid_order = [1, bzmesh[0], bzmesh[0] * bzmesh[1]]
num_triplets = len(triplets_at_q)
vertices = np.zeros((num_triplets, 2, 24, 4), dtype='uintp')
vertices = np.zeros((num_triplets, 2, 24, 4), dtype='int_')
for i, tp in enumerate(triplets_at_q):
for j, adrs_shift in enumerate(
(relative_address, -relative_address)):
@ -493,8 +495,8 @@ def _get_triplets_reciprocal_mesh_at_q(fixed_grid_number,
import phono3py._phono3py as phono3c
map_triplets = np.zeros(np.prod(mesh), dtype='uintp')
map_q = np.zeros(np.prod(mesh), dtype='uintp')
map_triplets = np.zeros(np.prod(mesh), dtype='int_')
map_q = np.zeros(np.prod(mesh), dtype='int_')
grid_address = np.zeros((np.prod(mesh), 3), dtype='intc')
phono3c.triplets_reciprocal_mesh_at_q(
@ -517,7 +519,7 @@ def _get_BZ_triplets_at_q(grid_point,
mesh):
import phono3py._phono3py as phono3c
weights = np.zeros(len(map_triplets), dtype='intc')
weights = np.zeros(len(map_triplets), dtype='int_')
for g in map_triplets:
weights[g] += 1
ir_weights = np.extract(weights > 0, weights)
@ -531,7 +533,7 @@ def _get_BZ_triplets_at_q(grid_point,
np.array(mesh, dtype='intc'))
assert num_ir_ret == len(ir_weights)
return triplets, np.array(ir_weights, dtype='intc')
return triplets, np.array(ir_weights, dtype='int_')
def _set_triplets_integration_weights_c(g,
@ -555,7 +557,7 @@ def _set_triplets_integration_weights_c(g,
len(unique_vertices) * len(triplets_at_q), dtype=bz_map.dtype)
phono3c.neighboring_grid_points(
neighboring_grid_points,
np.array(triplets_at_q[:, i], dtype='uintp').ravel(),
np.array(triplets_at_q[:, i], dtype='int_').ravel(),
np.array(j * unique_vertices, dtype='intc', order='C'),
mesh,
grid_address,