diff --git a/c/gridsys.h b/c/gridsys.h index fa87871f..9eb3daf8 100644 --- a/c/gridsys.h +++ b/c/gridsys.h @@ -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, diff --git a/ctest/test_gridsys.cpp b/ctest/test_gridsys.cpp index 633f4945..8e0bdc0b 100644 --- a/ctest/test_gridsys.cpp +++ b/ctest/test_gridsys.cpp @@ -4,6 +4,7 @@ extern "C" { #include #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); + } + } + } +} diff --git a/ctest/utils.c b/ctest/utils.c index f6efd3e8..14364293 100644 --- a/ctest/utils.c +++ b/ctest/utils.c @@ -3,52 +3,6 @@ #include #include -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; diff --git a/ctest/utils.h b/ctest/utils.h index f4a88afe..3b1386f9 100644 --- a/ctest/utils.h +++ b/ctest/utils.h @@ -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__ diff --git a/test/phonon3/test_triplets.py b/test/phonon3/test_triplets.py index c5109d26..5f04dd78 100644 --- a/test/phonon3/test_triplets.py +++ b/test/phonon3/test_triplets.py @@ -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())