diff --git a/anharmonic/phonon3/collision_matrix.py b/anharmonic/phonon3/collision_matrix.py index 7001cc32..8e52541f 100644 --- a/anharmonic/phonon3/collision_matrix.py +++ b/anharmonic/phonon3/collision_matrix.py @@ -107,15 +107,15 @@ class CollisionMatrix(ImagSelfEnergy): self._bz_map = self._interaction.get_bz_map() def _run_collision_matrix(self): - self._run_with_band_indices() # for Gamma + self._run_with_band_indices() if self._temperature > 0: - if self._no_kappa_stars: - self._run_py_collision_matrix() # for Omega - else: - if self._lang == 'C': - self._run_c_collision_matrix() # for Omega + if self._lang == 'C': + if self._no_kappa_stars: + self._run_c_collision_matrix_full() else: - self._run_py_collision_matrix() # for Omega + self._run_c_collision_matrix() + else: + self._run_py_collision_matrix() def _run_c_collision_matrix(self): import anharmonic._phono3py as phono3c @@ -133,6 +133,19 @@ class CollisionMatrix(ImagSelfEnergy): self._unit_conversion, self._cutoff_frequency) + def _run_c_collision_matrix_full(self): + import anharmonic._phono3py as phono3c + phono3c.collision_matrix_full(self._collision_matrix, + self._fc3_normal_squared, + self._frequencies, + self._g, + self._triplets_at_q, + self._triplets_map_at_q, + self._ir_map_at_q, + self._temperature, + self._unit_conversion, + self._cutoff_frequency) + def _run_py_collision_matrix(self): gp2tp_map = {} count = 0 diff --git a/anharmonic/phonon3/conductivity_LBTE.py b/anharmonic/phonon3/conductivity_LBTE.py index 6aad722b..9e980780 100644 --- a/anharmonic/phonon3/conductivity_LBTE.py +++ b/anharmonic/phonon3/conductivity_LBTE.py @@ -495,7 +495,8 @@ class Conductivity_LBTE(Conductivity): col_mat[:, :, :, :, :, i, j, :] = sum_col def _symmetrize_collision_matrix_no_kappa_stars(self): - self._py_symmetrize_collision_matrix_no_kappa_stars() + import anharmonic._phono3py as phono3c + phono3c.symmetrize_collision_matrix(self._collision_matrix) # Average matrix elements belonging to degenerate bands col_mat = self._collision_matrix diff --git a/c/_phono3py.c b/c/_phono3py.c index ef818100..ac311095 100644 --- a/c/_phono3py.c +++ b/c/_phono3py.c @@ -27,6 +27,7 @@ static PyObject * py_get_imag_self_energy_at_bands(PyObject *self, static PyObject * py_get_thm_imag_self_energy_at_bands(PyObject *self, PyObject *args); static PyObject * py_get_collision_matrix(PyObject *self, PyObject *args); +static PyObject * py_get_collision_matrix_full(PyObject *self, PyObject *args); static PyObject * py_symmetrize_collision_matrix(PyObject *self, PyObject *args); static PyObject * py_set_phonons_at_gridpoints(PyObject *self, PyObject *args); static PyObject * py_get_phonon(PyObject *self, PyObject *args); @@ -60,6 +61,7 @@ static PyMethodDef functions[] = { {"imag_self_energy_at_bands", py_get_imag_self_energy_at_bands, METH_VARARGS, "Imaginary part of self energy at phonon frequencies of bands"}, {"thm_imag_self_energy_at_bands", py_get_thm_imag_self_energy_at_bands, METH_VARARGS, "Imaginary part of self energy at phonon frequencies of bands for tetrahedron method"}, {"collision_matrix", py_get_collision_matrix, METH_VARARGS, "Collision matrix with g"}, + {"collision_matrix_full", py_get_collision_matrix_full, METH_VARARGS, "Collision matrix with g for reducible grid points"}, {"symmetrize_collision_matrix", py_symmetrize_collision_matrix, METH_VARARGS, "Symmetrize collision matrix"}, {"phonons_at_gridpoints", py_set_phonons_at_gridpoints, METH_VARARGS, "Set phonons at grid points"}, {"phonon", py_get_phonon, METH_VARARGS, "Get phonon"}, @@ -574,6 +576,56 @@ static PyObject * py_get_collision_matrix(PyObject *self, PyObject *args) Py_RETURN_NONE; } +static PyObject * py_get_collision_matrix_full(PyObject *self, PyObject *args) +{ + PyArrayObject* collision_matrix_py; + PyArrayObject* fc3_normal_squared_py; + PyArrayObject* frequencies_py; + PyArrayObject* triplets_py; + PyArrayObject* triplets_map_py; + PyArrayObject* stabilized_gp_map_py; + PyArrayObject* g_py; + double temperature, unit_conversion_factor, cutoff_frequency; + + if (!PyArg_ParseTuple(args, "OOOOOOOddd", + &collision_matrix_py, + &fc3_normal_squared_py, + &frequencies_py, + &g_py, + &triplets_py, + &triplets_map_py, + &stabilized_gp_map_py, + &temperature, + &unit_conversion_factor, + &cutoff_frequency)) { + return NULL; + } + + Darray* fc3_normal_squared = convert_to_darray(fc3_normal_squared_py); + double* collision_matrix = (double*)collision_matrix_py->data; + const double* g = (double*)g_py->data; + const double* frequencies = (double*)frequencies_py->data; + const int* triplets = (int*)triplets_py->data; + Iarray* triplets_map = convert_to_iarray(triplets_map_py); + const int* stabilized_gp_map = (int*)stabilized_gp_map_py->data; + + get_collision_matrix_full(collision_matrix, + fc3_normal_squared, + frequencies, + triplets, + triplets_map, + stabilized_gp_map, + g, + temperature, + unit_conversion_factor, + cutoff_frequency); + + free(fc3_normal_squared); + free(triplets_map); + + Py_RETURN_NONE; +} + static PyObject * py_symmetrize_collision_matrix(PyObject *self, PyObject *args) { PyArrayObject* collision_matrix_py; @@ -586,12 +638,16 @@ static PyObject * py_symmetrize_collision_matrix(PyObject *self, PyObject *args) double* collision_matrix = (double*)collision_matrix_py->data; const int num_sigma = (int)collision_matrix_py->dimensions[0]; const int num_temp = (int)collision_matrix_py->dimensions[1]; - const int num_ir_grid_points = (int)collision_matrix_py->dimensions[2]; + const int num_grid_points = (int)collision_matrix_py->dimensions[2]; const int num_band = (int)collision_matrix_py->dimensions[3]; int i, j, k, l, num_column, adrs_shift; double val; - num_column = num_ir_grid_points * num_band * 3; + if (collision_matrix_py->nd == 8) { + num_column = num_grid_points * num_band * 3; + } else { + num_column = num_grid_points * num_band; + } for (i = 0; i < num_sigma; i++) { for (j = 0; j < num_temp; j++) { diff --git a/c/anharmonic/phonon3/collision_matrix.c b/c/anharmonic/phonon3/collision_matrix.c index 049719e2..092adc67 100644 --- a/c/anharmonic/phonon3/collision_matrix.c +++ b/c/anharmonic/phonon3/collision_matrix.c @@ -86,6 +86,73 @@ void get_collision_matrix(double *collision_matrix, gp2tp_map = NULL; } +void get_collision_matrix_full(double *collision_matrix, + const Darray *fc3_normal_squared, + const double *frequencies, + const int *triplets, + const Iarray *triplets_map, + const int *stabilized_gp_map, + const double *g, + const double temperature, + const double unit_conversion_factor, + const double cutoff_frequency) +{ + int i, j, k, l, ti, gp2, num_triplets, num_band, num_gp; + int *gp2tp_map; + double f, collision; + double *inv_sinh; + + num_triplets = fc3_normal_squared->dims[0]; + num_band = fc3_normal_squared->dims[2]; + num_gp = triplets_map->dims[0]; + + gp2tp_map = create_gp2tp_map(triplets_map); + +#pragma omp parallel for private(j, k, l, ti, gp2, f, collision, inv_sinh) + for (i = 0; i < num_gp; i++) { + inv_sinh = (double*)malloc(sizeof(double) * num_band); + ti = gp2tp_map[triplets_map->data[i]]; + if (triplets_map->data[i] == stabilized_gp_map[i]) { + gp2 = triplets[ti * 3 + 2]; + } else { + gp2 = triplets[ti * 3 + 1]; + } + for (j = 0; j < num_band; j++) { + f = frequencies[gp2 * num_band + j]; + if (f > cutoff_frequency) { + inv_sinh[j] = inv_sinh_occupation(f, temperature); + } else { + inv_sinh[j] = 0; + } + } + + for (j = 0; j < num_band; j++) { + for (k = 0; k < num_band; k++) { + collision = 0; + for (l = 0; l < num_band; l++) { + collision += + fc3_normal_squared->data[ti * num_band * num_band * num_band + + j * num_band * num_band + + k * num_band + l] * + g[2 * num_triplets * num_band * num_band * num_band + + ti * num_band * num_band * num_band + + j * num_band * num_band + + k * num_band + l] * + inv_sinh[l] * unit_conversion_factor; + } + collision_matrix[j * num_gp * num_band + + i * num_band + k] += collision; + } + } + + free(inv_sinh); + inv_sinh = NULL; + } + + free(gp2tp_map); + gp2tp_map = NULL; +} + static int *create_gp2tp_map(const Iarray *triplets_map) { int i, max_i, count; diff --git a/c/anharmonic_h/phonon3_h/collision_matrix.h b/c/anharmonic_h/phonon3_h/collision_matrix.h index 33a3cdeb..15c2633c 100644 --- a/c/anharmonic_h/phonon3_h/collision_matrix.h +++ b/c/anharmonic_h/phonon3_h/collision_matrix.h @@ -17,4 +17,14 @@ void get_collision_matrix(double *collision_matrix, const double unit_conversion_factor, const double cutoff_frequency); +void get_collision_matrix_full(double *collision_matrix, + const Darray *fc3_normal_squared, + const double *frequencies, + const int *triplets, + const Iarray *triplets_map, + const int *stabilized_gp_map, + const double *g, + const double temperature, + const double unit_conversion_factor, + const double cutoff_frequency); #endif