diff --git a/anharmonic/phonon3/imag_self_energy.py b/anharmonic/phonon3/imag_self_energy.py index 6cc5c2f8..4977923b 100644 --- a/anharmonic/phonon3/imag_self_energy.py +++ b/anharmonic/phonon3/imag_self_energy.py @@ -2,9 +2,11 @@ import numpy as np from phonopy.units import THzToEv, Kb, VaspToTHz, Hbar, EV, Angstrom, THz, AMU from phonopy.phonon.group_velocity import degenerate_sets from phonopy.structure.tetrahedron_method import TetrahedronMethod +from phonopy.structure.spglib import get_neighboring_grid_points from anharmonic.phonon3.triplets import get_tetrahedra_vertices from anharmonic.phonon3.interaction import set_phonon_c + def gaussian(x, sigma): return 1.0 / np.sqrt(2 * np.pi) / sigma * np.exp(-x**2 / 2 / sigma**2) @@ -237,11 +239,10 @@ class ImagSelfEnergy: n2 = n[i, 0, j] n3 = n[i, 1, k] g1 = self._g1[i, j, k] - g2 = self._g2[i, j, k] - g3 = self._g3[i, j, k] + g2_g3 = self._g2_g3[i, j, k] # g2 - g3 self._imag_self_energy[:] += ( (n2 + n3 + 1) * g1 + - (n2 - n3) * (g2 - g3)) * interaction[:, j, k] * w + (n2 - n3) * (g2_g3)) * interaction[:, j, k] * w self._imag_self_energy *= self._unit_conversion @@ -314,12 +315,11 @@ class ImagSelfEnergy: n2 = occupation(f1, self._temperature) n3 = occupation(f2, self._temperature) g1 = self._g1[i, j, k] - g2 = self._g2[i, j, k] - g3 = self._g3[i, j, k] + g2_g3 = self._g2_g3[i, j, k] # g2 - g3 for l in range(len(interaction)): self._imag_self_energy[:, l] += ( (n2 + n3 + 1) * g1 + - (n2 - n3) * (g2 - g3)) * interaction[l, j, k] * w + (n2 - n3) * (g2_g3)) * interaction[l, j, k] * w self._imag_self_energy *= self._unit_conversion @@ -337,16 +337,27 @@ class ImagSelfEnergy: reciprocal_lattice = np.linalg.inv( self._interaction.get_primitive().get_cell()) thm = TetrahedronMethod(reciprocal_lattice) - grid_address = self._interaction.get_grid_address() bz_map = self._interaction.get_bz_map() + + unique_vertices = thm.get_unique_tetrahedra_vertices() + grid_points = np.array([self._triplets_at_q[0, 0]], dtype='intc') + for i, j in zip((1, 2), (1, -1)): + for gp in np.unique(self._triplets_at_q[:, i]): + grid_points = np.hstack( + (grid_points, + get_neighboring_grid_points(gp, + j * unique_vertices, + self._mesh, + grid_address, + bz_map))) + self._interaction.set_phonon(np.unique(grid_points)) tetrahedra_vertices = get_tetrahedra_vertices( thm.get_tetrahedra(), self._mesh, self._triplets_at_q, grid_address, bz_map) - self._interaction.set_phonon(tetrahedra_vertices.ravel()) if self._frequency_points is None: gp = self._grid_point frequency_points = self._frequencies[gp][self._band_indices] @@ -356,8 +367,7 @@ class ImagSelfEnergy: shape = self._fc3_normal_squared.shape self._g1 = np.zeros((shape[0], shape[2], shape[3], len(frequency_points)), dtype='double') - self._g2 = np.zeros_like(self._g1) - self._g3 = np.zeros_like(self._g1) + self._g2_g3 = np.zeros_like(self._g1) self._set_triplets_integration_weights_py(thm, tetrahedra_vertices, frequency_points) @@ -382,7 +392,7 @@ class ImagSelfEnergy: self._g1[i, j, k] = thm.get_integration_weight() thm.set_tetrahedra_omegas(-f1_v + f2_v) thm.run(frequency_points) - self._g2[i, j, k] = thm.get_integration_weight() + self._g2_g3[i, j, k] = thm.get_integration_weight() thm.set_tetrahedra_omegas(f1_v - f2_v) thm.run(frequency_points) - self._g3[i, j, k] = thm.get_integration_weight() + self._g2_g3[i, j, k] -= thm.get_integration_weight() diff --git a/c/_spglib.c b/c/_spglib.c index af787d91..a610195d 100644 --- a/c/_spglib.c +++ b/c/_spglib.c @@ -17,6 +17,7 @@ static PyObject * relocate_BZ_grid_address(PyObject *self, PyObject *args); static PyObject * get_triplets_reciprocal_mesh_at_q(PyObject *self, PyObject *args); static PyObject * get_BZ_triplets_at_q(PyObject *self, PyObject *args); +static PyObject *get_neighboring_grid_points(PyObject *self, PyObject *args); static PyObject * get_triplet_tetrahedra_vertices(PyObject *self, PyObject *args); static PyObject * @@ -47,6 +48,8 @@ static PyMethodDef functions[] = { METH_VARARGS, "Triplets on reciprocal mesh points at a specific q-point"}, {"BZ_triplets_at_q", get_BZ_triplets_at_q, METH_VARARGS, "Triplets in reciprocal primitive lattice are transformed to those in BZ."}, + {"neighboring_grid_points", get_neighboring_grid_points, + METH_VARARGS, "Neighboring grid points by relative grid addresses"}, {"triplet_tetrahedra_vertices", get_triplet_tetrahedra_vertices, METH_VARARGS, "Tetrahedra vertices of tetrahedron method for triplets"}, {"tetrahedra_relative_grid_address", get_tetrahedra_relative_grid_address, @@ -578,6 +581,42 @@ static PyObject * get_BZ_triplets_at_q(PyObject *self, PyObject *args) return PyInt_FromLong((long) num_ir); } +static PyObject *get_neighboring_grid_points(PyObject *self, PyObject *args) +{ + PyArrayObject* relative_grid_points_py; + PyArrayObject* relative_grid_address_py; + PyArrayObject* mesh_py; + PyArrayObject* bz_grid_address_py; + PyArrayObject* bz_map_py; + int grid_point; + if (!PyArg_ParseTuple(args, "OiOOOO", + &relative_grid_points_py, + &grid_point, + &relative_grid_address_py, + &mesh_py, + &bz_grid_address_py, + &bz_map_py)) { + return NULL; + } + + int* relative_grid_points = (int*)relative_grid_points_py->data; + SPGCONST 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; + SPGCONST int (*bz_grid_address)[3] = (int(*)[3])bz_grid_address_py->data; + const int *bz_map = (int*)bz_map_py->data; + + spg_get_neighboring_grid_points(relative_grid_points, + grid_point, + relative_grid_address, + num_relative_grid_address, + mesh, + bz_grid_address, + bz_map); + Py_RETURN_NONE; +} + static PyObject * get_triplet_tetrahedra_vertices(PyObject *self, PyObject *args) { diff --git a/c/spglib/kpoint.c b/c/spglib/kpoint.c index 1063e339..96636949 100644 --- a/c/spglib/kpoint.c +++ b/c/spglib/kpoint.c @@ -228,6 +228,40 @@ int kpt_get_BZ_triplets_at_q(int triplets[][3], mesh); } +void kpt_get_neighboring_grid_points(int relative_grid_points[], + const int grid_point, + SPGCONST int relative_grid_address[][3], + const int num_relative_grid_address, + const int mesh[3], + SPGCONST int bz_grid_address[][3], + const int bz_map[]) +{ + int mesh_double[3], bzmesh[3], bzmesh_double[3], + address_double[3], bz_address_double[3]; + int i, j, bz_gp; + + for (i = 0; i < 3; i++) { + mesh_double[i] = mesh[i] * 2; + bzmesh[i] = mesh[i] * 2; + bzmesh_double[i] = bzmesh[i] * 2; + } + for (i = 0; i < num_relative_grid_address; i++) { + for (j = 0; j < 3; j++) { + address_double[j] = (bz_grid_address[grid_point][j] + + relative_grid_address[i][j]) * 2; + bz_address_double[j] = address_double[j]; + } + get_vector_modulo(bz_address_double, bzmesh_double); + bz_gp = bz_map[get_grid_point(bz_address_double, bzmesh)]; + if (bz_gp == -1) { + get_vector_modulo(address_double, mesh_double); + relative_grid_points[i] = get_grid_point(address_double, mesh); + } else { + relative_grid_points[i] = bz_gp; + } + } +} + void kpt_get_triplet_tetrahedra_vertices (int vertices[2][24][4], SPGCONST int relative_grid_address[24][4][3], diff --git a/c/spglib/spglib.c b/c/spglib/spglib.c index 4708aeb1..b599fa17 100644 --- a/c/spglib/spglib.c +++ b/c/spglib/spglib.c @@ -550,6 +550,24 @@ int spg_get_BZ_triplets_at_q(int triplets[][3], mesh); } +void spg_get_neighboring_grid_points(int relative_grid_points[], + const int grid_point, + SPGCONST int relative_grid_address[][3], + const int num_relative_grid_address, + const int mesh[3], + SPGCONST int bz_grid_address[][3], + const int bz_map[]) +{ + kpt_get_neighboring_grid_points(relative_grid_points, + grid_point, + relative_grid_address, + num_relative_grid_address, + mesh, + bz_grid_address, + bz_map); +} + + void spg_get_triplet_tetrahedra_vertices (int vertices[2][24][4], SPGCONST int relative_grid_address[24][4][3], diff --git a/c/spglib_h/kpoint.h b/c/spglib_h/kpoint.h index 0a1e8f33..7f2ab6bc 100644 --- a/c/spglib_h/kpoint.h +++ b/c/spglib_h/kpoint.h @@ -39,6 +39,13 @@ int kpt_get_BZ_triplets_at_q(int triplets[][3], const int bz_map[], const int weights[], const int mesh[3]); +void kpt_get_neighboring_grid_points(int relative_grid_points[], + const int grid_point, + SPGCONST int relative_grid_address[][3], + const int num_relative_grid_address, + const int mesh[3], + SPGCONST int bz_grid_address[][3], + const int bz_map[]); void kpt_get_triplet_tetrahedra_vertices (int vertices[2][24][4], SPGCONST int relative_grid_address[24][4][3], diff --git a/c/spglib_h/spglib.h b/c/spglib_h/spglib.h index 5ba8a7ef..b416a3ee 100644 --- a/c/spglib_h/spglib.h +++ b/c/spglib_h/spglib.h @@ -339,6 +339,14 @@ int spg_get_BZ_triplets_at_q(int triplets[][3], const int triplet_weights[], const int mesh[3]); +void spg_get_neighboring_grid_points(int relative_grid_points[], + const int grid_point, + SPGCONST int relative_grid_address[][3], + const int num_relative_grid_address, + const int mesh[3], + SPGCONST int bz_grid_address[][3], + const int bz_map[]); + /* Tetrahedra vertices of tetrahedron method for a triplet */ void spg_get_triplet_tetrahedra_vertices (int vertices[2][24][4], diff --git a/phonopy/structure/spglib.py b/phonopy/structure/spglib.py index df0aec92..2fa6bb31 100644 --- a/phonopy/structure/spglib.py +++ b/phonopy/structure/spglib.py @@ -337,6 +337,21 @@ def get_BZ_triplets_at_q(grid_point, np.array(mesh, dtype='intc').copy()) return triplets +def get_neighboring_grid_points(grid_point, + relative_grid_address, + mesh, + bz_grid_address, + bz_map): + relative_grid_points = np.zeros(len(relative_grid_address), dtype='intc') + spg.neighboring_grid_points(relative_grid_points, + grid_point, + relative_grid_address, + mesh, + bz_grid_address, + bz_map) + return relative_grid_points + + def get_triplets_tetrahedra_vertices(relative_grid_address, mesh, triplets, diff --git a/phonopy/structure/tetrahedron_method.py b/phonopy/structure/tetrahedron_method.py index 87fa6d33..0286bf64 100644 --- a/phonopy/structure/tetrahedron_method.py +++ b/phonopy/structure/tetrahedron_method.py @@ -71,6 +71,18 @@ class TetrahedronMethod: """ return self._relative_grid_addresses + def get_unique_tetrahedra_vertices(self): + unique_vertices = [] + for adrs in self._relative_grid_addresses.reshape(-1, 3): + found = False + for uadrs in unique_vertices: + if (uadrs == adrs).all(): + found = True + break + if not found: + unique_vertices.append(adrs) + return np.array(unique_vertices, dtype='intc') + def set_tetrahedra_omegas(self, tetrahedra_omegas): """ tetrahedra_omegas: (24, 4) omegas at self._relative_grid_addresses