scnc
/
phono3py
mirror of https://github.com/phonopy/phono3py.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Fix distance calculation in BZ

This commit is contained in:
Atsushi Togo 2025-06-10 18:32:13 +09:00
parent 2d30f9baf0
commit 71bb34aa4a
10 changed files with 140 additions and 92 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
c/_phono3py.cpp
View File

@ -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;
}

3
c/gridsys.c
View File

@ -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];
}
}

3
c/gridsys.h
View File

@ -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

3
c/phono3py.c
View File

@ -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;

3
c/phono3py.h
View File

@ -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],

11
c/recgrid.c
View File

@ -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];
}

8
c/triplet_grid.c
View File

@ -355,6 +355,7 @@ static void get_BZ_triplets_at_q_type2(int64_t (*triplets)[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];
}
}

58
phono3py/phonon/grid.py
View File

@ -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],

14
phono3py/phonon3/triplets.py
View File

@ -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,
)

85
test/phonon3/test_triplets.py
View File

@ -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
# 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]:
# 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(", ".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("],")
# print("]")
# print("[")
# print(",".join(["%d" % x for x in ir_weights]))