Write C gridsys test_gridsys_get_BZ_triplets_at_q_wurtzite_force_SNF

This commit is contained in:
Atsushi Togo 2023-01-26 18:48:40 +09:00
parent 3f29c73c7f
commit 724427f255
5 changed files with 549 additions and 53 deletions

View File

@ -233,6 +233,23 @@ long gridsys_get_bz_grid_addresses(long (*bz_grid_addresses)[3], long *bz_map,
const long Q[3][3], const long PS[3],
const double rec_lattice[3][3],
const long type);
/**
* @brief Find independent q' of (q, q', q'') with given q.
*
* @param map_triplets Mapping table from all grid points to grid points of
* independent q'
* @param map_q Mapping table from all grid points to grid point indices of
* irreducible q-points under the stabilizer subgroup of q
* @param grid_point Grid point of q
* @param D_diag Diagonal elements of diagnoal matrix D of Smith normal form
* @param is_time_reversal With (1) or without (0) time reversal symmetry
* @param num_rot Number of rotation matrices
* @param rec_rotations Transformed rotation matrices in reciprocal space
* @param swappable With (1) or without (0) permutation symmetry between q' and
* q''
* @return long Number of unique element of map_triplets
*/
long gridsys_get_triplets_at_q(long *map_triplets, long *map_q,
const long grid_point, const long D_diag[3],
const long is_time_reversal, const long num_rot,

View File

@ -4,6 +4,7 @@ extern "C" {
#include <math.h>
#include "gridsys.h"
#include "lagrid.h"
#include "utils.h"
}
@ -981,6 +982,8 @@ TEST(test_gridsys, test_gridsys_get_triplets_at_q_AgNO2) {
* @brief gridsys_get_triplets_at_q by wurtzite rotations with and without
* force_SNF (i.e., transformed or not transformed rotations)
* @details Four patterns, is_time_reversal x swappable, are tested.
* The lattices generated with and without force_SNF are the same.
* Therefore numbers of unique triplets should agree, which is this test.
*/
TEST(test_gridsys, test_gridsys_get_triplets_at_q_wurtzite_force_SNF) {
long D_diag[2][3] = {{1, 5, 15}, {5, 5, 3}};
@ -1025,3 +1028,66 @@ TEST(test_gridsys, test_gridsys_get_triplets_at_q_wurtzite_force_SNF) {
}
}
}
/**
* @brief gridsys_get_triplets_at_q by wurtzite rotations with and without
* force_SNF (i.e., transformed or not transformed rotations)
* @details Four patterns, is_time_reversal x swappable, are tested.
* The lattices generated with and without force_SNF are the same.
* Therefore numbers of unique triplets should agree, which is this test.
*/
TEST(test_gridsys, test_gridsys_get_BZ_triplets_at_q_wurtzite_force_SNF) {
long D_diag[2][3] = {{1, 5, 15}, {5, 5, 3}};
long PS[3] = {0, 0, 0};
long Q[2][3][3] = {{{-1, 0, -6}, {0, -1, 0}, {-1, 0, -5}},
{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}};
double rec_lattice[3][3] = {{0.3214400514304082, 0.0, 0.0},
{0.1855835002216734, 0.3711670004433468, 0.0},
{0.0, 0.0, 0.20088388911209323}};
long grid_point = 1;
long map_triplets[75], map_q[75];
long i, j, k, ll, num_triplets_1, num_triplets_2, bz_size;
long ref_unique_elems[4][2] = {{18, 30}, {24, 45}, {30, 30}, {45, 45}};
long is_time_reversal[2] = {1, 0};
long swappable[2] = {1, 0};
long rec_rotations[2][12][3][3];
long triplets[75][3];
long bz_grid_addresses[108][3];
long bz_map[75];
long bzg2grg[108];
for (i = 0; i < 2; i++) {
for (j = 0; j < 12; j++) {
for (k = 0; k < 3; k++) {
for (ll = 0; ll < 3; ll++) {
if (i == 0) {
rec_rotations[i][j][k][ll] =
wurtzite_tilde_rec_rotations_without_time_reversal
[j][k][ll];
} else {
rec_rotations[i][j][k][ll] =
wurtzite_rec_rotations_without_time_reversal[j][k]
[ll];
}
}
}
}
}
for (i = 0; i < 2; i++) { // force_SNF True or False
bz_size =
gridsys_get_bz_grid_addresses(bz_grid_addresses, bz_map, bzg2grg,
D_diag[i], Q[i], PS, rec_lattice, 2);
for (j = 0; j < 2; j++) { // swappable
for (k = 0; k < 2; k++) { // is_time_reversal
num_triplets_1 = gridsys_get_triplets_at_q(
map_triplets, map_q, grid_point, D_diag[i],
is_time_reversal[k], 12, rec_rotations[i], swappable[j]);
num_triplets_2 = gridsys_get_BZ_triplets_at_q(
triplets, grid_point, bz_grid_addresses, bz_map,
map_triplets, 75, D_diag[i], Q[i], 2);
ASSERT_EQ(num_triplets_1, num_triplets_2);
}
}
}
}

View File

@ -3,52 +3,6 @@
#include <stdio.h>
#include <stdlib.h>
void lagmat_multiply_matrix_vector_l3(long v[3], const long a[3][3],
const long b[3]) {
long i;
long 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];
}
}
void lagmat_multiply_matrix_l3(long m[3][3], const long a[3][3],
const long b[3][3]) {
long i, j; /* a_ij */
long c[3][3];
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
c[i][j] = a[i][0] * b[0][j] + a[i][1] * b[1][j] + a[i][2] * b[2][j];
}
}
lagmat_copy_matrix_l3(m, c);
}
void lagmat_copy_matrix_l3(long a[3][3], const long b[3][3]) {
a[0][0] = b[0][0];
a[0][1] = b[0][1];
a[0][2] = b[0][2];
a[1][0] = b[1][0];
a[1][1] = b[1][1];
a[1][2] = b[1][2];
a[2][0] = b[2][0];
a[2][1] = b[2][1];
a[2][2] = b[2][2];
}
long lagmat_check_identity_matrix_l3(const long a[3][3], const long b[3][3]) {
if (a[0][0] - b[0][0] || a[0][1] - b[0][1] || a[0][2] - b[0][2] ||
a[1][0] - b[1][0] || a[1][1] - b[1][1] || a[1][2] - b[1][2] ||
a[2][0] - b[2][0] || a[2][1] - b[2][1] || a[2][2] - b[2][2]) {
return 0;
} else {
return 1;
}
}
long get_num_unique_elems(const long array[], const long array_size) {
long i, num_unique_elems;
long *unique_elems;

View File

@ -1,11 +1,5 @@
#ifndef __test_utils_H__
#define __test_utils_H__
void lagmat_multiply_matrix_vector_l3(long v[3], const long a[3][3],
const long b[3]);
void lagmat_multiply_matrix_l3(long m[3][3], const long a[3][3],
const long b[3][3]);
void lagmat_copy_matrix_l3(long a[3][3], const long b[3][3]);
long lagmat_check_identity_matrix_l3(const long a[3][3], const long b[3][3]);
long get_num_unique_elems(const long array[], const long array_size);
#endif // __test_utils_H__

View File

@ -6,6 +6,7 @@ from phonopy.structure.atoms import PhonopyAtoms
from phonopy.structure.symmetry import Symmetry
from phono3py.phonon3.triplets import (
_get_BZ_triplets_at_q,
_get_triplets_reciprocal_mesh_at_q,
get_triplets_at_q,
)
@ -824,7 +825,12 @@ def test_get_triplets_reciprocal_mesh_at_q_agno2(agno2_cell: PhonopyAtoms, param
def test_get_triplets_reciprocal_mesh_at_q_wurtzite_force(
aln_cell: PhonopyAtoms, params
):
"""Test BZGrid with shift using wurtzite with and without force_SNF."""
"""Test _get_triplets_reciprocal_mesh_at_q using wurtzite.
The lattices generated with and without force_SNF are the same.
Therefore numbers of unique triplets should agree, which is this test.
"""
grid_point = 1
mesh = 14
ph = Phonopy(aln_cell, supercell_matrix=[1, 1, 1], primitive_matrix="auto")
@ -865,3 +871,462 @@ def test_get_triplets_reciprocal_mesh_at_q_wurtzite_force(
ref_unique_elems[params[3] % 4],
[len(np.unique(map_triplets)), len(np.unique(map_q))],
)
@pytest.mark.parametrize(
"params",
[ # force_SNF, swappable, is_time_reversal
(True, True, True, 0),
(True, True, False, 1),
(True, False, True, 2),
(True, False, False, 3),
(False, True, True, 4),
(False, True, False, 5),
(False, False, True, 6),
(False, False, False, 7),
],
)
def test_get_BZ_triplets_at_q(aln_cell: PhonopyAtoms, params):
"""Test _get_BZ_triplets_at_q."""
ref_triplets = [
[
[1, 0, 4],
[1, 1, 3],
[1, 2, 2],
[1, 5, 91],
[1, 7, 90],
[1, 10, 87],
[1, 12, 85],
[1, 13, 84],
[1, 14, 83],
[1, 18, 79],
[1, 19, 77],
[1, 23, 74],
[1, 31, 66],
[1, 32, 65],
[1, 33, 64],
[1, 36, 60],
[1, 38, 59],
[1, 41, 56],
],
[
[1, 0, 4],
[1, 1, 3],
[1, 2, 2],
[1, 5, 91],
[1, 7, 90],
[1, 8, 88],
[1, 10, 87],
[1, 11, 86],
[1, 12, 85],
[1, 13, 84],
[1, 14, 83],
[1, 15, 81],
[1, 17, 80],
[1, 18, 79],
[1, 19, 77],
[1, 23, 74],
[1, 31, 66],
[1, 32, 65],
[1, 33, 64],
[1, 34, 63],
[1, 35, 62],
[1, 36, 60],
[1, 38, 59],
[1, 41, 56],
],
[
[1, 0, 4],
[1, 1, 3],
[1, 2, 2],
[1, 3, 1],
[1, 4, 0],
[1, 5, 91],
[1, 7, 90],
[1, 8, 88],
[1, 10, 87],
[1, 11, 86],
[1, 12, 85],
[1, 13, 84],
[1, 14, 83],
[1, 15, 81],
[1, 17, 80],
[1, 18, 79],
[1, 19, 77],
[1, 21, 76],
[1, 22, 75],
[1, 23, 74],
[1, 31, 66],
[1, 32, 65],
[1, 33, 64],
[1, 34, 63],
[1, 35, 62],
[1, 36, 60],
[1, 38, 59],
[1, 39, 57],
[1, 41, 56],
[1, 42, 55],
],
[
[1, 0, 4],
[1, 1, 3],
[1, 2, 2],
[1, 3, 1],
[1, 4, 0],
[1, 5, 91],
[1, 7, 90],
[1, 8, 88],
[1, 10, 87],
[1, 11, 86],
[1, 12, 85],
[1, 13, 84],
[1, 14, 83],
[1, 15, 81],
[1, 17, 80],
[1, 18, 79],
[1, 19, 77],
[1, 21, 76],
[1, 22, 75],
[1, 23, 74],
[1, 31, 66],
[1, 32, 65],
[1, 33, 64],
[1, 34, 63],
[1, 35, 62],
[1, 36, 60],
[1, 38, 59],
[1, 39, 57],
[1, 41, 56],
[1, 42, 55],
[1, 43, 54],
[1, 44, 53],
[1, 45, 52],
[1, 46, 50],
[1, 48, 49],
[1, 62, 35],
[1, 63, 34],
[1, 64, 33],
[1, 65, 32],
[1, 66, 31],
[1, 67, 29],
[1, 69, 28],
[1, 70, 26],
[1, 72, 25],
[1, 73, 24],
],
[
[1, 0, 4],
[1, 1, 3],
[1, 2, 2],
[1, 5, 30],
[1, 6, 28],
[1, 10, 25],
[1, 11, 23],
[1, 13, 21],
[1, 16, 19],
[1, 31, 66],
[1, 32, 65],
[1, 33, 64],
[1, 36, 92],
[1, 37, 90],
[1, 41, 87],
[1, 42, 85],
[1, 44, 83],
[1, 47, 81],
],
[
[1, 0, 4],
[1, 1, 3],
[1, 2, 2],
[1, 5, 30],
[1, 6, 28],
[1, 10, 25],
[1, 11, 23],
[1, 13, 21],
[1, 16, 19],
[1, 31, 66],
[1, 32, 65],
[1, 33, 64],
[1, 34, 63],
[1, 35, 62],
[1, 36, 92],
[1, 37, 90],
[1, 39, 89],
[1, 40, 88],
[1, 41, 87],
[1, 42, 85],
[1, 44, 83],
[1, 46, 82],
[1, 47, 81],
[1, 48, 80],
],
[
[1, 0, 4],
[1, 1, 3],
[1, 2, 2],
[1, 3, 1],
[1, 4, 0],
[1, 5, 30],
[1, 6, 28],
[1, 8, 27],
[1, 9, 26],
[1, 10, 25],
[1, 11, 23],
[1, 13, 21],
[1, 15, 20],
[1, 16, 19],
[1, 17, 18],
[1, 31, 66],
[1, 32, 65],
[1, 33, 64],
[1, 34, 63],
[1, 35, 62],
[1, 36, 92],
[1, 37, 90],
[1, 39, 89],
[1, 40, 88],
[1, 41, 87],
[1, 42, 85],
[1, 44, 83],
[1, 46, 82],
[1, 47, 81],
[1, 48, 80],
],
[
[1, 0, 4],
[1, 1, 3],
[1, 2, 2],
[1, 3, 1],
[1, 4, 0],
[1, 5, 30],
[1, 6, 28],
[1, 8, 27],
[1, 9, 26],
[1, 10, 25],
[1, 11, 23],
[1, 13, 21],
[1, 15, 20],
[1, 16, 19],
[1, 17, 18],
[1, 31, 66],
[1, 32, 65],
[1, 33, 64],
[1, 34, 63],
[1, 35, 62],
[1, 36, 92],
[1, 37, 90],
[1, 39, 89],
[1, 40, 88],
[1, 41, 87],
[1, 42, 85],
[1, 44, 83],
[1, 46, 82],
[1, 47, 81],
[1, 48, 80],
[1, 62, 35],
[1, 63, 34],
[1, 64, 33],
[1, 65, 32],
[1, 66, 31],
[1, 67, 61],
[1, 68, 59],
[1, 70, 58],
[1, 71, 57],
[1, 72, 56],
[1, 73, 54],
[1, 75, 52],
[1, 77, 51],
[1, 78, 50],
[1, 79, 49],
],
]
ref_ir_weights = [
[2, 2, 1, 8, 4, 8, 4, 8, 8, 4, 4, 2, 4, 4, 2, 4, 2, 4],
[2, 2, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 2, 2, 2, 2, 2, 2, 4, 2, 4],
[
1,
1,
1,
1,
1,
4,
4,
4,
4,
4,
4,
4,
4,
4,
4,
2,
2,
2,
2,
2,
2,
2,
2,
2,
2,
2,
2,
2,
2,
2,
],
[
1,
1,
1,
1,
1,
2,
2,
2,
2,
2,
2,
2,
2,
2,
2,
2,
2,
2,
2,
2,
1,
1,
1,
1,
1,
2,
2,
2,
2,
2,
2,
2,
2,
2,
2,
1,
1,
1,
1,
1,
2,
2,
2,
2,
2,
],
[2, 2, 1, 4, 4, 2, 4, 2, 4, 4, 4, 2, 8, 8, 4, 8, 4, 8],
[2, 2, 1, 4, 4, 2, 4, 2, 4, 2, 2, 2, 2, 2, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4],
[
1,
1,
1,
1,
1,
2,
2,
2,
2,
2,
2,
2,
2,
2,
2,
2,
2,
2,
2,
2,
4,
4,
4,
4,
4,
4,
4,
4,
4,
4,
],
[
1,
1,
1,
1,
1,
2,
2,
2,
2,
2,
2,
2,
2,
2,
2,
1,
1,
1,
1,
1,
2,
2,
2,
2,
2,
2,
2,
2,
2,
2,
1,
1,
1,
1,
1,
2,
2,
2,
2,
2,
2,
2,
2,
2,
2,
],
]
grid_point = 1
mesh = 14
ph = Phonopy(aln_cell, supercell_matrix=[1, 1, 1], primitive_matrix="auto")
bzgrid = BZGrid(
mesh,
lattice=ph.primitive.cell,
symmetry_dataset=ph.primitive_symmetry.dataset,
use_grg=True,
force_SNF=params[0],
is_time_reversal=False,
)
map_triplets, map_q = _get_triplets_reciprocal_mesh_at_q(
grid_point,
bzgrid.D_diag,
bzgrid.rotations,
is_time_reversal=params[2],
swappable=params[1],
)
triplets, ir_weights = _get_BZ_triplets_at_q(grid_point, bzgrid, map_triplets)
for tp in triplets:
print("{%d, %d, %d}," % tuple(tp))
# print(",".join(["%d" % x for x in ir_weights]))
np.testing.assert_equal(ref_triplets[params[3]], triplets)
np.testing.assert_equal(ref_ir_weights[params[3]], ir_weights)
print(np.linalg.inv(ph.primitive.cell).tolist())