diff --git a/anharmonic/phonon3/conductivity_RTA.py b/anharmonic/phonon3/conductivity_RTA.py index 6bcfafdd..d330fb92 100644 --- a/anharmonic/phonon3/conductivity_RTA.py +++ b/anharmonic/phonon3/conductivity_RTA.py @@ -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 @@ -450,10 +452,16 @@ class Conductivity_RTA(Conductivity): print("Number of triplets: %d" % len(self._pp.get_triplets_at_q()[0])) print("Calculating interaction...") - - self._collision.run_interaction() - self._averaged_pp_interaction[i] = ( - self._pp.get_averaged_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) if self._isotope is not None and not self._read_gamma_iso: @@ -506,9 +514,16 @@ class Conductivity_RTA(Conductivity): else: text += "sigma=%s" % sigma print(text) - self._collision.set_sigma(sigma) - if not sigma or self._run_with_g: - self._collision.set_integration_weights() + + 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 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) self._collision.run() @@ -554,8 +569,17 @@ class Conductivity_RTA(Conductivity): frequencies = self._frequencies[gp][self._pp.get_band_indices()] gv = self._gv[i] ave_pp = self._averaged_pp_interaction[i] - - text = "Frequency group velocity (x, y, z) |gv| Pqj" + + 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): - 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)) + 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): - 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)) + 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))) diff --git a/anharmonic/phonon3/imag_self_energy.py b/anharmonic/phonon3/imag_self_energy.py index 6548df85..b48033a5 100644 --- a/anharmonic/phonon3/imag_self_energy.py +++ b/anharmonic/phonon3/imag_self_energy.py @@ -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): - self._pp.run(lang=self._lang) + 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: diff --git a/anharmonic/phonon3/interaction.py b/anharmonic/phonon3/interaction.py index ddafc0c7..01516715 100644 --- a/anharmonic/phonon3/interaction.py +++ b/anharmonic/phonon3/interaction.py @@ -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, diff --git a/anharmonic/phonon3/joint_dos.py b/anharmonic/phonon3/joint_dos.py index 6d0cdf9b..e5decb3a 100644 --- a/anharmonic/phonon3/joint_dos.py +++ b/anharmonic/phonon3/joint_dos.py @@ -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, diff --git a/anharmonic/phonon3/triplets.py b/anharmonic/phonon3/triplets.py index c9e3918c..ed26df96 100644 --- a/anharmonic/phonon3/triplets.py +++ b/anharmonic/phonon3/triplets.py @@ -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() diff --git a/c/_phono3py.c b/c/_phono3py.c index ba5034e5..a98e14d3 100644 --- a/c/_phono3py.c +++ b/c/_phono3py.c @@ -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; diff --git a/c/_phonopy.c b/c/_phonopy.c index 61a2d147..4c57f37c 100644 --- a/c/_phonopy.c +++ b/c/_phonopy.c @@ -39,8 +39,7 @@ #include #include #include - -#define KB 8.6173382568083159E-05 +#include /* 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}; diff --git a/c/anharmonic/phonon3/interaction.c b/c/anharmonic/phonon3/interaction.c index 70e43158..7ef6e19d 100644 --- a/c/anharmonic/phonon3/interaction.c +++ b/c/anharmonic/phonon3/interaction.c @@ -36,6 +36,7 @@ #include #include #include +#include #include #include #include @@ -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]], diff --git a/c/anharmonic/phonon3/reciprocal_to_normal.c b/c/anharmonic/phonon3/reciprocal_to_normal.c index e4fb1334..460cc7c0 100644 --- a/c/anharmonic/phonon3/reciprocal_to_normal.c +++ b/c/anharmonic/phonon3/reciprocal_to_normal.c @@ -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, @@ -89,18 +89,18 @@ void reciprocal_to_normal_squared int i, j, k, bi, num_atom; num_atom = num_band / 3; - + for (i = 0; i < num_band0; i++) { bi = band_indices[i]; 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, @@ -108,7 +108,6 @@ void reciprocal_to_normal_squared num_atom, num_band, cutoff_frequency); - } } } } else { @@ -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 { - j = jk / num_band; - k = jk % num_band; - set_fc3_sum(fc3_normal_squared, - i, j, k, bi, + continue; + } + j = jk / num_band; + k = jk % num_band; + fc3_normal_squared[i * num_band * num_band + jk] = + get_fc3_sum(i, j, k, bi, freqs0, freqs1, freqs2, eigvecs0, eigvecs1, eigvecs2, fc3_reciprocal, @@ -169,7 +170,6 @@ void reciprocal_to_normal_squared_openmp num_atom, num_band, cutoff_frequency); - } } #ifdef MEASURE_R2N @@ -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; } } diff --git a/c/anharmonic/triplet/triplet.c b/c/anharmonic/triplet/triplet.c index f90a8716..9febb919 100644 --- a/c/anharmonic/triplet/triplet.c +++ b/c/anharmonic/triplet/triplet.c @@ -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, diff --git a/c/anharmonic/triplet/triplet_iw.c b/c/anharmonic/triplet/triplet_iw.c index 8808cb6c..bef16a6a 100644 --- a/c/anharmonic/triplet/triplet_iw.c +++ b/c/anharmonic/triplet/triplet_iw.c @@ -36,6 +36,15 @@ #include #include +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; diff --git a/c/anharmonic_h/phonon3_h/interaction.h b/c/anharmonic_h/phonon3_h/interaction.h index 71940806..0ad3eace 100644 --- a/c/anharmonic_h/phonon3_h/interaction.h +++ b/c/anharmonic_h/phonon3_h/interaction.h @@ -38,6 +38,7 @@ #include void get_interaction(Darray *fc3_normal_squared, + const char *g_zero, const Darray *frequencies, const Carray *eigenvectors, const Iarray *triplets, diff --git a/c/anharmonic_h/phonon3_h/reciprocal_to_normal.h b/c/anharmonic_h/phonon3_h/reciprocal_to_normal.h index 536077be..6ce3db7f 100644 --- a/c/anharmonic_h/phonon3_h/reciprocal_to_normal.h +++ b/c/anharmonic_h/phonon3_h/reciprocal_to_normal.h @@ -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, diff --git a/c/anharmonic_h/triplet_h/triplet.h b/c/anharmonic_h/triplet_h/triplet.h index 73fe40eb..fbdbbf37 100644 --- a/c/anharmonic_h/triplet_h/triplet.h +++ b/c/anharmonic_h/triplet_h/triplet.h @@ -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], diff --git a/c/anharmonic_h/triplet_h/triplet_iw.h b/c/anharmonic_h/triplet_h/triplet_iw.h index fca7ef36..f59fe74f 100644 --- a/c/anharmonic_h/triplet_h/triplet_iw.h +++ b/c/anharmonic_h/triplet_h/triplet_iw.h @@ -38,6 +38,7 @@ #include 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], diff --git a/c/harmonic_h/phonoc_const.h b/c/harmonic_h/phonoc_const.h index 7f543c2f..e6a1fffb 100644 --- a/c/harmonic_h/phonoc_const.h +++ b/c/harmonic_h/phonoc_const.h @@ -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.... */