An option is made: For grid points whose integration weights must be zero from calculation of tetrahedron method, ph-ph interaction calculation is omitted.

This commit is contained in:
Atsushi Togo 2016-03-09 19:45:50 +09:00
parent c545de41e6
commit c6b9d06963
16 changed files with 421 additions and 259 deletions

View File

@ -314,6 +314,7 @@ class Conductivity_RTA(Conductivity):
is_kappa_star=True,
gv_delta_q=None,
run_with_g=True,
with_g_zero=False,
log_level=0):
if sigmas is None:
@ -324,6 +325,7 @@ class Conductivity_RTA(Conductivity):
self._is_kappa_star = None
self._gv_delta_q = None
self._run_with_g = run_with_g
self._with_g_zero = with_g_zero
self._log_level = None
self._primitive = None
self._dm = None
@ -451,7 +453,13 @@ class Conductivity_RTA(Conductivity):
len(self._pp.get_triplets_at_q()[0]))
print("Calculating interaction...")
self._collision.run_interaction()
if (self._with_g_zero and
len(self._sigmas) == 1 and
self._sigmas[0] is None):
self._collision.set_sigma(None)
self._collision.set_integration_weights()
self._collision.run_interaction(with_g_zero=self._with_g_zero)
if not self._with_g_zero:
self._averaged_pp_interaction[i] = (
self._pp.get_averaged_interaction())
self._set_gamma_at_sigmas(i)
@ -506,8 +514,15 @@ class Conductivity_RTA(Conductivity):
else:
text += "sigma=%s" % sigma
print(text)
if (self._with_g_zero and
len(self._sigmas) == 1 and
self._sigmas[0] is None and
not self._use_ave_pp):
pass
else:
self._collision.set_sigma(sigma)
if not sigma or self._run_with_g:
if sigma is None or self._run_with_g:
self._collision.set_integration_weights()
for k, t in enumerate(self._temperatures):
self._collision.set_temperature(t)
@ -555,7 +570,16 @@ class Conductivity_RTA(Conductivity):
gv = self._gv[i]
ave_pp = self._averaged_pp_interaction[i]
show_Pqj = not (self._with_g_zero and
len(self._sigmas) == 1 and
self._sigmas[0] is None and
not self._use_ave_pp)
if show_Pqj:
text = "Frequency group velocity (x, y, z) |gv| Pqj"
else:
text = "Frequency group velocity (x, y, z) |gv|"
if self._gv_delta_q is None:
pass
else:
@ -578,10 +602,18 @@ class Conductivity_RTA(Conductivity):
for f, v, pp in zip(frequencies,
np.dot(rot_c, gv.T).T,
ave_pp):
if show_Pqj:
print("%8.3f (%8.3f %8.3f %8.3f) %8.3f %11.3e" %
(f, v[0], v[1], v[2], np.linalg.norm(v), pp))
else:
print("%8.3f (%8.3f %8.3f %8.3f) %8.3f" %
(f, v[0], v[1], v[2], np.linalg.norm(v)))
print('')
else:
for f, v, pp in zip(frequencies, gv, ave_pp):
if show_Pqj:
print("%8.3f (%8.3f %8.3f %8.3f) %8.3f %11.3e" %
(f, v[0], v[1], v[2], np.linalg.norm(v), pp))
else:
print("%8.3f (%8.3f %8.3f %8.3f) %8.3f" %
(f, v[0], v[1], v[2], np.linalg.norm(v)))

View File

@ -331,6 +331,7 @@ class ImagSelfEnergy:
self._cutoff_frequency = interaction.get_cutoff_frequency()
self._g = None # integration weights
self._g_zero = None # if calculate ph-ph interaction or not
self._mesh = self._pp.get_mesh_numbers()
self._band_indices = self._pp.get_band_indices()
self._is_collision_matrix = False
@ -363,11 +364,12 @@ class ImagSelfEnergy:
dtype='double')
self._run_with_frequency_points()
def run_interaction(self):
def run_interaction(self, with_g_zero=False):
if self._frequency_points is None and with_g_zero:
self._pp.run(g_zero=self._g_zero, lang=self._lang)
else:
self._pp.run(lang=self._lang)
self._pp_strength = self._pp.get_interaction_strength()
(self._frequencies,
self._eigenvectors) = self._pp.get_phonons()[:2]
def set_integration_weights(self, scattering_event_class=None):
if self._frequency_points is None:
@ -375,7 +377,7 @@ class ImagSelfEnergy:
else:
f_points = self._frequency_points
self._g = get_triplets_integration_weights(
self._g, self._g_zero = get_triplets_integration_weights(
self._pp,
np.array(f_points, dtype='double'),
self._sigma,
@ -419,6 +421,7 @@ class ImagSelfEnergy:
(self._triplets_at_q,
self._weights_at_q) = self._pp.get_triplets_at_q()[:2]
self._grid_point = grid_point
self._frequencies, self._eigenvectors, _ = self._pp.get_phonons()
def set_sigma(self, sigma):
if sigma is None:
@ -427,6 +430,7 @@ class ImagSelfEnergy:
self._sigma = float(sigma)
self._g = None
self._g_zero = None
def set_frequency_points(self, frequency_points):
if frequency_points is None:

View File

@ -72,7 +72,7 @@ class Interaction:
self._allocate_phonon()
def run(self, lang='C'):
def run(self, g_zero=None, lang='C'):
num_band = self._primitive.get_number_of_atoms() * 3
num_triplets = len(self._triplets_at_q)
@ -81,7 +81,7 @@ class Interaction:
(num_triplets, len(self._band_indices), num_band, num_band),
dtype='double')
if lang == 'C':
self._run_c()
self._run_c(g_zero=g_zero)
else:
self._run_py()
else:
@ -89,7 +89,6 @@ class Interaction:
self._interaction_strength = np.ones(
(num_triplets, len(self._band_indices), num_band, num_band),
dtype='double') * self._constant_averaged_interaction / num_grid
self.set_phonons(self._triplets_at_q.ravel())
def get_interaction_strength(self):
return self._interaction_strength
@ -185,6 +184,7 @@ class Interaction:
self._grid_address = grid_address
self._bz_map = bz_map
self._ir_map_at_q = ir_map_at_q
self.set_phonons(self._triplets_at_q.ravel())
def set_dynamical_matrix(self,
fc2,
@ -239,10 +239,9 @@ class Interaction:
num_band = self._primitive.get_number_of_atoms() * 3
return np.dot(w, v_sum) / num_band ** 2
def _run_c(self):
def _run_c(self, g_zero=None):
import anharmonic._phono3py as phono3c
self.set_phonons(self._triplets_at_q.ravel())
num_band = self._primitive.get_number_of_atoms() * 3
svecs, multiplicity = get_smallest_vectors(self._supercell,
self._primitive,
@ -251,7 +250,14 @@ class Interaction:
p2s = self._primitive.get_primitive_to_supercell_map()
s2p = self._primitive.get_supercell_to_primitive_map()
if g_zero is None:
_g_zero = np.zeros(self._interaction_strength.shape,
dtype='byte', order='C')
else:
_g_zero = g_zero
phono3c.interaction(self._interaction_strength,
_g_zero,
self._frequencies,
self._eigenvectors,
self._triplets_at_q,

View File

@ -167,7 +167,7 @@ class JointDos:
occupation(freqs, t), 0))
for i, freq_point in enumerate(self._frequency_points):
g = get_triplets_integration_weights(
g, _ = get_triplets_integration_weights(
self,
np.array([freq_point], dtype='double'),
self._sigma,

View File

@ -282,6 +282,7 @@ def get_triplets_integration_weights(interaction,
triplets = interaction.get_triplets_at_q()[0]
frequencies = interaction.get_phonons()[0]
num_band = frequencies.shape[1]
g_zero = None
if is_collision_matrix:
g = np.zeros(
@ -317,7 +318,7 @@ def get_triplets_integration_weights(interaction,
g[2, i, :, j, k] = g0 + g1 + g2
else:
if lang == 'C':
_set_triplets_integration_weights_c(
g_zero = _set_triplets_integration_weights_c(
g,
interaction,
frequency_points,
@ -326,7 +327,7 @@ def get_triplets_integration_weights(interaction,
_set_triplets_integration_weights_py(
g, interaction, frequency_points)
return g
return g, g_zero
def get_tetrahedra_vertices(relative_address,
mesh,
@ -420,8 +421,11 @@ def _set_triplets_integration_weights_c(g,
bz_map)
interaction.set_phonons(np.unique(neighboring_grid_points))
g_zero = np.zeros(g.shape[1:], dtype='byte', order='C')
phono3c.triplets_integration_weights(
g,
g_zero,
frequency_points,
thm.get_tetrahedra(),
mesh,
@ -430,6 +434,8 @@ def _set_triplets_integration_weights_c(g,
grid_address,
bz_map)
return g_zero
def _set_triplets_integration_weights_py(g, interaction, frequency_points):
reciprocal_lattice = np.linalg.inv(interaction.get_primitive().get_cell())
mesh = interaction.get_mesh_numbers()

View File

@ -91,7 +91,26 @@ static PyObject * py_inverse_collision_matrix(PyObject *self, PyObject *args);
static PyObject * py_inverse_collision_matrix_libflame(PyObject *self, PyObject *args);
#endif
static PyMethodDef functions[] = {
struct module_state {
PyObject *error;
};
#if PY_MAJOR_VERSION >= 3
#define GETSTATE(m) ((struct module_state*)PyModule_GetState(m))
#else
#define GETSTATE(m) (&_state)
static struct module_state _state;
#endif
static PyObject *
error_out(PyObject *m) {
struct module_state *st = GETSTATE(m);
PyErr_SetString(st->error, "something bad happened");
return NULL;
}
static PyMethodDef _phono3py_methods[] = {
{"error_out", (PyCFunction)error_out, METH_NOARGS, NULL},
{"interaction", py_get_interaction, METH_VARARGS, "Interaction of triplets"},
{"imag_self_energy", py_get_imag_self_energy, METH_VARARGS,
"Imaginary part of self energy at arbitrary frequency points"},
@ -141,15 +160,67 @@ static PyMethodDef functions[] = {
{NULL, NULL, 0, NULL}
};
PyMODINIT_FUNC init_phono3py(void)
#if PY_MAJOR_VERSION >= 3
static int _phono3py_traverse(PyObject *m, visitproc visit, void *arg) {
Py_VISIT(GETSTATE(m)->error);
return 0;
}
static int _phono3py_clear(PyObject *m) {
Py_CLEAR(GETSTATE(m)->error);
return 0;
}
static struct PyModuleDef moduledef = {
PyModuleDef_HEAD_INIT,
"_phono3py",
NULL,
sizeof(struct module_state),
_phono3py_methods,
NULL,
_phono3py_traverse,
_phono3py_clear,
NULL
};
#define INITERROR return NULL
PyObject *
PyInit__phono3py(void)
#else
#define INITERROR return
void
init_phono3py(void)
#endif
{
Py_InitModule3("_phono3py", functions, "C-extension for phono3py\n\n...\n");
return;
#if PY_MAJOR_VERSION >= 3
PyObject *module = PyModule_Create(&moduledef);
#else
PyObject *module = Py_InitModule("_phono3py", _phono3py_methods);
#endif
if (module == NULL)
INITERROR;
struct module_state *st = GETSTATE(module);
st->error = PyErr_NewException("_phono3py.Error", NULL, NULL);
if (st->error == NULL) {
Py_DECREF(module);
INITERROR;
}
#if PY_MAJOR_VERSION >= 3
return module;
#endif
}
static PyObject * py_get_interaction(PyObject *self, PyObject *args)
{
PyArrayObject* fc3_normal_squared_py;
PyArrayObject* g_zero_py;
PyArrayObject* frequencies;
PyArrayObject* eigenvectors;
PyArrayObject* grid_point_triplets;
@ -165,8 +236,9 @@ static PyObject * py_get_interaction(PyObject *self, PyObject *args)
double cutoff_frequency;
int symmetrize_fc3_q;
if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOid",
if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOOid",
&fc3_normal_squared_py,
&g_zero_py,
&frequencies,
&eigenvectors,
&grid_point_triplets,
@ -191,17 +263,19 @@ static PyObject * py_get_interaction(PyObject *self, PyObject *args)
/* So eigenvectors should not be used in Python side */
Carray* eigvecs = convert_to_carray(eigenvectors);
Iarray* triplets = convert_to_iarray(grid_point_triplets);
const int* grid_address = (int*)grid_address_py->data;
const int* mesh = (int*)mesh_py->data;
const char* g_zero = (char*)PyArray_DATA(g_zero_py);
const int* grid_address = (int*)PyArray_DATA(grid_address_py);
const int* mesh = (int*)PyArray_DATA(mesh_py);
Darray* fc3 = convert_to_darray(fc3_py);
Darray* svecs = convert_to_darray(shortest_vectors);
Iarray* multi = convert_to_iarray(multiplicity);
const double* masses = (double*)atomic_masses->data;
const int* p2s = (int*)p2s_map->data;
const int* s2p = (int*)s2p_map->data;
const int* band_indicies = (int*)band_indicies_py->data;
const double* masses = (double*)PyArray_DATA(atomic_masses);
const int* p2s = (int*)PyArray_DATA(p2s_map);
const int* s2p = (int*)PyArray_DATA(s2p_map);
const int* band_indicies = (int*)PyArray_DATA(band_indicies_py);
get_interaction(fc3_normal_squared,
g_zero,
freqs,
eigvecs,
triplets,
@ -254,10 +328,10 @@ static PyObject * py_get_imag_self_energy(PyObject *self, PyObject *args)
Darray* fc3_normal_squared = convert_to_darray(fc3_normal_squared_py);
double* gamma = (double*)gamma_py->data;
const double* frequencies = (double*)frequencies_py->data;
const int* grid_point_triplets = (int*)grid_point_triplets_py->data;
const int* triplet_weights = (int*)triplet_weights_py->data;
double* gamma = (double*)PyArray_DATA(gamma_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);
get_imag_self_energy(gamma,
fc3_normal_squared,
@ -302,11 +376,11 @@ static PyObject * py_get_imag_self_energy_at_bands(PyObject *self,
Darray* fc3_normal_squared = convert_to_darray(fc3_normal_squared_py);
double* gamma = (double*)gamma_py->data;
const double* frequencies = (double*)frequencies_py->data;
const int* band_indices = (int*)band_indices_py->data;
const int* grid_point_triplets = (int*)grid_point_triplets_py->data;
const int* triplet_weights = (int*)triplet_weights_py->data;
double* gamma = (double*)PyArray_DATA(gamma_py);
const double* frequencies = (double*)PyArray_DATA(frequencies_py);
const int* band_indices = (int*)PyArray_DATA(band_indices_py);
const int* grid_point_triplets = (int*)PyArray_DATA(grid_point_triplets_py);
const int* triplet_weights = (int*)PyArray_DATA(triplet_weights_py);
get_imag_self_energy_at_bands(gamma,
fc3_normal_squared,
@ -348,11 +422,11 @@ static PyObject * py_get_imag_self_energy_with_g(PyObject *self, PyObject *args)
}
Darray* fc3_normal_squared = convert_to_darray(fc3_normal_squared_py);
double* gamma = (double*)gamma_py->data;
const double* g = (double*)g_py->data;
const double* frequencies = (double*)frequencies_py->data;
const int* grid_point_triplets = (int*)grid_point_triplets_py->data;
const int* triplet_weights = (int*)triplet_weights_py->data;
double* gamma = (double*)PyArray_DATA(gamma_py);
const double* g = (double*)PyArray_DATA(g_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);
get_imag_self_energy_at_bands_with_g(gamma,
fc3_normal_squared,
@ -392,10 +466,10 @@ 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 = (double*)gamma_py->data;
const double* g = (double*)g_py->data;
const double* frequencies = (double*)frequencies_py->data;
const int* grid_point_triplets = (int*)grid_point_triplets_py->data;
double* gamma = (double*)PyArray_DATA(gamma_py);
const double* g = (double*)PyArray_DATA(g_py);
const double* frequencies = (double*)PyArray_DATA(frequencies_py);
const int* grid_point_triplets = (int*)PyArray_DATA(grid_point_triplets_py);
get_detailed_imag_self_energy_at_bands_with_g(gamma,
fc3_normal_squared,
@ -438,11 +512,11 @@ static PyObject * py_get_frequency_shift_at_bands(PyObject *self,
Darray* fc3_normal_squared = convert_to_darray(fc3_normal_squared_py);
double* shift = (double*)shift_py->data;
const double* frequencies = (double*)frequencies_py->data;
const int* band_indices = (int*)band_indices_py->data;
const int* grid_point_triplets = (int*)grid_point_triplets_py->data;
const int* triplet_weights = (int*)triplet_weights_py->data;
double* shift = (double*)PyArray_DATA(shift_py);
const double* frequencies = (double*)PyArray_DATA(frequencies_py);
const int* band_indices = (int*)PyArray_DATA(band_indices_py);
const int* grid_point_triplets = (int*)PyArray_DATA(grid_point_triplets_py);
const int* triplet_weights = (int*)PyArray_DATA(triplet_weights_py);
get_frequency_shift_at_bands(shift,
fc3_normal_squared,
@ -492,15 +566,16 @@ static PyObject * py_get_collision_matrix(PyObject *self, PyObject *args)
}
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;
double* collision_matrix = (double*)PyArray_DATA(collision_matrix_py);
const double* g = (double*)PyArray_DATA(g_py);
const double* frequencies = (double*)PyArray_DATA(frequencies_py);
const int* triplets = (int*)PyArray_DATA(triplets_py);
Iarray* triplets_map = convert_to_iarray(triplets_map_py);
const int* stabilized_gp_map = (int*)stabilized_gp_map_py->data;
const int* ir_grid_points = (int*)ir_grid_points_py->data;
const int* stabilized_gp_map = (int*)PyArray_DATA(stabilized_gp_map_py);
const int* ir_grid_points = (int*)PyArray_DATA(ir_grid_points_py);
Iarray* rotated_grid_points = convert_to_iarray(rotated_grid_points_py);
const double* rotations_cartesian = (double*)rotations_cartesian_py->data;
const double* rotations_cartesian =
(double*)PyArray_DATA(rotations_cartesian_py);
get_collision_matrix(collision_matrix,
fc3_normal_squared,
@ -549,12 +624,12 @@ static PyObject * py_get_reducible_collision_matrix(PyObject *self, PyObject *ar
}
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;
double* collision_matrix = (double*)PyArray_DATA(collision_matrix_py);
const double* g = (double*)PyArray_DATA(g_py);
const double* frequencies = (double*)PyArray_DATA(frequencies_py);
const int* triplets = (int*)PyArray_DATA(triplets_py);
Iarray* triplets_map = convert_to_iarray(triplets_map_py);
const int* stabilized_gp_map = (int*)stabilized_gp_map_py->data;
const int* stabilized_gp_map = (int*)PyArray_DATA(stabilized_gp_map_py);
get_reducible_collision_matrix(collision_matrix,
fc3_normal_squared,
@ -582,11 +657,11 @@ static PyObject * py_symmetrize_collision_matrix(PyObject *self, PyObject *args)
return NULL;
}
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_grid_points = (int)collision_matrix_py->dimensions[2];
const int num_band = (int)collision_matrix_py->dimensions[3];
double* collision_matrix = (double*)PyArray_DATA(collision_matrix_py);
const int num_sigma = PyArray_DIMS(collision_matrix_py)[0];
const int num_temp = PyArray_DIMS(collision_matrix_py)[1];
const int num_grid_points = PyArray_DIMS(collision_matrix_py)[2];
const int num_band = PyArray_DIMS(collision_matrix_py)[3];
int i, j, k, l, num_column, adrs_shift;
double val;
@ -641,14 +716,14 @@ static PyObject * py_get_isotope_strength(PyObject *self, PyObject *args)
}
double* gamma = (double*)gamma_py->data;
const double* frequencies = (double*)frequencies_py->data;
double* gamma = (double*)PyArray_DATA(gamma_py);
const double* frequencies = (double*)PyArray_DATA(frequencies_py);
const lapack_complex_double* eigenvectors =
(lapack_complex_double*)eigenvectors_py->data;
const int* band_indices = (int*)band_indices_py->data;
const double* mass_variances = (double*)mass_variances_py->data;
const int num_band = (int)frequencies_py->dimensions[1];
const int num_band0 = (int)band_indices_py->dimensions[0];
(lapack_complex_double*)PyArray_DATA(eigenvectors_py);
const int* band_indices = (int*)PyArray_DATA(band_indices_py);
const double* mass_variances = (double*)PyArray_DATA(mass_variances_py);
const int num_band = PyArray_DIMS(frequencies_py)[1];
const int num_band0 = PyArray_DIMS(band_indices_py)[0];
/* int i, j, k; */
/* double f, f0; */
@ -735,18 +810,19 @@ static PyObject * py_get_thm_isotope_strength(PyObject *self, PyObject *args)
}
double* gamma = (double*)gamma_py->data;
const double* frequencies = (double*)frequencies_py->data;
const int* ir_grid_points = (int*)ir_grid_points_py->data;
const int* weights = (int*)weights_py->data;
double* gamma = (double*)PyArray_DATA(gamma_py);
const double* frequencies = (double*)PyArray_DATA(frequencies_py);
const int* ir_grid_points = (int*)PyArray_DATA(ir_grid_points_py);
const int* weights = (int*)PyArray_DATA(weights_py);
const lapack_complex_double* eigenvectors =
(lapack_complex_double*)eigenvectors_py->data;
const int* band_indices = (int*)band_indices_py->data;
const double* mass_variances = (double*)mass_variances_py->data;
const int num_band = (int)frequencies_py->dimensions[1];
const int num_band0 = (int)band_indices_py->dimensions[0];
const double* integration_weights = (double*)integration_weights_py->data;
const int num_ir_grid_points = (int)ir_grid_points_py->dimensions[0];
(lapack_complex_double*)PyArray_DATA(eigenvectors_py);
const int* band_indices = (int*)PyArray_DATA(band_indices_py);
const double* mass_variances = (double*)PyArray_DATA(mass_variances_py);
const int num_band = PyArray_DIMS(frequencies_py)[1];
const int num_band0 = PyArray_DIMS(band_indices_py)[0];
const double* integration_weights =
(double*)PyArray_DATA(integration_weights_py);
const int num_ir_grid_points = PyArray_DIMS(ir_grid_points_py)[0];
get_thm_isotope_scattering_strength(gamma,
grid_point,
@ -782,11 +858,11 @@ static PyObject * py_distribute_fc3(PyObject *self, PyObject *args)
return NULL;
}
double* fc3_copy = (double*)force_constants_third_copy->data;
const double* fc3 = (double*)force_constants_third->data;
const double* rot_cart_inv = (double*)rotation_cart_inv->data;
const int* atom_mapping = (int*)atom_mapping_py->data;
const int num_atom = (int)atom_mapping_py->dimensions[0];
double* fc3_copy = (double*)PyArray_DATA(force_constants_third_copy);
const double* fc3 = (double*)PyArray_DATA(force_constants_third);
const double* rot_cart_inv = (double*)PyArray_DATA(rotation_cart_inv);
const int* atom_mapping = (int*)PyArray_DATA(atom_mapping_py);
const int num_atom = PyArray_DIMS(atom_mapping_py)[0];
distribute_fc3(fc3_copy,
fc3,
@ -807,8 +883,8 @@ static PyObject * py_set_permutation_symmetry_fc3(PyObject *self, PyObject *args
return NULL;
}
double* fc3 = (double*)fc3_py->data;
const int num_atom = (int)fc3_py->dimensions[0];
double* fc3 = (double*)PyArray_DATA(fc3_py);
const int num_atom = PyArray_DIMS(fc3_py)[0];
set_permutation_symmetry_fc3(fc3, num_atom);
@ -833,15 +909,17 @@ static PyObject * py_get_neighboring_gird_points(PyObject *self, PyObject *args)
return NULL;
}
int* relative_grid_points = (int*)relative_grid_points_py->data;
const int *grid_points = (int*)grid_points_py->data;
const int num_grid_points = (int)grid_points_py->dimensions[0];
int* relative_grid_points = (int*)PyArray_DATA(relative_grid_points_py);
const int *grid_points = (int*)PyArray_DATA(grid_points_py);
const int num_grid_points = PyArray_DIMS(grid_points_py)[0];
PHPYCONST int (*relative_grid_address)[3] =
(int(*)[3])relative_grid_address_py->data;
const int num_relative_grid_address = relative_grid_address_py->dimensions[0];
const int *mesh = (int*)mesh_py->data;
PHPYCONST int (*bz_grid_address)[3] = (int(*)[3])bz_grid_address_py->data;
const int *bz_map = (int*)bz_map_py->data;
(int(*)[3])PyArray_DATA(relative_grid_address_py);
const int num_relative_grid_address =
PyArray_DIMS(relative_grid_address_py)[0];
const int *mesh = (int*)PyArray_DATA(mesh_py);
PHPYCONST int (*bz_grid_address)[3] =
(int(*)[3])PyArray_DATA(bz_grid_address_py);
const int *bz_map = (int*)PyArray_DATA(bz_map_py);
int i;
#pragma omp parallel for
@ -881,18 +959,19 @@ static PyObject * py_set_integration_weights(PyObject *self, PyObject *args)
return NULL;
}
double *iw = (double*)iw_py->data;
const double *frequency_points = (double*)frequency_points_py->data;
const int num_band0 = frequency_points_py->dimensions[0];
double *iw = (double*)PyArray_DATA(iw_py);
const double *frequency_points = (double*)PyArray_DATA(frequency_points_py);
const int num_band0 = PyArray_DIMS(frequency_points_py)[0];
PHPYCONST int (*relative_grid_address)[4][3] =
(int(*)[4][3])relative_grid_address_py->data;
const int *mesh = (int*)mesh_py->data;
PHPYCONST int *grid_points = (int*)grid_points_py->data;
const int num_gp = (int)grid_points_py->dimensions[0];
PHPYCONST int (*bz_grid_address)[3] = (int(*)[3])bz_grid_address_py->data;
const int *bz_map = (int*)bz_map_py->data;
const double *frequencies = (double*)frequencies_py->data;
const int num_band = (int)frequencies_py->dimensions[1];
(int(*)[4][3])PyArray_DATA(relative_grid_address_py);
const int *mesh = (int*)PyArray_DATA(mesh_py);
PHPYCONST int *grid_points = (int*)PyArray_DATA(grid_points_py);
const int num_gp = PyArray_DIMS(grid_points_py)[0];
PHPYCONST int (*bz_grid_address)[3] =
(int(*)[3])PyArray_DATA(bz_grid_address_py);
const int *bz_map = (int*)PyArray_DATA(bz_map_py);
const double *frequencies = (double*)PyArray_DATA(frequencies_py);
const int num_band = PyArray_DIMS(frequencies_py)[1];
int i, j, k, bi;
int vertices[24][4];
@ -946,13 +1025,13 @@ py_tpl_get_triplets_reciprocal_mesh_at_q(PyObject *self, PyObject *args)
return NULL;
}
int (*grid_address)[3] = (int(*)[3])grid_address_py->data;
int *map_triplets_int = (int*)map_triplets->data;
int *map_q_int = (int*)map_q->data;
int (*grid_address)[3] = (int(*)[3])PyArray_DATA(grid_address_py);
int *map_triplets_int = (int*)PyArray_DATA(map_triplets);
int *map_q_int = (int*)PyArray_DATA(map_q);
const int* mesh_int = (int*)mesh->data;
PHPYCONST int (*rot)[3][3] = (int(*)[3][3])rotations->data;
const int num_rot = rotations->dimensions[0];
const int* mesh_int = (int*)PyArray_DATA(mesh);
PHPYCONST int (*rot)[3][3] = (int(*)[3][3])PyArray_DATA(rotations);
const int num_rot = PyArray_DIMS(rotations)[0];
const int num_ir =
tpl_get_triplets_reciprocal_mesh_at_q(map_triplets_int,
map_q_int,
@ -985,12 +1064,13 @@ static PyObject * py_tpl_get_BZ_triplets_at_q(PyObject *self, PyObject *args)
return NULL;
}
int (*triplets)[3] = (int(*)[3])triplets_py->data;
PHPYCONST int (*bz_grid_address)[3] = (int(*)[3])bz_grid_address_py->data;
const int *bz_map = (int*)bz_map_py->data;
const int *map_triplets = (int*)map_triplets_py->data;
const int num_map_triplets = (int)map_triplets_py->dimensions[0];
const int *mesh = (int*)mesh_py->data;
int (*triplets)[3] = (int(*)[3])PyArray_DATA(triplets_py);
PHPYCONST int (*bz_grid_address)[3] =
(int(*)[3])PyArray_DATA(bz_grid_address_py);
const int *bz_map = (int*)PyArray_DATA(bz_map_py);
const int *map_triplets = (int*)PyArray_DATA(map_triplets_py);
const int num_map_triplets = PyArray_DIMS(map_triplets_py)[0];
const int *mesh = (int*)PyArray_DATA(mesh_py);
int num_ir;
num_ir = tpl_get_BZ_triplets_at_q(triplets,
@ -1008,6 +1088,7 @@ static PyObject *
py_set_triplets_integration_weights(PyObject *self, PyObject *args)
{
PyArrayObject* iw_py;
PyArrayObject* iw_zero_py;
PyArrayObject* frequency_points_py;
PyArrayObject* relative_grid_address_py;
PyArrayObject* mesh_py;
@ -1015,8 +1096,9 @@ py_set_triplets_integration_weights(PyObject *self, PyObject *args)
PyArrayObject* frequencies_py;
PyArrayObject* bz_grid_address_py;
PyArrayObject* bz_map_py;
if (!PyArg_ParseTuple(args, "OOOOOOOO",
if (!PyArg_ParseTuple(args, "OOOOOOOOO",
&iw_py,
&iw_zero_py,
&frequency_points_py,
&relative_grid_address_py,
&mesh_py,
@ -1027,21 +1109,23 @@ py_set_triplets_integration_weights(PyObject *self, PyObject *args)
return NULL;
}
double *iw = (double*)iw_py->data;
const double *frequency_points = (double*)frequency_points_py->data;
const int num_band0 = frequency_points_py->dimensions[0];
double *iw = (double*)PyArray_DATA(iw_py);
char *iw_zero = (char*)PyArray_DATA(iw_zero_py);
const double *frequency_points = (double*)PyArray_DATA(frequency_points_py);
const int num_band0 = PyArray_DIMS(frequency_points_py)[0];
PHPYCONST int (*relative_grid_address)[4][3] =
(int(*)[4][3])relative_grid_address_py->data;
const int *mesh = (int*)mesh_py->data;
PHPYCONST int (*triplets)[3] = (int(*)[3])triplets_py->data;
const int num_triplets = (int)triplets_py->dimensions[0];
PHPYCONST int (*bz_grid_address)[3] = (int(*)[3])bz_grid_address_py->data;
const int *bz_map = (int*)bz_map_py->data;
const double *frequencies = (double*)frequencies_py->data;
const int num_band = (int)frequencies_py->dimensions[1];
const int num_iw = (int)iw_py->dimensions[0];
(int(*)[4][3])PyArray_DATA(relative_grid_address_py);
const int *mesh = (int*)PyArray_DATA(mesh_py);
PHPYCONST int (*triplets)[3] = (int(*)[3])PyArray_DATA(triplets_py);
const int num_triplets = PyArray_DIMS(triplets_py)[0];
PHPYCONST int (*bz_grid_address)[3] = (int(*)[3])PyArray_DATA(bz_grid_address_py);
const int *bz_map = (int*)PyArray_DATA(bz_map_py);
const double *frequencies = (double*)PyArray_DATA(frequencies_py);
const int num_band = PyArray_DIMS(frequencies_py)[1];
const int num_iw = PyArray_DIMS(iw_py)[0];
tpl_get_integration_weight(iw,
iw_zero,
frequency_points,
num_band0,
relative_grid_address,
@ -1075,14 +1159,14 @@ py_set_triplets_integration_weights_with_sigma(PyObject *self, PyObject *args)
return NULL;
}
double *iw = (double*)iw_py->data;
const double *frequency_points = (double*)frequency_points_py->data;
const int num_band0 = frequency_points_py->dimensions[0];
PHPYCONST int (*triplets)[3] = (int(*)[3])triplets_py->data;
const int num_triplets = (int)triplets_py->dimensions[0];
const double *frequencies = (double*)frequencies_py->data;
const int num_band = (int)frequencies_py->dimensions[1];
const int num_iw = (int)iw_py->dimensions[0];
double *iw = (double*)PyArray_DATA(iw_py);
const double *frequency_points = (double*)PyArray_DATA(frequency_points_py);
const int num_band0 = PyArray_DIMS(frequency_points_py)[0];
PHPYCONST int (*triplets)[3] = (int(*)[3])PyArray_DATA(triplets_py);
const int num_triplets = PyArray_DIMS(triplets_py)[0];
const double *frequencies = (double*)PyArray_DATA(frequencies_py);
const int num_band = PyArray_DIMS(frequencies_py)[1];
const int num_iw = PyArray_DIMS(iw_py)[0];
int i, j, k, l, adrs_shift;
double f0, f1, f2, g0, g1, g2;
@ -1133,11 +1217,11 @@ static PyObject * py_inverse_collision_matrix_libflame(PyObject *self, PyObject
}
double* collision_matrix = (double*)collision_matrix_py->data;
double* eigvals = (double*)eigenvalues_py->data;
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_band = (int)collision_matrix_py->dimensions[3];
double* collision_matrix = (double*)PyArray_DATA(collision_matrix_py);
double* eigvals = (double*)PyArray_DATA(eigenvalues_py);
const int num_temp = PyArray_DIMS(collision_matrix_py)[1];
const int num_ir_grid_points = PyArray_DIMS(collision_matrix_py)[2];
const int num_band = PyArray_DIMS(collision_matrix_py)[3];
int num_column, adrs_shift;
num_column = num_ir_grid_points * num_band * 3;
@ -1167,11 +1251,11 @@ static PyObject * py_inverse_collision_matrix(PyObject *self, PyObject *args)
return NULL;
}
double* collision_matrix = (double*)collision_matrix_py->data;
double* eigvals = (double*)eigenvalues_py->data;
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_band = (int)collision_matrix_py->dimensions[3];
double* collision_matrix = (double*)PyArray_DATA(collision_matrix_py);
double* eigvals = (double*)PyArray_DATA(eigenvalues_py);
const int num_temp = PyArray_DIMS(collision_matrix_py)[1];
const int num_ir_grid_points = PyArray_DIMS(collision_matrix_py)[2];
const int num_band = PyArray_DIMS(collision_matrix_py)[3];
int num_column, adrs_shift, info;
num_column = num_ir_grid_points * num_band * 3;

View File

@ -39,8 +39,7 @@
#include <derivative_dynmat.h>
#include <kgrid.h>
#include <tetrahedron_method.h>
#define KB 8.6173382568083159E-05
#include <phonoc_const.h>
/* Build dynamical matrix */
static PyObject * py_get_dynamical_matrix(PyObject *self, PyObject *args);
@ -591,11 +590,11 @@ static PyObject *py_thm_neighboring_grid_points(PyObject *self, PyObject *args)
}
int* relative_grid_points = (int*)PyArray_DATA(relative_grid_points_py);
THMCONST int (*relative_grid_address)[3] =
PHPYCONST int (*relative_grid_address)[3] =
(int(*)[3])PyArray_DATA(relative_grid_address_py);
const int num_relative_grid_address = PyArray_DIMS(relative_grid_address_py)[0];
const int *mesh = (int*)PyArray_DATA(mesh_py);
THMCONST int (*bz_grid_address)[3] = (int(*)[3])PyArray_DATA(bz_grid_address_py);
PHPYCONST int (*bz_grid_address)[3] = (int(*)[3])PyArray_DATA(bz_grid_address_py);
const int *bz_map = (int*)PyArray_DATA(bz_map_py);
thm_get_neighboring_grid_points(relative_grid_points,
@ -622,7 +621,7 @@ py_thm_relative_grid_address(PyObject *self, PyObject *args)
int (*relative_grid_address)[4][3] =
(int(*)[4][3])PyArray_DATA(relative_grid_address_py);
THMCONST double (*reciprocal_lattice)[3] =
PHPYCONST double (*reciprocal_lattice)[3] =
(double(*)[3])PyArray_DATA(reciprocal_lattice_py);
thm_get_relative_grid_address(relative_grid_address, reciprocal_lattice);
@ -661,7 +660,7 @@ py_thm_integration_weight(PyObject *self, PyObject *args)
return NULL;
}
THMCONST double (*tetrahedra_omegas)[4] =
PHPYCONST double (*tetrahedra_omegas)[4] =
(double(*)[4])PyArray_DATA(tetrahedra_omegas_py);
double iw = thm_get_integration_weight(omega,
@ -689,7 +688,7 @@ py_thm_integration_weight_at_omegas(PyObject *self, PyObject *args)
const double *omegas = (double*)PyArray_DATA(omegas_py);
double *iw = (double*)PyArray_DATA(integration_weights_py);
const int num_omegas = (int)PyArray_DIMS(omegas_py)[0];
THMCONST double (*tetrahedra_omegas)[4] =
PHPYCONST double (*tetrahedra_omegas)[4] =
(double(*)[4])PyArray_DATA(tetrahedra_omegas_py);
thm_get_integration_weight_at_omegas(iw,
@ -726,9 +725,9 @@ static PyObject * py_get_tetrahedra_frequenies(PyObject *self, PyObject *args)
const int* grid_points = (int*)PyArray_DATA(grid_points_py);
const int num_gp_in = (int)PyArray_DIMS(grid_points_py)[0];
const int* mesh = (int*)PyArray_DATA(mesh_py);
THMCONST int (*grid_address)[3] = (int(*)[3])PyArray_DATA(grid_address_py);
PHPYCONST int (*grid_address)[3] = (int(*)[3])PyArray_DATA(grid_address_py);
const int* gp_ir_index = (int*)PyArray_DATA(gp_ir_index_py);
THMCONST int (*relative_grid_address)[3] =
PHPYCONST int (*relative_grid_address)[3] =
(int(*)[3])PyArray_DATA(relative_grid_address_py);
const double* frequencies = (double*)PyArray_DATA(frequencies_py);
const int num_band = (int)PyArray_DIMS(frequencies_py)[1];
@ -796,12 +795,12 @@ static PyObject * py_run_tetrahedron_method(PyObject *self, PyObject *args)
const double* frequencies = (double*)PyArray_DATA(frequencies_py);
const int num_band = (int)PyArray_DIMS(frequencies_py)[1];
const int* weights = (int*)PyArray_DATA(weights_py);
THMCONST int (*grid_address)[3] = (int(*)[3])PyArray_DATA(grid_address_py);
PHPYCONST int (*grid_address)[3] = (int(*)[3])PyArray_DATA(grid_address_py);
const int num_gp = (int)PyArray_DIMS(grid_address_py)[0];
const int* grid_mapping_table = (int*)PyArray_DATA(grid_mapping_table_py);
const int* ir_gp = (int*)PyArray_DATA(ir_grid_points_py);
const int num_ir_gp = (int)PyArray_DIMS(ir_grid_points_py)[0];
THMCONST int (*relative_grid_address)[4][3] =
PHPYCONST int (*relative_grid_address)[4][3] =
(int(*)[4][3])PyArray_DATA(relative_grid_address_py);
const int is_shift[3] = {0, 0, 0};

View File

@ -36,6 +36,7 @@
#include <stdlib.h>
#include <lapacke.h>
#include <phonoc_array.h>
#include <phonoc_const.h>
#include <phonoc_utils.h>
#include <phonon3_h/interaction.h>
#include <phonon3_h/real_to_reciprocal.h>
@ -49,6 +50,7 @@ static const int index_exchange[6][3] = {{0, 1, 2},
{1, 0, 2}};
static void get_interaction_at_triplet(Darray *fc3_normal_squared,
const int i,
const char *g_zero,
const Darray *frequencies,
const Carray *eigenvectors,
const Iarray *triplets,
@ -66,6 +68,7 @@ static void get_interaction_at_triplet(Darray *fc3_normal_squared,
const int num_triplets,
const int openmp_at_bands);
static void real_to_normal(double *fc3_normal_squared,
const char *g_zero,
const double *freqs0,
const double *freqs1,
const double *freqs2,
@ -87,8 +90,9 @@ static void real_to_normal(double *fc3_normal_squared,
const int num_triplets,
const int openmp_at_bands);
static void real_to_normal_sym_q(double *fc3_normal_squared,
double *freqs[3],
lapack_complex_double *eigvecs[3],
const char *g_zero,
PHPYCONST double *freqs[3],
PHPYCONST lapack_complex_double *eigvecs[3],
const Darray *fc3,
const double q[9], /* q0, q1, q2 */
const Darray *shortest_vectors,
@ -106,6 +110,7 @@ static void real_to_normal_sym_q(double *fc3_normal_squared,
/* fc3_normal_squared[num_triplets, num_band0, num_band, num_band] */
void get_interaction(Darray *fc3_normal_squared,
const char *g_zero,
const Darray *frequencies,
const Carray *eigenvectors,
const Iarray *triplets,
@ -130,6 +135,7 @@ void get_interaction(Darray *fc3_normal_squared,
for (i = 0; i < triplets->dims[0]; i++) {
get_interaction_at_triplet(fc3_normal_squared,
i,
g_zero,
frequencies,
eigenvectors,
triplets,
@ -151,6 +157,7 @@ void get_interaction(Darray *fc3_normal_squared,
for (i = 0; i < triplets->dims[0]; i++) {
get_interaction_at_triplet(fc3_normal_squared,
i,
g_zero,
frequencies,
eigenvectors,
triplets,
@ -173,6 +180,7 @@ void get_interaction(Darray *fc3_normal_squared,
static void get_interaction_at_triplet(Darray *fc3_normal_squared,
const int i,
const char *g_zero,
const Darray *frequencies,
const Carray *eigenvectors,
const Iarray *triplets,
@ -210,6 +218,7 @@ static void get_interaction_at_triplet(Darray *fc3_normal_squared,
if (symmetrize_fc3_q) {
real_to_normal_sym_q((fc3_normal_squared->data +
i * num_band0 * num_band * num_band),
g_zero + i * num_band0 * num_band * num_band,
freqs,
eigvecs,
fc3,
@ -229,6 +238,7 @@ static void get_interaction_at_triplet(Darray *fc3_normal_squared,
} else {
real_to_normal((fc3_normal_squared->data +
i * num_band0 * num_band * num_band),
g_zero + i * num_band0 * num_band * num_band,
freqs[0],
freqs[1],
freqs[2],
@ -253,6 +263,7 @@ static void get_interaction_at_triplet(Darray *fc3_normal_squared,
}
static void real_to_normal(double *fc3_normal_squared,
const char* g_zero,
const double *freqs0,
const double *freqs1,
const double *freqs2,
@ -297,6 +308,7 @@ static void real_to_normal(double *fc3_normal_squared,
triplet_index, num_triplets, num_band0);
#endif
reciprocal_to_normal_squared_openmp(fc3_normal_squared,
g_zero,
fc3_reciprocal,
freqs0,
freqs1,
@ -311,6 +323,7 @@ static void real_to_normal(double *fc3_normal_squared,
cutoff_frequency);
} else {
reciprocal_to_normal_squared(fc3_normal_squared,
g_zero,
fc3_reciprocal,
freqs0,
freqs1,
@ -329,6 +342,7 @@ static void real_to_normal(double *fc3_normal_squared,
}
static void real_to_normal_sym_q(double *fc3_normal_squared,
const char *g_zero,
double *freqs[3],
lapack_complex_double *eigvecs[3],
const Darray *fc3,
@ -365,6 +379,7 @@ static void real_to_normal_sym_q(double *fc3_normal_squared,
}
}
real_to_normal(fc3_normal_squared_ex,
g_zero,
freqs[index_exchange[i][0]],
freqs[index_exchange[i][1]],
freqs[index_exchange[i][2]],

View File

@ -53,9 +53,8 @@ static lapack_complex_double fc3_sum_in_reciprocal_to_normal
const double *masses,
const int num_atom);
static void set_fc3_sum
(double *fc3_normal_squared,
const int i,
static double get_fc3_sum
(const int i,
const int j,
const int k,
const int bi,
@ -73,6 +72,7 @@ static void set_fc3_sum
void reciprocal_to_normal_squared
(double *fc3_normal_squared,
const char *g_zero,
const lapack_complex_double *fc3_reciprocal,
const double *freqs0,
const double *freqs1,
@ -95,12 +95,12 @@ void reciprocal_to_normal_squared
if (freqs0[bi] > cutoff_frequency) {
for (j = 0; j < num_band; j++) {
for (k = 0; k < num_band; k++) {
if (fc3_normal_squared
[i * num_band * num_band + j * num_band + k] < 0) {
if (g_zero[i * num_band * num_band + j * num_band + k]) {
fc3_normal_squared[i * num_band * num_band + j * num_band + k] = 0;
} else {
set_fc3_sum(fc3_normal_squared,
i, j, k, bi,
continue;
}
fc3_normal_squared[i * num_band * num_band + j * num_band + k] =
get_fc3_sum(i, j, k, bi,
freqs0, freqs1, freqs2,
eigvecs0, eigvecs1, eigvecs2,
fc3_reciprocal,
@ -110,7 +110,6 @@ void reciprocal_to_normal_squared
cutoff_frequency);
}
}
}
} else {
for (j = 0; j < num_band * num_band; j++) {
fc3_normal_squared[i * num_band * num_band + j] = 0;
@ -121,6 +120,7 @@ void reciprocal_to_normal_squared
void reciprocal_to_normal_squared_openmp
(double *fc3_normal_squared,
const char *g_zero,
const lapack_complex_double *fc3_reciprocal,
const double *freqs0,
const double *freqs1,
@ -155,13 +155,14 @@ void reciprocal_to_normal_squared_openmp
#pragma omp parallel for private(j, k)
for (jk = 0; jk < num_band * num_band; jk++) {
if (fc3_normal_squared[i * num_band * num_band + jk] < 0) {
if (g_zero[i * num_band * num_band + jk]) {
fc3_normal_squared[i * num_band * num_band + jk] = 0;
} else {
continue;
}
j = jk / num_band;
k = jk % num_band;
set_fc3_sum(fc3_normal_squared,
i, j, k, bi,
fc3_normal_squared[i * num_band * num_band + jk] =
get_fc3_sum(i, j, k, bi,
freqs0, freqs1, freqs2,
eigvecs0, eigvecs1, eigvecs2,
fc3_reciprocal,
@ -170,7 +171,6 @@ void reciprocal_to_normal_squared_openmp
num_band,
cutoff_frequency);
}
}
#ifdef MEASURE_R2N
loopTotalCPUTime = (double)(clock() - loopStartCPUTime) / CLOCKS_PER_SEC;
@ -201,9 +201,8 @@ void reciprocal_to_normal_squared_openmp
}
}
static void set_fc3_sum
(double *fc3_normal_squared,
const int i,
static double get_fc3_sum
(const int i,
const int j,
const int k,
const int bi,
@ -232,10 +231,9 @@ static void set_fc3_sum
num_atom);
sum_real = lapack_complex_double_real(fc3_sum);
sum_imag = lapack_complex_double_imag(fc3_sum);
fc3_normal_squared[i * num_band * num_band + j * num_band + k] =
(sum_real * sum_real + sum_imag * sum_imag) / fff;
return (sum_real * sum_real + sum_imag * sum_imag) / fff;
} else {
fc3_normal_squared[i * num_band * num_band + j * num_band + k] = 0;
return 0;
}
}

View File

@ -87,6 +87,7 @@ int tpl_get_triplets_reciprocal_mesh_at_q(int map_triplets[],
}
int tpl_get_integration_weight(double *iw,
char *iw_zero,
const double frequency_points[],
const int num_band0,
PHPYCONST int relative_grid_address[24][4][3],
@ -100,6 +101,7 @@ int tpl_get_integration_weight(double *iw,
const int num_iw)
{
return tpi_get_integration_weight(iw,
iw_zero,
frequency_points,
num_band0,
relative_grid_address,

View File

@ -36,6 +36,15 @@
#include <phonoc_const.h>
#include <tetrahedron_method.h>
static void set_freq_vertices(double freq_vertices[3][24][4],
const double frequencies[],
PHPYCONST int vertices[2][24][4],
const int num_band,
const int b1,
const int b2);
static int set_g(double g[3],
const double f0,
PHPYCONST double freq_vertices[3][24][4]);
static int in_tetrahedra(const double f0, PHPYCONST double freq_vertices[24][4]);
static void get_triplet_tetrahedra_vertices
(int vertices[2][24][4],
@ -46,6 +55,7 @@ static void get_triplet_tetrahedra_vertices
const int bz_map[]);
int tpi_get_integration_weight(double *iw,
char *iw_zero,
const double frequency_points[],
const int num_band0,
PHPYCONST int relative_grid_address[24][4][3],
@ -62,7 +72,7 @@ int tpi_get_integration_weight(double *iw,
int tp_relative_grid_address[2][24][4][3];
int vertices[2][24][4];
int adrs_shift;
double f0, f1, f2, g0, g1, g2;
double g[3];
double freq_vertices[3][24][4];
for (i = 0; i < 2; i++) {
@ -77,7 +87,7 @@ int tpi_get_integration_weight(double *iw,
}
}
#pragma omp parallel for private(j, k, b1, b2, vertices, adrs_shift, f0, f1, f2, g0, g1, g2, freq_vertices)
#pragma omp parallel for private(j, b1, b2, vertices, adrs_shift, g, freq_vertices)
for (i = 0; i < num_triplets; i++) {
get_triplet_tetrahedra_vertices(vertices,
tp_relative_grid_address,
@ -87,69 +97,18 @@ int tpi_get_integration_weight(double *iw,
bz_map);
for (b1 = 0; b1 < num_band; b1++) {
for (b2 = 0; b2 < num_band; b2++) {
for (j = 0; j < 24; j++) {
for (k = 0; k < 4; k++) {
f1 = frequencies[vertices[0][j][k] * num_band + b1];
f2 = frequencies[vertices[1][j][k] * num_band + b2];
freq_vertices[0][j][k] = f1 + f2;
freq_vertices[1][j][k] = -f1 + f2;
freq_vertices[2][j][k] = f1 - f2;
}
}
set_freq_vertices
(freq_vertices, frequencies, vertices, num_band, b1, b2);
for (j = 0; j < num_band0; j++) {
f0 = frequency_points[j];
if (in_tetrahedra(f0, freq_vertices[0])) {
g0 = thm_get_integration_weight(f0, freq_vertices[0], 'I');
} else {
g0 = -1;
}
if (in_tetrahedra(f0, freq_vertices[1])) {
g1 = thm_get_integration_weight(f0, freq_vertices[1], 'I');
} else {
g1 = -1;
}
if (in_tetrahedra(f0, freq_vertices[2])) {
g2 = thm_get_integration_weight(f0, freq_vertices[2], 'I');
} else {
g2 = -1;
}
adrs_shift = i * num_band0 * num_band * num_band +
j * num_band * num_band + b1 * num_band + b2;
if (g0 < 0) {
iw[adrs_shift] = 0;
} else {
iw[adrs_shift] = g0;
}
iw_zero[adrs_shift] = set_g(g, frequency_points[j], freq_vertices);
iw[adrs_shift] = g[0];
adrs_shift += num_triplets * num_band0 * num_band * num_band;
if (g1 < 0 && g2 < 0) {
iw[adrs_shift] = 0;
} else {
if (g1 < 0) {
iw[adrs_shift] = - g2;
} else {
if (g2 < 0) {
iw[adrs_shift] = g1;
} else {
iw[adrs_shift] = g1 - g2;
}
}
}
iw[adrs_shift] = g[1] - g[2];
if (num_iw == 3) {
adrs_shift += num_triplets * num_band0 * num_band * num_band;
if (g0 < 0 && g1 < 0 && g2 < 0) {
iw[adrs_shift] = 0;
} else {
if (g0 < 0) {
g0 = 0;
}
if (g1 < 0) {
g1 = 0;
}
if (g2 < 0) {
g2 = 0;
}
iw[adrs_shift] = g0 + g1 + g2;
}
iw[adrs_shift] = g[0] + g[1] + g[2];
}
}
}
@ -159,6 +118,57 @@ int tpi_get_integration_weight(double *iw,
return 0;
}
static void set_freq_vertices(double freq_vertices[3][24][4],
const double frequencies[],
PHPYCONST int vertices[2][24][4],
const int num_band,
const int b1,
const int b2)
{
int i, j;
double f1, f2;
for (i = 0; i < 24; i++) {
for (j = 0; j < 4; j++) {
f1 = frequencies[vertices[0][i][j] * num_band + b1];
f2 = frequencies[vertices[1][i][j] * num_band + b2];
freq_vertices[0][i][j] = f1 + f2;
freq_vertices[1][i][j] = -f1 + f2;
freq_vertices[2][i][j] = f1 - f2;
}
}
}
static int set_g(double g[3],
const double f0,
PHPYCONST double freq_vertices[3][24][4])
{
int iw_zero;
iw_zero = 1;
if (in_tetrahedra(f0, freq_vertices[0])) {
g[0] = thm_get_integration_weight(f0, freq_vertices[0], 'I');
iw_zero = 0;
} else {
g[0] = 0;
}
if (in_tetrahedra(f0, freq_vertices[1])) {
g[1] = thm_get_integration_weight(f0, freq_vertices[1], 'I');
iw_zero = 0;
} else {
g[1] = 0;
}
if (in_tetrahedra(f0, freq_vertices[2])) {
g[2] = thm_get_integration_weight(f0, freq_vertices[2], 'I');
iw_zero = 0;
} else {
g[2] = 0;
}
return iw_zero;
}
static int in_tetrahedra(const double f0, PHPYCONST double freq_vertices[24][4])
{
int i, j;

View File

@ -38,6 +38,7 @@
#include <phonoc_array.h>
void get_interaction(Darray *fc3_normal_squared,
const char *g_zero,
const Darray *frequencies,
const Carray *eigenvectors,
const Iarray *triplets,

View File

@ -40,6 +40,7 @@
void reciprocal_to_normal_squared
(double *fc3_normal_squared,
const char *g_zero,
const lapack_complex_double *fc3_reciprocal,
const double *freqs0,
const double *freqs1,
@ -55,6 +56,7 @@ void reciprocal_to_normal_squared
void reciprocal_to_normal_squared_openmp
(double *fc3_normal_squared,
const char *g_zero,
const lapack_complex_double *fc3_reciprocal,
const double *freqs0,
const double *freqs1,

View File

@ -70,6 +70,7 @@ int tpl_get_BZ_triplets_at_q(int triplets[][3],
int tpl_get_integration_weight(double *iw,
char *iw_zero,
const double frequency_points[],
const int num_band0,
PHPYCONST int relative_grid_address[24][4][3],

View File

@ -38,6 +38,7 @@
#include <phonoc_const.h>
int tpi_get_integration_weight(double *iw,
char *iw_zero,
const double frequency_points[],
const int num_band0,
PHPYCONST int relative_grid_address[24][4][3],

View File

@ -36,6 +36,7 @@
#define __phonoc_const_H__
#define M_2PI 6.283185307179586
#define KB 8.6173382568083159E-05
/* PHPYCONST is used instead of 'const' so to avoid gcc warning. */
/* However there should be better way than this way.... */