mirror of https://github.com/phonopy/phono3py.git
Put relocate_BZ_grid_address in phono3py
This commit is contained in:
parent
5a5ba5a32d
commit
6bbb489213
|
@ -82,6 +82,8 @@ static PyObject *
|
|||
py_set_triplets_integration_weights(PyObject *self, PyObject *args);
|
||||
static PyObject *
|
||||
py_set_triplets_integration_weights_with_sigma(PyObject *self, PyObject *args);
|
||||
static PyObject * py_relocate_BZ_grid_address(PyObject *self, PyObject *args);
|
||||
|
||||
static PyObject *
|
||||
py_diagonalize_collision_matrix(PyObject *self, PyObject *args);
|
||||
static PyObject * py_pinv_from_eigensolution(PyObject *self, PyObject *args);
|
||||
|
@ -216,6 +218,9 @@ static PyMethodDef _phono3py_methods[] = {
|
|||
(PyCFunction)py_set_triplets_integration_weights_with_sigma,
|
||||
METH_VARARGS,
|
||||
"Integration weights of smearing method for triplets"},
|
||||
{"BZ_grid_address",
|
||||
(PyCFunction)py_relocate_BZ_grid_address, METH_VARARGS,
|
||||
"Relocate grid addresses inside Brillouin zone"},
|
||||
{"diagonalize_collision_matrix",
|
||||
(PyCFunction)py_diagonalize_collision_matrix,
|
||||
METH_VARARGS,
|
||||
|
@ -1773,6 +1778,51 @@ py_set_triplets_integration_weights_with_sigma(PyObject *self, PyObject *args)
|
|||
Py_RETURN_NONE;
|
||||
}
|
||||
|
||||
static PyObject * py_relocate_BZ_grid_address(PyObject *self, PyObject *args)
|
||||
{
|
||||
PyArrayObject* py_bz_grid_address;
|
||||
PyArrayObject* py_bz_map;
|
||||
PyArrayObject* py_grid_address;
|
||||
PyArrayObject* py_mesh;
|
||||
PyArrayObject* py_is_shift;
|
||||
PyArrayObject* py_reciprocal_lattice;
|
||||
|
||||
long (*bz_grid_address)[3];
|
||||
long *bz_map;
|
||||
long (*grid_address)[3];
|
||||
long* mesh;
|
||||
long* is_shift;
|
||||
double (*reciprocal_lattice)[3];
|
||||
long num_ir_gp;
|
||||
|
||||
if (!PyArg_ParseTuple(args, "OOOOOO",
|
||||
&py_bz_grid_address,
|
||||
&py_bz_map,
|
||||
&py_grid_address,
|
||||
&py_mesh,
|
||||
&py_reciprocal_lattice,
|
||||
&py_is_shift)) {
|
||||
return NULL;
|
||||
}
|
||||
|
||||
bz_grid_address = (long(*)[3])PyArray_DATA(py_bz_grid_address);
|
||||
bz_map = (long*)PyArray_DATA(py_bz_map);
|
||||
grid_address = (long(*)[3])PyArray_DATA(py_grid_address);
|
||||
mesh = (long*)PyArray_DATA(py_mesh);
|
||||
is_shift = (long*)PyArray_DATA(py_is_shift);
|
||||
reciprocal_lattice = (double(*)[3])PyArray_DATA(py_reciprocal_lattice);
|
||||
|
||||
num_ir_gp = ph3py_relocate_BZ_grid_address(bz_grid_address,
|
||||
bz_map,
|
||||
grid_address,
|
||||
mesh,
|
||||
reciprocal_lattice,
|
||||
is_shift);
|
||||
|
||||
return PyLong_FromLong(num_ir_gp);
|
||||
}
|
||||
|
||||
|
||||
static PyObject *
|
||||
py_diagonalize_collision_matrix(PyObject *self, PyObject *args)
|
||||
{
|
||||
|
|
274
c/kpoint.c
274
c/kpoint.c
|
@ -45,6 +45,135 @@
|
|||
#define warning_print(...)
|
||||
#endif
|
||||
|
||||
#define KPT_NUM_BZ_SEARCH_SPACE 125
|
||||
static long bz_search_space[KPT_NUM_BZ_SEARCH_SPACE][3] = {
|
||||
{ 0, 0, 0},
|
||||
{ 0, 0, 1},
|
||||
{ 0, 0, 2},
|
||||
{ 0, 0, -2},
|
||||
{ 0, 0, -1},
|
||||
{ 0, 1, 0},
|
||||
{ 0, 1, 1},
|
||||
{ 0, 1, 2},
|
||||
{ 0, 1, -2},
|
||||
{ 0, 1, -1},
|
||||
{ 0, 2, 0},
|
||||
{ 0, 2, 1},
|
||||
{ 0, 2, 2},
|
||||
{ 0, 2, -2},
|
||||
{ 0, 2, -1},
|
||||
{ 0, -2, 0},
|
||||
{ 0, -2, 1},
|
||||
{ 0, -2, 2},
|
||||
{ 0, -2, -2},
|
||||
{ 0, -2, -1},
|
||||
{ 0, -1, 0},
|
||||
{ 0, -1, 1},
|
||||
{ 0, -1, 2},
|
||||
{ 0, -1, -2},
|
||||
{ 0, -1, -1},
|
||||
{ 1, 0, 0},
|
||||
{ 1, 0, 1},
|
||||
{ 1, 0, 2},
|
||||
{ 1, 0, -2},
|
||||
{ 1, 0, -1},
|
||||
{ 1, 1, 0},
|
||||
{ 1, 1, 1},
|
||||
{ 1, 1, 2},
|
||||
{ 1, 1, -2},
|
||||
{ 1, 1, -1},
|
||||
{ 1, 2, 0},
|
||||
{ 1, 2, 1},
|
||||
{ 1, 2, 2},
|
||||
{ 1, 2, -2},
|
||||
{ 1, 2, -1},
|
||||
{ 1, -2, 0},
|
||||
{ 1, -2, 1},
|
||||
{ 1, -2, 2},
|
||||
{ 1, -2, -2},
|
||||
{ 1, -2, -1},
|
||||
{ 1, -1, 0},
|
||||
{ 1, -1, 1},
|
||||
{ 1, -1, 2},
|
||||
{ 1, -1, -2},
|
||||
{ 1, -1, -1},
|
||||
{ 2, 0, 0},
|
||||
{ 2, 0, 1},
|
||||
{ 2, 0, 2},
|
||||
{ 2, 0, -2},
|
||||
{ 2, 0, -1},
|
||||
{ 2, 1, 0},
|
||||
{ 2, 1, 1},
|
||||
{ 2, 1, 2},
|
||||
{ 2, 1, -2},
|
||||
{ 2, 1, -1},
|
||||
{ 2, 2, 0},
|
||||
{ 2, 2, 1},
|
||||
{ 2, 2, 2},
|
||||
{ 2, 2, -2},
|
||||
{ 2, 2, -1},
|
||||
{ 2, -2, 0},
|
||||
{ 2, -2, 1},
|
||||
{ 2, -2, 2},
|
||||
{ 2, -2, -2},
|
||||
{ 2, -2, -1},
|
||||
{ 2, -1, 0},
|
||||
{ 2, -1, 1},
|
||||
{ 2, -1, 2},
|
||||
{ 2, -1, -2},
|
||||
{ 2, -1, -1},
|
||||
{-2, 0, 0},
|
||||
{-2, 0, 1},
|
||||
{-2, 0, 2},
|
||||
{-2, 0, -2},
|
||||
{-2, 0, -1},
|
||||
{-2, 1, 0},
|
||||
{-2, 1, 1},
|
||||
{-2, 1, 2},
|
||||
{-2, 1, -2},
|
||||
{-2, 1, -1},
|
||||
{-2, 2, 0},
|
||||
{-2, 2, 1},
|
||||
{-2, 2, 2},
|
||||
{-2, 2, -2},
|
||||
{-2, 2, -1},
|
||||
{-2, -2, 0},
|
||||
{-2, -2, 1},
|
||||
{-2, -2, 2},
|
||||
{-2, -2, -2},
|
||||
{-2, -2, -1},
|
||||
{-2, -1, 0},
|
||||
{-2, -1, 1},
|
||||
{-2, -1, 2},
|
||||
{-2, -1, -2},
|
||||
{-2, -1, -1},
|
||||
{-1, 0, 0},
|
||||
{-1, 0, 1},
|
||||
{-1, 0, 2},
|
||||
{-1, 0, -2},
|
||||
{-1, 0, -1},
|
||||
{-1, 1, 0},
|
||||
{-1, 1, 1},
|
||||
{-1, 1, 2},
|
||||
{-1, 1, -2},
|
||||
{-1, 1, -1},
|
||||
{-1, 2, 0},
|
||||
{-1, 2, 1},
|
||||
{-1, 2, 2},
|
||||
{-1, 2, -2},
|
||||
{-1, 2, -1},
|
||||
{-1, -2, 0},
|
||||
{-1, -2, 1},
|
||||
{-1, -2, 2},
|
||||
{-1, -2, -2},
|
||||
{-1, -2, -1},
|
||||
{-1, -1, 0},
|
||||
{-1, -1, 1},
|
||||
{-1, -1, 2},
|
||||
{-1, -1, -2},
|
||||
{-1, -1, -1}
|
||||
};
|
||||
|
||||
static MatLONG *get_point_group_reciprocal(const MatLONG * rotations,
|
||||
const long is_time_reversal);
|
||||
static MatLONG *get_point_group_reciprocal_with_q(const MatLONG * rot_reciprocal,
|
||||
|
@ -66,6 +195,14 @@ static long get_ir_reciprocal_mesh_distortion(long grid_address[][3],
|
|||
const long mesh[3],
|
||||
const long is_shift[3],
|
||||
const MatLONG *rot_reciprocal);
|
||||
static long relocate_BZ_grid_address(long bz_grid_address[][3],
|
||||
long bz_map[],
|
||||
KPTCONST long grid_address[][3],
|
||||
const long mesh[3],
|
||||
KPTCONST double rec_lattice[3][3],
|
||||
const long is_shift[3]);
|
||||
static double get_tolerance_for_BZ_reduction(KPTCONST double rec_lattice[3][3],
|
||||
const long mesh[3]);
|
||||
static long get_num_ir(long ir_mapping_table[], const long mesh[3]);
|
||||
static long check_mesh_symmetry(const long mesh[3],
|
||||
const long is_shift[3],
|
||||
|
@ -83,6 +220,10 @@ static void multiply_matrix_vector_ld3(double v[3],
|
|||
static void multiply_matrix_vector_l3(long v[3],
|
||||
KPTCONST long a[3][3],
|
||||
const long b[3]);
|
||||
static void multiply_matrix_vector_d3(double v[3],
|
||||
KPTCONST double a[3][3],
|
||||
const double b[3]);
|
||||
static double norm_squared_d3(const double a[3]);
|
||||
|
||||
|
||||
long kpt_get_irreducible_reciprocal_mesh(long grid_address[][3],
|
||||
|
@ -119,6 +260,21 @@ MatLONG *kpt_get_point_group_reciprocal_with_q(const MatLONG * rot_reciprocal,
|
|||
qpoints);
|
||||
}
|
||||
|
||||
long kpt_relocate_BZ_grid_address(long bz_grid_address[][3],
|
||||
long bz_map[],
|
||||
KPTCONST long grid_address[][3],
|
||||
const long mesh[3],
|
||||
KPTCONST double rec_lattice[3][3],
|
||||
const long is_shift[3])
|
||||
{
|
||||
return relocate_BZ_grid_address(bz_grid_address,
|
||||
bz_map,
|
||||
grid_address,
|
||||
mesh,
|
||||
rec_lattice,
|
||||
is_shift);
|
||||
}
|
||||
|
||||
void kpt_copy_matrix_l3(long a[3][3], KPTCONST long b[3][3])
|
||||
{
|
||||
a[0][0] = b[0][0];
|
||||
|
@ -446,6 +602,103 @@ get_ir_reciprocal_mesh_distortion(long grid_address[][3],
|
|||
return get_num_ir(ir_mapping_table, mesh);
|
||||
}
|
||||
|
||||
static long relocate_BZ_grid_address(long bz_grid_address[][3],
|
||||
long bz_map[],
|
||||
KPTCONST long grid_address[][3],
|
||||
const long mesh[3],
|
||||
KPTCONST double rec_lattice[3][3],
|
||||
const long is_shift[3])
|
||||
{
|
||||
double tolerance, min_distance;
|
||||
double q_vector[3], distance[KPT_NUM_BZ_SEARCH_SPACE];
|
||||
long bzmesh[3], bz_address_double[3];
|
||||
long i, boundary_num_gp, total_num_gp, bzgp, gp, num_bzmesh;
|
||||
long j, k, min_index;
|
||||
|
||||
tolerance = get_tolerance_for_BZ_reduction(rec_lattice, mesh);
|
||||
for (j = 0; j < 3; j++) {
|
||||
bzmesh[j] = mesh[j] * 2;
|
||||
}
|
||||
|
||||
num_bzmesh = bzmesh[0] * bzmesh[1] * bzmesh[2];
|
||||
for (i = 0; i < num_bzmesh; i++) {
|
||||
bz_map[i] = num_bzmesh;
|
||||
}
|
||||
|
||||
boundary_num_gp = 0;
|
||||
total_num_gp = mesh[0] * mesh[1] * mesh[2];
|
||||
|
||||
/* Multithreading doesn't work for this loop since gp calculated */
|
||||
/* with boundary_num_gp is unstable to store bz_grid_address. */
|
||||
for (i = 0; i < total_num_gp; i++) {
|
||||
for (j = 0; j < KPT_NUM_BZ_SEARCH_SPACE; j++) {
|
||||
for (k = 0; k < 3; k++) {
|
||||
q_vector[k] =
|
||||
((grid_address[i][k] + bz_search_space[j][k] * mesh[k]) * 2 +
|
||||
is_shift[k]) / ((double)mesh[k]) / 2;
|
||||
}
|
||||
multiply_matrix_vector_d3(q_vector, rec_lattice, q_vector);
|
||||
distance[j] = norm_squared_d3(q_vector);
|
||||
}
|
||||
min_distance = distance[0];
|
||||
min_index = 0;
|
||||
for (j = 1; j < KPT_NUM_BZ_SEARCH_SPACE; j++) {
|
||||
if (distance[j] < min_distance) {
|
||||
min_distance = distance[j];
|
||||
min_index = j;
|
||||
}
|
||||
}
|
||||
|
||||
for (j = 0; j < KPT_NUM_BZ_SEARCH_SPACE; j++) {
|
||||
if (distance[j] < min_distance + tolerance) {
|
||||
if (j == min_index) {
|
||||
gp = i;
|
||||
} else {
|
||||
gp = boundary_num_gp + total_num_gp;
|
||||
}
|
||||
|
||||
for (k = 0; k < 3; k++) {
|
||||
bz_grid_address[gp][k] =
|
||||
grid_address[i][k] + bz_search_space[j][k] * mesh[k];
|
||||
bz_address_double[k] = bz_grid_address[gp][k] * 2 + is_shift[k];
|
||||
}
|
||||
bzgp = rgd_get_double_grid_index(bz_address_double, bzmesh);
|
||||
bz_map[bzgp] = gp;
|
||||
if (j != min_index) {
|
||||
boundary_num_gp++;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return boundary_num_gp + total_num_gp;
|
||||
}
|
||||
|
||||
static double get_tolerance_for_BZ_reduction(KPTCONST double rec_lattice[3][3],
|
||||
const long mesh[3])
|
||||
{
|
||||
long i, j;
|
||||
double tolerance;
|
||||
double length[3];
|
||||
|
||||
for (i = 0; i < 3; i++) {
|
||||
length[i] = 0;
|
||||
for (j = 0; j < 3; j++) {
|
||||
length[i] += rec_lattice[j][i] * rec_lattice[j][i];
|
||||
}
|
||||
length[i] /= mesh[i] * mesh[i];
|
||||
}
|
||||
tolerance = length[0];
|
||||
for (i = 1; i < 3; i++) {
|
||||
if (tolerance < length[i]) {
|
||||
tolerance = length[i];
|
||||
}
|
||||
}
|
||||
tolerance *= 0.01;
|
||||
|
||||
return tolerance;
|
||||
}
|
||||
|
||||
static long get_num_ir(long ir_mapping_table[], const long mesh[3])
|
||||
{
|
||||
long i, num_ir;
|
||||
|
@ -453,14 +706,14 @@ static long get_num_ir(long ir_mapping_table[], const long mesh[3])
|
|||
num_ir = 0;
|
||||
|
||||
#pragma omp parallel for reduction(+:num_ir)
|
||||
for (i = 0; i < mesh[0] * mesh[1] * (long)(mesh[2]); i++) {
|
||||
for (i = 0; i < mesh[0] * mesh[1] * mesh[2]; i++) {
|
||||
if (ir_mapping_table[i] == i) {
|
||||
num_ir++;
|
||||
}
|
||||
}
|
||||
|
||||
#ifdef _OPENMP
|
||||
for (i = 0; i < mesh[0] * mesh[1] * (long)(mesh[2]); i++) {
|
||||
for (i = 0; i < mesh[0] * mesh[1] * mesh[2]; i++) {
|
||||
ir_mapping_table[i] = ir_mapping_table[ir_mapping_table[i]];
|
||||
}
|
||||
#endif
|
||||
|
@ -599,3 +852,20 @@ static void multiply_matrix_vector_l3(long v[3],
|
|||
for (i = 0; i < 3; i++)
|
||||
v[i] = c[i];
|
||||
}
|
||||
|
||||
static void multiply_matrix_vector_d3(double v[3],
|
||||
KPTCONST double a[3][3],
|
||||
const double b[3])
|
||||
{
|
||||
long i;
|
||||
double c[3];
|
||||
for (i = 0; i < 3; i++)
|
||||
c[i] = a[i][0] * b[0] + a[i][1] * b[1] + a[i][2] * b[2];
|
||||
for (i = 0; i < 3; i++)
|
||||
v[i] = c[i];
|
||||
}
|
||||
|
||||
static double norm_squared_d3(const double a[3])
|
||||
{
|
||||
return a[0] * a[0] + a[1] * a[1] + a[2] * a[2];
|
||||
}
|
||||
|
|
|
@ -55,6 +55,12 @@ MatLONG *kpt_get_point_group_reciprocal_with_q(const MatLONG * rot_reciprocal,
|
|||
const double symprec,
|
||||
const long num_q,
|
||||
KPTCONST double qpoints[][3]);
|
||||
long kpt_relocate_BZ_grid_address(long bz_grid_address[][3],
|
||||
long bz_map[],
|
||||
KPTCONST long grid_address[][3],
|
||||
const long mesh[3],
|
||||
KPTCONST double rec_lattice[3][3],
|
||||
const long is_shift[3]);
|
||||
void kpt_copy_matrix_l3(long a[3][3], KPTCONST long b[3][3]);
|
||||
MatLONG * kpt_alloc_MatLONG(const long size);
|
||||
void kpt_free_MatLONG(MatLONG * matlong);
|
||||
|
|
17
c/phono3py.c
17
c/phono3py.c
|
@ -46,6 +46,7 @@
|
|||
#include "tetrahedron_method.h"
|
||||
#include "triplet.h"
|
||||
#include "triplet_iw.h"
|
||||
#include "kpoint.h"
|
||||
|
||||
#include <stdio.h>
|
||||
|
||||
|
@ -573,7 +574,6 @@ void ph3py_get_integration_weight(double *iw,
|
|||
openmp_per_bands);
|
||||
}
|
||||
|
||||
|
||||
void ph3py_get_integration_weight_with_sigma(double *iw,
|
||||
char *iw_zero,
|
||||
const double sigma,
|
||||
|
@ -599,6 +599,21 @@ void ph3py_get_integration_weight_with_sigma(double *iw,
|
|||
tp_type);
|
||||
}
|
||||
|
||||
long ph3py_relocate_BZ_grid_address(long bz_grid_address[][3],
|
||||
long bz_map[],
|
||||
PHPYCONST long grid_address[][3],
|
||||
const long mesh[3],
|
||||
PHPYCONST double rec_lattice[3][3],
|
||||
const long is_shift[3])
|
||||
{
|
||||
return kpt_relocate_BZ_grid_address(bz_grid_address,
|
||||
bz_map,
|
||||
grid_address,
|
||||
mesh,
|
||||
rec_lattice,
|
||||
is_shift);
|
||||
}
|
||||
|
||||
void ph3py_symmetrize_collision_matrix(double *collision_matrix,
|
||||
const long num_column,
|
||||
const long num_temp,
|
||||
|
|
|
@ -281,6 +281,12 @@ void ph3py_get_integration_weight_with_sigma(double *iw,
|
|||
const double *frequencies,
|
||||
const long num_band,
|
||||
const long tp_type);
|
||||
long ph3py_relocate_BZ_grid_address(long bz_grid_address[][3],
|
||||
long bz_map[],
|
||||
PHPYCONST long grid_address[][3],
|
||||
const long mesh[3],
|
||||
PHPYCONST double rec_lattice[3][3],
|
||||
const long is_shift[3]);
|
||||
|
||||
|
||||
void ph3py_symmetrize_collision_matrix(double *collision_matrix,
|
||||
|
|
|
@ -108,7 +108,7 @@ class Isotope(object):
|
|||
if self._grid_address is None:
|
||||
primitive_lattice = np.linalg.inv(self._primitive.cell)
|
||||
self._grid_address, self._bz_map = get_bz_grid_address(
|
||||
self._mesh, primitive_lattice, with_boundary=True)
|
||||
self._mesh, primitive_lattice)
|
||||
|
||||
if self._phonon_done is None:
|
||||
self._allocate_phonon()
|
||||
|
|
|
@ -549,7 +549,7 @@ class Interaction(object):
|
|||
def _allocate_phonon(self):
|
||||
primitive_lattice = np.linalg.inv(self._primitive.cell)
|
||||
self._grid_address, self._bz_map = get_bz_grid_address(
|
||||
self._mesh, primitive_lattice, with_boundary=True)
|
||||
self._mesh, primitive_lattice)
|
||||
num_band = len(self._primitive) * 3
|
||||
num_grid = len(self._grid_address)
|
||||
self._phonon_done = np.zeros(num_grid, dtype='byte')
|
||||
|
|
|
@ -108,13 +108,12 @@ def get_triplets_at_q(grid_point,
|
|||
point_group,
|
||||
is_time_reversal=is_time_reversal,
|
||||
swappable=swappable)
|
||||
bz_grid_address, bz_map = spglib.relocate_BZ_grid_address(
|
||||
np.array(grid_address, dtype='intc', order='C'),
|
||||
mesh,
|
||||
reciprocal_lattice,
|
||||
is_dense=True)
|
||||
bz_map = np.array(bz_map, dtype='int_')
|
||||
bz_grid_address = np.array(bz_grid_address, dtype='int_', order='C')
|
||||
|
||||
bz_grid_address, bz_map = _relocate_BZ_grid_address(grid_address,
|
||||
mesh,
|
||||
reciprocal_lattice)
|
||||
# bz_map = np.array(bz_map, dtype='int_')
|
||||
# bz_grid_address = np.array(bz_grid_address, dtype='int_', order='C')
|
||||
triplets_at_q, weights = _get_BZ_triplets_at_q(
|
||||
grid_point,
|
||||
bz_grid_address,
|
||||
|
@ -152,13 +151,10 @@ def get_nosym_triplets_at_q(grid_point,
|
|||
mesh,
|
||||
reciprocal_lattice,
|
||||
stores_triplets_map=False):
|
||||
bz_grid_address, bz_map = spglib.relocate_BZ_grid_address(
|
||||
get_grid_address(mesh),
|
||||
bz_grid_address, bz_map = _relocate_BZ_grid_address(
|
||||
np.array(get_grid_address(mesh), dtype='int_', order='C'),
|
||||
mesh,
|
||||
reciprocal_lattice,
|
||||
is_dense=True)
|
||||
bz_map = np.array(bz_map, dtype='int_')
|
||||
bz_grid_address = np.array(bz_grid_address, dtype='int_', order='C')
|
||||
reciprocal_lattice)
|
||||
map_triplets = np.arange(np.prod(mesh), dtype=bz_map.dtype)
|
||||
triplets_at_q, weights = _get_BZ_triplets_at_q(
|
||||
grid_point,
|
||||
|
@ -187,19 +183,12 @@ def get_grid_address(mesh):
|
|||
return grid_address
|
||||
|
||||
|
||||
def get_bz_grid_address(mesh, reciprocal_lattice, with_boundary=False):
|
||||
grid_address = get_grid_address(mesh)
|
||||
bz_grid_address, bz_map = spglib.relocate_BZ_grid_address(
|
||||
grid_address,
|
||||
def get_bz_grid_address(mesh, reciprocal_lattice):
|
||||
bz_grid_address, bz_map = _relocate_BZ_grid_address(
|
||||
np.array(get_grid_address(mesh), dtype='int_', order='C'),
|
||||
mesh,
|
||||
reciprocal_lattice,
|
||||
is_dense=True)
|
||||
bz_map = np.array(bz_map, dtype='int_')
|
||||
bz_grid_address = np.array(bz_grid_address, dtype='int_', order='C')
|
||||
if with_boundary:
|
||||
return bz_grid_address, bz_map
|
||||
else:
|
||||
return bz_grid_address[:np.prod(mesh)]
|
||||
reciprocal_lattice)
|
||||
return bz_grid_address, bz_map
|
||||
|
||||
|
||||
def get_grid_point_from_address_py(address, mesh):
|
||||
|
@ -350,11 +339,11 @@ def get_coarse_ir_grid_points(primitive,
|
|||
ir_grid_weights = ir_grid_weights
|
||||
|
||||
reciprocal_lattice = np.linalg.inv(primitive.cell)
|
||||
bz_grid_address, bz_map = spglib.relocate_BZ_grid_address(
|
||||
np.array(grid_address, dtype='intc', order='C'),
|
||||
|
||||
bz_grid_address, bz_map = _relocate_BZ_grid_address(
|
||||
grid_address,
|
||||
mesh,
|
||||
reciprocal_lattice,
|
||||
is_dense=True)
|
||||
reciprocal_lattice)
|
||||
|
||||
return (ir_grid_points,
|
||||
ir_grid_weights,
|
||||
|
@ -470,6 +459,62 @@ def get_tetrahedra_vertices(relative_address,
|
|||
return vertices
|
||||
|
||||
|
||||
def _relocate_BZ_grid_address(grid_address,
|
||||
mesh,
|
||||
reciprocal_lattice, # column vectors
|
||||
is_shift=None):
|
||||
"""Grid addresses are relocated to be inside first Brillouin zone.
|
||||
|
||||
Number of ir-grid-points inside Brillouin zone is returned.
|
||||
It is assumed that the following arrays have the shapes of
|
||||
bz_grid_address : (num_grid_points_in_FBZ, 3)
|
||||
bz_map (prod(mesh * 2), )
|
||||
|
||||
Note that the shape of grid_address is (prod(mesh), 3) and the
|
||||
addresses in grid_address are arranged to be in parallelepiped
|
||||
made of reciprocal basis vectors. The addresses in bz_grid_address
|
||||
are inside the first Brillouin zone or on its surface. Each
|
||||
address in grid_address is mapped to one of those in
|
||||
bz_grid_address by a reciprocal lattice vector (including zero
|
||||
vector) with keeping element order. For those inside first
|
||||
Brillouin zone, the mapping is one-to-one. For those on the first
|
||||
Brillouin zone surface, more than one addresses in bz_grid_address
|
||||
that are equivalent by the reciprocal lattice translations are
|
||||
mapped to one address in grid_address. In this case, those grid
|
||||
points except for one of them are appended to the tail of this array,
|
||||
for which bz_grid_address has the following data storing:
|
||||
|
||||
|------------------array size of bz_grid_address-------------------------|
|
||||
|--those equivalent to grid_address--|--those on surface except for one--|
|
||||
|-----array size of grid_address-----|
|
||||
|
||||
Number of grid points stored in bz_grid_address is returned.
|
||||
bz_map is used to recover grid point index expanded to include BZ
|
||||
surface from grid address. The grid point indices are mapped to
|
||||
(mesh[0] * 2) x (mesh[1] * 2) x (mesh[2] * 2) space (bz_map).
|
||||
|
||||
"""
|
||||
|
||||
import phono3py._phono3py as phono3c
|
||||
|
||||
if is_shift is None:
|
||||
_is_shift = np.zeros(3, dtype='int_')
|
||||
else:
|
||||
_is_shift = np.array(is_shift, dtype='int_')
|
||||
bz_grid_address = np.zeros((np.prod(np.add(mesh, 1)), 3), dtype='int_')
|
||||
bz_map = np.zeros(np.prod(np.multiply(mesh, 2)), dtype='int_')
|
||||
num_bz_ir = phono3c.BZ_grid_address(
|
||||
bz_grid_address,
|
||||
bz_map,
|
||||
grid_address,
|
||||
np.array(mesh, dtype='int_'),
|
||||
np.array(reciprocal_lattice, dtype='double', order='C'),
|
||||
_is_shift)
|
||||
|
||||
bz_grid_address = np.array(bz_grid_address[:num_bz_ir], dtype='int_')
|
||||
return bz_grid_address, bz_map
|
||||
|
||||
|
||||
def _get_triplets_reciprocal_mesh_at_q(fixed_grid_number,
|
||||
mesh,
|
||||
rotations,
|
||||
|
|
|
@ -23,6 +23,28 @@ def si_pbesol():
|
|||
log_level=1)
|
||||
|
||||
|
||||
@pytest.fixture(scope='session')
|
||||
def si_pbesol_nosym():
|
||||
yaml_filename = os.path.join(current_dir, "phono3py_si_pbesol.yaml")
|
||||
forces_fc3_filename = os.path.join(current_dir, "FORCES_FC3_si_pbesol")
|
||||
return phono3py.load(yaml_filename,
|
||||
forces_fc3_filename=forces_fc3_filename,
|
||||
is_symmetry=False,
|
||||
produce_fc=False,
|
||||
log_level=1)
|
||||
|
||||
|
||||
@pytest.fixture(scope='session')
|
||||
def si_pbesol_nomeshsym():
|
||||
yaml_filename = os.path.join(current_dir, "phono3py_si_pbesol.yaml")
|
||||
forces_fc3_filename = os.path.join(current_dir, "FORCES_FC3_si_pbesol")
|
||||
return phono3py.load(yaml_filename,
|
||||
forces_fc3_filename=forces_fc3_filename,
|
||||
is_mesh_symmetry=False,
|
||||
produce_fc=False,
|
||||
log_level=1)
|
||||
|
||||
|
||||
@pytest.fixture(scope='session')
|
||||
def si_pbesol_compact_fc():
|
||||
yaml_filename = os.path.join(current_dir, "phono3py_si_pbesol.yaml")
|
||||
|
|
|
@ -4,6 +4,9 @@ si_pbesol_kappa_RTA = [107.991, 107.991, 107.991, 0, 0, 0]
|
|||
si_pbesol_kappa_RTA_with_sigmas = [109.6985, 109.6985, 109.6985, 0, 0, 0]
|
||||
si_pbesol_kappa_RTA_iso = [96.92419, 96.92419, 96.92419, 0, 0, 0]
|
||||
si_pbesol_kappa_RTA_with_sigmas_iso = [96.03248, 96.03248, 96.03248, 0, 0, 0]
|
||||
si_pbesol_kappa_RTA_si_nosym = [38.242347, 38.700219, 39.198018,
|
||||
0.3216, 0.207731, 0.283]
|
||||
si_pbesol_kappa_RTA_si_nomeshsym = [38.15450, 38.15450, 38.15450, 0, 0, 0]
|
||||
nacl_pbe_kappa_RTA = [7.72798252, 7.72798252, 7.72798252, 0, 0, 0]
|
||||
nacl_pbe_kappa_RTA_with_sigma = [7.71913708, 7.71913708, 7.71913708, 0, 0, 0]
|
||||
|
||||
|
@ -11,6 +14,7 @@ nacl_pbe_kappa_RTA_with_sigma = [7.71913708, 7.71913708, 7.71913708, 0, 0, 0]
|
|||
aln_lda_kappa_RTA = [203.304059, 203.304059, 213.003125, 0, 0, 0]
|
||||
aln_lda_kappa_RTA_with_sigmas = [213.820000, 213.820000, 224.800121, 0, 0, 0]
|
||||
|
||||
|
||||
def test_kappa_RTA_si(si_pbesol):
|
||||
kappa = _get_kappa(si_pbesol, [9, 9, 9]).ravel()
|
||||
np.testing.assert_allclose(si_pbesol_kappa_RTA, kappa, atol=0.5)
|
||||
|
@ -37,6 +41,7 @@ def test_kappa_RTA_si_with_sigma(si_pbesol):
|
|||
def test_kappa_RTA_si_with_sigma_full_pp(si_pbesol):
|
||||
si_pbesol.sigmas = [0.1, ]
|
||||
kappa = _get_kappa(si_pbesol, [9, 9, 9], is_full_pp=True).ravel()
|
||||
print(kappa)
|
||||
np.testing.assert_allclose(
|
||||
si_pbesol_kappa_RTA_with_sigmas, kappa, atol=0.5)
|
||||
si_pbesol.sigmas = None
|
||||
|
@ -55,6 +60,22 @@ def test_kappa_RTA_si_compact_fc(si_pbesol_compact_fc):
|
|||
np.testing.assert_allclose(si_pbesol_kappa_RTA, kappa, atol=0.5)
|
||||
|
||||
|
||||
def test_kappa_RTA_si_nosym(si_pbesol, si_pbesol_nosym):
|
||||
si_pbesol_nosym.fc2 = si_pbesol.fc2
|
||||
si_pbesol_nosym.fc3 = si_pbesol.fc3
|
||||
kappa = _get_kappa(si_pbesol_nosym, [4, 4, 4]).reshape(-1, 3).sum(axis=1)
|
||||
kappa_ref = np.reshape(si_pbesol_kappa_RTA_si_nosym, (-1, 3)).sum(axis=1)
|
||||
np.testing.assert_allclose(kappa_ref, kappa, atol=0.5)
|
||||
|
||||
|
||||
def test_kappa_RTA_si_nomeshsym(si_pbesol, si_pbesol_nomeshsym):
|
||||
si_pbesol_nomeshsym.fc2 = si_pbesol.fc2
|
||||
si_pbesol_nomeshsym.fc3 = si_pbesol.fc3
|
||||
kappa = _get_kappa(si_pbesol_nomeshsym, [4, 4, 4]).ravel()
|
||||
kappa_ref = si_pbesol_kappa_RTA_si_nomeshsym
|
||||
np.testing.assert_allclose(kappa_ref, kappa, atol=0.5)
|
||||
|
||||
|
||||
def test_kappa_RTA_nacl(nacl_pbe):
|
||||
kappa = _get_kappa(nacl_pbe, [9, 9, 9]).ravel()
|
||||
np.testing.assert_allclose(nacl_pbe_kappa_RTA, kappa, atol=0.5)
|
||||
|
|
Loading…
Reference in New Issue