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

Considering BZ boundary for triplets

This commit is contained in:
Atsushi Togo 2014-05-10 14:04:24 +09:00
parent a833b510ac
commit a6045d5f14
13 changed files with 74 additions and 82 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
anharmonic/phonon3/collision_matrix.py
View File

@ -124,8 +124,11 @@ class CollisionMatrix(ImagSelfEnergy):
def _run_py_collision_matrix(self):
gp2tp_map = {}
for i, j in enumerate(self._triplets_at_q[:, 1]):
gp2tp_map[j] = i
count = 0
for i, j in enumerate(self._triplets_map_at_q):
if i == j:
gp2tp_map[i] = count
count += 1
num_band = self._fc3_normal_squared.shape[1]
for i, ir_gp in enumerate(self._ir_grid_points):

15
anharmonic/phonon3/conductivity_LBTE.py
View File

@ -346,8 +346,7 @@ class Conductivity_LBTE(Conductivity):
self._rot_grid_points[i] = get_grid_points_by_rotations(
self._grid_address[ir_gp],
self._point_operations,
self._mesh,
self._pp.get_bz_map())
self._mesh)
self._collision = CollisionMatrix(self._pp,
self._point_operations,
self._ir_grid_points,
@ -419,11 +418,9 @@ class Conductivity_LBTE(Conductivity):
for j, k in list(np.ndindex((len(self._sigmas),
len(self._temperatures)))):
for i, ir_gp in enumerate(self._ir_grid_points):
r_gps = self._rot_grid_points[i]
diffs = (self._grid_address[r_gps] -
self._grid_address[ir_gp]) % self._mesh
for r, diff in zip(self._rotations_cartesian, diffs):
if (diff != [0, 0, 0]).any():
for r, r_gp in zip(self._rotations_cartesian,
self._rot_grid_points[i]):
if ir_gp != r_gp:
continue
main_diagonal = self._gamma[j, k, i].copy()
@ -437,9 +434,7 @@ class Conductivity_LBTE(Conductivity):
def _get_weights(self):
weights = []
for r_gps in self._rot_grid_points:
multi = ((self._grid_address[r_gps] - self._grid_address[r_gps[0]])
% self._mesh == [0, 0, 0]).all(axis=1).sum()
weights.append(np.sqrt(multi) / np.sqrt(len(r_gps)))
weights.append(np.sqrt(len(np.unique(r_gps))) / np.sqrt(len(r_gps)))
return weights
def _symmetrize_collision_matrix(self, write_to_hdf5=False):

6
anharmonic/phonon3/conductivity_RTA.py
View File

@ -380,8 +380,7 @@ class Conductivity_RTA(Conductivity):
rotation_map = get_grid_points_by_rotations(
self._grid_address[self._grid_points[i]],
self._point_operations,
self._mesh,
self._pp.get_bz_map())
self._mesh)
gv_by_gv = np.zeros((len(self._gv[i]), 3, 3), dtype='double')
if self._no_kappa_stars:
@ -441,8 +440,7 @@ class Conductivity_RTA(Conductivity):
rotation_map = get_grid_points_by_rotations(
self._grid_address[gp],
self._point_operations,
self._mesh,
self._pp.get_bz_map())
self._mesh)
for i, j in enumerate(np.unique(rotation_map)):
for k, (rot, rot_c) in enumerate(zip(self._point_operations,
self._rotations_cartesian)):

2
anharmonic/phonon3/triplets.py
View File

@ -138,13 +138,11 @@ def get_ir_grid_points(mesh, rotations, mesh_shifts=[False, False, False]):
def get_grid_points_by_rotations(grid_point,
reciprocal_rotations,
mesh,
bz_map,
mesh_shifts=[False, False, False]):
return spg.get_grid_points_by_rotations(
grid_point,
reciprocal_rotations,
mesh,
bz_map,
is_shift=np.where(mesh_shifts, 1, 0))
def reduce_grid_points(mesh_divisors,

3
c/_phono3py.c
View File

@ -541,7 +541,7 @@ static PyObject * py_get_collision_matrix(PyObject *self, PyObject *args)
const double* g = (double*)g_py->data;
const double* frequencies = (double*)frequencies_py->data;
const int* triplets = (int*)triplets_py->data;
const int* triplets_map = (int*)triplets_map_py->data;
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;
Iarray* rotated_grid_points = convert_to_iarray(rotated_grid_points_py);
@ -562,6 +562,7 @@ static PyObject * py_get_collision_matrix(PyObject *self, PyObject *args)
cutoff_frequency);
free(fc3_normal_squared);
free(triplets_map);
free(rotated_grid_points);
Py_RETURN_NONE;

10
c/_spglib.c
View File

@ -480,14 +480,12 @@ static PyObject * get_grid_points_by_rotations(PyObject *self, PyObject *args)
PyArrayObject* rot_reciprocal_py;
PyArrayObject* mesh_py;
PyArrayObject* is_shift_py;
PyArrayObject* bz_map_py;
if (!PyArg_ParseTuple(args, "OOOOOO",
if (!PyArg_ParseTuple(args, "OOOOO",
&rot_grid_points_py,
&address_orig_py,
&rot_reciprocal_py,
&mesh_py,
&is_shift_py,
&bz_map_py)) {
&is_shift_py)) {
return NULL;
}
@ -497,15 +495,13 @@ static PyObject * get_grid_points_by_rotations(PyObject *self, PyObject *args)
const int num_rot = rot_reciprocal_py->dimensions[0];
const int* mesh = (int*)mesh_py->data;
const int* is_shift = (int*)is_shift_py->data;
const int* bz_map = (int*)bz_map_py->data;
spg_get_grid_points_by_rotations(rot_grid_points,
address_orig,
num_rot,
rot_reciprocal,
mesh,
is_shift,
bz_map);
is_shift);
Py_RETURN_NONE;
}

33
c/anharmonic/phonon3/collision_matrix.c
View File

@ -5,14 +5,13 @@
#include "phonoc_utils.h"
#include "phonon3_h/collision_matrix.h"
static int *create_gp2tp_map(const int num_triplets,
const int *triplets);
static int *create_gp2tp_map(const Iarray *triplets);
void get_collision_matrix(double *collision_matrix,
const Darray *fc3_normal_squared,
const double *frequencies,
const int *triplets,
const int *triplets_map,
const Iarray *triplets_map,
const int *stabilized_gp_map,
const int *ir_grid_points,
const Iarray *rotated_grid_points,
@ -32,15 +31,15 @@ void get_collision_matrix(double *collision_matrix,
num_ir_gp = rotated_grid_points->dims[0];
num_rot = rotated_grid_points->dims[1];
gp2tp_map = create_gp2tp_map(num_triplets, triplets);
gp2tp_map = create_gp2tp_map(triplets_map);
#pragma omp parallel for private(j, k, l, m, n, ti, gp2, r_gp, f, collision, inv_sinh)
for (i = 0; i < num_ir_gp; i++) {
inv_sinh = (double*)malloc(sizeof(double) * num_band);
for (j = 0; j < num_rot; j++) {
r_gp = rotated_grid_points->data[i * num_rot + j];
ti = gp2tp_map[triplets_map[r_gp]];
if (triplets_map[r_gp] == stabilized_gp_map[r_gp]) {
ti = gp2tp_map[triplets_map->data[r_gp]];
if (triplets_map->data[r_gp] == stabilized_gp_map[r_gp]) {
gp2 = triplets[ti * 3 + 2];
} else {
gp2 = triplets[ti * 3 + 1];
@ -87,23 +86,29 @@ void get_collision_matrix(double *collision_matrix,
gp2tp_map = NULL;
}
static int *create_gp2tp_map(const int num_triplets,
const int *triplets)
static int *create_gp2tp_map(const Iarray *triplets_map)
{
int i, max_i;
int i, max_i, count;
int *gp2tp_map;
max_i = 0;
for (i = 0; i < num_triplets; i++) {
if (max_i < triplets[3 * i + 1]) {
max_i = triplets[3 * i + 1];
for (i = 0; i < triplets_map->dims[0]; i++) {
if (max_i < triplets_map->data[i]) {
max_i = triplets_map->data[i];
}
}
gp2tp_map = (int*)malloc(sizeof(int) * (max_i + 1));
for (i = 0; i < max_i + 1; i++) {
gp2tp_map[i] = 0;
}
for (i = 0; i < num_triplets; i++) {
gp2tp_map[triplets[3 * i + 1]] = i;
count = 0;
for (i = 0; i < triplets_map->dims[0]; i++) {
if (triplets_map->data[i] == i) {
gp2tp_map[i] = count;
count++;
}
}
return gp2tp_map;

2
c/anharmonic_h/phonon3_h/collision_matrix.h
View File

@ -7,7 +7,7 @@ void get_collision_matrix(double *collision_matrix,
const Darray *fc3_normal_squared,
const double *frequencies,
const int *triplets,
const int *triplets_map,
const Iarray *triplets_map,
const int *stabilized_gp_map,
const int *ir_grid_points,
const Iarray *rotated_grid_points,

60
c/spglib/kpoint.c
View File

@ -75,11 +75,12 @@ static int get_BZ_triplets_at_q(int triplets[][3],
const int map_triplets[],
const int num_map_triplets,
const int mesh[3]);
static void get_third_q_of_triplets_at_q(int address[3][3],
const int bz_map[],
const int mesh[3],
const int bzmesh[3],
const int bzmesh_double[3]);
static int get_third_q_of_triplets_at_q(int address[3][3],
const int q_index,
const int bz_map[],
const int mesh[3],
const int bzmesh[3],
const int bzmesh_double[3]);
static int get_grid_point(const int grid_double[3],
const int mesh[3]);
static void grid_point_to_grid_double(int grid_double[3],
@ -176,15 +177,13 @@ void kpt_get_grid_points_by_rotations(int rot_grid_points[],
const int address_orig[3],
const MatINT * rot_reciprocal,
const int mesh[3],
const int is_shift[3],
const int bz_map[])
const int is_shift[3])
{
int i, j;
int address_double_orig[3], address_double[3], mesh_double[3], bzmesh_double[3];
int i;
int address_double_orig[3], address_double[3], mesh_double[3];
for (i = 0; i < 3; i++) {
mesh_double[i] = mesh[i] * 2;
bzmesh_double[i] = mesh[i] * 4;
address_double_orig[i] = address_orig[i] * 2 + is_shift[i];
}
for (i = 0; i < rot_reciprocal->size; i++) {
@ -193,12 +192,6 @@ void kpt_get_grid_points_by_rotations(int rot_grid_points[],
address_double_orig);
get_vector_modulo(address_double, mesh_double);
rot_grid_points[i] = get_grid_point(address_double, mesh);
/* for (j = 0; j < 3; j++) { */
/* if (address_double[j] < 0) { */
/* address_double[j] += bzmesh_double[j]; */
/* } */
/* } */
/* rot_grid_points[i] = bz_map[get_grid_point(address_double, mesh_double)]; */
}
}
@ -762,11 +755,16 @@ static int get_BZ_triplets_at_q(int triplets[][3],
address[1][j] = bz_grid_address[ir_grid_points[i]][j];
address[2][j] = - address[0][j] - address[1][j];
}
get_third_q_of_triplets_at_q(address,
bz_map,
mesh,
bzmesh,
bzmesh_double);
for (j = 2; j > -1; j--) {
if (get_third_q_of_triplets_at_q(address,
j,
bz_map,
mesh,
bzmesh,
bzmesh_double) == 0) {
break;
}
}
for (j = 0; j < 3; j++) {
for (k = 0; k < 3; k++) {
address_double[k] = address[j][k] * 2;
@ -783,16 +781,17 @@ static int get_BZ_triplets_at_q(int triplets[][3],
return num_ir;
}
static void get_third_q_of_triplets_at_q(int address[3][3],
const int bz_map[],
const int mesh[3],
const int bzmesh[3],
const int bzmesh_double[3])
static int get_third_q_of_triplets_at_q(int address[3][3],
const int q_index,
const int bz_map[],
const int mesh[3],
const int bzmesh[3],
const int bzmesh_double[3])
{
int i, j, smallest_g, smallest_index, sum_g, delta_g[3];
int bzgp[27], address_double[3];
get_vector_modulo(address[2], mesh);
get_vector_modulo(address[q_index], mesh);
for (i = 0; i < 3; i++) {
delta_g[i] = 0;
for (j = 0; j < 3; j++) {
@ -803,7 +802,8 @@ static void get_third_q_of_triplets_at_q(int address[3][3],
for (i = 0; i < 27; i++) {
for (j = 0; j < 3; j++) {
address_double[j] = (address[2][j] + search_space[i][j] * mesh[j]) * 2;
address_double[j] = (address[q_index][j] +
search_space[i][j] * mesh[j]) * 2;
}
if ((address_double[0] < bzmesh[0]) &&
(address_double[1] < bzmesh[1]) &&
@ -848,8 +848,10 @@ static void get_third_q_of_triplets_at_q(int address[3][3],
}
for (i = 0; i < 3; i++) {
address[2][i] += search_space[smallest_index][i] * mesh[i];
address[q_index][i] += search_space[smallest_index][i] * mesh[i];
}
return smallest_g;
}
static int get_grid_point(const int grid_double[3],

6
c/spglib/spglib.c
View File

@ -522,8 +522,7 @@ void spg_get_grid_points_by_rotations(int rot_grid_points[],
const int num_rot,
SPGCONST int rot_reciprocal[][3][3],
const int mesh[3],
const int is_shift[3],
const int bz_map[])
const int is_shift[3])
{
int i;
MatINT *rot;
@ -536,8 +535,7 @@ void spg_get_grid_points_by_rotations(int rot_grid_points[],
address_orig,
rot,
mesh,
is_shift,
bz_map);
is_shift);
mat_free_MatINT(rot);
}

3
c/spglib_h/kpoint.h
View File

@ -29,8 +29,7 @@ void kpt_get_grid_points_by_rotations(int rot_grid_points[],
const int address_orig[3],
const MatINT * rot_reciprocal,
const int mesh[3],
const int is_shift[3],
const int bz_map[]);
const int is_shift[3]);
int kpt_relocate_BZ_grid_address(int bz_grid_address[][3],
int bz_map[],
SPGCONST int grid_address[][3],

5
c/spglib_h/spglib.h
View File

@ -298,15 +298,14 @@ int spg_get_stabilized_reciprocal_mesh(int grid_address[][3],
SPGCONST double qpoints[][3]);
/* Rotation operations in reciprocal space ``rot_reciprocal`` are applied */
/* to a grid point ``grid_point`` and resulting grid points are stored in */
/* to a grid address ``address_orig`` and resulting grid points are stored in */
/* ``rot_grid_points``. */
void spg_get_grid_points_by_rotations(int rot_grid_points[],
const int address_orig[3],
const int num_rot,
SPGCONST int rot_reciprocal[][3][3],
const int mesh[3],
const int is_shift[3],
const int bz_map[]);
const int is_shift[3]);
/* Grid addresses are relocated inside Brillouin zone. */
/* Number of ir-grid-points inside Brillouin zone is returned. */

4
phonopy/structure/spglib.py
View File

@ -233,7 +233,6 @@ def get_ir_reciprocal_mesh(mesh,
def get_grid_points_by_rotations(address_orig,
reciprocal_rotations,
mesh,
bz_map,
is_shift=np.zeros(3, dtype='intc')):
"""
Rotation operations in reciprocal space ``reciprocal_rotations`` are applied
@ -246,8 +245,7 @@ def get_grid_points_by_rotations(address_orig,
np.array(address_orig, dtype='intc'),
np.array(reciprocal_rotations, dtype='intc', order='C'),
np.array(mesh, dtype='intc'),
np.array(is_shift, dtype='intc'),
bz_map)
np.array(is_shift, dtype='intc'))
return rot_grid_points