From 71bb34aa4a95211f02af95514eb520a8ebff3ab7 Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Tue, 10 Jun 2025 18:32:13 +0900 Subject: [PATCH] Fix distance calculation in BZ --- c/_phono3py.cpp | 20 ++++---- c/gridsys.c | 7 +-- c/gridsys.h | 3 +- c/phono3py.c | 3 +- c/phono3py.h | 3 +- c/recgrid.c | 15 ++---- c/triplet_grid.c | 12 +++-- phono3py/phonon/grid.py | 58 +++++++++++++++++---- phono3py/phonon3/triplets.py | 14 +++-- test/phonon3/test_triplets.py | 97 ++++++++++++++++++----------------- 10 files changed, 140 insertions(+), 92 deletions(-) diff --git a/c/_phono3py.cpp b/c/_phono3py.cpp index 032d2393..c986732b 100644 --- a/c/_phono3py.cpp +++ b/c/_phono3py.cpp @@ -804,13 +804,11 @@ int64_t py_tpl_get_triplets_reciprocal_mesh_at_q( return num_ir; } -int64_t py_tpl_get_BZ_triplets_at_q(nb::ndarray<> py_triplets, - int64_t grid_point, - nb::ndarray<> py_bz_grid_address, - nb::ndarray<> py_bz_map, - nb::ndarray<> py_map_triplets, - nb::ndarray<> py_D_diag, nb::ndarray<> py_Q, - int64_t bz_grid_type) { +int64_t py_tpl_get_BZ_triplets_at_q( + nb::ndarray<> py_triplets, int64_t grid_point, + nb::ndarray<> py_bz_grid_address, nb::ndarray<> py_bz_map, + nb::ndarray<> py_map_triplets, nb::ndarray<> py_D_diag, nb::ndarray<> py_Q, + nb::ndarray<> py_reciprocal_lattice, int64_t bz_grid_type) { int64_t (*triplets)[3]; int64_t (*bz_grid_address)[3]; int64_t *bz_map; @@ -818,6 +816,7 @@ int64_t py_tpl_get_BZ_triplets_at_q(nb::ndarray<> py_triplets, int64_t num_map_triplets; int64_t *D_diag; int64_t (*Q)[3]; + double (*reciprocal_lattice)[3]; int64_t num_ir; triplets = (int64_t (*)[3])py_triplets.data(); @@ -827,10 +826,11 @@ int64_t py_tpl_get_BZ_triplets_at_q(nb::ndarray<> py_triplets, num_map_triplets = (int64_t)py_map_triplets.shape(0); D_diag = (int64_t *)py_D_diag.data(); Q = (int64_t (*)[3])py_Q.data(); + reciprocal_lattice = (double (*)[3])py_reciprocal_lattice.data(); - num_ir = ph3py_get_BZ_triplets_at_q(triplets, grid_point, bz_grid_address, - bz_map, map_triplets, num_map_triplets, - D_diag, Q, bz_grid_type); + num_ir = ph3py_get_BZ_triplets_at_q( + triplets, grid_point, bz_grid_address, bz_map, map_triplets, + num_map_triplets, D_diag, Q, reciprocal_lattice, bz_grid_type); return num_ir; } diff --git a/c/gridsys.c b/c/gridsys.c index 067ebf66..b7100c4c 100644 --- a/c/gridsys.c +++ b/c/gridsys.c @@ -138,7 +138,7 @@ int64_t gridsys_get_bz_grid_addresses( if (!niggli_reduce(niggli_lattice, GRIDSYS_NIGGLI_TOLERANCE)) { return 0; } - if (!lagmat_inverse_matrix_d3(inv_Lr, (double(*)[3])niggli_lattice, + if (!lagmat_inverse_matrix_d3(inv_Lr, (double (*)[3])niggli_lattice, GRIDSYS_NIGGLI_TOLERANCE)) { return 0; } @@ -156,7 +156,7 @@ int64_t gridsys_get_bz_grid_addresses( bzgrid->bzg2grg = bzg2grg; bzgrid->type = bz_grid_type; lagmat_multiply_matrix_l3(bzgrid->Q, inv_Mpr_int, Q); - lagmat_copy_matrix_d3(bzgrid->reclat, (double(*)[3])niggli_lattice); + lagmat_copy_matrix_d3(bzgrid->reclat, (double (*)[3])niggli_lattice); for (i = 0; i < 3; i++) { bzgrid->D_diag[i] = D_diag[i]; bzgrid->PS[i] = PS[i]; @@ -220,7 +220,7 @@ int64_t gridsys_get_bz_triplets_at_q( const int64_t (*bz_grid_addresses)[3], const int64_t *bz_map, const int64_t *map_triplets, const int64_t num_map_triplets, const int64_t D_diag[3], const int64_t Q[3][3], - const int64_t bz_grid_type) { + const double reciprocal_lattice[3][3], const int64_t bz_grid_type) { RecgridConstBZGrid *bzgrid; int64_t i, j, num_ir; @@ -237,6 +237,7 @@ int64_t gridsys_get_bz_triplets_at_q( bzgrid->D_diag[i] = D_diag[i]; bzgrid->PS[i] = 0; for (j = 0; j < 3; j++) { + bzgrid->reclat[i][j] = reciprocal_lattice[i][j]; bzgrid->Q[i][j] = Q[i][j]; } } diff --git a/c/gridsys.h b/c/gridsys.h index 6c427688..3ba91cdf 100644 --- a/c/gridsys.h +++ b/c/gridsys.h @@ -272,7 +272,8 @@ int64_t gridsys_get_bz_triplets_at_q( int64_t (*ir_triplets)[3], const int64_t bz_grid_index, const int64_t (*bz_grid_addresses)[3], const int64_t *bz_map, const int64_t *map_triplets, const int64_t num_map_triplets, - const int64_t D_diag[3], const int64_t Q[3][3], const int64_t bz_grid_type); + const int64_t D_diag[3], const int64_t Q[3][3], + const double reciprocal_lattice[3][3], const int64_t bz_grid_type); /** * @brief Return integration weight of linear tetrahedron method diff --git a/c/phono3py.c b/c/phono3py.c index 47029370..8ea3d1b5 100644 --- a/c/phono3py.c +++ b/c/phono3py.c @@ -411,7 +411,7 @@ int64_t ph3py_get_BZ_triplets_at_q( const int64_t (*bz_grid_addresses)[3], const int64_t *bz_map, const int64_t *map_triplets, const int64_t num_map_triplets, const int64_t D_diag[3], const int64_t Q[3][3], - const int64_t bz_grid_type) { + const double reciprocal_lattice[3][3], const int64_t bz_grid_type) { RecgridConstBZGrid *bzgrid; int64_t i, j, num_ir; @@ -429,6 +429,7 @@ int64_t ph3py_get_BZ_triplets_at_q( bzgrid->PS[i] = 0; for (j = 0; j < 3; j++) { bzgrid->Q[i][j] = Q[i][j]; + bzgrid->reclat[i][j] = reciprocal_lattice[i][j]; } } bzgrid->size = num_map_triplets; diff --git a/c/phono3py.h b/c/phono3py.h index e5751d79..9cd02fcd 100644 --- a/c/phono3py.h +++ b/c/phono3py.h @@ -173,7 +173,8 @@ int64_t ph3py_get_BZ_triplets_at_q( int64_t (*triplets)[3], const int64_t grid_point, const int64_t (*bz_grid_addresses)[3], const int64_t *bz_map, const int64_t *map_triplets, const int64_t num_map_triplets, - const int64_t D_diag[3], const int64_t Q[3][3], const int64_t bz_grid_type); + const int64_t D_diag[3], const int64_t Q[3][3], + const double reciprocal_lattice[3][3], const int64_t bz_grid_type); int64_t ph3py_get_integration_weight( double *iw, char *iw_zero, const double *frequency_points, const int64_t num_band0, const int64_t relative_grid_address[24][4][3], diff --git a/c/recgrid.c b/c/recgrid.c index 1fe6f554..0f0d65da 100644 --- a/c/recgrid.c +++ b/c/recgrid.c @@ -193,20 +193,11 @@ double recgrid_get_tolerance_for_BZ_reduction(const RecgridBZGrid *bzgrid) { int64_t i, j; double tolerance; double length[3]; - double reclatQ[3][3]; - - for (i = 0; i < 3; i++) { - for (j = 0; j < 3; j++) { - reclatQ[i][j] = bzgrid->reclat[i][0] * bzgrid->Q[0][j] + - bzgrid->reclat[i][1] * bzgrid->Q[1][j] + - bzgrid->reclat[i][2] * bzgrid->Q[2][j]; - } - } for (i = 0; i < 3; i++) { length[i] = 0; for (j = 0; j < 3; j++) { - length[i] += reclatQ[j][i] * reclatQ[j][i]; + length[i] += bzgrid->reclat[j][i] * bzgrid->reclat[j][i]; } length[i] /= bzgrid->D_diag[i] * bzgrid->D_diag[i]; } @@ -233,8 +224,8 @@ RecgridMats *recgrid_alloc_RotMats(const int64_t size) { rotmats->size = size; if (size > 0) { - if ((rotmats->mat = (int64_t(*)[3][3])malloc(sizeof(int64_t[3][3]) * - size)) == NULL) { + if ((rotmats->mat = (int64_t (*)[3][3])malloc(sizeof(int64_t[3][3]) * + size)) == NULL) { warning_print("Memory could not be allocated "); warning_print("(RecgridMats, line %d, %s).\n", __LINE__, __FILE__); free(rotmats); diff --git a/c/triplet_grid.c b/c/triplet_grid.c index c562e8d6..63213014 100644 --- a/c/triplet_grid.c +++ b/c/triplet_grid.c @@ -262,7 +262,7 @@ static void get_BZ_triplets_at_q_type1(int64_t (*triplets)[3], int64_t bzgp[3], G[3]; int64_t bz_adrs0[3], bz_adrs1[3], bz_adrs2[3]; const int64_t *gp_map; - const int64_t(*bz_adrs)[3]; + const int64_t (*bz_adrs)[3]; double d2, min_d2, tolerance; double LQD_inv[3][3]; @@ -348,13 +348,14 @@ static void get_BZ_triplets_at_q_type2(int64_t (*triplets)[3], int64_t bzgp[3], G[3]; int64_t bz_adrs0[3], bz_adrs1[3], bz_adrs2[3]; const int64_t *gp_map; - const int64_t(*bz_adrs)[3]; + const int64_t (*bz_adrs)[3]; double d2, min_d2, tolerance; double LQD_inv[3][3]; gp_map = bzgrid->gp_map; bz_adrs = bzgrid->addresses; get_LQD_inv(LQD_inv, bzgrid); + /* This tolerance is used to be consistent to BZ reduction in bzgrid. */ tolerance = recgrid_get_tolerance_for_BZ_reduction((RecgridBZGrid *)bzgrid); @@ -423,10 +424,15 @@ static void get_LQD_inv(double LQD_inv[3][3], int64_t i, j, k; /* LQD^-1 */ + for (i = 0; i < 3; i++) { + for (j = 0; j < 3; j++) { + LQD_inv[i][j] = 0; + } + } for (i = 0; i < 3; i++) { for (j = 0; j < 3; j++) { for (k = 0; k < 3; k++) { - LQD_inv[i][k] = + LQD_inv[i][k] += bzgrid->reclat[i][j] * bzgrid->Q[j][k] / bzgrid->D_diag[k]; } } diff --git a/phono3py/phonon/grid.py b/phono3py/phonon/grid.py index f448d1df..6fcab8b7 100644 --- a/phono3py/phonon/grid.py +++ b/phono3py/phonon/grid.py @@ -39,6 +39,8 @@ from __future__ import annotations import warnings from collections.abc import Sequence +from numpy.typing import ArrayLike, NDArray + try: from spglib import SpglibDataset # type: ignore except ImportError: @@ -356,6 +358,16 @@ class BZGrid: """ return self._microzone_lattice + @property + def reciprocal_lattice(self) -> np.ndarray: + """Reciprocal basis vectors of primitive cell. + + Reciprocal basis vectors of primitive cell in column vectors. + shape=(3, 3), dtype='double', order='C'. + + """ + return self._reciprocal_lattice + @property def store_dense_gp_map(self): """Return gp_map type. @@ -1172,22 +1184,15 @@ def _relocate_BZ_grid_address( else: bz_map = np.zeros(np.prod(D_diag) * 9 + 1, dtype="int64") - # Mpr^-1 = Lr^-1 Lp - reclat_T = np.array(np.transpose(reciprocal_lattice), dtype="double", order="C") - reduced_basis = get_reduced_bases(reclat_T) - assert reduced_basis is not None, "Reduced basis is not found." - tmat_inv = np.dot(np.linalg.inv(reduced_basis.T), reclat_T.T) - tmat_inv_int = np.rint(tmat_inv).astype("int64") - assert (np.abs(tmat_inv - tmat_inv_int) < 1e-5).all() - + reduced_basis, tmat_inv_int = get_reduced_bases_and_tmat_inv(reciprocal_lattice) num_gp = recgrid.bz_grid_addresses( bz_grid_addresses, bz_map, bzg2grg, np.array(D_diag, dtype="int64"), - np.array(np.dot(tmat_inv_int, Q), dtype="int64", order="C"), + np.array(tmat_inv_int @ Q, dtype="int64", order="C"), _PS, - np.array(reduced_basis.T, dtype="double", order="C"), + reduced_basis, store_dense_gp_map * 1 + 1, ) @@ -1196,6 +1201,39 @@ def _relocate_BZ_grid_address( return bz_grid_addresses, bz_map, bzg2grg +def get_reduced_bases_and_tmat_inv( + reciprocal_lattice: ArrayLike, +) -> tuple[NDArray, NDArray]: + """Return reduced bases and inverse transformation matrix. + + Parameters + ---------- + reciprocal_lattice : ArrayLike + Reciprocal lattice vectors in column vectors. + shape=(3, 3), dtype='double' + + Returns + ------- + reduced_basis : ndarray + Reduced basis vectors in column vectors. + shape=(3, 3), dtype='double', order='C' + tmat_inv_int : ndarray + Inverse transformation matrix in integer. + This is used to transform reciprocal lattice vectors to + conventional lattice vectors. + shape=(3, 3), dtype='int64' + + """ + # Mpr^-1 = Lr^-1 Lp + reclat_T = np.array(np.transpose(reciprocal_lattice), dtype="double", order="C") + reduced_basis = get_reduced_bases(reclat_T) + assert reduced_basis is not None, "Reduced basis is not found." + tmat_inv = np.linalg.inv(reduced_basis.T) @ reclat_T.T + tmat_inv_int = np.rint(tmat_inv).astype("int64") + assert (np.abs(tmat_inv - tmat_inv_int) < 1e-5).all() + return np.array(reduced_basis.T, dtype="double", order="C"), tmat_inv_int + + def _get_ir_grid_map( D_diag: Union[np.ndarray, Sequence], grg_rotations: Union[np.ndarray, Sequence], diff --git a/phono3py/phonon3/triplets.py b/phono3py/phonon3/triplets.py index b1c36216..30d731a8 100644 --- a/phono3py/phonon3/triplets.py +++ b/phono3py/phonon3/triplets.py @@ -42,7 +42,11 @@ from phonopy.structure.tetrahedron_method import TetrahedronMethod from phono3py.other.tetrahedron_method import get_tetrahedra_relative_grid_address from phono3py.phonon.func import gaussian -from phono3py.phonon.grid import BZGrid, get_grid_point_from_address_py +from phono3py.phonon.grid import ( + BZGrid, + get_grid_point_from_address_py, + get_reduced_bases_and_tmat_inv, +) if TYPE_CHECKING: from phono3py.phonon3.interaction import Interaction @@ -327,14 +331,18 @@ def _get_BZ_triplets_at_q(bz_grid_index, bz_grid: BZGrid, map_triplets): weights[g] += 1 ir_weights = np.extract(weights > 0, weights) triplets = -np.ones((len(ir_weights), 3), dtype="int64", order="C") + reduced_basis, tmat_inv_int = get_reduced_bases_and_tmat_inv( + bz_grid.reciprocal_lattice + ) num_ir_ret = phono3c.BZ_triplets_at_q( triplets, bz_grid_index, bz_grid.addresses, bz_grid.gp_map, map_triplets, - np.array(bz_grid.D_diag, dtype="int64"), - bz_grid.Q, + bz_grid.D_diag, + np.array(tmat_inv_int @ bz_grid.Q, dtype="int64", order="C"), + reduced_basis, bz_grid.store_dense_gp_map * 1 + 1, ) diff --git a/test/phonon3/test_triplets.py b/test/phonon3/test_triplets.py index ef5ddf5d..f3bdfedd 100644 --- a/test/phonon3/test_triplets.py +++ b/test/phonon3/test_triplets.py @@ -911,7 +911,7 @@ def test_get_BZ_triplets_at_q(aln_cell: PhonopyAtoms, params): [1, 36, 60], [1, 38, 59], [1, 41, 56], - ], + ], # 0 [ [1, 0, 4], [1, 1, 3], @@ -937,7 +937,7 @@ def test_get_BZ_triplets_at_q(aln_cell: PhonopyAtoms, params): [1, 36, 60], [1, 38, 59], [1, 41, 56], - ], + ], # 1 [ [1, 0, 4], [1, 1, 3], @@ -969,7 +969,7 @@ def test_get_BZ_triplets_at_q(aln_cell: PhonopyAtoms, params): [1, 39, 57], [1, 41, 56], [1, 42, 55], - ], + ], # 2 [ [1, 0, 4], [1, 1, 3], @@ -1016,7 +1016,7 @@ def test_get_BZ_triplets_at_q(aln_cell: PhonopyAtoms, params): [1, 70, 26], [1, 72, 25], [1, 73, 24], - ], + ], # 3 [ [1, 0, 4], [1, 1, 3], @@ -1150,11 +1150,11 @@ def test_get_BZ_triplets_at_q(aln_cell: PhonopyAtoms, params): [8, 2, 87], [8, 3, 86], [9, 4, 92], - [8, 5, 84], - [8, 6, 82], - [8, 8, 81], - [8, 10, 80], - [8, 11, 85], + [9, 5, 84], + [8, 7, 82], + [8, 9, 81], + [9, 10, 80], + [9, 11, 85], [8, 12, 78], [8, 13, 76], [8, 14, 75], @@ -1164,23 +1164,23 @@ def test_get_BZ_triplets_at_q(aln_cell: PhonopyAtoms, params): [9, 22, 67], [8, 24, 65], [8, 27, 62], - [8, 29, 66], + [9, 30, 66], [8, 31, 58], [8, 32, 57], [8, 40, 50], [8, 48, 48], - ], + ], # 0 [ [8, 0, 89], [8, 1, 88], [8, 2, 87], [8, 3, 86], [9, 4, 92], - [8, 5, 84], - [8, 6, 82], - [8, 8, 81], - [8, 10, 80], - [8, 11, 85], + [9, 5, 84], + [8, 7, 82], + [8, 9, 81], + [9, 10, 80], + [9, 11, 85], [8, 12, 78], [8, 13, 76], [8, 14, 75], @@ -1190,23 +1190,23 @@ def test_get_BZ_triplets_at_q(aln_cell: PhonopyAtoms, params): [9, 22, 67], [8, 24, 65], [8, 27, 62], - [8, 29, 66], + [9, 30, 66], [8, 31, 58], [8, 32, 57], [8, 40, 50], [8, 48, 48], - ], + ], # 1 [ [8, 0, 89], [8, 1, 88], [8, 2, 87], [8, 3, 86], [9, 4, 92], - [8, 5, 84], - [8, 6, 82], - [8, 8, 81], - [8, 10, 80], - [8, 11, 85], + [9, 5, 84], + [8, 7, 82], + [8, 9, 81], + [9, 10, 80], + [9, 11, 85], [8, 12, 78], [8, 13, 76], [8, 14, 75], @@ -1219,7 +1219,7 @@ def test_get_BZ_triplets_at_q(aln_cell: PhonopyAtoms, params): [8, 24, 65], [8, 25, 64], [8, 27, 62], - [8, 29, 66], + [9, 30, 66], [8, 31, 58], [8, 32, 57], [8, 33, 56], @@ -1227,7 +1227,7 @@ def test_get_BZ_triplets_at_q(aln_cell: PhonopyAtoms, params): [8, 40, 50], [9, 41, 49], [9, 42, 54], - [8, 43, 46], + [8, 43, 47], [8, 44, 45], [8, 48, 48], [8, 50, 40], @@ -1240,20 +1240,20 @@ def test_get_BZ_triplets_at_q(aln_cell: PhonopyAtoms, params): [8, 71, 19], [9, 72, 18], [8, 79, 17], - [8, 81, 8], + [8, 81, 9], [8, 89, 0], - ], + ], # 2 [ [8, 0, 89], [8, 1, 88], [8, 2, 87], [8, 3, 86], [9, 4, 92], - [8, 5, 84], - [8, 6, 82], - [8, 8, 81], - [8, 10, 80], - [8, 11, 85], + [9, 5, 84], + [8, 7, 82], + [8, 9, 81], + [9, 10, 80], + [9, 11, 85], [8, 12, 78], [8, 13, 76], [8, 14, 75], @@ -1266,7 +1266,7 @@ def test_get_BZ_triplets_at_q(aln_cell: PhonopyAtoms, params): [8, 24, 65], [8, 25, 64], [8, 27, 62], - [8, 29, 66], + [9, 30, 66], [8, 31, 58], [8, 32, 57], [8, 33, 56], @@ -1274,7 +1274,7 @@ def test_get_BZ_triplets_at_q(aln_cell: PhonopyAtoms, params): [8, 40, 50], [9, 41, 49], [9, 42, 54], - [8, 43, 46], + [8, 43, 47], [8, 44, 45], [8, 48, 48], [8, 50, 40], @@ -1287,9 +1287,9 @@ def test_get_BZ_triplets_at_q(aln_cell: PhonopyAtoms, params): [8, 71, 19], [9, 72, 18], [8, 79, 17], - [8, 81, 8], + [8, 81, 9], [8, 89, 0], - ], + ], # 3 [ [7, 0, 28], [8, 1, 27], @@ -1794,22 +1794,23 @@ def test_get_BZ_triplets_at_q(aln_cell: PhonopyAtoms, params): ir_weights_with_zeros = np.zeros(45, dtype=int) ir_weights_with_zeros[: len(ir_weights)] = ir_weights - np.testing.assert_equal(ref_triplets[i][params[3]], ir_triplets) - np.testing.assert_equal(ref_ir_weights[i][params[3]], ir_weights) - if i == 1 and params[0]: - # print("{") - # for j, tp in enumerate(triplets_with_zeros): - # print("%d, %d, %d, " % tuple(tp), end="") - # if (j + 1) % 5 == 0: - # print("&") - print(", ".join(["%d" % x for x in ir_weights_with_zeros])) - # print("},") - # print(len(ir_triplets)) - # print("[") # print("[") # for tp in ir_triplets: # print("[%d, %d, %d]," % tuple(tp)) # print("],") + + np.testing.assert_equal(ref_triplets[i][params[3]], ir_triplets) + np.testing.assert_equal(ref_ir_weights[i][params[3]], ir_weights) + # if i == 1 and params[0]: + # print("{") + # for j, tp in enumerate(triplets_with_zeros): + # print("%d, %d, %d, " % tuple(tp), end="") + # if (j + 1) % 5 == 0: + # print("&") + # print(", ".join(["%d" % x for x in ir_weights_with_zeros])) + # print("},") + # print(len(ir_triplets)) + # print("[") # print("]") # print("[") # print(",".join(["%d" % x for x in ir_weights]))