scnc
/
phonopy
mirror of https://github.com/phonopy/phonopy.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Update tetrahedron method implementing in C

This commit is contained in:
Atsushi Togo 2014-01-21 00:54:07 +09:00
parent 06e45f82fb
commit b3be2f54d9
8 changed files with 155 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

34
anharmonic/phonon3/imag_self_energy.py
View File

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

39
c/_spglib.c
View File

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

34
c/spglib/kpoint.c
View File

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

18
c/spglib/spglib.c
View File

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

7
c/spglib_h/kpoint.h
View File

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

8
c/spglib_h/spglib.h
View File

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

15
phonopy/structure/spglib.py
View File

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

12
phonopy/structure/tetrahedron_method.py
View File

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