Add fortran interface and tests to gridsys_rotate_bz_grid_index

c/gridsys.c

@@ -100,22 +100,6 @@ long gridsys_transform_rotations(long (*transformed_rots)[3][3],
     return succeeded;
 }
 
-double gridsys_get_thm_integration_weight(const double omega,
-                                          const double tetrahedra_omegas[24][4],
-                                          const char function) {
-    return thm_get_integration_weight(omega, tetrahedra_omegas, function);
-}
-
-void gridsys_get_thm_all_relative_grid_address(
-    long relative_grid_address[4][24][4][3]) {
-    thm_get_all_relative_grid_address(relative_grid_address);
-}
-
-long gridsys_get_thm_relative_grid_address(
-    long relative_grid_addresses[24][4][3], const double rec_lattice[3][3]) {
-    return thm_get_relative_grid_address(relative_grid_addresses, rec_lattice);
-}
-
 void gridsys_get_ir_grid_map(long *ir_grid_map, const long (*rotations)[3][3],
                              const long num_rot, const long D_diag[3],
                              const long PS[3]) {
@@ -182,6 +166,35 @@ long gridsys_get_bz_grid_addresses(long (*bz_grid_addresses)[3], long *bz_map,
     return size;
 }
 
+long gridsys_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] = PS[i];
+    }
+
+    rot_bz_gp = bzg_rotate_grid_index(bz_grid_index, rotation, bzgrid);
+
+    free(bzgrid);
+    bzgrid = NULL;
+
+    return rot_bz_gp;
+}
+
 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,
@@ -226,6 +239,22 @@ long gridsys_get_bz_triplets_at_q(long (*triplets)[3], const long grid_point,
     return num_ir;
 }
 
+double gridsys_get_thm_integration_weight(const double omega,
+                                          const double tetrahedra_omegas[24][4],
+                                          const char function) {
+    return thm_get_integration_weight(omega, tetrahedra_omegas, function);
+}
+
+void gridsys_get_thm_all_relative_grid_address(
+    long relative_grid_address[4][24][4][3]) {
+    thm_get_all_relative_grid_address(relative_grid_address);
+}
+
+long gridsys_get_thm_relative_grid_address(
+    long relative_grid_addresses[24][4][3], const double rec_lattice[3][3]) {
+    return thm_get_relative_grid_address(relative_grid_addresses, rec_lattice);
+}
+
 /* relative_grid_addresses are given as P multipled with those from */
 /* dataset, i.e., */
 /*     np.dot(relative_grid_addresses, P.T) */

c/gridsys.h

@@ -112,7 +112,7 @@ long gridsys_get_grid_index_from_address(const long address[3],
  * @param rotation Transformed rotation in reciprocal space tilde-R^T
  * @param D_diag Numbers of divisions along a, b, c directions of GR-grid
  * @param PS Shift in GR-grid
- * @return * long
+ * @return long
  */
 long gridsys_rotate_grid_index(const long grid_index, const long rotation[3][3],
                                const long D_diag[3], const long PS[3]);
@@ -125,7 +125,7 @@ long gridsys_rotate_grid_index(const long grid_index, const long rotation[3][3],
  * @param rotations Rotations in direct space {R}
  * @param num_rot Number of given rotations |{R}|
  * @param is_time_reversal With (1) or without (0) time reversal symmetry
- * @return * long
+ * @return long
  */
 long gridsys_get_reciprocal_point_group(long rec_rotations[48][3][3],
                                         const long (*rotations)[3][3],
@@ -161,40 +161,6 @@ long gridsys_transform_rotations(long (*transformed_rots)[3][3],
                                  const long num_rot, const long D_diag[3],
                                  const long Q[3][3]);
 
-/**
- * @brief Return integration weight of linear tetrahedron method
- *
- * @param omega A frequency point where integration weight is computed
- * @param tetrahedra_omegas Frequencies at vertices of 6 tetrahedra
- * @param function I (delta function) or J (theta function)
- * @return double
- */
-double gridsys_get_thm_integration_weight(const double omega,
-                                          const double tetrahedra_omegas[24][4],
-                                          const char function);
-
-/**
- * @brief Return predefined relative grid addresses of 4 different
- * main diagonals for linear tetrahedron method
- *
- * @param relative_grid_address predefined relative grid addresses for linear
- * tetrahedron method
- */
-void gridsys_get_thm_all_relative_grid_address(
-    long relative_grid_address[4][24][4][3]);
-
-/**
- * @brief Return predefined relative grid addresses of main diagonal determined
- * from reciprocal basis vectors for linear tetrahedron method
- *
- * @param relative_grid_addresses predefined relative grid addresses of given
- * reciprocal basis vectors
- * @param rec_lattice Reciprocal basis vectors in column vectors
- * @return * long
- */
-long gridsys_get_thm_relative_grid_address(
-    long relative_grid_addresses[24][4][3], const double rec_lattice[3][3]);
-
 /**
  * @brief Return mapping table from GR-grid points to GR-ir-grid points
  *
@@ -234,6 +200,27 @@ long gridsys_get_bz_grid_addresses(long (*bz_grid_addresses)[3], long *bz_map,
                                    const double rec_lattice[3][3],
                                    const long bz_grid_type);
 
+/**
+ * @brief Return index of rotated bz grid point
+ *
+ * @param bz_grid_index BZ grid point index
+ * @param rotation Transformed rotation in reciprocal space tilde-R^T
+ * @param bz_grid_addresses BZ grid point adddresses
+ * @param bz_map List of accumulated numbers of BZ grid points from the
+ * first GR grid point to the last grid point. In type-II, [0, 1, 3, 4, ...]
+ * means multiplicities of [1, 2, 1, ...], with len(bz_map)=product(D_diag) + 1.
+ * @param D_diag Numbers of divisions along a, b, c directions of GR-grid
+ * @param PS Shift in GR-grid
+ * @param bz_grid_type Data structure type I (old and sparse) or II (new and
+ * dense, recommended) of bz_map
+ * @return long
+ */
+long gridsys_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);
+
 /**
  * @brief Find independent q' of (q, q', q'') with given q.
  *
@@ -246,8 +233,8 @@ long gridsys_get_bz_grid_addresses(long (*bz_grid_addresses)[3], long *bz_map,
  * @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
+ * @param swappable With (1) or without (0) permutation symmetry between q'
- * q''
+ * and q''
  * @return long Number of unique element of map_triplets
  */
 long gridsys_get_triplets_at_q(long *map_triplets, long *map_q,
@@ -281,6 +268,40 @@ long gridsys_get_bz_triplets_at_q(long (*triplets)[3], const long grid_point,
                                   const long D_diag[3], const long Q[3][3],
                                   const long bz_grid_type);
 
+/**
+ * @brief Return integration weight of linear tetrahedron method
+ *
+ * @param omega A frequency point where integration weight is computed
+ * @param tetrahedra_omegas Frequencies at vertices of 6 tetrahedra
+ * @param function I (delta function) or J (theta function)
+ * @return double
+ */
+double gridsys_get_thm_integration_weight(const double omega,
+                                          const double tetrahedra_omegas[24][4],
+                                          const char function);
+
+/**
+ * @brief Return predefined relative grid addresses of 4 different
+ * main diagonals for linear tetrahedron method
+ *
+ * @param relative_grid_address predefined relative grid addresses for linear
+ * tetrahedron method
+ */
+void gridsys_get_thm_all_relative_grid_address(
+    long relative_grid_address[4][24][4][3]);
+
+/**
+ * @brief Return predefined relative grid addresses of main diagonal determined
+ * from reciprocal basis vectors for linear tetrahedron method
+ *
+ * @param relative_grid_addresses predefined relative grid addresses of given
+ * reciprocal basis vectors
+ * @param rec_lattice Reciprocal basis vectors in column vectors
+ * @return * long
+ */
+long gridsys_get_thm_relative_grid_address(
+    long relative_grid_addresses[24][4][3], const double rec_lattice[3][3]);
+
 long gridsys_get_integration_weight(
     double *iw, char *iw_zero, const double *frequency_points,
     const long num_band0, const long relative_grid_address[24][4][3],

ctest/fortran/gridsysf/test_gridsysf.f90

@@ -1,7 +1,13 @@
 program test_gridsysf
     use, intrinsic :: iso_c_binding
     use gridsysf, only: &
+        gridsys_get_grid_index_from_address, &
+        gridsys_rotate_grid_index, &
+        gridsys_get_double_grid_address, &
+        gridsys_get_grid_address_from_index, &
+        gridsys_get_double_grid_index, &
         gridsys_get_bz_grid_addresses, &
+        gridsys_rotate_bz_grid_index, &
         gridsys_get_triplets_at_q, &
         gridsys_get_bz_triplets_at_q
     implicit none
@@ -31,6 +37,12 @@ program test_gridsysf
                  11, 1, 4, 0, 1, 0, -30, -3, -11, &
                  11, 0, 4, 5, 1, 2, -30, 0, -11], [3, 3, 12])
 
+    write (*, '("[test_gridsys_get_grid_index_from_address]")')
+    call test_gridsys_get_grid_index_from_address()
+    write (*, '("[test_gridsys_rotate_grid_index]")')
+    call test_gridsys_rotate_grid_index()
+    write (*, '("[test_gridsys_rotate_bz_grid_index]")')
+    call test_gridsys_rotate_bz_grid_index()
     write (*, '("[test_gridsys_get_bz_grid_addresses_wurtzite]")')
     call test_gridsys_get_bz_grid_addresses_wurtzite()
     write (*, '("[test_gridsys_get_triplets_at_q_wurtzite]")')
@@ -39,6 +51,130 @@ program test_gridsysf
     call test_gridsys_get_bz_triplets_at_q_wurtzite_force_SNF()
 
 contains
+    subroutine test_gridsys_get_grid_index_from_address() bind(C)
+        integer(c_long) :: address(3)
+        integer(c_long) :: D_diag(3)
+        integer :: i, j, k, grid_index, ret_grid_index
+
+        D_diag(:) = [3, 4, 5]
+        grid_index = 0
+        do k = 0, D_diag(3) - 1
+            address(3) = k
+            do j = 0, D_diag(2) - 1
+                address(2) = j
+                do i = 0, D_diag(1) - 1
+                    address(1) = i
+                    ret_grid_index = gridsys_get_grid_index_from_address(address, D_diag)
+                    call assert_int(ret_grid_index, grid_index)
+                    grid_index = grid_index + 1
+                end do
+            end do
+        end do
+        write (*, '("  OK")')
+    end subroutine test_gridsys_get_grid_index_from_address
+
+    subroutine test_gridsys_rotate_grid_index() bind(C)
+        integer(c_long) :: address(3), d_address(3), rot_address(3)
+        integer(c_long) :: D_diag(3, 2)
+        integer(c_long) :: PS(3, 2, 2)
+        integer(c_long) :: grid_index
+        integer(c_long) :: rotation(3, 3)
+        integer :: i, j, k, rot_grid_index, ref_rot_grid_index, i_tilde, i_ps, i_rot
+        integer :: rot_grid_indices(60, 8)
+
+        rotation(:, :) = reshape([0, 1, 0, -1, 0, 0, 0, 0, -1], [3, 3])
+        D_diag(:, :) = reshape([1, 5, 15, 5, 5, 3], [3, 2])
+        PS(:, :, :) = reshape([ &
+                              0, 0, 0, -2, 0, 5, 0, 0, 0, 0, 0, 1], &
+                              [3, 2, 2])
+
+        do i_tilde = 1, 2
+        do i_ps = 1, 2
+        do i_rot = 1, 12
+            if (i_tilde == 1) then
+                rotation(:, :) = &
+                    wurtzite_tilde_rec_rotations_without_time_reversal(:, :, i_rot)
+            else
+                rotation(:, :) = &
+                    wurtzite_rec_rotations_without_time_reversal(:, :, i_rot)
+            end if
+            do grid_index = 0, 74
+                call gridsys_get_grid_address_from_index(address, grid_index, &
+                                                         D_diag(:, i_tilde))
+                call gridsys_get_double_grid_address(d_address, address, &
+                                                     PS(:, i_ps, i_tilde))
+                rot_address = matmul(transpose(rotation), d_address)
+                ref_rot_grid_index = gridsys_get_double_grid_index( &
+                                     rot_address, D_diag(:, i_tilde), &
+                                     PS(:, i_ps, i_tilde))
+                rot_grid_index = gridsys_rotate_grid_index( &
+                                 grid_index, rotation, D_diag(:, i_tilde), &
+                                 PS(:, i_ps, i_tilde))
+                call assert_int(rot_grid_index, ref_rot_grid_index)
+            end do
+        end do
+        end do
+        end do
+        write (*, '("  OK")')
+    end subroutine test_gridsys_rotate_grid_index
+
+    subroutine test_gridsys_rotate_bz_grid_index() bind(C)
+        integer(c_long) :: address(3), d_address(3), rot_address(3), ref_d_address(3)
+        integer(c_long) :: D_diag(3, 2)
+        integer(c_long) :: PS(3, 2, 2)
+        integer(c_long) :: grid_index
+        integer(c_long) :: rotation(3, 3)
+        integer :: i, j, k, i_tilde, i_ps, i_rot, rot_bz_gp, bz_size
+        integer :: rot_grid_indices(60, 8)
+        integer(c_long) :: Q(3, 3, 2)
+        real(c_double) :: rec_lattice(3, 3)
+        integer(c_long) :: bz_grid_addresses(3, 144)
+        integer(c_long) :: bz_map(76)
+        integer(c_long) :: bzg2grg(144)
+
+        Q(:, :, :) = reshape([-1, 0, -6, 0, -1, 0, -1, 0, -5, &
+                              1, 0, 0, 0, 1, 0, 0, 0, 1], [3, 3, 2])
+        rec_lattice(:, :) = reshape([0.3214400514304082, 0.0, 0.0, &
+                                     0.1855835002216734, 0.3711670004433468, 0.0, &
+                                     0.0, 0.0, 0.20088388911209323], [3, 3])
+        D_diag(:, :) = reshape([1, 5, 15, 5, 5, 3], [3, 2])
+        PS(:, :, :) = reshape([ &
+                              0, 0, 0, -2, 0, 5, 0, 0, 0, 0, 0, 1], &
+                              [3, 2, 2])
+
+        do i_tilde = 1, 2
+        do i_ps = 1, 2
+            bz_size = gridsys_get_bz_grid_addresses(bz_grid_addresses, bz_map, bzg2grg, &
+                                                    D_diag(:, i_tilde), Q(:, :, i_tilde), PS(:, i_ps, i_tilde), &
+                                                    rec_lattice, int(2, c_long))
+            do i_rot = 1, 12
+                if (i_tilde == 1) then
+                    rotation(:, :) = &
+                        wurtzite_tilde_rec_rotations_without_time_reversal(:, :, i_rot)
+                else
+                    rotation(:, :) = &
+                        wurtzite_rec_rotations_without_time_reversal(:, :, i_rot)
+                end if
+                do grid_index = 0, 74
+                    call gridsys_get_double_grid_address( &
+                        d_address, bz_grid_addresses(:, grid_index + 1), &
+                        PS(:, i_ps, i_tilde))
+                    rot_address = matmul(transpose(rotation), d_address)
+                    rot_bz_gp = gridsys_rotate_bz_grid_index( &
+                                grid_index, rotation, &
+                                bz_grid_addresses, bz_map, D_diag(:, i_tilde), &
+                                PS(:, i_ps, i_tilde), int(2, c_long));
+                    call gridsys_get_double_grid_address( &
+                        ref_d_address, &
+                        bz_grid_addresses(:, rot_bz_gp + 1), PS(:, i_ps, i_tilde))
+                    call assert_1D_array_c_long(ref_d_address, rot_address, 3)
+                end do
+            end do
+        end do
+        end do
+        write (*, '("  OK")')
+    end subroutine test_gridsys_rotate_bz_grid_index
+
     subroutine test_gridsys_get_bz_grid_addresses_wurtzite() bind(C)
         integer(c_long) :: bz_size
         integer(c_long) :: PS(3), D_diag(3), Q(3, 3), bz_grid_addresses(3, 144)
@@ -268,11 +404,13 @@ contains
                 num_triplets_2 = gridsys_get_bz_triplets_at_q( &
                                  triplets, grid_point, bz_grid_addresses, bz_map, &
                                  map_triplets, num_gp, D_diag, Q, int(2, c_long))
+                write (*, '("swappable:", i0, ", is_time_reversal:", i0)', advance='no') swappable, is_time_reversal
                 call assert_int(num_triplets_1, num_triplets_2)
                 call assert_int(num_triplets_1, ref_num_triplets(i))
                 shape_of_array(:) = [3, num_triplets_2]
                 call assert_2D_array_c_long( &
                     triplets, ref_triplets(:, :, i), shape_of_array)
+                write (*, '("  OK")')
                 i = i + 1
             end do
         end do

ctest/test_gridsys.cpp

@@ -64,6 +64,7 @@ const long wurtzite_rec_rotations_without_time_reversal[12][3][3] = {
 
 // Transformed point group operations of wurtzite {R^T} (without time reversal)
 // D_diag=[1, 5, 15], Q=[[-1, 0, -6], [0, -1, 0], [-1, 0, -5]]
+// P=[[1, 0, -2], [0, -1, 0], [-3, 0, 5]]
 const long wurtzite_tilde_rec_rotations_without_time_reversal[12][3][3] = {
     {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}},
     {{1, -1, 0}, {-5, 0, -2}, {0, 3, 1}},
@@ -257,72 +258,118 @@ TEST(test_gridsys, test_gridsys_get_grid_index_from_address) {
  * Return grid point index of rotated address of given grid point index.
  */
 TEST(test_gridsys, test_gridsys_rotate_grid_index) {
-    long address[3], rot_address[3];
+    long D_diag[2][3] = {{1, 5, 15}, {5, 5, 3}};
-    long D_diag[3] = {3, 4, 5};
+    long PS[2][2][3] = {{{0, 0, 0}, {-2, 0, 5}}, {{0, 0, 0}, {0, 0, 1}}};
-    long PS[8][3] = {{0, 0, 0}, {0, 0, 3}, {0, 3, 0}, {0, 3, 3},
+    long address[3], rot_address[3], d_address[3];
-                     {3, 0, 0}, {3, 0, 3}, {3, 3, 0}, {3, 3, 3}};
+    long i, j, k, ll, i_rot, i_tilde, i_ps, grid_index;
-    long i, j, k, l;
+    long rec_rotations[2][12][3][3];
-    long grid_index = 0;
-    long rot_grid_index;
-    long rot_grid_indices[8][60] = {
-        {0,  9,  6,  1,  10, 7,  2,  11, 8,  0,  9,  6,  48, 57, 54,
-         49, 58, 55, 50, 59, 56, 48, 57, 54, 36, 45, 42, 37, 46, 43,
-         38, 47, 44, 36, 45, 42, 24, 33, 30, 25, 34, 31, 26, 35, 32,
-         24, 33, 30, 12, 21, 18, 13, 22, 19, 14, 23, 20, 12, 21, 18},
-        {24, 33, 30, 25, 34, 31, 26, 35, 32, 24, 33, 30, 12, 21, 18,
-         13, 22, 19, 14, 23, 20, 12, 21, 18, 0,  9,  6,  1,  10, 7,
-         2,  11, 8,  0,  9,  6,  48, 57, 54, 49, 58, 55, 50, 59, 56,
-         48, 57, 54, 36, 45, 42, 37, 46, 43, 38, 47, 44, 36, 45, 42},
-        {10, 7,  4,  11, 8,  5,  9,  6,  3,  10, 7,  4,  58, 55, 52,
-         59, 56, 53, 57, 54, 51, 58, 55, 52, 46, 43, 40, 47, 44, 41,
-         45, 42, 39, 46, 43, 40, 34, 31, 28, 35, 32, 29, 33, 30, 27,
-         34, 31, 28, 22, 19, 16, 23, 20, 17, 21, 18, 15, 22, 19, 16},
-        {34, 31, 28, 35, 32, 29, 33, 30, 27, 34, 31, 28, 22, 19, 16,
-         23, 20, 17, 21, 18, 15, 22, 19, 16, 10, 7,  4,  11, 8,  5,
-         9,  6,  3,  10, 7,  4,  58, 55, 52, 59, 56, 53, 57, 54, 51,
-         58, 55, 52, 46, 43, 40, 47, 44, 41, 45, 42, 39, 46, 43, 40},
-        {11, 8,  5,  9,  6,  3,  9,  6,  3,  10, 7,  4,  59, 56, 53,
-         57, 54, 51, 57, 54, 51, 58, 55, 52, 47, 44, 41, 45, 42, 39,
-         45, 42, 39, 46, 43, 40, 35, 32, 29, 33, 30, 27, 33, 30, 27,
-         34, 31, 28, 23, 20, 17, 21, 18, 15, 21, 18, 15, 22, 19, 16},
-        {35, 32, 29, 33, 30, 27, 33, 30, 27, 34, 31, 28, 23, 20, 17,
-         21, 18, 15, 21, 18, 15, 22, 19, 16, 11, 8,  5,  9,  6,  3,
-         9,  6,  3,  10, 7,  4,  59, 56, 53, 57, 54, 51, 57, 54, 51,
-         58, 55, 52, 47, 44, 41, 45, 42, 39, 45, 42, 39, 46, 43, 40},
-        {3,  0,  9,  4,  1,  10, 5,  2,  11, 3,  0,  9,  51, 48, 57,
-         52, 49, 58, 53, 50, 59, 51, 48, 57, 39, 36, 45, 40, 37, 46,
-         41, 38, 47, 39, 36, 45, 27, 24, 33, 28, 25, 34, 29, 26, 35,
-         27, 24, 33, 15, 12, 21, 16, 13, 22, 17, 14, 23, 15, 12, 21},
-        {27, 24, 33, 28, 25, 34, 29, 26, 35, 27, 24, 33, 15, 12, 21,
-         16, 13, 22, 17, 14, 23, 15, 12, 21, 3,  0,  9,  4,  1,  10,
-         5,  2,  11, 3,  0,  9,  51, 48, 57, 52, 49, 58, 53, 50, 59,
-         51, 48, 57, 39, 36, 45, 40, 37, 46, 41, 38, 47, 39, 36, 45}};
 
-    // Rutile R^T of a screw operation.
+    for (i = 0; i < 2; i++) {
-    long rotation[3][3] = {{0, 1, 0}, {-1, 0, 0}, {0, 0, -1}};
+        for (j = 0; j < 12; j++) {
-
+            for (k = 0; k < 3; k++) {
-    grid_index = 0;
+                for (ll = 0; ll < 3; ll++) {
-    for (k = 0; k < D_diag[2]; k++) {
+                    if (i == 0) {
-        address[2] = k;
+                        rec_rotations[i][j][k][ll] =
-        for (j = 0; j < D_diag[1]; j++) {
+                            wurtzite_tilde_rec_rotations_without_time_reversal
-            address[1] = j;
+                                [j][k][ll];
-            for (i = 0; i < D_diag[0]; i++) {
+                    } else {
-                address[0] = i;
+                        rec_rotations[i][j][k][ll] =
-                lagmat_multiply_matrix_vector_l3(rot_address, rotation,
+                            wurtzite_rec_rotations_without_time_reversal[j][k]
-                                                 address);
+                                                                        [ll];
-                ASSERT_EQ(
+                    }
-                    (gridsys_get_grid_index_from_address(rot_address, D_diag)),
+                }
-                    (gridsys_rotate_grid_index(grid_index, rotation, D_diag,
-                                               PS[0])));
-                grid_index++;
             }
         }
     }
-    for (l = 0; l < 8; l++) {
+
-        for (grid_index = 0; grid_index < 60; grid_index++) {
+    for (i_tilde = 0; i_tilde < 1; i_tilde++) {
-            ASSERT_EQ(
+        for (i_ps = 0; i_ps < 2; i_ps++) {
-                rot_grid_indices[l][grid_index],
+            for (i_rot = 0; i_rot < 12; i_rot++) {
-                gridsys_rotate_grid_index(grid_index, rotation, D_diag, PS[l]));
+                for (grid_index = 0; grid_index < 75; grid_index++) {
+                    gridsys_get_grid_address_from_index(address, grid_index,
+                                                        D_diag[i_tilde]);
+                    gridsys_get_double_grid_address(d_address, address,
+                                                    PS[i_tilde][i_ps]);
+                    lagmat_multiply_matrix_vector_l3(
+                        rot_address, rec_rotations[i_tilde][i_rot], d_address);
+                    ASSERT_EQ(
+                        (gridsys_get_double_grid_index(
+                            rot_address, D_diag[i_tilde], PS[i_tilde][i_ps])),
+                        (gridsys_rotate_grid_index(
+                            grid_index, rec_rotations[i_tilde][i_rot],
+                            D_diag[i_tilde], PS[i_tilde][i_ps])));
+                }
+            }
+        }
+    }
+}
+
+/**
+ * @brief gridsys_rotate_bz_grid_index
+ * Return bz grid point index of rotated bz address of given bz grid point
+ * index.
+ */
+TEST(test_gridsys, test_gridsys_rotate_bz_grid_index) {
+    long D_diag[2][3] = {{1, 5, 15}, {5, 5, 3}};
+    long PS[2][2][3] = {{{0, 0, 0}, {-2, 0, 5}}, {{0, 0, 0}, {0, 0, 1}}};
+    long address[3], rot_address[3], d_address[3], ref_d_address[3];
+    long i, j, k, ll, i_rot, i_tilde, i_ps, grid_index, rot_bz_gp, bz_size;
+    long rec_rotations[2][12][3][3];
+    long bz_grid_addresses[144][3];
+    long bz_map[76];
+    long bzg2grg[144];
+    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.18558350022167336, 0.37116700044334666, 0.0},
+                                {0.0, 0.0, 0.20088388911209318}};
+
+    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_tilde = 0; i_tilde < 1; i_tilde++) {
+        for (i_ps = 0; i_ps < 2; i_ps++) {
+            bz_size = gridsys_get_bz_grid_addresses(
+                bz_grid_addresses, bz_map, bzg2grg, D_diag[i_tilde], Q[i_tilde],
+                PS[i_tilde][i_ps], rec_lattice, 2);
+            for (i_rot = 0; i_rot < 12; i_rot++) {
+                for (grid_index = 0; grid_index < bz_size; grid_index++) {
+                    gridsys_get_double_grid_address(
+                        d_address, bz_grid_addresses[grid_index],
+                        PS[i_tilde][i_ps]);
+                    lagmat_multiply_matrix_vector_l3(
+                        rot_address, rec_rotations[i_tilde][i_rot], d_address);
+                    rot_bz_gp = gridsys_rotate_bz_grid_index(
+                        grid_index, rec_rotations[i_tilde][i_rot],
+                        bz_grid_addresses, bz_map, D_diag[i_tilde],
+                        PS[i_tilde][i_ps], 2);
+                    gridsys_get_double_grid_address(
+                        ref_d_address, bz_grid_addresses[rot_bz_gp],
+                        PS[i_tilde][i_ps]);
+                    printf("[%d-%d-%d-%d]\n", i_tilde, i_ps, i_rot, grid_index);
+                    for (i = 0; i < 3; i++) {
+                        printf("%d | %d | %d\n", d_address[i], ref_d_address[i],
+                               rot_address[i]);
+                        ASSERT_EQ(ref_d_address[i], rot_address[i]);
+                    }
+                    printf("------\n");
+                }
+            }
         }
     }
 }
@@ -442,8 +489,8 @@ TEST(test_gridsys, test_gridsys_get_snf3x3_AgNO2) {
 
 /**
  * @brief gridsys_transform_rotations
- * Transform {R^T} to {R^T} with respect to transformed microzone basis vectors
+ * Transform {R^T} to {R^T} with respect to transformed microzone basis
- * in GR-grid
+ * vectors in GR-grid
  */
 TEST(test_gridsys, test_gridsys_transform_rotations) {
     long transformed_rotations[8][3][3];
@@ -977,8 +1024,8 @@ TEST(test_gridsys, test_gridsys_get_triplets_at_q_wurtzite_force_SNF) {
 }
 
 /**
- * @brief gridsys_get_BZ_triplets_at_q by wurtzite rotations with and without
+ * @brief gridsys_get_BZ_triplets_at_q by wurtzite rotations with and
- * force_SNF (i.e., transformed or not transformed rotations)
+ * without force_SNF (i.e., transformed or not transformed rotations)
  * @details Four patterns, is_time_reversal x swappable, are tested.
  */
 TEST(test_gridsys, test_gridsys_get_bz_triplets_at_q_wurtzite_force_SNF) {
@@ -1131,8 +1178,8 @@ TEST(test_gridsys, test_gridsys_get_bz_triplets_at_q_wurtzite_force_SNF) {
                 ASSERT_EQ(num_triplets_1, num_triplets_2);
                 ASSERT_EQ(num_triplets_1, ref_num_triplets[count]);
                 for (ll = 0; ll < num_triplets_2; ll++) {
-                    // printf("%ld %ld %ld %ld [%ld %ld %ld] [%ld %ld %ld]\n",
+                    // printf("%ld %ld %ld %ld [%ld %ld %ld] [%ld %ld
-                    // i,
+                    // %ld]\n", i,
                     //        j, k, ll, ref_triplets[count][ll][0],
                     //        ref_triplets[count][ll][1],
                     //        ref_triplets[count][ll][2], triplets[ll][0],

fortran/gridsysf.f90

@@ -64,8 +64,14 @@ module gridsysf
             integer(c_long), intent(in) :: D_diag(3)
         end subroutine gridsys_get_grid_address_from_index
 
-        function gridsys_rotate_grid_index(grid_index, rotation, D_diag, PS) &
+        function gridsys_get_grid_index_from_address(address, D_diag) bind(c)
-            bind(c)
+            import c_long
+            integer(c_long), intent(in) :: address(3)
+            integer(c_long), intent(in) :: D_diag(3)
+            integer(c_long) :: gridsys_get_grid_index_from_address
+        end function gridsys_get_grid_index_from_address
+
+        function gridsys_rotate_grid_index(grid_index, rotation, D_diag, PS) bind(c)
             import c_long
             integer(c_long), intent(in), value :: grid_index
             integer(c_long), intent(in) :: rotation(3, 3)
@@ -112,22 +118,6 @@ module gridsysf
             integer(c_long) :: gridsys_transform_rotations
         end function gridsys_transform_rotations
 
-        function gridsys_get_thm_integration_weight(omega, &
-                                                    tetrahedra_omegas, function_char) bind(c)
-            import c_char, c_double
-            real(c_double), intent(in), value :: omega
-            real(c_double), intent(in) :: tetrahedra_omegas(4, 24)
-            character(kind=c_char), intent(in), value :: function_char
-            real(c_double) :: gridsys_get_thm_integration_weight
-        end function gridsys_get_thm_integration_weight
-
-        subroutine gridsys_get_thm_relative_grid_address(relative_grid_addresses, &
-                                                         rec_lattice) bind(c)
-            import c_long, c_double
-            integer(c_long), intent(inout) :: relative_grid_addresses(3, 4, 24)
-            real(c_double), intent(in) :: rec_lattice(3, 3)
-        end subroutine gridsys_get_thm_relative_grid_address
-
         subroutine gridsys_get_ir_grid_map(ir_grid_map, rotations, num_rot, &
                                            D_diag, PS) bind(c)
             import c_long
@@ -139,7 +129,8 @@ module gridsysf
         end subroutine gridsys_get_ir_grid_map
 
         function gridsys_get_bz_grid_addresses(bz_grid_addresses, bz_map, &
-                                               bzg2grg, D_diag, Q, PS, rec_lattice, type) bind(c)
+                                               bzg2grg, D_diag, Q, PS, &
+                                               rec_lattice, bz_grid_type) bind(c)
             import c_long, c_double
             integer(c_long), intent(inout) :: bz_grid_addresses(3, *)
             integer(c_long), intent(inout) :: bz_map(*)
@@ -148,12 +139,26 @@ module gridsysf
             integer(c_long), intent(in) :: Q(3, 3)
             integer(c_long), intent(in) :: PS(3)
             real(c_double), intent(in) :: rec_lattice(3, 3)
-            integer(c_long), intent(in), value :: type
+            integer(c_long), intent(in), value :: bz_grid_type
             integer(c_long) :: gridsys_get_bz_grid_addresses
         end function gridsys_get_bz_grid_addresses
 
+        function gridsys_rotate_bz_grid_index(bz_grid_index, rotation, bz_grid_addresses, &
+                                              bz_map, D_diag, PS, bz_grid_type) bind(c)
+            import c_long
+            integer(c_long), intent(in), value :: bz_grid_index
+            integer(c_long), intent(in) :: rotation(3, 3)
+            integer(c_long), intent(inout) :: bz_grid_addresses(3, *)
+            integer(c_long), intent(inout) :: bz_map(*)
+            integer(c_long), intent(in) :: D_diag(3)
+            integer(c_long), intent(in) :: PS(3)
+            integer(c_long), intent(in), value :: bz_grid_type
+            integer(c_long) :: gridsys_rotate_bz_grid_index
+        end function gridsys_rotate_bz_grid_index
+
         function gridsys_get_triplets_at_q(map_triplets, map_q, grid_point, &
-                                           D_diag, is_time_reversal, num_rot, rec_rotations, swappable) bind(c)
+                                           D_diag, is_time_reversal, num_rot, &
+                                           rec_rotations, swappable) bind(c)
             import c_long
             integer(c_long), intent(inout) :: map_triplets(*)
             integer(c_long), intent(inout) :: map_q(*)
@@ -181,21 +186,39 @@ module gridsysf
             integer(c_long) :: gridsys_get_bz_triplets_at_q
         end function gridsys_get_bz_triplets_at_q
 
+        function gridsys_get_thm_integration_weight(omega, &
+                                                    tetrahedra_omegas, function_char) bind(c)
+            import c_char, c_double
+            real(c_double), intent(in), value :: omega
+            real(c_double), intent(in) :: tetrahedra_omegas(4, 24)
+            character(kind=c_char), intent(in), value :: function_char
+            real(c_double) :: gridsys_get_thm_integration_weight
+        end function gridsys_get_thm_integration_weight
+
+        subroutine gridsys_get_thm_relative_grid_address(relative_grid_addresses, &
+                                                         rec_lattice) bind(c)
+            import c_long, c_double
+            integer(c_long), intent(inout) :: relative_grid_addresses(3, 4, 24)
+            real(c_double), intent(in) :: rec_lattice(3, 3)
+        end subroutine gridsys_get_thm_relative_grid_address
+
     end interface
 
     public :: gridsys_get_all_grid_addresses, &
               gridsys_get_double_grid_address, &
               gridsys_get_grid_address_from_index, &
+              gridsys_get_grid_index_from_address, &
               gridsys_rotate_grid_index, &
-              gridsys_get_thm_relative_grid_address, &
               gridsys_get_double_grid_index, &
               gridsys_get_reciprocal_point_group, &
               gridsys_get_snf3x3, &
               gridsys_transform_rotations, &
-              gridsys_get_thm_integration_weight, &
               gridsys_get_ir_grid_map, &
               gridsys_get_bz_grid_addresses, &
+              gridsys_rotate_bz_grid_index, &
               gridsys_get_triplets_at_q, &
-              gridsys_get_bz_triplets_at_q
+              gridsys_get_bz_triplets_at_q, &
+              gridsys_get_thm_integration_weight, &
+              gridsys_get_thm_relative_grid_address
 
 end module gridsysf

test/phonon/test_grid.py

@@ -1714,7 +1714,7 @@ def test_relocate_BZ_grid_address_FCC():
 
 
 def test_relocate_BZ_grid_address_aln_grg():
-    """Test of _relocate_BZ_grid_address by zincblende in GR-grid."""
+    """Test of _relocate_BZ_grid_address by wurtzite in GR-grid."""
     D_diag = [1, 5, 15]
     Q = [[-1, 0, -6], [0, -1, 0], [-1, 0, -5]]
     reciprocal_lattice = np.array(
@@ -2011,7 +2011,7 @@ def test_relocate_BZ_grid_address_aln_grg():
 
 
 def test_relocate_BZ_grid_address_aln_553():
-    """Test of _relocate_BZ_grid_address by zincblende."""
+    """Test of _relocate_BZ_grid_address by wurtzite."""
     D_diag = [5, 5, 3]
     Q = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
     reciprocal_lattice = np.array(
@@ -2352,3 +2352,35 @@ def test_relocate_BZ_grid_address_aln_compare():
         assert len(indices[0]) == 1
         check[indices[0][0]] = True
     assert all(check)
+
+
+def test_get_aln_bzgrid(aln_cell: PhonopyAtoms):
+    """Return BZGrid of wurtzite AlN."""
+    mesh = 14
+    symmetry = Symmetry(aln_cell)
+    bzgrid = BZGrid(mesh, lattice=aln_cell.cell, symmetry_dataset=symmetry.dataset)
+    bzgrid = BZGrid(
+        mesh,
+        lattice=aln_cell.cell,
+        symmetry_dataset=symmetry.dataset,
+        use_grg=True,
+        force_SNF=True,
+    )
+    print(bzgrid.D_diag)
+    print(bzgrid.P)
+    print(bzgrid.Q)
+    shifts = np.array(
+        [
+            [0, 0, 0],
+            [0, 0, 1],
+            [0, 1, 0],
+            [0, 1, 1],
+            [1, 0, 0],
+            [1, 0, 1],
+            [1, 1, 0],
+            [1, 1, 1],
+        ]
+    )
+    print(np.dot(bzgrid.P, shifts.T).T)
+
+    return bzgrid