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

Formatted.

This commit is contained in:
Atsushi Togo 2021-10-12 18:12:56 +09:00
parent 117fe05c24
commit 55df7fb796
3 changed files with 1211 additions and 972 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

306
phono3py/phonon3/conductivity.py
View File

@ -43,11 +43,11 @@ from phonopy.units import THzToEv, EV, THz, Angstrom
from phono3py.file_IO import write_pp_to_hdf5
from phono3py.phonon3.triplets import get_all_triplets
from phono3py.other.isotope import Isotope
from phono3py.phonon.grid import (get_ir_grid_points,
get_grid_points_by_rotations)
from phono3py.phonon.grid import get_ir_grid_points, get_grid_points_by_rotations
unit_to_WmK = ((THz * Angstrom) ** 2 / (Angstrom ** 3) * EV / THz /
(2 * np.pi)) # 2pi comes from definition of lifetime.
unit_to_WmK = (
(THz * Angstrom) ** 2 / (Angstrom ** 3) * EV / THz / (2 * np.pi)
) # 2pi comes from definition of lifetime.
def all_bands_exist(interaction):
@ -60,11 +60,7 @@ def all_bands_exist(interaction):
return False
def write_pp(conductivity,
pp,
i,
filename=None,
compression="gzip"):
def write_pp(conductivity, pp, i, filename=None, compression="gzip"):
"""Write ph-ph interaction strength in hdf5 file."""
grid_point = conductivity.grid_points[i]
sigmas = conductivity.sigmas
@ -77,7 +73,8 @@ def write_pp(conductivity,
print("Multiple smearing parameters were given. The last one in ")
print("ph-ph interaction calculations was written in the file.")
write_pp_to_hdf5(mesh,
write_pp_to_hdf5(
mesh,
pp=pp.interaction_strength,
g_zero=pp.zero_value_positions,
grid_point=grid_point,
@ -87,17 +84,16 @@ def write_pp(conductivity,
sigma=sigmas[-1],
sigma_cutoff=sigma_cutoff,
filename=filename,
compression=compression)
compression=compression,
)
class ConductivityBase:
"""Base class of Conductivity classes."""
def __init__(self,
interaction,
grid_points=None,
is_kappa_star=True,
gv_delta_q=None):
def __init__(
self, interaction, grid_points=None, is_kappa_star=True, gv_delta_q=None
):
"""Init method.
interaction : Interaction
@ -130,36 +126,43 @@ class ConductivityBase:
self._set_grid_properties(grid_points)
self._qpoints = np.array(
np.dot(self._pp.bz_grid.addresses[self._grid_points],
self._pp.bz_grid.QDinv.T), dtype='double', order='C')
np.dot(
self._pp.bz_grid.addresses[self._grid_points], self._pp.bz_grid.QDinv.T
),
dtype="double",
order="C",
)
self._gv_obj = GroupVelocity(
self._pp.dynamical_matrix,
q_length=gv_delta_q,
symmetry=self._pp.primitive_symmetry,
frequency_factor_to_THz=self._pp.frequency_factor_to_THz)
frequency_factor_to_THz=self._pp.frequency_factor_to_THz,
)
def _set_point_operations(self):
if not self._is_kappa_star:
self._point_operations = np.array(
[np.eye(3, dtype='int_')], dtype='int_', order='C')
[np.eye(3, dtype="int_")], dtype="int_", order="C"
)
else:
self._point_operations = np.array(
self._pp.bz_grid.reciprocal_operations,
dtype='int_', order='C')
self._pp.bz_grid.reciprocal_operations, dtype="int_", order="C"
)
rec_lat = np.linalg.inv(self._pp.primitive.cell)
self._rotations_cartesian = np.array(
[similarity_transformation(rec_lat, r)
for r in self._point_operations], dtype='double', order='C')
[similarity_transformation(rec_lat, r) for r in self._point_operations],
dtype="double",
order="C",
)
def _set_grid_properties(self, grid_points):
if grid_points is not None: # Specify grid points
self._grid_points = grid_points
(self._ir_grid_points,
self._ir_grid_weights) = self._get_ir_grid_points()
(self._ir_grid_points, self._ir_grid_weights) = self._get_ir_grid_points()
elif not self._is_kappa_star: # All grid points
self._grid_points = self._pp.bz_grid.grg2bzg
self._grid_weights = np.ones(len(self._grid_points), dtype='int_')
self._grid_weights = np.ones(len(self._grid_points), dtype="int_")
self._ir_grid_points = self._grid_points
self._ir_grid_weights = self._grid_weights
else: # Automatic sampling
@ -169,17 +172,18 @@ class ConductivityBase:
def _get_ir_grid_points(self):
"""Find irreducible grid points."""
ir_grid_points, ir_grid_weights, _ = get_ir_grid_points(
self._pp.bz_grid)
ir_grid_points, ir_grid_weights, _ = get_ir_grid_points(self._pp.bz_grid)
ir_grid_points = np.array(
self._pp.bz_grid.grg2bzg[ir_grid_points], dtype='int_')
self._pp.bz_grid.grg2bzg[ir_grid_points], dtype="int_"
)
return ir_grid_points, ir_grid_weights
class Conductivity(ConductivityBase):
"""Thermal conductivity base class."""
def __init__(self,
def __init__(
self,
interaction,
grid_points=None,
temperatures=None,
@ -191,7 +195,8 @@ class Conductivity(ConductivityBase):
is_kappa_star=True,
gv_delta_q=None, # finite difference for group veolocity
is_full_pp=False,
log_level=0):
log_level=0,
):
"""Init method."""
self._pp = None
self._is_kappa_star = None
@ -203,10 +208,12 @@ class Conductivity(ConductivityBase):
self._point_operations = None
self._gv_obj = None
super().__init__(interaction,
super().__init__(
interaction,
grid_points=grid_points,
is_kappa_star=is_kappa_star,
gv_delta_q=gv_delta_q)
gv_delta_q=gv_delta_q,
)
self._is_full_pp = is_full_pp
self._log_level = log_level
@ -221,7 +228,7 @@ class Conductivity(ConductivityBase):
if temperatures is None:
self._temperatures = None
else:
self._temperatures = np.array(temperatures, dtype='double')
self._temperatures = np.array(temperatures, dtype="double")
self._boundary_mfp = boundary_mfp
self._read_gamma = False
@ -241,9 +248,11 @@ class Conductivity(ConductivityBase):
self._conversion_factor = unit_to_WmK / volume
self._pp.nac_q_direction = None
(self._frequencies,
(
self._frequencies,
self._eigenvectors,
self._phonon_done) = self._pp.get_phonons()
self._phonon_done,
) = self._pp.get_phonons()
if (self._phonon_done == 0).any():
self._pp.run_phonon_solver()
@ -263,8 +272,10 @@ class Conductivity(ConductivityBase):
"""Return grid point count for mode kappa."""
if self._grid_point_count == len(self._grid_points):
if self._log_level:
print("=================== End of collection of collisions "
"===================")
print(
"=================== End of collection of collisions "
"==================="
)
raise StopIteration
else:
self._run_at_grid_point()
@ -282,9 +293,11 @@ class Conductivity(ConductivityBase):
def get_mesh_numbers(self):
"""Return mesh numbers of GR-grid."""
warnings.warn("Use attribute, Conductivity.mesh_numbers "
warnings.warn(
"Use attribute, Conductivity.mesh_numbers "
"instead of Conductivity.get_mesh_numbers().",
DeprecationWarning)
DeprecationWarning,
)
return self.mesh_numbers
@property
@ -302,9 +315,11 @@ class Conductivity(ConductivityBase):
Grid points are those at mode kappa are calculated.
"""
warnings.warn("Use attribute, Conductivity.mode_heat_capacities "
warnings.warn(
"Use attribute, Conductivity.mode_heat_capacities "
"instead of Conductivity.get_mode_heat_capacities().",
DeprecationWarning)
DeprecationWarning,
)
return self.mode_heat_capacities
@property
@ -322,9 +337,11 @@ class Conductivity(ConductivityBase):
Grid points are those at mode kappa are calculated.
"""
warnings.warn("Use attribute, Conductivity.group_velocities "
warnings.warn(
"Use attribute, Conductivity.group_velocities "
"instead of Conductivity.get_group_velocities().",
DeprecationWarning)
DeprecationWarning,
)
return self.group_velocities
@property
@ -334,9 +351,11 @@ class Conductivity(ConductivityBase):
def get_gv_by_gv(self):
"""Return gv_by_gv at grid points where mode kappa are calculated."""
warnings.warn("Use attribute, Conductivity.gv_by_gv "
warnings.warn(
"Use attribute, Conductivity.gv_by_gv "
"instead of Conductivity.get_gv_by_gv().",
DeprecationWarning)
DeprecationWarning,
)
return self.gv_by_gv
@property
@ -354,9 +373,11 @@ class Conductivity(ConductivityBase):
Grid points are those at mode kappa are calculated.
"""
warnings.warn("Use attribute, Conductivity.frequencies "
warnings.warn(
"Use attribute, Conductivity.frequencies "
"instead of Conductivity.get_frequencies().",
DeprecationWarning)
DeprecationWarning,
)
return self.frequencies
@property
@ -366,9 +387,11 @@ class Conductivity(ConductivityBase):
def get_qpoints(self):
"""Return q-points where mode kappa are calculated."""
warnings.warn("Use attribute, Conductivity.qpoints "
warnings.warn(
"Use attribute, Conductivity.qpoints "
"instead of Conductivity.get_qpoints().",
DeprecationWarning)
DeprecationWarning,
)
return self.qpoints
@property
@ -386,9 +409,11 @@ class Conductivity(ConductivityBase):
Grid point indices are given in BZ-grid.
"""
warnings.warn("Use attribute, Conductivity.grid_points "
warnings.warn(
"Use attribute, Conductivity.grid_points "
"instead of Conductivity.get_grid_points().",
DeprecationWarning)
DeprecationWarning,
)
return self.grid_points
@property
@ -398,9 +423,11 @@ class Conductivity(ConductivityBase):
def get_grid_weights(self):
"""Return grid point weights where mode kappa are calculated."""
warnings.warn("Use attribute, Conductivity.grid_weights "
warnings.warn(
"Use attribute, Conductivity.grid_weights "
"instead of Conductivity.get_grid_weights().",
DeprecationWarning)
DeprecationWarning,
)
return self.grid_weights
@property
@ -415,16 +442,20 @@ class Conductivity(ConductivityBase):
def get_temperatures(self):
"""Return temperatures."""
warnings.warn("Use attribute, Conductivity.temperatures "
warnings.warn(
"Use attribute, Conductivity.temperatures "
"instead of Conductivity.get_temperatures().",
DeprecationWarning)
DeprecationWarning,
)
return self.temperatures
def set_temperatures(self, temperatures):
"""Set temperatures."""
warnings.warn("Use attribute, Conductivity.temperatures "
warnings.warn(
"Use attribute, Conductivity.temperatures "
"instead of Conductivity.set_temperatures().",
DeprecationWarning)
DeprecationWarning,
)
self.temperatures = temperatures
@property
@ -439,16 +470,18 @@ class Conductivity(ConductivityBase):
def get_gamma(self):
"""Return gamma."""
warnings.warn("Use attribute, Conductivity.gamma "
"instead of Conductivity.get_gamma().",
DeprecationWarning)
warnings.warn(
"Use attribute, Conductivity.gamma " "instead of Conductivity.get_gamma().",
DeprecationWarning,
)
return self.gamma
def set_gamma(self, gamma):
"""Set gamma."""
warnings.warn("Use attribute, Conductivity.gamma "
"instead of Conductivity.set_gamma().",
DeprecationWarning)
warnings.warn(
"Use attribute, Conductivity.gamma " "instead of Conductivity.set_gamma().",
DeprecationWarning,
)
self.gamma = gamma
@property
@ -463,16 +496,20 @@ class Conductivity(ConductivityBase):
def get_gamma_isotope(self):
"""Return gamma from isotope."""
warnings.warn("Use attribute, Conductivity.gamma_isotope "
warnings.warn(
"Use attribute, Conductivity.gamma_isotope "
"instead of Conductivity.get_gamma_isotope().",
DeprecationWarning)
DeprecationWarning,
)
return self.gamma_isotope
def set_gamma_isotope(self, gamma_iso):
"""Set gamma from isotope."""
warnings.warn("Use attribute, Conductivity.gamma_isotope "
warnings.warn(
"Use attribute, Conductivity.gamma_isotope "
"instead of Conductivity.set_gamma_isotope().",
DeprecationWarning)
DeprecationWarning,
)
self.gamma_isotope = gamma_iso
@property
@ -482,9 +519,10 @@ class Conductivity(ConductivityBase):
def get_kappa(self):
"""Return kappa."""
warnings.warn("Use attribute, Conductivity.kappa "
"instead of Conductivity.get_kappa().",
DeprecationWarning)
warnings.warn(
"Use attribute, Conductivity.kappa " "instead of Conductivity.get_kappa().",
DeprecationWarning,
)
return self.kappa
@property
@ -494,9 +532,11 @@ class Conductivity(ConductivityBase):
def get_mode_kappa(self):
"""Return mode_kappa."""
warnings.warn("Use attribute, Conductivity.mode_kappa "
warnings.warn(
"Use attribute, Conductivity.mode_kappa "
"instead of Conductivity.get_mode_kappa().",
DeprecationWarning)
DeprecationWarning,
)
return self.mode_kappa
@property
@ -506,9 +546,11 @@ class Conductivity(ConductivityBase):
def get_sigmas(self):
"""Return sigmas."""
warnings.warn("Use attribute, Conductivity.sigmas "
warnings.warn(
"Use attribute, Conductivity.sigmas "
"instead of Conductivity.get_sigmas().",
DeprecationWarning)
DeprecationWarning,
)
return self.sigmas
@property
@ -518,9 +560,11 @@ class Conductivity(ConductivityBase):
def get_sigma_cutoff_width(self):
"""Return smearing width cutoff."""
warnings.warn("Use attribute, Conductivity.sigma_cutoff_width "
warnings.warn(
"Use attribute, Conductivity.sigma_cutoff_width "
"instead of Conductivity.get_sigma_cutoff_width().",
DeprecationWarning)
DeprecationWarning,
)
return self.sigma_cutoff_width
@property
@ -530,9 +574,11 @@ class Conductivity(ConductivityBase):
def get_grid_point_count(self):
"""Return interator count of self."""
warnings.warn("Use attribute, Conductivity.grid_point_count "
warnings.warn(
"Use attribute, Conductivity.grid_point_count "
"instead of Conductivity.get_grid_point_count().",
DeprecationWarning)
DeprecationWarning,
)
return self.grid_point_count
@property
@ -542,9 +588,11 @@ class Conductivity(ConductivityBase):
def get_averaged_pp_interaction(self):
"""Return averaged pp interaction strength."""
warnings.warn("Use attribute, Conductivity.averaged_pp_interaction "
warnings.warn(
"Use attribute, Conductivity.averaged_pp_interaction "
"instead of Conductivity.get_averaged_pp_interaction().",
DeprecationWarning)
DeprecationWarning,
)
return self.averaged_pp_interaction
def get_number_of_sampling_grid_points(self):
@ -577,16 +625,18 @@ class Conductivity(ConductivityBase):
print(text)
self._isotope.sigma = sigma
self._isotope.set_phonons(self._frequencies,
self._isotope.set_phonons(
self._frequencies,
self._eigenvectors,
self._phonon_done,
dm=self._pp.dynamical_matrix)
dm=self._pp.dynamical_matrix,
)
gp = self._grid_points[i]
self._isotope.set_grid_point(gp)
self._isotope.run()
gamma_iso.append(self._isotope.gamma)
return np.array(gamma_iso, dtype='double', order='C')
return np.array(gamma_iso, dtype="double", order="C")
def _set_isotope(self, mass_variances):
if mass_variances is True:
@ -601,7 +651,8 @@ class Conductivity(ConductivityBase):
frequency_factor_to_THz=self._pp.frequency_factor_to_THz,
symprec=self._pp.primitive_symmetry.tolerance,
cutoff_frequency=self._pp.cutoff_frequency,
lapack_zheev_uplo=self._pp.lapack_zheev_uplo)
lapack_zheev_uplo=self._pp.lapack_zheev_uplo,
)
self._mass_variances = self._isotope.mass_variances
def _set_harmonic_properties(self, i_irgp, i_data):
@ -609,21 +660,28 @@ class Conductivity(ConductivityBase):
grid_point = self._grid_points[i_irgp]
freqs = self._frequencies[grid_point][self._pp.band_indices]
self._cv[:, i_data, :] = self._get_cv(freqs)
self._gv_obj.run([self._qpoints[i_irgp], ])
self._gv_obj.run(
[
self._qpoints[i_irgp],
]
)
gv = self._gv_obj.group_velocities[0, self._pp.band_indices, :]
self._gv[i_data] = gv
def _get_cv(self, freqs):
cv = np.zeros((len(self._temperatures), len(freqs)), dtype='double')
cv = np.zeros((len(self._temperatures), len(freqs)), dtype="double")
# T/freq has to be large enough to avoid divergence.
# Otherwise just set 0.
for i, f in enumerate(freqs):
finite_t = (self._temperatures > f / 100)
finite_t = self._temperatures > f / 100
if f > self._pp.cutoff_frequency:
cv[:, i] = np.where(
finite_t, get_mode_cv(
np.where(finite_t, self._temperatures, 10000),
f * THzToEv), 0)
finite_t,
get_mode_cv(
np.where(finite_t, self._temperatures, 10000), f * THzToEv
),
0,
)
return cv
def _set_gv_by_gv(self, i_irgp, i_data):
@ -636,23 +694,23 @@ class Conductivity(ConductivityBase):
self._num_sampling_grid_points += order_kstar
# Sum all vxv at k*
for j, vxv in enumerate(
([0, 0], [1, 1], [2, 2], [1, 2], [0, 2], [0, 1])):
for j, vxv in enumerate(([0, 0], [1, 1], [2, 2], [1, 2], [0, 2], [0, 1])):
self._gv_sum2[i_data, :, j] = gv_by_gv_tensor[:, vxv[0], vxv[1]]
def _get_gv_by_gv(self, i_irgp, i_data):
if self._is_kappa_star:
rotation_map = get_grid_points_by_rotations(
self._grid_points[i_irgp],
self._pp.bz_grid)
self._grid_points[i_irgp], self._pp.bz_grid
)
else:
rotation_map = get_grid_points_by_rotations(
self._grid_points[i_irgp],
self._pp.bz_grid,
reciprocal_rotations=self._point_operations)
reciprocal_rotations=self._point_operations,
)
gv = self._gv[i_data]
gv_by_gv = np.zeros((len(gv), 3, 3), dtype='double')
gv_by_gv = np.zeros((len(gv), 3, 3), dtype="double")
for r in self._rotations_cartesian:
gvs_rot = np.dot(gv, r.T)
@ -663,15 +721,16 @@ class Conductivity(ConductivityBase):
if self._grid_weights is not None:
if order_kstar != self._grid_weights[i_irgp]:
if self._log_level:
text = ("Number of elements in k* is unequal "
text = (
"Number of elements in k* is unequal "
"to number of equivalent grid-points. "
"This means that the mesh sampling grids break "
"symmetry. Please check carefully "
"the convergence over grid point densities.")
msg = textwrap.fill(text,
initial_indent=" ",
subsequent_indent=" ",
width=70)
"the convergence over grid point densities."
)
msg = textwrap.fill(
text, initial_indent=" ", subsequent_indent=" ", width=70
)
print("*" * 30 + "Warning" + "*" * 30)
print(msg)
print("*" * 67)
@ -688,27 +747,40 @@ class Conductivity(ConductivityBase):
def _get_boundary_scattering(self, i):
num_band = len(self._pp.primitive) * 3
g_boundary = np.zeros(num_band, dtype='double')
g_boundary = np.zeros(num_band, dtype="double")
for ll in range(num_band):
g_boundary[ll] = (np.linalg.norm(self._gv[i, ll]) * Angstrom * 1e6
/ (4 * np.pi * self._boundary_mfp))
g_boundary[ll] = (
np.linalg.norm(self._gv[i, ll])
* Angstrom
* 1e6
/ (4 * np.pi * self._boundary_mfp)
)
return g_boundary
def _show_log_header(self, i):
if self._log_level:
gp = self._grid_points[i]
print("======================= Grid point %d (%d/%d) "
"=======================" %
(gp, i + 1, len(self._grid_points)))
print(
"======================= Grid point %d (%d/%d) "
"=======================" % (gp, i + 1, len(self._grid_points))
)
print("q-point: (%5.2f %5.2f %5.2f)" % tuple(self._qpoints[i]))
if self._boundary_mfp is not None:
if self._boundary_mfp > 1000:
print("Boundary mean free path (millimetre): %.3f" %
(self._boundary_mfp / 1000.0))
print(
"Boundary mean free path (millimetre): %.3f"
% (self._boundary_mfp / 1000.0)
)
else:
print("Boundary mean free path (micrometre): %.5f" %
self._boundary_mfp)
print(
"Boundary mean free path (micrometre): %.5f"
% self._boundary_mfp
)
if self._is_isotope:
print(("Mass variance parameters: " +
"%5.2e " * len(self._mass_variances)) %
tuple(self._mass_variances))
print(
(
"Mass variance parameters: "
+ "%5.2e " * len(self._mass_variances)
)
% tuple(self._mass_variances)
)

729
phono3py/phonon3/conductivity_LBTE.py
View File

File diff suppressed because it is too large Load Diff

402
phono3py/phonon3/conductivity_RTA.py
View File

@ -36,13 +36,15 @@
import sys
import numpy as np
from phonopy.structure.tetrahedron_method import TetrahedronMethod
from phono3py.file_IO import (write_kappa_to_hdf5, read_gamma_from_hdf5,
write_gamma_detail_to_hdf5, read_pp_from_hdf5)
from phono3py.phonon3.conductivity import (Conductivity, all_bands_exist,
unit_to_WmK)
from phono3py.file_IO import (
write_kappa_to_hdf5,
read_gamma_from_hdf5,
write_gamma_detail_to_hdf5,
read_pp_from_hdf5,
)
from phono3py.phonon3.conductivity import Conductivity, all_bands_exist, unit_to_WmK
from phono3py.phonon3.conductivity import write_pp as _write_pp
from phono3py.phonon3.imag_self_energy import (ImagSelfEnergy,
average_by_degeneracy)
from phono3py.phonon3.imag_self_energy import ImagSelfEnergy, average_by_degeneracy
from phono3py.phonon3.triplets import get_all_triplets
from phono3py.phonon.grid import get_grid_points_by_rotations
@ -71,16 +73,19 @@ def get_thermal_conductivity_RTA(
compression="gzip",
input_filename=None,
output_filename=None,
log_level=0):
log_level=0,
):
"""Run RTA thermal conductivity calculation."""
if temperatures is None:
_temperatures = np.arange(0, 1001, 10, dtype='double')
_temperatures = np.arange(0, 1001, 10, dtype="double")
else:
_temperatures = temperatures
if log_level:
print("-------------------- Lattice thermal conducitivity (RTA) "
"--------------------")
print(
"-------------------- Lattice thermal conducitivity (RTA) "
"--------------------"
)
br = Conductivity_RTA(
interaction,
grid_points=grid_points,
@ -100,7 +105,8 @@ def get_thermal_conductivity_RTA(
pp_filename=input_filename,
is_N_U=is_N_U,
is_gamma_detail=write_gamma_detail,
log_level=log_level)
log_level=log_level,
)
if read_gamma:
if not _set_gamma_from_file(br, filename=input_filename):
@ -109,42 +115,47 @@ def get_thermal_conductivity_RTA(
for i in br:
if write_pp:
_write_pp(br,
interaction,
i,
compression=compression,
filename=output_filename)
_write_pp(
br, interaction, i, compression=compression, filename=output_filename
)
if write_gamma:
_write_gamma(br,
_write_gamma(
br,
interaction,
i,
compression=compression,
filename=output_filename,
verbose=log_level)
verbose=log_level,
)
if write_gamma_detail:
_write_gamma_detail(br,
_write_gamma_detail(
br,
interaction,
i,
compression=compression,
filename=output_filename,
verbose=log_level)
verbose=log_level,
)
if grid_points is None and all_bands_exist(interaction):
br.set_kappa_at_sigmas()
if log_level:
_show_kappa(br, log_level)
if write_kappa:
_write_kappa(br,
_write_kappa(
br,
interaction.primitive.volume,
compression=compression,
filename=output_filename,
log_level=log_level)
log_level=log_level,
)
return br
def _write_gamma_detail(br, interaction, i, compression="gzip", filename=None,
verbose=True):
def _write_gamma_detail(
br, interaction, i, compression="gzip", filename=None, verbose=True
):
gamma_detail = br.get_gamma_detail_at_q()
temperatures = br.get_temperatures()
mesh = br.get_mesh_numbers()
@ -169,7 +180,8 @@ def _write_gamma_detail(br, interaction, i, compression="gzip", filename=None,
sigma_cutoff=sigma_cutoff,
compression=compression,
filename=filename,
verbose=verbose)
verbose=verbose,
)
else:
for j, sigma in enumerate(sigmas):
for k, bi in enumerate(interaction.get_band_indices()):
@ -185,11 +197,11 @@ def _write_gamma_detail(br, interaction, i, compression="gzip", filename=None,
sigma_cutoff=sigma_cutoff,
compression=compression,
filename=filename,
verbose=verbose)
verbose=verbose,
)
def _write_gamma(br, interaction, i, compression="gzip", filename=None,
verbose=True):
def _write_gamma(br, interaction, i, compression="gzip", filename=None, verbose=True):
"""Write mode kappa related properties into a hdf5 file."""
grid_points = br.grid_points
group_velocities = br.group_velocities
@ -226,7 +238,8 @@ def _write_gamma(br, interaction, i, compression="gzip", filename=None,
else:
gamma_U_at_sigma = gamma_U[j, :, i]
write_kappa_to_hdf5(temperatures,
write_kappa_to_hdf5(
temperatures,
mesh,
frequency=frequencies,
group_velocity=group_velocities[i],
@ -243,7 +256,8 @@ def _write_gamma(br, interaction, i, compression="gzip", filename=None,
kappa_unit_conversion=unit_to_WmK / volume,
compression=compression,
filename=filename,
verbose=verbose)
verbose=verbose,
)
else:
for j, sigma in enumerate(sigmas):
for k, bi in enumerate(interaction.band_indices):
@ -283,7 +297,8 @@ def _write_gamma(br, interaction, i, compression="gzip", filename=None,
kappa_unit_conversion=unit_to_WmK / volume,
compression=compression,
filename=filename,
verbose=verbose)
verbose=verbose,
)
def _show_kappa(br, log_level):
@ -301,18 +316,27 @@ def _show_kappa(br, log_level):
text += "with tetrahedron method -----------"
print(text)
if log_level > 1:
print(("#%6s " + " %-10s" * 6 + "#ipm") %
("T(K)", "xx", "yy", "zz", "yz", "xz", "xy"))
print(
("#%6s " + " %-10s" * 6 + "#ipm")
% ("T(K)", "xx", "yy", "zz", "yz", "xz", "xy")
)
for j, (t, k) in enumerate(zip(temperatures, kappa[i])):
print(("%7.1f" + " %10.3f" * 6 + " %d/%d") %
((t,) + tuple(k) +
(num_ignored_phonon_modes[i, j], num_phonon_modes)))
print(
("%7.1f" + " %10.3f" * 6 + " %d/%d")
% (
(t,)
+ tuple(k)
+ (num_ignored_phonon_modes[i, j], num_phonon_modes)
)
)
else:
print(("#%6s " + " %-10s" * 6) %
("T(K)", "xx", "yy", "zz", "yz", "xz", "xy"))
print(
("#%6s " + " %-10s" * 6)
% ("T(K)", "xx", "yy", "zz", "yz", "xz", "xy")
)
for j, (t, k) in enumerate(zip(temperatures, kappa[i])):
print(("%7.1f " + " %10.3f" * 6) % ((t,) + tuple(k)))
print('')
print("")
def _write_kappa(br, volume, compression="gzip", filename=None, log_level=0):
@ -348,7 +372,8 @@ def _write_kappa(br, volume, compression="gzip", filename=None, log_level=0):
else:
gamma_U_at_sigma = gamma_U[i]
write_kappa_to_hdf5(temperatures,
write_kappa_to_hdf5(
temperatures,
mesh,
frequency=frequencies,
group_velocity=gv,
@ -368,7 +393,8 @@ def _write_kappa(br, volume, compression="gzip", filename=None, log_level=0):
kappa_unit_conversion=unit_to_WmK / volume,
compression=compression,
filename=filename,
verbose=log_level)
verbose=log_level,
)
def _set_gamma_from_file(br, filename=None, verbose=True):
@ -380,16 +406,13 @@ def _set_gamma_from_file(br, filename=None, verbose=True):
temperatures = br.get_temperatures()
num_band = br.get_frequencies().shape[1]
gamma = np.zeros((len(sigmas),
len(temperatures),
len(grid_points),
num_band), dtype='double')
gamma = np.zeros(
(len(sigmas), len(temperatures), len(grid_points), num_band), dtype="double"
)
gamma_N = np.zeros_like(gamma)
gamma_U = np.zeros_like(gamma)
gamma_iso = np.zeros((len(sigmas),
len(grid_points),
num_band), dtype='double')
ave_pp = np.zeros((len(grid_points), num_band), dtype='double')
gamma_iso = np.zeros((len(sigmas), len(grid_points), num_band), dtype="double")
ave_pp = np.zeros((len(grid_points), num_band), dtype="double")
is_gamma_N_U_in = False
is_ave_pp_in = False
@ -401,18 +424,19 @@ def _set_gamma_from_file(br, filename=None, verbose=True):
sigma=sigma,
sigma_cutoff=sigma_cutoff,
filename=filename,
verbose=verbose)
verbose=verbose,
)
if data:
gamma[j] = data['gamma']
if 'gamma_isotope' in data:
gamma_iso[j] = data['gamma_isotope']
if 'gamma_N' in data:
gamma[j] = data["gamma"]
if "gamma_isotope" in data:
gamma_iso[j] = data["gamma_isotope"]
if "gamma_N" in data:
is_gamma_N_U_in = True
gamma_N[j] = data['gamma_N']
gamma_U[j] = data['gamma_U']
if 'ave_pp' in data:
gamma_N[j] = data["gamma_N"]
gamma_U[j] = data["gamma_U"]
if "ave_pp" in data:
is_ave_pp_in = True
ave_pp[:] = data['ave_pp']
ave_pp[:] = data["ave_pp"]
else:
for i, gp in enumerate(grid_points):
data_gp = read_gamma_from_hdf5(
@ -421,18 +445,19 @@ def _set_gamma_from_file(br, filename=None, verbose=True):
sigma=sigma,
sigma_cutoff=sigma_cutoff,
filename=filename,
verbose=verbose)
verbose=verbose,
)
if data_gp:
gamma[j, :, i] = data_gp['gamma']
if 'gamma_iso' in data_gp:
gamma_iso[j, i] = data_gp['gamma_iso']
if 'gamma_N' in data_gp:
gamma[j, :, i] = data_gp["gamma"]
if "gamma_iso" in data_gp:
gamma_iso[j, i] = data_gp["gamma_iso"]
if "gamma_N" in data_gp:
is_gamma_N_U_in = True
gamma_N[j, :, i] = data_gp['gamma_N']
gamma_U[j, :, i] = data_gp['gamma_U']
if 'ave_pp' in data_gp:
gamma_N[j, :, i] = data_gp["gamma_N"]
gamma_U[j, :, i] = data_gp["gamma_U"]
if "ave_pp" in data_gp:
is_ave_pp_in = True
ave_pp[i] = data_gp['ave_pp']
ave_pp[i] = data_gp["ave_pp"]
else:
for bi in range(num_band):
data_band = read_gamma_from_hdf5(
@ -442,18 +467,19 @@ def _set_gamma_from_file(br, filename=None, verbose=True):
sigma=sigma,
sigma_cutoff=sigma_cutoff,
filename=filename,
verbose=verbose)
verbose=verbose,
)
if data_band:
gamma[j, :, i, bi] = data_band['gamma']
if 'gamma_iso' in data_band:
gamma_iso[j, i, bi] = data_band['gamma_iso']
if 'gamma_N' in data_band:
gamma[j, :, i, bi] = data_band["gamma"]
if "gamma_iso" in data_band:
gamma_iso[j, i, bi] = data_band["gamma_iso"]
if "gamma_N" in data_band:
is_gamma_N_U_in = True
gamma_N[j, :, i, bi] = data_band['gamma_N']
gamma_U[j, :, i, bi] = data_band['gamma_U']
if 'ave_pp' in data_band:
gamma_N[j, :, i, bi] = data_band["gamma_N"]
gamma_U[j, :, i, bi] = data_band["gamma_U"]
if "ave_pp" in data_band:
is_ave_pp_in = True
ave_pp[i, bi] = data_band['ave_pp']
ave_pp[i, bi] = data_band["ave_pp"]
else:
read_succeeded = False
@ -471,7 +497,8 @@ def _set_gamma_from_file(br, filename=None, verbose=True):
class Conductivity_RTA(Conductivity):
"""Lattice thermal conductivity calculation with RTA."""
def __init__(self,
def __init__(
self,
interaction,
grid_points=None,
temperatures=None,
@ -491,7 +518,8 @@ class Conductivity_RTA(Conductivity):
is_N_U=False,
is_gamma_detail=False,
is_frequency_shift_by_bubble=False,
log_level=0):
log_level=0,
):
"""Init method."""
self._pp = None
self._gv_obj = None
@ -538,7 +566,8 @@ class Conductivity_RTA(Conductivity):
self._mass_variances = None
self._grid_point_count = None
Conductivity.__init__(self,
Conductivity.__init__(
self,
interaction,
grid_points=grid_points,
temperatures=temperatures,
@ -550,7 +579,8 @@ class Conductivity_RTA(Conductivity):
is_kappa_star=is_kappa_star,
gv_delta_q=gv_delta_q,
is_full_pp=is_full_pp,
log_level=log_level)
log_level=log_level,
)
self._use_const_ave_pp = self._pp.get_constant_averaged_interaction()
self._read_pp = read_pp
@ -577,20 +607,27 @@ class Conductivity_RTA(Conductivity):
self._num_ignored_phonon_modes[j, k] += 1
continue
old_settings = np.seterr(all='raise')
old_settings = np.seterr(all="raise")
try:
self._mode_kappa[j, k, i, ll] = (
self._gv_sum2[i, ll] * cv[k, ll] /
(g_sum[ll] * 2) * self._conversion_factor)
self._gv_sum2[i, ll]
* cv[k, ll]
/ (g_sum[ll] * 2)
* self._conversion_factor
)
except FloatingPointError:
# supposed that g is almost 0 and |gv|=0
pass
except Exception:
print("=" * 26 + " Warning " + "=" * 26)
print(" Unexpected physical condition of ph-ph "
"interaction calculation was found.")
print(" g=%f at gp=%d, band=%d, freq=%f" %
(g_sum[ll], gp, ll + 1, frequencies[ll]))
print(
" Unexpected physical condition of ph-ph "
"interaction calculation was found."
)
print(
" g=%f at gp=%d, band=%d, freq=%f"
% (g_sum[ll], gp, ll + 1, frequencies[ll])
)
print("=" * 61)
np.seterr(**old_settings)
@ -635,12 +672,14 @@ class Conductivity_RTA(Conductivity):
if self._log_level:
print("Number of triplets: %d" % num_triplets)
if (self._is_full_pp or
self._read_pp or
self._store_pp or
self._use_ave_pp or
self._use_const_ave_pp or
self._is_gamma_detail): # noqa E129
if (
self._is_full_pp
or self._read_pp
or self._store_pp
or self._use_ave_pp
or self._use_const_ave_pp
or self._is_gamma_detail
):
self._set_gamma_at_sigmas(i)
else: # can save memory space
self._set_gamma_at_sigmas_lowmem(i)
@ -656,43 +695,48 @@ class Conductivity_RTA(Conductivity):
num_band0 = len(self._pp.band_indices)
num_grid_points = len(self._grid_points)
num_temp = len(self._temperatures)
self._kappa = np.zeros((len(self._sigmas), num_temp, 6),
order='C', dtype='double')
self._mode_kappa = np.zeros((len(self._sigmas),
num_temp,
num_grid_points,
num_band0,
6),
order='C', dtype='double')
self._kappa = np.zeros(
(len(self._sigmas), num_temp, 6), order="C", dtype="double"
)
self._mode_kappa = np.zeros(
(len(self._sigmas), num_temp, num_grid_points, num_band0, 6),
order="C",
dtype="double",
)
if not self._read_gamma:
self._gamma = np.zeros((len(self._sigmas),
num_temp,
num_grid_points,
num_band0),
order='C', dtype='double')
self._gamma = np.zeros(
(len(self._sigmas), num_temp, num_grid_points, num_band0),
order="C",
dtype="double",
)
if self._is_gamma_detail or self._is_N_U:
self._gamma_N = np.zeros_like(self._gamma)
self._gamma_U = np.zeros_like(self._gamma)
self._gv = np.zeros((num_grid_points, num_band0, 3),
order='C', dtype='double')
self._gv_sum2 = np.zeros((num_grid_points, num_band0, 6),
order='C', dtype='double')
self._cv = np.zeros((num_temp, num_grid_points, num_band0),
order='C', dtype='double')
self._gv = np.zeros((num_grid_points, num_band0, 3), order="C", dtype="double")
self._gv_sum2 = np.zeros(
(num_grid_points, num_band0, 6), order="C", dtype="double"
)
self._cv = np.zeros(
(num_temp, num_grid_points, num_band0), order="C", dtype="double"
)
if self._isotope is not None:
self._gamma_iso = np.zeros(
(len(self._sigmas), num_grid_points, num_band0),
order='C', dtype='double')
order="C",
dtype="double",
)
if self._is_full_pp or self._use_ave_pp or self._use_const_ave_pp:
self._averaged_pp_interaction = np.zeros(
(num_grid_points, num_band0),
order='C', dtype='double')
(num_grid_points, num_band0), order="C", dtype="double"
)
self._num_ignored_phonon_modes = np.zeros(
(len(self._sigmas), num_temp), order='C', dtype='intc')
(len(self._sigmas), num_temp), order="C", dtype="intc"
)
self._collision = ImagSelfEnergy(
self._pp,
with_detail=(self._is_gamma_detail or self._is_N_U),
unit_conversion=self._gamma_unit_conversion)
unit_conversion=self._gamma_unit_conversion,
)
def _set_gamma_at_sigmas(self, i):
for j, sigma in enumerate(self._sigmas):
@ -718,25 +762,30 @@ class Conductivity_RTA(Conductivity):
sigma=sigma,
sigma_cutoff=self._sigma_cutoff,
filename=self._pp_filename,
verbose=(self._log_level > 0))
verbose=(self._log_level > 0),
)
_, g_zero = self._collision.get_integration_weights()
if self._log_level:
if len(self._sigmas) > 1:
print("Multiple sigmas or mixing smearing and "
"tetrahedron method is not supported.")
print(
"Multiple sigmas or mixing smearing and "
"tetrahedron method is not supported."
)
if _g_zero is not None and (_g_zero != g_zero).any():
raise ValueError("Inconsistency found in g_zero.")
self._collision.set_interaction_strength(pp)
elif self._use_ave_pp:
self._collision.set_averaged_pp_interaction(
self._averaged_pp_interaction[i])
self._averaged_pp_interaction[i]
)
elif self._use_const_ave_pp:
if self._log_level:
print("Constant ph-ph interaction of %6.3e is used." %
self._pp.get_constant_averaged_interaction())
print(
"Constant ph-ph interaction of %6.3e is used."
% self._pp.get_constant_averaged_interaction()
)
self._collision.run_interaction()
self._averaged_pp_interaction[i] = (
self._pp.get_averaged_interaction())
self._averaged_pp_interaction[i] = self._pp.get_averaged_interaction()
elif j != 0 and (self._is_full_pp or self._sigma_cutoff is None):
if self._log_level:
print("Existing ph-ph interaction is used.")
@ -745,8 +794,9 @@ class Conductivity_RTA(Conductivity):
print("Calculating ph-ph interaction...")
self._collision.run_interaction(is_full_pp=self._is_full_pp)
if self._is_full_pp:
self._averaged_pp_interaction[i] = (
self._pp.get_averaged_interaction())
self._averaged_pp_interaction[
i
] = self._pp.get_averaged_interaction()
# Number of triplets depends on q-point.
# So this is allocated each time.
@ -754,7 +804,9 @@ class Conductivity_RTA(Conductivity):
num_temp = len(self._temperatures)
self._gamma_detail_at_q = np.empty(
((num_temp,) + self._pp.get_interaction_strength().shape),
dtype='double', order='C')
dtype="double",
order="C",
)
self._gamma_detail_at_q[:] = 0
if self._log_level:
@ -768,8 +820,9 @@ class Conductivity_RTA(Conductivity):
self._gamma_N[j, k, i] = g_N
self._gamma_U[j, k, i] = g_U
if self._is_gamma_detail:
self._gamma_detail_at_q[k] = (
self._collision.get_detailed_imag_self_energy())
self._gamma_detail_at_q[
k
] = self._collision.get_detailed_imag_self_energy()
def _set_gamma_at_sigmas_lowmem(self, i):
"""Calculate gamma without storing ph-ph interaction strength.
@ -782,11 +835,13 @@ class Conductivity_RTA(Conductivity):
"""
band_indices = self._pp.band_indices
(svecs,
(
svecs,
multi,
p2s,
s2p,
masses) = self._pp.get_primitive_and_supercell_correspondence()
masses,
) = self._pp.get_primitive_and_supercell_correspondence()
fc3 = self._pp.fc3
triplets_at_q, weights_at_q, _, _ = self._pp.get_triplets_at_q()
symmetrize_fc3_q = 0
@ -798,20 +853,27 @@ class Conductivity_RTA(Conductivity):
for j, sigma in enumerate(self._sigmas):
self._collision.set_sigma(sigma)
if self._is_N_U:
collisions = np.zeros((2, len(self._temperatures),
len(band_indices)),
dtype='double', order='C')
collisions = np.zeros(
(2, len(self._temperatures), len(band_indices)),
dtype="double",
order="C",
)
else:
collisions = np.zeros((len(self._temperatures),
len(band_indices)),
dtype='double', order='C')
collisions = np.zeros(
(len(self._temperatures), len(band_indices)),
dtype="double",
order="C",
)
import phono3py._phono3py as phono3c
if sigma is None:
phono3c.pp_collision(
collisions,
np.array(
np.dot(thm.get_tetrahedra(), self._pp.bz_grid.P.T),
dtype='int_', order='C'),
dtype="int_",
order="C",
),
self._frequencies,
self._eigenvectors,
triplets_at_q,
@ -831,13 +893,15 @@ class Conductivity_RTA(Conductivity):
self._temperatures,
self._is_N_U * 1,
symmetrize_fc3_q,
self._pp.cutoff_frequency)
self._pp.cutoff_frequency,
)
else:
if self._sigma_cutoff is None:
sigma_cutoff = -1
else:
sigma_cutoff = float(self._sigma_cutoff)
phono3c.pp_collision_with_sigma(collisions,
phono3c.pp_collision_with_sigma(
collisions,
sigma,
sigma_cutoff,
self._frequencies,
@ -857,7 +921,8 @@ class Conductivity_RTA(Conductivity):
self._temperatures,
self._is_N_U * 1,
symmetrize_fc3_q,
self._pp.cutoff_frequency)
self._pp.cutoff_frequency,
)
col_unit_conv = self._collision.get_unit_conversion_factor()
pp_unit_conv = self._pp.get_unit_conversion_factor()
if self._is_N_U:
@ -870,16 +935,19 @@ class Conductivity_RTA(Conductivity):
self._gamma[j, k, i, :] = average_by_degeneracy(
col[k] * col_unit_conv * pp_unit_conv,
band_indices,
self._frequencies[self._grid_points[i]])
self._frequencies[self._grid_points[i]],
)
if self._is_N_U:
self._gamma_N[j, k, i, :] = average_by_degeneracy(
col_N[k] * col_unit_conv * pp_unit_conv,
band_indices,
self._frequencies[self._grid_points[i]])
self._frequencies[self._grid_points[i]],
)
self._gamma_U[j, k, i, :] = average_by_degeneracy(
col_U[k] * col_unit_conv * pp_unit_conv,
band_indices,
self._frequencies[self._grid_points[i]])
self._frequencies[self._grid_points[i]],
)
def _show_log(self, q, i):
gp = self._grid_points[i]
@ -901,36 +969,42 @@ class Conductivity_RTA(Conductivity):
def _show_log_values(self, frequencies, gv, ave_pp):
if self._is_full_pp or self._use_ave_pp or self._use_const_ave_pp:
for f, v, pp in zip(frequencies, gv, ave_pp):
print("%8.3f (%8.3f %8.3f %8.3f) %8.3f %11.3e" %
(f, v[0], v[1], v[2], np.linalg.norm(v), pp))
print(
"%8.3f (%8.3f %8.3f %8.3f) %8.3f %11.3e"
% (f, v[0], v[1], v[2], np.linalg.norm(v), pp)
)
else:
for f, v in zip(frequencies, gv):
print("%8.3f (%8.3f %8.3f %8.3f) %8.3f" %
(f, v[0], v[1], v[2], np.linalg.norm(v)))
print(
"%8.3f (%8.3f %8.3f %8.3f) %8.3f"
% (f, v[0], v[1], v[2], np.linalg.norm(v))
)
def _show_log_values_on_kstar(self, frequencies, gv, ave_pp, gp, q):
rotation_map = get_grid_points_by_rotations(gp, self._pp.bz_grid)
for i, j in enumerate(np.unique(rotation_map)):
for k, (rot, rot_c) in enumerate(zip(
self._point_operations, self._rotations_cartesian)):
for k, (rot, rot_c) in enumerate(
zip(self._point_operations, self._rotations_cartesian)
):
if rotation_map[k] != j:
continue
print(" k*%-2d (%5.2f %5.2f %5.2f)" %
((i + 1,) + tuple(np.dot(rot, q))))
if (self._is_full_pp or
self._use_ave_pp or
self._use_const_ave_pp): # noqa E129
for f, v, pp in zip(frequencies,
np.dot(rot_c, gv.T).T,
ave_pp):
print("%8.3f (%8.3f %8.3f %8.3f) %8.3f %11.3e" %
(f, v[0], v[1], v[2], np.linalg.norm(v), pp))
print(
" k*%-2d (%5.2f %5.2f %5.2f)" % ((i + 1,) + tuple(np.dot(rot, q)))
)
if self._is_full_pp or self._use_ave_pp or self._use_const_ave_pp:
for f, v, pp in zip(frequencies, np.dot(rot_c, gv.T).T, ave_pp):
print(
"%8.3f (%8.3f %8.3f %8.3f) %8.3f %11.3e"
% (f, v[0], v[1], v[2], np.linalg.norm(v), pp)
)
else:
for f, v in zip(frequencies, np.dot(rot_c, gv.T).T):
print("%8.3f (%8.3f %8.3f %8.3f) %8.3f" %
(f, v[0], v[1], v[2], np.linalg.norm(v)))
print('')
print(
"%8.3f (%8.3f %8.3f %8.3f) %8.3f"
% (f, v[0], v[1], v[2], np.linalg.norm(v))
)
print("")
def _show_log_value_names(self):
if self._is_full_pp or self._use_ave_pp or self._use_const_ave_pp: