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

Bug fix of ph-isotope scattering 

Bug in handling BZ-grid was fixed. Grid points on BZ surface were all included
whwere as only translationally independent ones should be included in the
summation of Tamura's formula. The impact is considered small. Some usually
unused functions that were not migrated to BZ-grid and GRG system were fixed
to work with those grid system. These functions exist for comparisons with
those in C. The tests were added for those python functions.
This commit is contained in:
Atsushi Togo 2022-12-23 13:48:44 +09:00
parent 2e90642ad6
commit c4c54c73a7
14 changed files with 119 additions and 139 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

25
c/_phono3py.c
View File

@ -912,37 +912,44 @@ static PyObject *py_get_isotope_strength(PyObject *self, PyObject *args) {
PyArrayObject *py_eigenvectors;
PyArrayObject *py_band_indices;
PyArrayObject *py_mass_variances;
PyArrayObject *py_ir_grid_points;
PyArrayObject *py_weights;
long grid_point;
long num_grid_points;
double cutoff_frequency;
double sigma;
double *gamma;
double *frequencies;
long *ir_grid_points;
long *weights;
_lapack_complex_double *eigenvectors;
long *band_indices;
double *mass_variances;
long num_band, num_band0;
long num_band, num_band0, num_ir_grid_points;
if (!PyArg_ParseTuple(args, "OlOOOOldd", &py_gamma, &grid_point,
&py_mass_variances, &py_frequencies, &py_eigenvectors,
&py_band_indices, &num_grid_points, &sigma,
&cutoff_frequency)) {
if (!PyArg_ParseTuple(args, "OlOOOOOOdd", &py_gamma, &grid_point,
&py_ir_grid_points, &py_weights, &py_mass_variances,
&py_frequencies, &py_eigenvectors, &py_band_indices,
&sigma, &cutoff_frequency)) {
return NULL;
}
gamma = (double *)PyArray_DATA(py_gamma);
frequencies = (double *)PyArray_DATA(py_frequencies);
eigenvectors = (_lapack_complex_double *)PyArray_DATA(py_eigenvectors);
ir_grid_points = (long *)PyArray_DATA(py_ir_grid_points);
weights = (long *)PyArray_DATA(py_weights);
band_indices = (long *)PyArray_DATA(py_band_indices);
mass_variances = (double *)PyArray_DATA(py_mass_variances);
num_band = (long)PyArray_DIMS(py_frequencies)[1];
num_band0 = (long)PyArray_DIMS(py_band_indices)[0];
num_ir_grid_points = (long)PyArray_DIMS(py_ir_grid_points)[0];
ph3py_get_isotope_scattering_strength(
gamma, grid_point, mass_variances, frequencies, eigenvectors,
num_grid_points, band_indices, num_band, num_band0, sigma,
cutoff_frequency);
gamma, grid_point, ir_grid_points, weights, mass_variances, frequencies,
eigenvectors, num_ir_grid_points, band_indices, num_band, num_band0,
sigma, cutoff_frequency);
Py_RETURN_NONE;
}

16
c/isotope.c
View File

@ -41,11 +41,12 @@
#include "phonoc_utils.h"
void iso_get_isotope_scattering_strength(
double *gamma, const long grid_point, const double *mass_variances,
double *gamma, const long grid_point, const long *ir_grid_points,
const long *weights, const double *mass_variances,
const double *frequencies, const lapack_complex_double *eigenvectors,
const long num_grid_points, const long *band_indices, const long num_band,
const long num_band0, const double sigma, const double cutoff_frequency) {
long i, j, k, l, m;
long i, j, k, l, m, gp;
double *e0_r, *e0_i, e1_r, e1_i, a, b, f, *f0, dist, sum_g, sum_g_k;
e0_r = (double *)malloc(sizeof(double) * num_band * num_band0);
@ -74,13 +75,14 @@ void iso_get_isotope_scattering_strength(
}
sum_g = 0;
#ifdef _OPENMP
#pragma omp parallel for private(k, l, m, f, e1_r, e1_i, a, b, dist, sum_g_k) reduction(+ \
#pragma omp parallel for private(gp, k, l, m, f, e1_r, e1_i, a, b, dist, sum_g_k) reduction(+ \
: sum_g)
#endif
for (j = 0; j < num_grid_points; j++) {
gp = ir_grid_points[j];
sum_g_k = 0;
for (k = 0; k < num_band; k++) { /* band index */
f = frequencies[j * num_band + k];
f = frequencies[gp * num_band + k];
if (f < cutoff_frequency) {
continue;
}
@ -90,10 +92,10 @@ void iso_get_isotope_scattering_strength(
b = 0;
for (m = 0; m < 3; m++) {
e1_r = lapack_complex_double_real(
eigenvectors[j * num_band * num_band +
eigenvectors[gp * num_band * num_band +
(l * 3 + m) * num_band + k]);
e1_i = lapack_complex_double_imag(
eigenvectors[j * num_band * num_band +
eigenvectors[gp * num_band * num_band +
(l * 3 + m) * num_band + k]);
a += (e0_r[i * num_band + l * 3 + m] * e1_r +
e0_i[i * num_band + l * 3 + m] * e1_i);
@ -103,7 +105,7 @@ void iso_get_isotope_scattering_strength(
sum_g_k += (a * a + b * b) * mass_variances[l] * dist;
}
}
sum_g += sum_g_k;
sum_g += sum_g_k * weights[gp];
}
gamma[i] = sum_g;
}

3
c/isotope.h
View File

@ -38,7 +38,8 @@
#include "lapack_wrapper.h"
void iso_get_isotope_scattering_strength(
double *gamma, const long grid_point, const double *mass_variances,
double *gamma, const long grid_point, const long *ir_grid_points,
const long *weights, const double *mass_variances,
const double *frequencies, const lapack_complex_double *eigenvectors,
const long num_grid_points, const long *band_indices, const long num_band,
const long num_band0, const double sigma, const double cutoff_frequency);

12
c/phono3py.c
View File

@ -261,13 +261,15 @@ void ph3py_get_reducible_collision_matrix(
}
void ph3py_get_isotope_scattering_strength(
double *gamma, const long grid_point, const double *mass_variances,
double *gamma, const long grid_point, const long *ir_grid_points,
const long *weights, const double *mass_variances,
const double *frequencies, const _lapack_complex_double *eigenvectors,
const long num_grid_points, const long *band_indices, const long num_band,
const long num_band0, const double sigma, const double cutoff_frequency) {
const long num_ir_grid_points, const long *band_indices,
const long num_band, const long num_band0, const double sigma,
const double cutoff_frequency) {
iso_get_isotope_scattering_strength(
gamma, grid_point, mass_variances, frequencies,
(lapack_complex_double *)eigenvectors, num_grid_points, band_indices,
gamma, grid_point, ir_grid_points, weights, mass_variances, frequencies,
(lapack_complex_double *)eigenvectors, num_ir_grid_points, band_indices,
num_band, num_band0, sigma, cutoff_frequency);
}

8
c/phono3py.h
View File

@ -116,10 +116,12 @@ void ph3py_get_reducible_collision_matrix(
const long num_gp, const double temperature,
const double unit_conversion_factor, const double cutoff_frequency);
void ph3py_get_isotope_scattering_strength(
double *gamma, const long grid_point, const double *mass_variances,
double *gamma, const long grid_point, const long *ir_grid_points,
const long *weights, const double *mass_variances,
const double *frequencies, const _lapack_complex_double *eigenvectors,
const long num_grid_points, const long *band_indices, const long num_band,
const long num_band0, const double sigma, const double cutoff_frequency);
const long num_ir_grid_points, const long *band_indices,
const long num_band, const long num_band0, const double sigma,
const double cutoff_frequency);
void ph3py_get_thm_isotope_scattering_strength(
double *gamma, const long grid_point, const long *ir_grid_points,
const long *weights, const double *mass_variances,

11
c/triplet_iw.c
View File

@ -189,6 +189,17 @@ void tpi_get_integration_weight_with_sigma(
}
}
/**
* @brief Return grid points of relative grid adddresses in BZ-grid
*
* @param neighboring_grid_points Grid points of relative grid addresses in
* BZ-grid.
* @param grid_point Grid point of interest.
* @param relative_grid_address Relative grid address wrt grid point of
* interest.
* @param num_relative_grid_address Number of relative grid addresses.
* @param bzgrid
*/
void tpi_get_neighboring_grid_points(long *neighboring_grid_points,
const long grid_point,
const long (*relative_grid_address)[3],

4
phono3py/api_isotope.py
View File

@ -101,7 +101,7 @@ class Phono3pyIsotope:
"""Return grid points in BZ-grid."""
return self._grid_points
def run(self, grid_points):
def run(self, grid_points, lang="C"):
"""Calculate isotope scattering."""
gamma = np.zeros(
(len(self._sigmas), len(grid_points), len(self._iso.band_indices)),
@ -127,7 +127,7 @@ class Phono3pyIsotope:
else:
print("Sigma: %s" % sigma)
self._iso.sigma = sigma
self._iso.run()
self._iso.run(lang=lang)
gamma[i, j] = self._iso.gamma
frequencies = self._iso.get_phonons()[0]
print("")

4
phono3py/conductivity/base.py
View File

@ -674,10 +674,10 @@ class ConductivityBase(ABC):
Returns
-------
grid_points : ndarray
Grid point indices to be iterated over.
Grid point indices in BZ-grid to be iterated over.
shape=(len(grid_points),), dtype='int_'
ir_grid_points : ndarray
Irreducible grid points on regular grid.
Irreducible grid points in BZ-grid on regular grid.
shape=(len(ir_grid_points),), dtype='int_'
grid_weights : ndarray
Grid weights of `grid_points`. If grid symmetry is not broken,

50
phono3py/other/isotope.py
View File

@ -157,9 +157,12 @@ class Isotope:
if self._phonon_done is None:
self._allocate_phonon()
def run(self):
def run(self, lang="C"):
"""Run isotope scattering calculation."""
self._run_c()
if lang == "C":
self._run_c()
else:
self._run_py()
@property
def sigma(self):
@ -266,14 +269,14 @@ class Isotope:
import phono3py._phono3py as phono3c
gamma = np.zeros(len(self._band_indices), dtype="double")
weights_in_bzgp = np.ones(len(self._grid_points), dtype="int_")
if self._sigma is None:
self._set_integration_weights()
weights = np.ones(len(self._grid_points), dtype="int_")
phono3c.thm_isotope_strength(
gamma,
self._grid_point,
self._grid_points,
weights,
self._bz_grid.grg2bzg,
weights_in_bzgp,
self._mass_variances,
self._frequencies,
self._eigenvectors,
@ -285,21 +288,31 @@ class Isotope:
phono3c.isotope_strength(
gamma,
self._grid_point,
self._bz_grid.grg2bzg,
weights_in_bzgp,
self._mass_variances,
self._frequencies,
self._eigenvectors,
self._band_indices,
np.prod(self._bz_grid.D_diag),
self._sigma,
self._cutoff_frequency,
)
self._gamma = gamma / np.prod(self._bz_grid.D_diag)
def _set_integration_weights(self):
self._set_integration_weights_c()
def _set_integration_weights(self, lang="C"):
if lang == "C":
self._set_integration_weights_c()
else:
self._set_integration_weights_py()
def _set_integration_weights_c(self):
"""Set tetrahedron method integration weights.
self._frequencies are those on all BZ-grid. So all those grid points in
BZ-grid, i.e., self._grid_points, are passed to get_integration_weights.
"""
unique_grid_points = get_unique_grid_points(self._grid_points, self._bz_grid)
self._run_phonon_solver_c(unique_grid_points)
freq_points = np.array(
@ -312,10 +325,11 @@ class Isotope:
)
def _set_integration_weights_py(self):
if self._bz_grid.store_dense_gp_map:
raise NotImplementedError("Only for type-I bz_map.")
if self._bz_grid.grid_matrix is not None:
raise NotImplementedError("Generalized regular grid is not supported.")
"""Set tetrahedron method integration weights.
Python implementation corresponding to _set_integration_weights_c.
"""
thm = TetrahedronMethod(self._bz_grid.microzone_lattice)
num_grid_points = len(self._grid_points)
num_band = len(self._primitive) * 3
@ -325,7 +339,7 @@ class Isotope:
for i, gp in enumerate(self._grid_points):
tfreqs = get_tetrahedra_frequencies(
gp,
gp, # In BZ-grid used only to access self._bz_grid.addresses.
self._bz_grid.D_diag,
self._bz_grid.addresses,
np.array(
@ -333,7 +347,7 @@ class Isotope:
dtype="int_",
order="C",
),
self._grid_points,
self._bz_grid.grg2bzg,
self._frequencies,
grid_order=[
1,
@ -354,16 +368,16 @@ class Isotope:
self._run_phonon_solver_py(gp)
if self._sigma is None:
self._set_integration_weights()
self._set_integration_weights(lang="Py")
t_inv = []
for bi in self._band_indices:
vec0 = self._eigenvectors[self._grid_point][:, bi].conj()
f0 = self._frequencies[self._grid_point][bi]
ti_sum = 0.0
for i, gp in enumerate(self._grid_points):
for gp in self._bz_grid.grg2bzg:
for j, (f, vec) in enumerate(
zip(self._frequencies[i], self._eigenvectors[i].T)
zip(self._frequencies[gp], self._eigenvectors[gp].T)
):
if f < self._cutoff_frequency:
continue
@ -372,7 +386,7 @@ class Isotope:
* self._mass_variances
)
if self._sigma is None:
ti_sum += ti_sum_band * self._integration_weights[i, bi, j]
ti_sum += ti_sum_band * self._integration_weights[gp, bi, j]
else:
ti_sum += ti_sum_band * gaussian(f0 - f, self._sigma)
t_inv.append(np.pi / 2 / np.prod(self._bz_grid.D_diag) * f0**2 * ti_sum)

17
phono3py/other/tetrahedron_method.py
View File

@ -101,26 +101,25 @@ def get_integration_weights(
Values at which the integration weights are computed.
shape=(sampling_points, ), dtype='double'
grid_values : array_like
Values of tetrahedron vertices. Usually they are phonon frequencies,
but the same shape array can be used instead of frequencies.
Values of tetrahedron vertices. Usually they are phonon frequencies, but
the same shape array can be used instead of frequencies.
shape=(regular_grid_points, num_band), dtype='double'
bz_grid : BZGrid
Grid information in reciprocal space.
grid_points : array_like, optional, default=None
Grid point indices in BZ-grid. If None, all regular grid points.
shape=(grid_points, ), dtype='int_'
Grid point indices in BZ-grid. If None, all regular grid points in
BZ-grid. shape=(grid_points, ), dtype='int_'
bzgp2irgp_map : array_like, optional, default=None
Grid point index mapping from bz_grid to index of the first
dimension of `grid_values` array, i.e., usually irreducible
grid point count.
Grid point index mapping from bz_grid to index of the first dimension of
`grid_values` array, i.e., usually irreducible grid point count.
function : str, 'I' or 'J', optional, default='I'
'J' is for intetration and 'I' is for its derivative.
Returns
-------
integration_weights : ndarray
shape=(grid_points, sampling_points, num_band),
dtype='double', order='C'
shape=(grid_points, sampling_points, num_band), dtype='double',
order='C'
"""
import phono3py._phono3py as phono3c

2
phono3py/phonon/grid.py
View File

@ -444,7 +444,7 @@ class BZGrid:
store_dense_gp_map=self._store_dense_gp_map,
)
if self._store_dense_gp_map:
self._grg2bzg = self._gp_map[:-1]
self._grg2bzg = np.array(self._gp_map[:-1], dtype="int_")
else:
self._grg2bzg = np.arange(np.prod(self._D_diag), dtype="int_")

4
test/conductivity/test_kappa_RTA.py
View File

@ -18,8 +18,8 @@ si_pbesol_kappa_RTA_si_nosym = [
]
si_pbesol_kappa_RTA_si_nomeshsym = [81.31304, 81.31304, 81.31304, 0, 0, 0]
si_pbesol_kappa_RTA_grg = [93.99526, 93.99526, 93.99526, 0, 0, 0]
si_pbesol_kappa_RTA_grg_iso = [103.2607, 103.2607, 103.2607, 0, 0, 0]
si_pbesol_kappa_RTA_grg_sigma_iso = [107.9677, 107.9677, 107.9677, 0, 0, 0]
si_pbesol_kappa_RTA_grg_iso = [104.281556, 104.281556, 104.281556, 0, 0, 0]
si_pbesol_kappa_RTA_grg_sigma_iso = [107.2834, 107.2834, 107.2834, 0, 0, 0]
nacl_pbe_kappa_RTA = [7.72798252, 7.72798252, 7.72798252, 0, 0, 0]
nacl_pbe_kappa_RTA_with_sigma = [7.71913708, 7.71913708, 7.71913708, 0, 0, 0]

81
test/conductivity/test_kappa_RTA_Wigner.py
View File

@ -2,98 +2,35 @@
import numpy as np
# first list is k_P, second list is k_C
si_pbesol_kappa_P_RTA = [
108.723,
108.723,
108.723,
0.000,
0.000,
0.000,
] # old value 107.991
si_pbesol_kappa_P_RTA = [108.723, 108.723, 108.723, 0.000, 0.000, 0.000]
si_pbesol_kappa_C = [0.167, 0.167, 0.167, 0.000, 0.000, 0.000]
si_pbesol_kappa_P_RTA_iso = [
97.494,
97.494,
97.494,
0.000,
0.000,
0.000,
] # old value 96.92419
si_pbesol_kappa_P_RTA_iso = [98.008, 98.008, 98.008, 0.000, 0.000, 0.000]
si_pbesol_kappa_C_iso = [0.177, 0.177, 0.177, 0.000, 0.000, 0.000]
si_pbesol_kappa_P_RTA_with_sigmas = [
110.534,
110.534,
110.534,
0,
0,
0,
] # old value 109.6985
si_pbesol_kappa_P_RTA_with_sigmas = [110.534, 110.534, 110.534, 0, 0, 0]
si_pbesol_kappa_C_with_sigmas = [0.163, 0.163, 0.163, 0.000, 0.000, 0.000]
si_pbesol_kappa_P_RTA_with_sigmas_iso = [
97.268,
97.268,
97.268,
0,
0,
0,
] # old value 96.03248
si_pbesol_kappa_P_RTA_with_sigmas_iso = [97.268, 97.268, 97.268, 0, 0, 0]
si_pbesol_kappa_C_with_sigmas_iso = [0.179, 0.179, 0.179, 0.000, 0.000, 0.000]
si_pbesol_kappa_P_RTA_si_nosym = [
39.325, # old value 38.242347
39.323, # old value 38.700219
39.496, # old value 39.198018
-0.004, # old value 0.3216,
0.020, # old value 0.207731,
0.018, # old value 0.283,
]
si_pbesol_kappa_P_RTA_si_nosym = [39.325, 39.323, 39.496, -0.004, 0.020, 0.018]
si_pbesol_kappa_C_si_nosym = [0.009, 0.009, 0.009, 0.000, 0.000, 0.000]
si_pbesol_kappa_P_RTA_si_nomeshsym = [
39.411,
39.411,
39.411,
0,
0,
0,
] # old value 38.90918
si_pbesol_kappa_P_RTA_si_nomeshsym = [39.411, 39.411, 39.411, 0, 0, 0]
si_pbesol_kappa_C_si_nomeshsym = [0.009, 0.009, 0.009, 0.000, 0.000, 0.000]
nacl_pbe_kappa_P_RTA = [
7.753,
7.753,
7.753,
0.000,
0.000,
0.000,
] # old value 7.72798252
nacl_pbe_kappa_P_RTA = [7.753, 7.753, 7.753, 0.000, 0.000, 0.000]
nacl_pbe_kappa_C = [0.081, 0.081, 0.081, 0.000, 0.000, 0.000]
nacl_pbe_kappa_RTA_with_sigma = [7.742, 7.742, 7.742, 0, 0, 0] # old value 7.71913708
nacl_pbe_kappa_C_with_sigma = [0.081, 0.081, 0.081, 0.000, 0.000, 0.000]
aln_lda_kappa_P_RTA = [
203.304,
203.304,
213.003,
0,
0,
0,
] # old value [203.304059, 203.304059, 213.003125, 0, 0, 0]
aln_lda_kappa_P_RTA = [203.304, 203.304, 213.003, 0, 0, 0]
aln_lda_kappa_C = [0.084, 0.084, 0.037, 0, 0, 0]
aln_lda_kappa_P_RTA_with_sigmas = [
213.820000,
213.820000,
224.800000,
0,
0,
0,
] # old value [213.820000, 213.820000, 224.800121, 0, 0, 0]
aln_lda_kappa_P_RTA_with_sigmas = [213.820000, 213.820000, 224.800000, 0, 0, 0]
aln_lda_kappa_C_with_sigmas = [0.084, 0.084, 0.036, 0, 0, 0]

21
test/other/test_isotope.py
View File

@ -1,5 +1,6 @@
"""Tests for isotope scatterings."""
import numpy as np
import pytest
from phono3py import Phono3pyIsotope
@ -49,7 +50,8 @@ si_pbesol_grg_iso_sigma = [
]
def test_Phono3pyIsotope(si_pbesol):
@pytest.mark.parametrize("lang", ["C", "Py"])
def test_Phono3pyIsotope(si_pbesol, lang):
"""Phono3pyIsotope with tetrahedron method."""
si_pbesol.mesh_numbers = [21, 21, 21]
iso = Phono3pyIsotope(
@ -63,12 +65,13 @@ def test_Phono3pyIsotope(si_pbesol):
si_pbesol.phonon_primitive,
nac_params=si_pbesol.nac_params,
)
iso.run([23, 103])
iso.run([23, 103], lang=lang)
# print(iso.gamma[0])
np.testing.assert_allclose(si_pbesol_iso, iso.gamma[0], atol=3e-4)
def test_Phono3pyIsotope_with_sigma(si_pbesol):
@pytest.mark.parametrize("lang", ["C", "Py"])
def test_Phono3pyIsotope_with_sigma(si_pbesol, lang):
"""Phono3pyIsotope with smearing method."""
si_pbesol.mesh_numbers = [21, 21, 21]
iso = Phono3pyIsotope(
@ -85,12 +88,13 @@ def test_Phono3pyIsotope_with_sigma(si_pbesol):
si_pbesol.phonon_primitive,
nac_params=si_pbesol.nac_params,
)
iso.run([23, 103])
iso.run([23, 103], lang=lang)
# print(iso.gamma[0])
np.testing.assert_allclose(si_pbesol_iso_sigma, iso.gamma[0], atol=3e-4)
def test_Phono3pyIsotope_grg(si_pbesol_grg):
@pytest.mark.parametrize("lang", ["C", "Py"])
def test_Phono3pyIsotope_grg(si_pbesol_grg, lang):
"""Phono3pyIsotope with tetrahedron method and GR-grid."""
ph3 = si_pbesol_grg
iso = Phono3pyIsotope(
@ -105,14 +109,15 @@ def test_Phono3pyIsotope_grg(si_pbesol_grg):
ph3.phonon_primitive,
nac_params=ph3.nac_params,
)
iso.run([23, 103])
np.testing.assert_equal(
iso.grid.grid_matrix, [[-15, 15, 15], [15, -15, 15], [15, 15, -15]]
)
iso.run([23, 103], lang=lang)
np.testing.assert_allclose(si_pbesol_grg_iso, iso.gamma[0], atol=2e-3)
def test_Phono3pyIsotope_grg_with_sigma(si_pbesol_grg):
@pytest.mark.parametrize("lang", ["C", "Py"])
def test_Phono3pyIsotope_grg_with_sigma(si_pbesol_grg, lang):
"""Phono3pyIsotope with smearing method and GR-grid."""
ph3 = si_pbesol_grg
iso = Phono3pyIsotope(
@ -130,7 +135,7 @@ def test_Phono3pyIsotope_grg_with_sigma(si_pbesol_grg):
ph3.phonon_primitive,
nac_params=ph3.nac_params,
)
iso.run([23, 103])
iso.run([23, 103], lang=lang)
np.testing.assert_equal(
iso.grid.grid_matrix, [[-15, 15, 15], [15, -15, 15], [15, 15, -15]]
)