diff --git a/c/gridsys.c b/c/gridsys.c index 1261e387..e15b7512 100644 --- a/c/gridsys.c +++ b/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) */ diff --git a/c/gridsys.h b/c/gridsys.h index 90a2e8a1..92222de8 100644 --- a/c/gridsys.h +++ b/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 - * q'' + * @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, @@ -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], diff --git a/ctest/fortran/gridsysf/test_gridsysf.f90 b/ctest/fortran/gridsysf/test_gridsysf.f90 index c6dd2271..fb881c17 100644 --- a/ctest/fortran/gridsysf/test_gridsysf.f90 +++ b/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 diff --git a/ctest/test_gridsys.cpp b/ctest/test_gridsys.cpp index e5c3d847..811500b5 100644 --- a/ctest/test_gridsys.cpp +++ b/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[3] = {3, 4, 5}; - long PS[8][3] = {{0, 0, 0}, {0, 0, 3}, {0, 3, 0}, {0, 3, 3}, - {3, 0, 0}, {3, 0, 3}, {3, 3, 0}, {3, 3, 3}}; - long i, j, k, l; - 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}}; + 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]; + long i, j, k, ll, i_rot, i_tilde, i_ps, grid_index; + long rec_rotations[2][12][3][3]; - // Rutile R^T of a screw operation. - long rotation[3][3] = {{0, 1, 0}, {-1, 0, 0}, {0, 0, -1}}; - - grid_index = 0; - for (k = 0; k < D_diag[2]; k++) { - address[2] = k; - for (j = 0; j < D_diag[1]; j++) { - address[1] = j; - for (i = 0; i < D_diag[0]; i++) { - address[0] = i; - lagmat_multiply_matrix_vector_l3(rot_address, rotation, - address); - 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 (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 (l = 0; l < 8; l++) { - for (grid_index = 0; grid_index < 60; grid_index++) { - ASSERT_EQ( - rot_grid_indices[l][grid_index], - gridsys_rotate_grid_index(grid_index, rotation, D_diag, PS[l])); + + for (i_tilde = 0; i_tilde < 1; i_tilde++) { + for (i_ps = 0; i_ps < 2; i_ps++) { + for (i_rot = 0; i_rot < 12; i_rot++) { + 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 - * in GR-grid + * Transform {R^T} to {R^T} with respect to transformed microzone basis + * 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 - * force_SNF (i.e., transformed or not transformed rotations) + * @brief gridsys_get_BZ_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. */ 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", - // i, + // printf("%ld %ld %ld %ld [%ld %ld %ld] [%ld %ld + // %ld]\n", i, // j, k, ll, ref_triplets[count][ll][0], // ref_triplets[count][ll][1], // ref_triplets[count][ll][2], triplets[ll][0], diff --git a/fortran/gridsysf.f90 b/fortran/gridsysf.f90 index 0b30c32e..3d85e9a6 100644 --- a/fortran/gridsysf.f90 +++ b/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) & - bind(c) + function gridsys_get_grid_index_from_address(address, D_diag) 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 diff --git a/test/phonon/test_grid.py b/test/phonon/test_grid.py index a9228bc3..e648d359 100644 --- a/test/phonon/test_grid.py +++ b/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