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

Rotation of grid points considering BZ-surface

This commit is contained in:
Atsushi Togo 2021-04-25 23:28:08 +09:00
parent d2bf885781
commit 19b2b596bb
13 changed files with 429 additions and 178 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

92
c/_phono3py.c
View File

@ -92,7 +92,7 @@ py_get_snf3x3(PyObject *self, PyObject *args);
static PyObject *
py_get_ir_grid_map(PyObject *self, PyObject *args);
static PyObject * py_get_bz_grid_addresses(PyObject *self, PyObject *args);
static PyObject * py_rotate_bz_grid_addresses(PyObject *self, PyObject *args);
static PyObject *
py_diagonalize_collision_matrix(PyObject *self, PyObject *args);
static PyObject * py_pinv_from_eigensolution(PyObject *self, PyObject *args);
@ -251,6 +251,10 @@ static PyMethodDef _phono3py_methods[] = {
(PyCFunction)py_get_bz_grid_addresses,
METH_VARARGS,
"Get grid addresses including Brillouin zone surface"},
{"rotate_bz_grid_index",
(PyCFunction)py_rotate_bz_grid_addresses,
METH_VARARGS,
"Rotate grid point considering Brillouin zone surface"},
{"diagonalize_collision_matrix",
(PyCFunction)py_diagonalize_collision_matrix,
METH_VARARGS,
@ -324,6 +328,7 @@ PyInit__phono3py(void)
#endif
}
static PyObject * py_get_interaction(PyObject *self, PyObject *args)
{
PyArrayObject *py_fc3_normal_squared;
@ -442,6 +447,7 @@ static PyObject * py_get_interaction(PyObject *self, PyObject *args)
Py_RETURN_NONE;
}
static PyObject * py_get_pp_collision(PyObject *self, PyObject *args)
{
PyArrayObject *py_gamma;
@ -578,6 +584,7 @@ static PyObject * py_get_pp_collision(PyObject *self, PyObject *args)
Py_RETURN_NONE;
}
static PyObject * py_get_pp_collision_with_sigma(PyObject *self, PyObject *args)
{
PyArrayObject *py_gamma;
@ -707,6 +714,7 @@ static PyObject * py_get_pp_collision_with_sigma(PyObject *self, PyObject *args)
Py_RETURN_NONE;
}
static PyObject * py_get_imag_self_energy_with_g(PyObject *self, PyObject *args)
{
PyArrayObject *py_gamma;
@ -769,6 +777,7 @@ static PyObject * py_get_imag_self_energy_with_g(PyObject *self, PyObject *args)
Py_RETURN_NONE;
}
static PyObject *
py_get_detailed_imag_self_energy_with_g(PyObject *self, PyObject *args)
{
@ -841,6 +850,7 @@ py_get_detailed_imag_self_energy_with_g(PyObject *self, PyObject *args)
Py_RETURN_NONE;
}
static PyObject * py_get_real_self_energy_at_bands(PyObject *self,
PyObject *args)
{
@ -898,6 +908,7 @@ static PyObject * py_get_real_self_energy_at_bands(PyObject *self,
Py_RETURN_NONE;
}
static PyObject * py_get_real_self_energy_at_frequency_point(PyObject *self,
PyObject *args)
{
@ -958,6 +969,7 @@ static PyObject * py_get_real_self_energy_at_frequency_point(PyObject *self,
Py_RETURN_NONE;
}
static PyObject * py_get_collision_matrix(PyObject *self, PyObject *args)
{
PyArrayObject *py_collision_matrix;
@ -1036,6 +1048,7 @@ static PyObject * py_get_collision_matrix(PyObject *self, PyObject *args)
Py_RETURN_NONE;
}
static PyObject * py_get_reducible_collision_matrix(PyObject *self, PyObject *args)
{
PyArrayObject *py_collision_matrix;
@ -1097,6 +1110,7 @@ static PyObject * py_get_reducible_collision_matrix(PyObject *self, PyObject *ar
Py_RETURN_NONE;
}
static PyObject * py_symmetrize_collision_matrix(PyObject *self, PyObject *args)
{
PyArrayObject *py_collision_matrix;
@ -1130,6 +1144,7 @@ static PyObject * py_symmetrize_collision_matrix(PyObject *self, PyObject *args)
Py_RETURN_NONE;
}
static PyObject * py_expand_collision_matrix(PyObject *self, PyObject *args)
{
PyArrayObject *py_collision_matrix;
@ -1171,6 +1186,7 @@ static PyObject * py_expand_collision_matrix(PyObject *self, PyObject *args)
Py_RETURN_NONE;
}
static PyObject * py_get_isotope_strength(PyObject *self, PyObject *args)
{
PyArrayObject *py_gamma;
@ -1227,6 +1243,7 @@ static PyObject * py_get_isotope_strength(PyObject *self, PyObject *args)
Py_RETURN_NONE;
}
static PyObject * py_get_thm_isotope_strength(PyObject *self, PyObject *args)
{
PyArrayObject *py_gamma;
@ -1294,6 +1311,7 @@ static PyObject * py_get_thm_isotope_strength(PyObject *self, PyObject *args)
Py_RETURN_NONE;
}
static PyObject * py_distribute_fc3(PyObject *self, PyObject *args)
{
PyArrayObject *force_constants_third;
@ -1331,6 +1349,7 @@ static PyObject * py_distribute_fc3(PyObject *self, PyObject *args)
Py_RETURN_NONE;
}
static PyObject * py_rotate_delta_fc2s(PyObject *self, PyObject *args)
{
PyArrayObject *py_fc3;
@ -1382,6 +1401,7 @@ static PyObject * py_rotate_delta_fc2s(PyObject *self, PyObject *args)
Py_RETURN_NONE;
}
static PyObject *
py_get_permutation_symmetry_fc3(PyObject *self, PyObject *args)
{
@ -1402,6 +1422,7 @@ py_get_permutation_symmetry_fc3(PyObject *self, PyObject *args)
Py_RETURN_NONE;
}
static PyObject *
py_get_permutation_symmetry_compact_fc3(PyObject *self, PyObject *args)
{
@ -1446,6 +1467,7 @@ py_get_permutation_symmetry_compact_fc3(PyObject *self, PyObject *args)
Py_RETURN_NONE;
}
static PyObject * py_transpose_compact_fc3(PyObject *self, PyObject *args)
{
PyArrayObject* py_fc3;
@ -1492,6 +1514,7 @@ static PyObject * py_transpose_compact_fc3(PyObject *self, PyObject *args)
Py_RETURN_NONE;
}
static PyObject * py_get_neighboring_grid_points(PyObject *self, PyObject *args)
{
PyArrayObject *py_relative_grid_points;
@ -1543,6 +1566,7 @@ static PyObject * py_get_neighboring_grid_points(PyObject *self, PyObject *args)
Py_RETURN_NONE;
}
static PyObject *
py_get_thm_integration_weights_at_grid_points(PyObject *self, PyObject *args)
{
@ -1608,6 +1632,7 @@ py_get_thm_integration_weights_at_grid_points(PyObject *self, PyObject *args)
Py_RETURN_NONE;
}
static PyObject *
py_tpl_get_triplets_reciprocal_mesh_at_q(PyObject *self, PyObject *args)
{
@ -1654,6 +1679,7 @@ py_tpl_get_triplets_reciprocal_mesh_at_q(PyObject *self, PyObject *args)
return PyLong_FromLong(num_ir);
}
static PyObject * py_tpl_get_BZ_triplets_at_q(PyObject *self, PyObject *args)
{
PyArrayObject *py_triplets;
@ -1707,6 +1733,7 @@ static PyObject * py_tpl_get_BZ_triplets_at_q(PyObject *self, PyObject *args)
return PyLong_FromLong(num_ir);
}
static PyObject *
py_get_triplets_integration_weights(PyObject *self, PyObject *args)
{
@ -1787,6 +1814,7 @@ py_get_triplets_integration_weights(PyObject *self, PyObject *args)
Py_RETURN_NONE;
}
static PyObject *
py_get_triplets_integration_weights_with_sigma(PyObject *self, PyObject *args)
{
@ -1840,6 +1868,7 @@ py_get_triplets_integration_weights_with_sigma(PyObject *self, PyObject *args)
Py_RETURN_NONE;
}
static PyObject *
py_get_grid_index_from_address(PyObject *self, PyObject *args)
{
@ -1887,6 +1916,7 @@ py_get_gr_grid_addresses(PyObject *self, PyObject *args)
Py_RETURN_NONE;
}
static PyObject *
py_transform_rotations(PyObject *self, PyObject *args)
{
@ -1927,8 +1957,8 @@ py_transform_rotations(PyObject *self, PyObject *args)
}
}
static PyObject *
py_get_snf3x3(PyObject *self, PyObject *args)
static PyObject *py_get_snf3x3(PyObject *self, PyObject *args)
{
PyArrayObject* py_D_diag;
PyArrayObject* py_P;
@ -1962,8 +1992,8 @@ py_get_snf3x3(PyObject *self, PyObject *args)
}
}
static PyObject *
py_get_ir_grid_map(PyObject *self, PyObject *args)
static PyObject *py_get_ir_grid_map(PyObject *self, PyObject *args)
{
PyArrayObject* py_grid_mapping_table;
PyArrayObject* py_D_diag;
@ -2000,6 +2030,7 @@ py_get_ir_grid_map(PyObject *self, PyObject *args)
return PyLong_FromLong(num_ir);
}
static PyObject * py_get_bz_grid_addresses(PyObject *self, PyObject *args)
{
PyArrayObject *py_bz_grid_address;
@ -2057,6 +2088,53 @@ static PyObject * py_get_bz_grid_addresses(PyObject *self, PyObject *args)
return PyLong_FromLong(num_total_gp);
}
static PyObject * py_rotate_bz_grid_addresses(PyObject *self, PyObject *args)
{
PyArrayObject *py_bz_grid_addresses;
PyArrayObject *py_rotation;
PyArrayObject *py_bz_map;
PyArrayObject *py_D_diag;
PyArrayObject *py_PS;
long bz_grid_index;
long type;
long (*bz_grid_addresses)[3];
long (*rotation)[3];
long *bz_map;
long *D_diag;
long *PS;
long ret_bz_gp;
if (!PyArg_ParseTuple(args, "lOOOOOl",
&bz_grid_index,
&py_rotation,
&py_bz_grid_addresses,
&py_bz_map,
&py_D_diag,
&py_PS,
&type)) {
return NULL;
}
bz_grid_addresses = (long(*)[3])PyArray_DATA(py_bz_grid_addresses);
rotation = (long(*)[3])PyArray_DATA(py_rotation);
bz_map = (long*)PyArray_DATA(py_bz_map);
D_diag = (long*)PyArray_DATA(py_D_diag);
PS = (long*)PyArray_DATA(py_PS);
ret_bz_gp = ph3py_rotate_bz_grid_index(bz_grid_index,
rotation,
bz_grid_addresses,
bz_map,
D_diag,
PS,
type);
return PyLong_FromLong(ret_bz_gp);
}
static PyObject *
py_diagonalize_collision_matrix(PyObject *self, PyObject *args)
{
@ -2113,6 +2191,7 @@ py_diagonalize_collision_matrix(PyObject *self, PyObject *args)
return PyLong_FromLong(info);
}
static PyObject * py_pinv_from_eigensolution(PyObject *self, PyObject *args)
{
PyArrayObject *py_collision_matrix;
@ -2157,6 +2236,7 @@ static PyObject * py_pinv_from_eigensolution(PyObject *self, PyObject *args)
Py_RETURN_NONE;
}
static PyObject * py_get_default_colmat_solver(PyObject *self, PyObject *args)
{
if (!PyArg_ParseTuple(args, "")) {
@ -2171,6 +2251,7 @@ static PyObject * py_get_default_colmat_solver(PyObject *self, PyObject *args)
}
static void pinv_from_eigensolution(double *data,
const double *eigvals,
const long size,
@ -2259,6 +2340,7 @@ static void pinv_from_eigensolution(double *data,
tmp_data = NULL;
}
static void show_colmat_info(const PyArrayObject *py_collision_matrix,
const long i_sigma,
const long i_temp,

173
c/bzgrid.c
View File

@ -40,6 +40,7 @@
#include "lagrid.h"
#define BZG_NUM_BZ_SEARCH_SPACE 125
#define GRID_TOLERANCE_FACTOR 0.01
static long bz_search_space[BZG_NUM_BZ_SEARCH_SPACE][3] = {
{ 0, 0, 0},
{ 0, 0, 1},
@ -169,10 +170,6 @@ static long bz_search_space[BZG_NUM_BZ_SEARCH_SPACE][3] = {
};
static long get_ir_grid_map(long ir_mapping_table[],
const long D_diag[3],
const long PS[3],
const RotMats *rot_reciprocal);
static void get_bz_grid_addresses_type1(BZGrid *bzgrid,
const long (*gr_grid_addresses)[3],
const long Qinv[3][3]);
@ -193,64 +190,53 @@ static double get_bz_distances(long nint[3],
static void multiply_matrix_vector_d3(double v[3],
const double a[3][3],
const double b[3]);
static long get_inverse_unimodular_matrix_l3(long m[3][3],
const long a[3][3]);
static double norm_squared_d3(const double a[3]);
RotMats *bzg_get_reciprocal_point_group(const long (*rec_rotations_in)[3][3],
const long num_rot,
const long is_time_reversal,
const long is_transpose)
long bzg_rotate_grid_index(const long bz_grid_index,
const long rotation[3][3],
const ConstBZGrid *bzgrid)
{
long i, num_rot_out;
long rec_rotations_out[48][3][3];
RotMats *rec_rotations;
long i, gp, num_bzgp, num_grgp;
long dadrs[3], dadrs_rot[3], adrs_rot[3];
num_rot_out = grg_get_reciprocal_point_group(rec_rotations_out,
rec_rotations_in,
num_rot,
is_time_reversal,
is_transpose);
if (num_rot_out == 0) {
return NULL;
grg_get_double_grid_address(dadrs, bzgrid->addresses[bz_grid_index], bzgrid->PS);
lagmat_multiply_matrix_vector_l3(dadrs_rot, rotation, dadrs);
grg_get_grid_address(adrs_rot, dadrs_rot, bzgrid->PS);
gp = grg_get_grid_index(adrs_rot, bzgrid->D_diag);
if (bzgrid->type == 1) {
if (bzgrid->addresses[gp][0] == adrs_rot[0] &&
bzgrid->addresses[gp][1] == adrs_rot[1] &&
bzgrid->addresses[gp][2] == adrs_rot[2]) {
return gp;
}
num_grgp = bzgrid->D_diag[0] * bzgrid->D_diag[1] * bzgrid->D_diag[2];
num_bzgp = num_grgp * 8;
for (i = bzgrid->gp_map[num_bzgp + gp] + num_grgp;
i < bzgrid->gp_map[num_bzgp + gp + 1] + num_grgp; i++) {
if (bzgrid->addresses[i][0] == adrs_rot[0] &&
bzgrid->addresses[i][1] == adrs_rot[1] &&
bzgrid->addresses[i][2] == adrs_rot[2]) {
return i;
}
}
} else {
for (i = bzgrid->gp_map[gp]; i < bzgrid->gp_map[gp + 1]; i++) {
if (bzgrid->addresses[i][0] == adrs_rot[0] &&
bzgrid->addresses[i][1] == adrs_rot[1] &&
bzgrid->addresses[i][2] == adrs_rot[2]) {
return i;
}
}
}
rec_rotations = bzg_alloc_RotMats(num_rot_out);
for (i = 0; i < num_rot_out; i++) {
lagmat_copy_matrix_l3(rec_rotations->mat[i], rec_rotations_out[i]);
}
return rec_rotations;
/* This should not happen, but possible when bzgrid is ill-defined. */
return bzgrid->gp_map[gp];
}
long bzg_get_ir_grid_map(long *ir_mapping_table,
const long D_diag[3],
const long PS[3],
const long (*rotations_in)[3][3],
const long num_rot,
const long is_time_reversal,
const long is_transpose)
{
long num_ir;
RotMats *rot_reciprocal;
rot_reciprocal = bzg_get_reciprocal_point_group(rotations_in,
num_rot,
is_time_reversal,
is_transpose);
if (rot_reciprocal == NULL) {
return 0;
}
num_ir = get_ir_grid_map(ir_mapping_table,
D_diag,
PS,
rot_reciprocal);
bzg_free_RotMats(rot_reciprocal);
rot_reciprocal = NULL;
return num_ir;
}
long bzg_get_bz_grid_addresses(BZGrid *bzgrid,
const long (*grid_address)[3])
@ -258,7 +244,7 @@ long bzg_get_bz_grid_addresses(BZGrid *bzgrid,
long det;
long Qinv[3][3];
det = bzg_inverse_unimodular_matrix_l3(Qinv, bzgrid->Q);
det = get_inverse_unimodular_matrix_l3(Qinv, bzgrid->Q);
if (det == 0) {
return 0;
}
@ -303,7 +289,7 @@ double bzg_get_tolerance_for_BZ_reduction(const BZGrid *bzgrid)
tolerance = length[i];
}
}
tolerance *= 0.01;
tolerance *= GRID_TOLERANCE_FACTOR;
return tolerance;
}
@ -356,56 +342,6 @@ void bzg_multiply_matrix_vector_ld3(double v[3],
}
}
long bzg_inverse_unimodular_matrix_l3(long m[3][3],
const long a[3][3])
{
long det;
long c[3][3];
det = lagmat_get_determinant_l3(a);
if (labs(det) != 1) {
return 0;
}
c[0][0] = (a[1][1] * a[2][2] - a[1][2] * a[2][1]) / det;
c[1][0] = (a[1][2] * a[2][0] - a[1][0] * a[2][2]) / det;
c[2][0] = (a[1][0] * a[2][1] - a[1][1] * a[2][0]) / det;
c[0][1] = (a[2][1] * a[0][2] - a[2][2] * a[0][1]) / det;
c[1][1] = (a[2][2] * a[0][0] - a[2][0] * a[0][2]) / det;
c[2][1] = (a[2][0] * a[0][1] - a[2][1] * a[0][0]) / det;
c[0][2] = (a[0][1] * a[1][2] - a[0][2] * a[1][1]) / det;
c[1][2] = (a[0][2] * a[1][0] - a[0][0] * a[1][2]) / det;
c[2][2] = (a[0][0] * a[1][1] - a[0][1] * a[1][0]) / det;
lagmat_copy_matrix_l3(m, c);
return det;
}
/* It is assumed that the rotations have been examined by
* grg_transform_rotations, i.e., no broken symmetry of grid is ensured. */
static long get_ir_grid_map(long ir_mapping_table[],
const long D_diag[3],
const long PS[3],
const RotMats *rot_reciprocal)
{
long i, num_ir;
grg_get_ir_grid_map(ir_mapping_table,
rot_reciprocal->mat,
rot_reciprocal->size,
D_diag,
PS);
num_ir = 0;
for (i = 0; i < D_diag[0] * D_diag[1] * D_diag[2]; i++) {
if (ir_mapping_table[i] == i) {
num_ir++;
}
}
return num_ir;
}
static void get_bz_grid_addresses_type1(BZGrid *bzgrid,
const long (*gr_grid_addresses)[3],
const long Qinv[3][3])
@ -462,6 +398,8 @@ static void get_bz_grid_addresses_type1(BZGrid *bzgrid,
}
}
/* This is used in get_BZ_triplets_at_q_type1. */
/* The first one among those found is treated specially, so */
/* excluded from the gp_map address shift by -1. */
id_shift += count - 1;
bzgrid->gp_map[num_bzmesh + i + 1] = id_shift;
}
@ -579,6 +517,31 @@ static void multiply_matrix_vector_d3(double v[3],
}
}
static long get_inverse_unimodular_matrix_l3(long m[3][3],
const long a[3][3])
{
long det;
long c[3][3];
det = lagmat_get_determinant_l3(a);
if (labs(det) != 1) {
return 0;
}
c[0][0] = (a[1][1] * a[2][2] - a[1][2] * a[2][1]) / det;
c[1][0] = (a[1][2] * a[2][0] - a[1][0] * a[2][2]) / det;
c[2][0] = (a[1][0] * a[2][1] - a[1][1] * a[2][0]) / det;
c[0][1] = (a[2][1] * a[0][2] - a[2][2] * a[0][1]) / det;
c[1][1] = (a[2][2] * a[0][0] - a[2][0] * a[0][2]) / det;
c[2][1] = (a[2][0] * a[0][1] - a[2][1] * a[0][0]) / det;
c[0][2] = (a[0][1] * a[1][2] - a[0][2] * a[1][1]) / det;
c[1][2] = (a[0][2] * a[1][0] - a[0][0] * a[1][2]) / det;
c[2][2] = (a[0][0] * a[1][1] - a[0][1] * a[1][0]) / det;
lagmat_copy_matrix_l3(m, c);
return det;
}
static double norm_squared_d3(const double a[3])
{
return a[0] * a[0] + a[1] * a[1] + a[2] * a[2];

16
c/bzgrid.h
View File

@ -92,17 +92,9 @@ typedef struct {
} ConstBZGrid;
RotMats *bzg_get_reciprocal_point_group(const long (*rec_rotations_in)[3][3],
const long num_rot,
const long is_time_reversal,
const long is_transpose);
long bzg_get_ir_grid_map(long *ir_mapping_table,
const long D_diag[3],
const long PS[3],
const long (*rotations_in)[3][3],
const long num_rot,
const long is_time_reversal,
const long is_transpose);
long bzg_rotate_grid_index(const long grid_index,
const long rotation[3][3],
const ConstBZGrid *bzgrid);
long bzg_get_bz_grid_addresses(BZGrid *bzgrid,
const long grid_address[][3]);
double bzg_get_tolerance_for_BZ_reduction(const BZGrid *bzgrid);
@ -111,7 +103,5 @@ void bzg_free_RotMats(RotMats * rotmats);
void bzg_multiply_matrix_vector_ld3(double v[3],
const long a[3][3],
const double b[3]);
long bzg_inverse_unimodular_matrix_l3(long m[3][3],
const long a[3][3]);
#endif

27
c/grgrid.c
View File

@ -42,8 +42,6 @@
#define IDENTITY_TOL 1e-5
static void reduce_grid_address(long address[3], const long D_diag[3]);
static void reduce_double_grid_address(long address_double[3],
const long D_diag[3]);
static long get_double_grid_index(const long address_double[3],
const long D_diag[3],
const long PS[3]);
@ -154,33 +152,31 @@ void grg_get_all_grid_addresses(long (*grid_address)[3], const long D_diag[3])
/* -------------------------------------------------------*/
/* Get address in double grid from address in single grid */
/* -------------------------------------------------------*/
/* This function doubles single-grid address and shifts it by PS. */
/* No modulo operation is applied to returned double-grid address. */
/* address_double : Double grid address. */
/* address : Single grid address. */
/* D_diag : Diagnal elements of D. */
/* PS : Shifts transformed by P. s_i is 0 or 1. */
void grg_get_double_grid_address(long address_double[3],
const long address[3],
const long D_diag[3],
const long PS[3])
{
get_double_grid_address(address_double, address, PS);
reduce_double_grid_address(address_double, D_diag);
}
/* -------------------------------------------------------*/
/* Get address in single grid from address in double grid */
/* -------------------------------------------------------*/
/* This function shifts double-grid adress by PS and divides it by 2. */
/* No modulo operation is applied to returned single-grid address. */
/* address : Single grid address. */
/* address_double : Double grid address. */
/* D_diag : Diagnal elements of D. */
/* PS : Shifts transformed by P. s_i is 0 or 1. */
void grg_get_grid_address(long address[3],
const long address_double[3],
const long D_diag[3],
const long PS[3])
{
get_grid_address(address, address_double, PS);
reduce_grid_address(address, D_diag);
}
/* -------------------------------------------------*/
@ -332,15 +328,6 @@ static void reduce_grid_address(long address[3], const long D_diag[3])
}
}
static void reduce_double_grid_address(long address_double[3],
const long D_diag[3])
{
long i;
for (i = 0; i < 3; i++) {
address_double[i] = lagmat_modulo_l(address_double[i], 2 * D_diag[i]);
}
}
static long get_double_grid_index(const long address_double[3],
const long D_diag[3],
@ -348,10 +335,8 @@ static long get_double_grid_index(const long address_double[3],
{
long address[3];
grg_get_grid_address(address,
address_double,
D_diag,
PS);
get_grid_address(address, address_double, PS);
reduce_grid_address(address, D_diag);
return get_grid_index_from_address(address, D_diag);
}

2
c/grgrid.h
View File

@ -49,11 +49,9 @@ long grg_transform_rotations(long (*transformed_rots)[3][3],
void grg_get_all_grid_addresses(long (*grid_address)[3], const long D_diag[3]);
void grg_get_double_grid_address(long address_double[3],
const long address[3],
const long D_diag[3],
const long PS[3]);
void grg_get_grid_address(long address[3],
const long address_double[3],
const long D_diag[3],
const long PS[3]);
long grg_get_grid_index(const long address[3],
const long D_diag[3]);

3
c/gridsys.c
View File

@ -53,10 +53,9 @@ void gridsys_get_all_grid_addresses(long (*gr_grid_addresses)[3],
void gridsys_get_double_grid_address(long address_double[3],
const long address[3],
const long D_diag[3],
const long PS[3])
{
grg_get_double_grid_address(address_double, address, D_diag, PS);
grg_get_double_grid_address(address_double, address, PS);
}

1
c/gridsys.h
View File

@ -43,7 +43,6 @@ void gridsys_get_all_grid_addresses(long (*gr_grid_addresses)[3],
const long D_diag[3]);
void gridsys_get_double_grid_address(long address_double[3],
const long address[3],
const long D_diag[3],
const long PS[3]);
void gridsys_get_grid_address_from_index(long address[3],
const long grid_index,

33
c/phono3py.c
View File

@ -813,6 +813,39 @@ long ph3py_get_bz_grid_address(long (*bz_grid_addresses)[3],
return size;
}
long ph3py_rotate_bz_grid_index(const long bz_grid_index,
const long rotation[3][3],
const long (*bz_grid_addresses)[3],
const long *bz_map,
const long D_diag[3],
const long PS[3],
const long bz_grid_type)
{
ConstBZGrid *bzgrid;
long i, rot_bz_gp;
if ((bzgrid = (ConstBZGrid*) malloc(sizeof(ConstBZGrid))) == NULL) {
warning_print("Memory could not be allocated.");
return 0;
}
bzgrid->addresses = bz_grid_addresses;
bzgrid->gp_map = bz_map;
bzgrid->type = bz_grid_type;
for (i = 0; i < 3; i++) {
bzgrid->D_diag[i] = D_diag[i];
bzgrid->PS[i] = 0;
}
rot_bz_gp = bzg_rotate_grid_index(bz_grid_index, rotation, bzgrid);
free(bzgrid);
bzgrid = NULL;
return rot_bz_gp;
}
void ph3py_symmetrize_collision_matrix(double *collision_matrix,
const long num_column,
const long num_temp,

8
c/phono3py.h
View File

@ -319,7 +319,13 @@ long ph3py_get_bz_grid_address(long (*bz_grid_addresses)[3],
const long PS[3],
const double rec_lattice[3][3],
const long type);
long ph3py_rotate_bz_grid_index(const long bz_grid_index,
const long rotation[3][3],
const long (*bz_grid_addresses)[3],
const long *bz_map,
const long D_diag[3],
const long PS[3],
const long bz_grid_type);
void ph3py_symmetrize_collision_matrix(double *collision_matrix,
const long num_column,
const long num_temp,

45
c/triplet_grid.c
View File

@ -76,9 +76,13 @@ static void get_BZ_triplets_at_q_type2(long (*triplets)[3],
static double get_squared_distance(const long G[3],
const double LQD_inv[3][3]);
static void get_LQD_inv(double LQD_inv[3][3], const ConstBZGrid *bzgrid);
static RotMats *get_point_group_reciprocal_with_q(const RotMats * rot_reciprocal,
static RotMats *get_reciprocal_point_group_with_q(const RotMats * rot_reciprocal,
const long D_diag[3],
const long grid_point);
static RotMats *get_reciprocal_point_group(const long (*rec_rotations_in)[3][3],
const long num_rot,
const long is_time_reversal,
const long is_transpose);
long tpk_get_ir_triplets_at_q(long *map_triplets,
long *map_q,
@ -92,10 +96,10 @@ long tpk_get_ir_triplets_at_q(long *map_triplets,
long num_ir;
RotMats *rotations;
rotations = bzg_get_reciprocal_point_group(rec_rotations_in,
num_rot,
is_time_reversal,
0);
rotations = get_reciprocal_point_group(rec_rotations_in,
num_rot,
is_time_reversal,
0);
if (rotations == NULL) {
return 0;
}
@ -141,7 +145,7 @@ static long get_ir_triplets_at_q(long *map_triplets,
}
/* Search irreducible q-points (map_q) with a stabilizer. */
rot_reciprocal_q = get_point_group_reciprocal_with_q(rot_reciprocal,
rot_reciprocal_q = get_reciprocal_point_group_with_q(rot_reciprocal,
D_diag,
grid_point);
@ -508,7 +512,7 @@ static void get_LQD_inv(double LQD_inv[3][3], const ConstBZGrid *bzgrid)
}
/* Return NULL if failed */
static RotMats *get_point_group_reciprocal_with_q(const RotMats * rot_reciprocal,
static RotMats *get_reciprocal_point_group_with_q(const RotMats * rot_reciprocal,
const long D_diag[3],
const long grid_point)
{
@ -553,3 +557,30 @@ static RotMats *get_point_group_reciprocal_with_q(const RotMats * rot_reciprocal
return rot_reciprocal_q;
}
static RotMats *get_reciprocal_point_group(const long (*rec_rotations_in)[3][3],
const long num_rot,
const long is_time_reversal,
const long is_transpose)
{
long i, num_rot_out;
long rec_rotations_out[48][3][3];
RotMats *rec_rotations;
num_rot_out = grg_get_reciprocal_point_group(rec_rotations_out,
rec_rotations_in,
num_rot,
is_time_reversal,
is_transpose);
if (num_rot_out == 0) {
return NULL;
}
rec_rotations = bzg_alloc_RotMats(num_rot_out);
for (i = 0; i < num_rot_out; i++) {
lagmat_copy_matrix_l3(rec_rotations->mat[i], rec_rotations_out[i]);
}
return rec_rotations;
}

3
phono3py/api_phono3py.py
View File

@ -134,6 +134,9 @@ class Phono3py(object):
if masses is not None:
self._set_masses(masses)
# Grid
self._bz_grid = None
# Set supercell, primitive, and phonon supercell symmetries
self._symmetry = None
self._primitive_symmetry = None

78
phono3py/phonon/grid.py
View File

@ -514,21 +514,44 @@ def get_grid_points_by_rotations(bz_gp,
else:
rec_rots = bz_grid.rotations
rot_adrs = np.dot(rec_rots, bz_grid.addresses[bz_gp])
grgps = get_grid_point_from_address(rot_adrs, bz_grid.D_diag)
if with_surface:
if bz_grid.is_dense_gp_map:
return _get_bz_grid_points_by_rotations_py(
grgps, rot_adrs, bz_grid)
else:
msg = "with_surface can be used only with is_dense_gp_map=True."
raise RuntimeError(msg)
return _get_grid_points_by_bz_rotations(bz_gp, bz_grid, rec_rots)
else:
return bz_grid.grg2bzg[grgps]
return _get_grid_points_by_rotations(bz_gp, bz_grid, rec_rots)
def _get_bz_grid_points_by_rotations_py(grgps, rot_adrs, bz_grid):
def _get_grid_points_by_rotations(bz_gp, bz_grid, rotations):
"""Grid point rotations without surface treatment"""
rot_adrs = np.dot(rotations, bz_grid.addresses[bz_gp])
grgps = get_grid_point_from_address(rot_adrs, bz_grid.D_diag)
return bz_grid.grg2bzg[grgps]
def _get_grid_points_by_bz_rotations(bz_gp, bz_grid, rotations, lang='C'):
"""Grid point rotations with surface treatment"""
if lang == 'C':
return _get_grid_points_by_bz_rotations_c(bz_gp, bz_grid, rotations)
else:
return _get_grid_points_by_bz_rotations_py(bz_gp, bz_grid, rotations)
def _get_grid_points_by_bz_rotations_c(bz_gp, bz_grid, rotations):
import phono3py._phono3py as phono3c
bzgps = np.zeros(len(rotations), dtype='int_')
for i, r in enumerate(rotations):
bzgps[i] = phono3c.rotate_bz_grid_index(
bz_gp,
r,
bz_grid.addresses,
bz_grid.gp_map,
bz_grid.D_diag,
bz_grid.PS,
bz_grid.is_dense_gp_map * 1 + 1)
return bzgps
def _get_grid_points_by_bz_rotations_py(bz_gp, bz_grid, rotations):
"""Return BZ-grid point indices generated by rotations
Rotated BZ-grid addresses are compared with translationally
@ -536,15 +559,32 @@ def _get_bz_grid_points_by_rotations_py(grgps, rot_adrs, bz_grid):
indices.
"""
rot_adrs = np.dot(rotations, bz_grid.addresses[bz_gp])
grgps = get_grid_point_from_address(rot_adrs, bz_grid.D_diag)
bzgps = np.zeros(len(grgps), dtype='int_')
for i, (gp, adrs) in enumerate(zip(grgps, rot_adrs)):
indices = np.where(
(bz_grid.addresses[bz_grid.gp_map[gp] : bz_grid.gp_map[gp + 1]]
== adrs).all(axis=1))[0]
if len(indices) == 0:
msg = "with_surface did not work properly."
raise RuntimeError(msg)
bzgps[i] = bz_grid.gp_map[gp] + indices[0]
if bz_grid.is_dense_gp_map:
for i, (gp, adrs) in enumerate(zip(grgps, rot_adrs)):
indices = np.where(
(bz_grid.addresses[bz_grid.gp_map[gp] : bz_grid.gp_map[gp + 1]]
== adrs).all(axis=1))[0]
if len(indices) == 0:
msg = "with_surface did not work properly."
raise RuntimeError(msg)
bzgps[i] = bz_grid.gp_map[gp] + indices[0]
else:
num_grgp = np.prod(bz_grid.D_diag)
num_bzgp = num_grgp * 8
for i, (gp, adrs) in enumerate(zip(grgps, rot_adrs)):
gps = (np.arange(
bz_grid.gp_map[num_bzgp + gp],
bz_grid.gp_map[num_bzgp + gp + 1]) + num_grgp).tolist()
gps.insert(0, gp)
indices = np.where((bz_grid.addresses[gps] == adrs).all(axis=1))[0]
if len(indices) == 0:
msg = "with_surface did not work properly."
raise RuntimeError(msg)
bzgps[i] = gps[indices[0]]
return bzgps

126
test/phonon/test_grid.py
View File

@ -1,7 +1,9 @@
import numpy as np
from phonopy.structure.tetrahedron_method import TetrahedronMethod
from phono3py.phonon.grid import (get_grid_point_from_address,
get_grid_point_from_address_py, BZGrid)
from phono3py.phonon.grid import (
get_grid_point_from_address, get_grid_point_from_address_py, BZGrid,
_get_grid_points_by_rotations, _get_grid_points_by_bz_rotations_c,
_get_grid_points_by_bz_rotations_py)
def test_get_grid_point_from_address(agno2_cell):
@ -179,3 +181,123 @@ def test_SNF_tetrahedra_relative_grid(aln_lda):
atol=1e-8)
# print(np.dot(mtet, mlat.T))
# print(np.dot(np.dot(ptet, bz_grid.QDinv.T), plat.T))
def test_get_grid_points_by_bz_rotations(si_pbesol_111):
"""Rotate grid point by rotations with and without considering BZ surface
The following three methods are tested between type-1 and type-2.
_get_grid_points_by_rotations
_get_grid_points_by_bz_rotations_c
_get_grid_points_by_bz_rotations_py
"""
ref10_type1 = [10, 26, 10, 26, 26, 10, 26, 10, 88, 80, 200, 208,
200, 208, 88, 80, 208, 88, 80, 200, 208, 88, 80, 200,
26, 10, 26, 10, 10, 26, 10, 26, 200, 208, 88, 80,
88, 80, 200, 208, 80, 200, 208, 88, 80, 200, 208, 88]
ref12_type2 = [12, 39, 12, 39, 39, 12, 39, 12, 122, 109, 265, 278,
265, 278, 122, 109, 278, 122, 109, 265, 278, 122, 109, 265,
39, 12, 39, 12, 12, 39, 12, 39, 265, 278, 122, 109,
122, 109, 265, 278, 109, 265, 278, 122, 109, 265, 278, 122]
ref10_bz_type1 = [
10, 26, 260, 270, 269, 258, 271, 259, 88, 285, 200, 328,
322, 208, 291, 286, 327, 292, 287, 321, 326, 290, 80, 323,
269, 258, 271, 259, 10, 26, 260, 270, 200, 328, 88, 285,
291, 286, 322, 208, 80, 323, 326, 290, 287, 321, 327, 292]
ref12_bz_type2 = [
12, 39, 15, 41, 40, 13, 42, 14, 122, 110, 265, 281,
267, 278, 124, 111, 280, 125, 112, 266, 279, 123, 109, 268,
40, 13, 42, 14, 12, 39, 15, 41, 265, 281, 122, 110,
124, 111, 267, 278, 109, 268, 279, 123, 112, 266, 280, 125]
lat = si_pbesol_111.primitive.cell
mesh = 20
bz_grid_type1 = BZGrid(
mesh,
lattice=lat,
primitive_symmetry=si_pbesol_111.primitive_symmetry,
is_dense_gp_map=False)
bz_grid_type2 = BZGrid(
mesh,
lattice=lat,
primitive_symmetry=si_pbesol_111.primitive_symmetry,
is_dense_gp_map=True)
# Check data consistency by reducing to GR-grid.
# Grid point 10 in type-1 and 12 in type-2 are the same points in GR-grid.
assert bz_grid_type1.bzg2grg[10] == bz_grid_type2.bzg2grg[12]
np.testing.assert_equal(bz_grid_type1.bzg2grg[ref10_type1],
bz_grid_type2.bzg2grg[ref12_type2])
np.testing.assert_equal(bz_grid_type1.bzg2grg[ref10_type1],
bz_grid_type1.bzg2grg[ref10_bz_type1])
np.testing.assert_equal(bz_grid_type1.bzg2grg[ref10_type1],
bz_grid_type2.bzg2grg[ref12_bz_type2])
bzgps = _get_grid_points_by_rotations(
10, bz_grid_type1, bz_grid_type1.rotations)
np.testing.assert_equal(bzgps, ref10_type1)
bzgps = _get_grid_points_by_rotations(
12, bz_grid_type2, bz_grid_type2.rotations)
np.testing.assert_equal(bzgps, ref12_type2)
bzgps = _get_grid_points_by_bz_rotations_c(
10, bz_grid_type1, bz_grid_type1.rotations)
np.testing.assert_equal(bzgps, ref10_bz_type1)
bzgps = _get_grid_points_by_bz_rotations_c(
12, bz_grid_type2, bz_grid_type2.rotations)
np.testing.assert_equal(bzgps, ref12_bz_type2)
bzgps = _get_grid_points_by_bz_rotations_py(
10, bz_grid_type1, bz_grid_type1.rotations)
np.testing.assert_equal(bzgps, ref10_bz_type1)
bzgps = _get_grid_points_by_bz_rotations_py(
12, bz_grid_type2, bz_grid_type2.rotations)
np.testing.assert_equal(bzgps, ref12_bz_type2)
# Exhaustive consistency check among methods
for grgp in range(np.prod(bz_grid_type1.D_diag)):
bzgp_type1 = bz_grid_type1.grg2bzg[grgp]
bzgp_type2 = bz_grid_type2.grg2bzg[grgp]
rot_grgps = bz_grid_type1.bzg2grg[
_get_grid_points_by_rotations(
bzgp_type1, bz_grid_type1, bz_grid_type1.rotations)]
np.testing.assert_equal(
rot_grgps,
bz_grid_type2.bzg2grg[_get_grid_points_by_rotations(
bzgp_type2, bz_grid_type2, bz_grid_type2.rotations)])
np.testing.assert_equal(
_get_grid_points_by_bz_rotations_c(
bzgp_type1, bz_grid_type1, bz_grid_type1.rotations),
_get_grid_points_by_bz_rotations_py(
bzgp_type1, bz_grid_type1, bz_grid_type1.rotations))
np.testing.assert_equal(
_get_grid_points_by_bz_rotations_c(
bzgp_type2, bz_grid_type2, bz_grid_type2.rotations),
_get_grid_points_by_bz_rotations_py(
bzgp_type2, bz_grid_type2, bz_grid_type2.rotations))
np.testing.assert_equal(
rot_grgps,
bz_grid_type1.bzg2grg[_get_grid_points_by_bz_rotations_c(
bzgp_type1, bz_grid_type1, bz_grid_type1.rotations)])
np.testing.assert_equal(
rot_grgps,
bz_grid_type2.bzg2grg[_get_grid_points_by_bz_rotations_c(
bzgp_type2, bz_grid_type2, bz_grid_type2.rotations)])
# for gps in bzgps.reshape(-1, 12):
# print("".join(["%d, " % gp for gp in gps]))