Separate N and U processes in detailed_imag_self_energy_at_bands

This commit is contained in:
Atsushi Togo 2016-11-08 18:27:57 +09:00
parent cbed08e7b6
commit 5224e1628e
4 changed files with 118 additions and 28 deletions

View File

@ -451,21 +451,25 @@ static PyObject *
py_get_detailed_imag_self_energy_with_g(PyObject *self, PyObject *args)
{
PyArrayObject* gamma_detail_py;
PyArrayObject* gamma_py;
PyArrayObject* gamma_N_py;
PyArrayObject* gamma_U_py;
PyArrayObject* fc3_normal_squared_py;
PyArrayObject* frequencies_py;
PyArrayObject* grid_point_triplets_py;
PyArrayObject* triplet_weights_py;
PyArrayObject* grid_address_py;
PyArrayObject* g_py;
PyArrayObject* g_zero_py;
double unit_conversion_factor, cutoff_frequency, temperature;
if (!PyArg_ParseTuple(args, "OOOOOOdOOdd",
if (!PyArg_ParseTuple(args, "OOOOOOOOdOOdd",
&gamma_detail_py,
&gamma_py,
&gamma_N_py,
&gamma_U_py,
&fc3_normal_squared_py,
&grid_point_triplets_py,
&triplet_weights_py,
&grid_address_py,
&frequencies_py,
&temperature,
&g_py,
@ -477,19 +481,23 @@ py_get_detailed_imag_self_energy_with_g(PyObject *self, PyObject *args)
Darray* fc3_normal_squared = convert_to_darray(fc3_normal_squared_py);
double* gamma_detail = (double*)PyArray_DATA(gamma_detail_py);
double* gamma = (double*)PyArray_DATA(gamma_py);
double* gamma_N = (double*)PyArray_DATA(gamma_N_py);
double* gamma_U = (double*)PyArray_DATA(gamma_U_py);
const double* g = (double*)PyArray_DATA(g_py);
const char* g_zero = (char*)PyArray_DATA(g_zero_py);
const double* frequencies = (double*)PyArray_DATA(frequencies_py);
const int* grid_point_triplets = (int*)PyArray_DATA(grid_point_triplets_py);
const int* triplet_weights = (int*)PyArray_DATA(triplet_weights_py);
const int* grid_address = (int*)PyArray_DATA(grid_address_py);
get_detailed_imag_self_energy_at_bands_with_g(gamma_detail,
gamma,
gamma_N,
gamma_U,
fc3_normal_squared,
frequencies,
grid_point_triplets,
triplet_weights,
grid_address,
g,
g_zero,
temperature,

View File

@ -96,6 +96,15 @@ static void sum_imag_self_energy_along_triplets
const int num_band0,
const int num_triplets,
const double unit_conversion_factor);
static void sum_imag_self_energy_N_and_U_along_triplets
(double *imag_self_energy_N,
double *imag_self_energy_U,
const double *ise,
const int *triplets,
const int *weights,
const int *grid_address,
const int num_band0,
const int num_triplets);
static void set_tmp_values(double *n1,
double *n2,
const int i,
@ -144,11 +153,13 @@ void get_imag_self_energy_at_bands_with_g(double *imag_self_energy,
void get_detailed_imag_self_energy_at_bands_with_g
(double *detailed_imag_self_energy,
double *imag_self_energy,
double *imag_self_energy_N,
double *imag_self_energy_U,
const Darray *fc3_normal_squared,
const double *frequencies,
const int *triplets,
const int *weights,
const int *grid_address,
const double *g,
const char *g_zero,
const double temperature,
@ -174,12 +185,15 @@ void get_detailed_imag_self_energy_at_bands_with_g
unit_conversion_factor,
cutoff_frequency);
sum_imag_self_energy_along_triplets(imag_self_energy,
ise,
weights,
num_band0,
num_triplets,
1);
sum_imag_self_energy_N_and_U_along_triplets
(imag_self_energy_N,
imag_self_energy_U,
ise,
triplets,
weights,
grid_address,
num_band0,
num_triplets);
free(ise);
}
@ -430,6 +444,54 @@ static void sum_imag_self_energy_along_triplets
}
}
static void sum_imag_self_energy_N_and_U_along_triplets
(double *imag_self_energy_N,
double *imag_self_energy_U,
const double *ise,
const int *triplets,
const int *weights,
const int *grid_address,
const int num_band0,
const int num_triplets)
{
int i, j, k, sum_q;
char *is_N;
double sum_g_N, sum_g_U, g;
is_N = (char*)malloc(sizeof(char) * num_triplets);
for (i = 0; i < num_triplets; i++) {
is_N[i] = 1;
for (j = 0; j < 3; j++) {
sum_q = 0;
for (k = 0; k < 3; k++) { /* 1st, 2nd, 3rd triplet */
sum_q += grid_address[triplets[i * 3 + k] * 3 + j];
}
if (sum_q) {
is_N[i] = 0;
break;
}
}
}
for (i = 0; i < num_band0; i++) {
sum_g_N = 0;
sum_g_U = 0;
for (j = 0; j < num_triplets; j++) {
g = ise[j * num_band0 + i] * weights[j];
if (is_N[j]) {
sum_g_N += g;
} else {
sum_g_U += g;
}
}
imag_self_energy_N[i] = sum_g_N;
imag_self_energy_U[i] = sum_g_U;
}
free(is_N);
}
static void set_tmp_values(double *n1,
double *n2,
const int i,

View File

@ -49,11 +49,13 @@ void get_imag_self_energy_at_bands_with_g(double *imag_self_energy,
const double cutoff_frequency);
void get_detailed_imag_self_energy_at_bands_with_g
(double *detailed_imag_self_energy,
double *imag_self_energy,
double *imag_self_energy_N,
double *imag_self_energy_U,
const Darray *fc3_normal_squared,
const double *frequencies,
const int *triplets,
const int *weights,
const int *grid_address,
const double *g,
const char *g_zero,
const double temperature,

View File

@ -115,7 +115,7 @@ def get_imag_self_energy(interaction,
num_band = frequencies.shape[1]
detailed_gamma = np.zeros(
(len(temperatures), len(frequency_points_at_sigma),
num_band0, num_band, num_band, len(weights)),
len(weights), num_band0, num_band, num_band),
dtype='double')
for k, freq_point in enumerate(frequency_points_at_sigma):
@ -129,9 +129,8 @@ def get_imag_self_energy(interaction,
ise.run()
gamma[l, k] = ise.get_imag_self_energy()[0]
if write_detail:
detailed_gamma[l, k] = np.transpose(
ise.get_detailed_imag_self_energy()[0],
axes=(1, 2, 3, 0))
detailed_gamma[l, k] = (
ise.get_detailed_imag_self_energy()[0])
gamma_sigmas.append(gamma)
@ -220,8 +219,9 @@ def get_linewidth(interaction,
num_band0 = len(interaction.get_band_indices())
num_band = frequencies.shape[1]
num_temp = len(temperatures)
num_triplets = len(weights)
detailed_gamma = np.zeros(
(num_temp, num_band0, num_band, num_band, len(weights)),
(num_temp, num_triplets, num_band0, num_band, num_band),
dtype='double')
for k, t in enumerate(temperatures):
@ -229,9 +229,7 @@ def get_linewidth(interaction,
ise.run()
gamma[i, j, k] = ise.get_imag_self_energy()
if write_detail:
detailed_gamma[k] = np.transpose(
ise.get_detailed_imag_self_energy(),
axes=(1, 2, 3, 0))
detailed_gamma[k] = ise.get_detailed_imag_self_energy()
if write_detail:
filename = write_gamma_detail_to_hdf5(
@ -469,13 +467,12 @@ class ImagSelfEnergy(object):
if self._g is not None:
if self._lang == 'C':
if self._with_detail:
# (num_triplets, num_band0, num_band, num_band)
# self._detailed_imag_self_energy.shape =
# (num_triplets, num_band0, num_band, num_band)
# self._imag_self_energy is also set.
self._run_c_detailed_with_band_indices_with_g()
# self._imag_self_energy[:] = np.dot(
# self._weights_at_q,
# self._detailed_imag_self_energy.sum(axis=2).sum(axis=2))
else:
# (num_band0,)
# self._imag_self_energy.shape = (num_band0,)
self._run_c_with_band_indices_with_g()
else:
print("Running into _run_py_with_band_indices_with_g()")
@ -492,15 +489,20 @@ class ImagSelfEnergy(object):
def _run_with_frequency_points(self):
if self._g is not None:
if self._lang == 'C':
self._run_c_with_frequency_points_with_g()
if self._with_detail:
self._run_c_detailed_with_frequency_points_with_g()
else:
self._run_c_with_frequency_points_with_g()
else:
print("Running into _run_py_with_frequency_points_with_g()")
print("This routine is super slow and only for the test.")
self._run_py_with_frequency_points_with_g()
else:
if self._lang == 'C':
self._run_c_with_frequency_points()
else:
print("Running into _run_py_with_frequency_points()")
print("This routine is super slow and only for the test.")
self._run_py_with_frequency_points()
def _run_c_with_band_indices(self):
@ -545,12 +547,17 @@ class ImagSelfEnergy(object):
else:
_g_zero = self._g_zero
imag_self_energy_N = np.zeros_like(self._imag_self_energy)
imag_self_energy_U = np.zeros_like(self._imag_self_energy)
phono3c.detailed_imag_self_energy_with_g(
self._detailed_imag_self_energy,
self._imag_self_energy,
imag_self_energy_N, # Normal
imag_self_energy_U, # Umklapp
self._pp_strength,
self._triplets_at_q,
self._weights_at_q,
self._pp.get_grid_address(),
self._frequencies,
self._temperature,
self._g,
@ -558,6 +565,8 @@ class ImagSelfEnergy(object):
self._unit_conversion,
self._cutoff_frequency)
self._imag_self_energy = imag_self_energy_N + imag_self_energy_U;
def _run_c_with_frequency_points(self):
import phono3py._phono3py as phono3c
ise_at_f = np.zeros(self._imag_self_energy.shape[1], dtype='double')
@ -600,22 +609,31 @@ class ImagSelfEnergy(object):
def _run_c_detailed_with_frequency_points_with_g(self):
import phono3py._phono3py as phono3c
num_band0 = self._pp_strength.shape[1]
g_shape = list(self._g.shape)
g_shape[2] = num_band0
g = np.zeros((2,) + self._pp_strength.shape, dtype='double')
detailed_ise_at_f = np.zeros(self._detailed_imag_self_energy.shape[1:5],
dtype='double')
ise_at_f = np.zeros(num_band0, dtype='double')
_g_zero = np.zeros(g_shape, dtype='byte', order='C')
for i in range(len(self._frequency_points)):
for j in range(g.shape[2]):
g[:, :, j, :, :] = self._g[:, :, i, :, :]
phono3c.detailed_imag_self_energy_with_g(detailed_ise_at_f,
ise_at_f,
self._pp_strength,
self._triplets_at_q,
self._weights_at_q,
self._frequencies,
self._temperature,
g,
_g_zero,
self._unit_conversion,
self._cutoff_frequency)
self._detailed_imag_self_energy[i] = detailed_ise_at_f
self._imag_self_energy[i] = ise_at_f
def _run_py_with_band_indices(self):
for i, (triplet, w, interaction) in enumerate(