Some implementation in C for no-kappa-star option of LBTE

This commit is contained in:
Atsushi Togo 2014-05-15 15:54:51 +09:00
parent b8c6cdeecf
commit 06d9fb5e2c
5 changed files with 157 additions and 10 deletions

View File

@ -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

View File

@ -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

View File

@ -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++) {

View File

@ -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;

View File

@ -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