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

Merge branch 'develop'

This commit is contained in:
Atsushi Togo 2025-02-01 16:15:06 +09:00
parent 30ecb1b697 cacc4428d5
commit 64432ef141
7 changed files with 196 additions and 125 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

228
c/collision_matrix.c
View File

@ -51,6 +51,15 @@ static void get_collision_matrix(
const double *rotations_cartesian, const double *g,
const double temperature, const double unit_conversion_factor,
const double cutoff_frequency);
static void get_collision_matrix_at_gp(
double *collision_matrix, const double *fc3_normal_squared,
const int64_t num_band0, const int64_t num_band, const double *frequencies,
const int64_t (*triplets)[3], const int64_t *triplets_map,
const int64_t *map_q, const int64_t *rot_grid_points,
const int64_t num_ir_gp, const int64_t num_rot,
const double *rotations_cartesian, const double *g,
const double temperature, const double unit_conversion_factor,
const double cutoff_frequency, const int64_t *gp2tp_map, const int64_t i);
static void get_reducible_collision_matrix(
double *collision_matrix, const double *fc3_normal_squared,
const int64_t num_band0, const int64_t num_band, const double *frequencies,
@ -58,6 +67,13 @@ static void get_reducible_collision_matrix(
const int64_t num_gp, const int64_t *map_q, const double *g,
const double temperature, const double unit_conversion_factor,
const double cutoff_frequency);
static void get_reducible_collision_matrix_at_gp(
double *collision_matrix, const double *fc3_normal_squared,
const int64_t num_band0, const int64_t num_band, const double *frequencies,
const int64_t (*triplets)[3], const int64_t *triplets_map,
const int64_t num_gp, const int64_t *map_q, const double *g,
const double temperature, const double unit_conversion_factor,
const double cutoff_frequency, const int64_t *gp2tp_map, const int64_t i);
static int64_t get_inv_sinh(double *inv_sinh, const int64_t gp,
const double temperature, const double *frequencies,
const int64_t triplet[3],
@ -117,69 +133,87 @@ static void get_collision_matrix(
const double *rotations_cartesian, const double *g,
const double temperature, const double unit_conversion_factor,
const double cutoff_frequency) {
int64_t i, j, k, l, m, n, ti, r_gp, swapped;
int64_t i;
int64_t *gp2tp_map;
double collision;
double *inv_sinh;
gp2tp_map = create_gp2tp_map(triplets_map, num_gp);
#ifdef _OPENMP
#pragma omp parallel for private(j, k, l, m, n, ti, r_gp, collision, inv_sinh)
#pragma omp parallel for
#endif
for (i = 0; i < num_ir_gp; i++) {
inv_sinh = (double *)malloc(sizeof(double) * num_band);
for (j = 0; j < num_rot; j++) {
r_gp = rot_grid_points[i * num_rot + j];
ti = gp2tp_map[triplets_map[r_gp]];
swapped = get_inv_sinh(inv_sinh, r_gp, temperature, frequencies,
triplets[ti], triplets_map, map_q, num_band,
cutoff_frequency);
for (k = 0; k < num_band0; k++) {
for (l = 0; l < num_band; l++) {
collision = 0;
for (m = 0; m < num_band; m++) {
if (swapped) {
collision +=
fc3_normal_squared[ti * num_band0 * num_band *
num_band +
k * num_band * num_band +
m * num_band + l] *
g[ti * num_band0 * num_band * num_band +
k * num_band * num_band + m * num_band + l] *
inv_sinh[m] * unit_conversion_factor;
} else {
collision +=
fc3_normal_squared[ti * num_band0 * num_band *
num_band +
k * num_band * num_band +
l * num_band + m] *
g[ti * num_band0 * num_band * num_band +
k * num_band * num_band + l * num_band + m] *
inv_sinh[m] * unit_conversion_factor;
}
}
for (m = 0; m < 3; m++) {
for (n = 0; n < 3; n++) {
collision_matrix[k * 3 * num_ir_gp * num_band * 3 +
m * num_ir_gp * num_band * 3 +
i * num_band * 3 + l * 3 + n] +=
collision *
rotations_cartesian[j * 9 + m * 3 + n];
}
}
}
}
}
free(inv_sinh);
inv_sinh = NULL;
get_collision_matrix_at_gp(
collision_matrix, fc3_normal_squared, num_band0, num_band,
frequencies, triplets, triplets_map, map_q, rot_grid_points,
num_ir_gp, num_rot, rotations_cartesian, g, temperature,
unit_conversion_factor, cutoff_frequency, gp2tp_map, i);
}
free(gp2tp_map);
gp2tp_map = NULL;
}
static void get_collision_matrix_at_gp(
double *collision_matrix, const double *fc3_normal_squared,
const int64_t num_band0, const int64_t num_band, const double *frequencies,
const int64_t (*triplets)[3], const int64_t *triplets_map,
const int64_t *map_q, const int64_t *rot_grid_points,
const int64_t num_ir_gp, const int64_t num_rot,
const double *rotations_cartesian, const double *g,
const double temperature, const double unit_conversion_factor,
const double cutoff_frequency, const int64_t *gp2tp_map, const int64_t i) {
int64_t j, k, l, m, n, ti, r_gp, swapped;
double collision;
double *inv_sinh;
inv_sinh = (double *)malloc(sizeof(double) * num_band);
for (j = 0; j < num_rot; j++) {
r_gp = rot_grid_points[i * num_rot + j];
ti = gp2tp_map[triplets_map[r_gp]];
swapped =
get_inv_sinh(inv_sinh, r_gp, temperature, frequencies, triplets[ti],
triplets_map, map_q, num_band, cutoff_frequency);
for (k = 0; k < num_band0; k++) {
for (l = 0; l < num_band; l++) {
collision = 0;
for (m = 0; m < num_band; m++) {
if (swapped) {
collision +=
fc3_normal_squared[ti * num_band0 * num_band *
num_band +
k * num_band * num_band +
m * num_band + l] *
g[ti * num_band0 * num_band * num_band +
k * num_band * num_band + m * num_band + l] *
inv_sinh[m] * unit_conversion_factor;
} else {
collision +=
fc3_normal_squared[ti * num_band0 * num_band *
num_band +
k * num_band * num_band +
l * num_band + m] *
g[ti * num_band0 * num_band * num_band +
k * num_band * num_band + l * num_band + m] *
inv_sinh[m] * unit_conversion_factor;
}
}
for (m = 0; m < 3; m++) {
for (n = 0; n < 3; n++) {
collision_matrix[k * 3 * num_ir_gp * num_band * 3 +
m * num_ir_gp * num_band * 3 +
i * num_band * 3 + l * 3 + n] +=
collision * rotations_cartesian[j * 9 + m * 3 + n];
}
}
}
}
}
free(inv_sinh);
inv_sinh = NULL;
}
static void get_reducible_collision_matrix(
double *collision_matrix, const double *fc3_normal_squared,
const int64_t num_band0, const int64_t num_band, const double *frequencies,
@ -187,60 +221,72 @@ static void get_reducible_collision_matrix(
const int64_t num_gp, const int64_t *map_q, const double *g,
const double temperature, const double unit_conversion_factor,
const double cutoff_frequency) {
int64_t i, j, k, l, ti, swapped;
int64_t i;
int64_t *gp2tp_map;
double collision;
double *inv_sinh;
gp2tp_map = create_gp2tp_map(triplets_map, num_gp);
#ifdef _OPENMP
#pragma omp parallel for private(j, k, l, ti, collision, inv_sinh)
#pragma omp parallel for
#endif
for (i = 0; i < num_gp; i++) {
inv_sinh = (double *)malloc(sizeof(double) * num_band);
ti = gp2tp_map[triplets_map[i]];
swapped =
get_inv_sinh(inv_sinh, i, temperature, frequencies, triplets[ti],
triplets_map, map_q, num_band, cutoff_frequency);
for (j = 0; j < num_band0; j++) {
for (k = 0; k < num_band; k++) {
collision = 0;
for (l = 0; l < num_band; l++) {
if (swapped) {
collision +=
fc3_normal_squared[ti * num_band0 * num_band *
num_band +
j * num_band * num_band +
l * num_band + k] *
g[ti * num_band0 * num_band * num_band +
j * num_band * num_band + l * num_band + k] *
inv_sinh[l] * unit_conversion_factor;
} else {
collision +=
fc3_normal_squared[ti * num_band0 * num_band *
num_band +
j * num_band * num_band +
k * num_band + l] *
g[ti * num_band0 * num_band * num_band +
j * num_band * num_band + k * num_band + l] *
inv_sinh[l] * unit_conversion_factor;
}
}
collision_matrix[j * num_gp * num_band + i * num_band + k] +=
collision;
}
}
free(inv_sinh);
inv_sinh = NULL;
get_reducible_collision_matrix_at_gp(
collision_matrix, fc3_normal_squared, num_band0, num_band,
frequencies, triplets, triplets_map, num_gp, map_q, g, temperature,
unit_conversion_factor, cutoff_frequency, gp2tp_map, i);
}
free(gp2tp_map);
gp2tp_map = NULL;
}
static void get_reducible_collision_matrix_at_gp(
double *collision_matrix, const double *fc3_normal_squared,
const int64_t num_band0, const int64_t num_band, const double *frequencies,
const int64_t (*triplets)[3], const int64_t *triplets_map,
const int64_t num_gp, const int64_t *map_q, const double *g,
const double temperature, const double unit_conversion_factor,
const double cutoff_frequency, const int64_t *gp2tp_map, const int64_t i) {
int64_t j, k, l, ti, swapped;
double collision;
double *inv_sinh;
inv_sinh = (double *)malloc(sizeof(double) * num_band);
ti = gp2tp_map[triplets_map[i]];
swapped = get_inv_sinh(inv_sinh, i, temperature, frequencies, triplets[ti],
triplets_map, map_q, num_band, cutoff_frequency);
for (j = 0; j < num_band0; j++) {
for (k = 0; k < num_band; k++) {
collision = 0;
for (l = 0; l < num_band; l++) {
if (swapped) {
collision += fc3_normal_squared[ti * num_band0 * num_band *
num_band +
j * num_band * num_band +
l * num_band + k] *
g[ti * num_band0 * num_band * num_band +
j * num_band * num_band + l * num_band + k] *
inv_sinh[l] * unit_conversion_factor;
} else {
collision += fc3_normal_squared[ti * num_band0 * num_band *
num_band +
j * num_band * num_band +
k * num_band + l] *
g[ti * num_band0 * num_band * num_band +
j * num_band * num_band + k * num_band + l] *
inv_sinh[l] * unit_conversion_factor;
}
}
collision_matrix[j * num_gp * num_band + i * num_band + k] +=
collision;
}
}
free(inv_sinh);
inv_sinh = NULL;
}
static int64_t get_inv_sinh(double *inv_sinh, const int64_t gp,
const double temperature, const double *frequencies,
const int64_t triplet[3],

14
doc/changelog.md
View File

@ -2,27 +2,31 @@
# Change Log
## Jan-28-2024: Version 3.12.1
## Feb-1-2025: Version 3.12.2
- Fix an openmp related bug in computing collision matrix in C
## Jan-28-2025: Version 3.12.1
- Update `pyproject.toml`.
## Jan-28-2024: Version 3.12.0
## Jan-28-2025: Version 3.12.0
- `dtype="long"` was replaced by `dtype="int64"` aiming making Windows build. In
C, `long` was replaced by `int64_t`.
- Fix `phono3py-kaccum`.
## Jan-18-2024: Version 3.11.2
## Jan-18-2025: Version 3.11.2
- Maintenance release.
## Jan-12-2024: Version 3.11.1
## Jan-12-2025: Version 3.11.1
- `-i`, `-o`, `--io` options have been deprecated.
- The `--amplitude` option can now be used to specify the displacement distance
for `phono3py-load --pypolymlp`.
## Jan-2-2024: Version 3.11.0
## Jan-2-2025: Version 3.11.0
- Release to follow the change of phonopy
- Add `--rd auto` and `--rd-fc2 auto` options

2
doc/conf.py
View File

@ -60,7 +60,7 @@ copyright = "2015, Atsushi Togo"
# The short X.Y version.
version = "3.12"
# The full version, including alpha/beta/rc tags.
release = "3.12.1"
release = "3.12.2"
# The language for content autogenerated by Sphinx. Refer to documentation
# for a list of supported languages.

61
phono3py/conductivity/direct_solution.py
View File

@ -38,7 +38,7 @@ import sys
import time
import warnings
from abc import abstractmethod
from typing import Type, Union
from typing import Optional, Type, Union
import numpy as np
from phonopy.phonon.degeneracy import degenerate_sets
@ -133,7 +133,9 @@ class ConductivityLBTEBase(ConductivityBase):
if self._is_reducible_collision_matrix:
self._collision = CollisionMatrix(
self._pp, is_reducible_collision_matrix=True, log_level=self._log_level
self._pp,
is_reducible_collision_matrix=True,
log_level=self._log_level,
)
else:
self._rot_grid_points = self._get_rot_grid_points()
@ -817,6 +819,8 @@ class ConductivityLBTEBase(ConductivityBase):
def _get_Y(self, i_sigma, i_temp, weights, X):
r"""Calculate Y = (\Omega^-1, X)."""
solver = select_colmat_solver(self._pinv_solver)
if self._pinv_solver == 6:
solver = 6
num_band = len(self._pp.primitive) * 3
if self._is_reducible_collision_matrix:
@ -832,7 +836,7 @@ class ConductivityLBTEBase(ConductivityBase):
start = time.time()
if self._log_level:
if self._log_level and solver != 7:
if self._pinv_method == 0:
eig_str = "abs(eig)"
else:
@ -856,12 +860,11 @@ class ConductivityLBTEBase(ConductivityBase):
Y = np.dot(v, X1)
else:
Y = np.dot(v, e * np.dot(v.T, X.ravel())).reshape(-1, 3)
else: # solver=6 This is slower as far as tested.
elif solver == 6: # solver=6 This is slower as far as tested.
import phono3py._phono3py as phono3c
if self._log_level:
print(" (built-in) ", end="")
sys.stdout.flush()
print(" (built-in-pinv) ", end="", flush=True)
w = self._collision_eigenvalues[i_sigma, i_temp]
phono3c.pinv_from_eigensolution(
@ -876,11 +879,18 @@ class ConductivityLBTEBase(ConductivityBase):
Y = np.dot(v, X)
else:
Y = np.dot(v, X.ravel()).reshape(-1, 3)
elif solver == 7:
if self._is_reducible_collision_matrix:
Y = np.dot(v, X)
else:
Y = np.dot(v, X.ravel()).reshape(-1, 3)
else:
raise ValueError(f"Unknown collision matrix solver {solver}")
self._set_f_vectors(Y, num_grid_points, weights)
if self._log_level:
print("[%.3fs]" % (time.time() - start))
if self._log_level and solver != 7:
print("[%.3fs]" % (time.time() - start), flush=True)
sys.stdout.flush()
return Y
@ -1359,7 +1369,8 @@ class ConductivityLBTE(ConductivityMixIn, ConductivityLBTEBase):
pinv_solver=self._pinv_solver,
log_level=self._log_level,
)
self._collision_eigenvalues[j, k] = w
if w is not None:
self._collision_eigenvalues[j, k] = w
self._set_kappa(j, k, weights)
@ -1829,7 +1840,7 @@ def get_thermal_conductivity_LBTE(
def diagonalize_collision_matrix(
collision_matrices, i_sigma=None, i_temp=None, pinv_solver=0, log_level=0
):
) -> Optional[np.ndarray]:
"""Diagonalize collision matrices.
Note
@ -1862,6 +1873,7 @@ def diagonalize_collision_matrix(
w : ndarray, optional
Eigenvalues.
shape=(size_of_collision_matrix,), dtype='double'
When pinv_solve==7, None is returned.
"""
start = time.time()
@ -1890,8 +1902,7 @@ def diagonalize_collision_matrix(
if solver in [1, 2]:
if log_level:
routine = ["dsyev", "dsyevd"][solver - 1]
sys.stdout.write("Diagonalizing by lapacke %s ... " % routine)
sys.stdout.flush()
print("Diagonalizing by lapacke %s ... " % routine, end="", flush=True)
import phono3py._phono3py as phono3c
w = np.zeros(size, dtype="double")
@ -1908,31 +1919,39 @@ def diagonalize_collision_matrix(
) # only diagonalization
elif solver == 3: # np.linalg.eigh depends on dsyevd.
if log_level:
sys.stdout.write("Diagonalize by np.linalg.eigh ")
sys.stdout.flush()
print("Diagonalize by np.linalg.eigh ", end="", flush=True)
col_mat = collision_matrices[i_sigma, i_temp].reshape(size, size)
w, col_mat[:] = np.linalg.eigh(col_mat)
elif solver == 4: # fully scipy dsyev
if log_level:
sys.stdout.write("Diagonalize by scipy.linalg.lapack.dsyev ")
sys.stdout.flush()
print("Diagonalize by scipy.linalg.lapack.dsyev ", end="", flush=True)
import scipy.linalg
col_mat = collision_matrices[i_sigma, i_temp].reshape(size, size)
w, _, info = scipy.linalg.lapack.dsyev(col_mat.T, overwrite_a=1)
elif solver == 5: # fully scipy dsyevd
if log_level:
sys.stdout.write("Diagnalize by scipy.linalg.lapack.dsyevd ")
sys.stdout.flush()
print("Diagnalize by scipy.linalg.lapack.dsyevd ", end="", flush=True)
import scipy.linalg
col_mat = collision_matrices[i_sigma, i_temp].reshape(size, size)
w, _, info = scipy.linalg.lapack.dsyevd(col_mat.T, overwrite_a=1)
elif solver == 7:
if log_level:
print(
"Pseudo inversion using np.linalg.pinv(a, hermitian=False) ",
end="",
flush=True,
)
col_mat = collision_matrices[i_sigma, i_temp].reshape(size, size)
# hermitian=True calls eigh, which is not what we want.
col_mat[:, :] = np.linalg.pinv(col_mat, hermitian=False)
w = None
if log_level:
print(f"sum={w.sum():<.1e} d={trace - w.sum():<.1e} ", end="")
print("[%.3fs]" % (time.time() - start))
sys.stdout.flush()
if w is not None:
print(f"sum={w.sum():<.1e} d={trace - w.sum():<.1e} ", end="")
print("[%.3fs]" % (time.time() - start), flush=True)
return w

6
phono3py/conductivity/utils.py
View File

@ -693,11 +693,13 @@ def select_colmat_solver(pinv_solver):
"phono3py is not compiled with LAPACKE."
)
solver_numbers = (1, 2, 3, 4, 5, 6)
solver_numbers = (1, 2, 3, 4, 5, 6, 7)
solver = pinv_solver
if solver == 6: # 6 must return 3 for not transposing unitary matrix.
solver = 3
if solver == 0: # default solver
if default_solver in (4, 5, 6):
if default_solver in (3, 4, 5):
try:
import scipy.linalg # noqa F401
except ImportError:

2
phono3py/version.py
View File

@ -34,4 +34,4 @@
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
__version__ = "3.12.1"
__version__ = "3.12.2"

8
test/conductivity/test_kappa_LBTE.py
View File

@ -10,14 +10,14 @@ from phono3py.api_phono3py import Phono3py
@pytest.mark.skipif(
not phono3c.include_lapacke(), reason="test for phono3py compliled with lapacke"
)
@pytest.mark.parametrize("pinv_solver", [1, 2])
def test_kappa_LBTE_12(si_pbesol: Phono3py, pinv_solver: int):
@pytest.mark.parametrize("pinv_solver", [1, 2, 6])
def test_kappa_LBTE_126(si_pbesol: Phono3py, pinv_solver: int):
"""Test for symmetry reduced collision matrix."""
_test_kappa_LBTE(si_pbesol, pinv_solver)
@pytest.mark.parametrize("pinv_solver", [3, 4, 5])
def test_kappa_LBTE_345(si_pbesol: Phono3py, pinv_solver: int):
@pytest.mark.parametrize("pinv_solver", [3, 4, 5, 7])
def test_kappa_LBTE_3457(si_pbesol: Phono3py, pinv_solver: int):
"""Test for symmetry reduced collision matrix."""
_test_kappa_LBTE(si_pbesol, pinv_solver)