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

Minor refactorings

This commit is contained in:
Atsushi Togo 2021-10-22 10:04:47 +09:00
parent fa922b30a2
commit 923487c2bc
7 changed files with 1714 additions and 1683 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
phono3py/api_phono3py.py
View File

@ -2375,8 +2375,8 @@ class Phono3py(object):
self._primitive_symmetry = Symmetry(
self._primitive, self._symprec, self._is_symmetry
)
if len(self._symmetry.get_pointgroup_operations()) != len(
self._primitive_symmetry.get_pointgroup_operations()
if len(self._symmetry.pointgroup_operations) != len(
self._primitive_symmetry.pointgroup_operations
): # noqa E129
print(
"Warning: point group symmetries of supercell and primitive"

105
phono3py/phonon3/conductivity.py
View File

@ -35,6 +35,7 @@
import warnings
import textwrap
from typing import Optional, List
import numpy as np
from phonopy.phonon.group_velocity import GroupVelocity
from phonopy.harmonic.force_constants import similarity_transformation
@ -42,6 +43,7 @@ from phonopy.phonon.thermal_properties import mode_cv as get_mode_cv
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.phonon3.interaction import Interaction
from phono3py.other.isotope import Isotope
from phono3py.phonon.grid import get_ir_grid_points, get_grid_points_by_rotations
@ -50,49 +52,15 @@ unit_to_WmK = (
) # 2pi comes from definition of lifetime.
def all_bands_exist(interaction):
"""Return if all bands are selected or not."""
band_indices = interaction.band_indices
num_band = len(interaction.primitive) * 3
if len(band_indices) == num_band:
if (band_indices - np.arange(num_band) == 0).all():
return True
return False
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
sigma_cutoff = conductivity.sigma_cutoff_width
mesh = conductivity.mesh_numbers
triplets, weights, _, _ = pp.get_triplets_at_q()
all_triplets = get_all_triplets(grid_point, pp.bz_grid)
if len(sigmas) > 1:
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,
pp=pp.interaction_strength,
g_zero=pp.zero_value_positions,
grid_point=grid_point,
triplet=triplets,
weight=weights,
triplet_all=all_triplets,
sigma=sigmas[-1],
sigma_cutoff=sigma_cutoff,
filename=filename,
compression=compression,
)
class ConductivityBase:
"""Base class of Conductivity classes."""
def __init__(
self, interaction, grid_points=None, is_kappa_star=True, gv_delta_q=None
self,
interaction: Interaction,
grid_points=None,
is_kappa_star=True,
gv_delta_q=None,
):
"""Init method.
@ -112,7 +80,7 @@ class ConductivityBase:
See phonopy's GroupVelocity class. Default is None.
"""
self._pp = interaction
self._pp: Interaction = interaction
self._is_kappa_star = is_kappa_star
self._rotations_cartesian = None
@ -184,10 +152,10 @@ class Conductivity(ConductivityBase):
def __init__(
self,
interaction,
interaction: Interaction,
grid_points=None,
temperatures=None,
sigmas=None,
sigmas: Optional[List] = None,
sigma_cutoff=None,
is_isotope=False,
mass_variances=None,
@ -198,15 +166,15 @@ class Conductivity(ConductivityBase):
log_level=0,
):
"""Init method."""
self._pp = None
self._pp: Interaction
self._qpoints: np.ndarray
self._gv_obj: GroupVelocity
self._is_kappa_star = None
self._grid_points = None
self._grid_weights = None
self._ir_grid_points = None
self._ir_grid_weights = None
self._qpoints = None
self._point_operations = None
self._gv_obj = None
super().__init__(
interaction,
@ -218,7 +186,8 @@ class Conductivity(ConductivityBase):
self._is_full_pp = is_full_pp
self._log_level = log_level
self._grid_point_count = 0
self._grid_point_count: int = 0
self._sigmas: Optional[List]
if sigmas is None:
self._sigmas = []
else:
@ -784,3 +753,47 @@ class Conductivity(ConductivityBase):
)
% tuple(self._mass_variances)
)
def all_bands_exist(interaction: Interaction):
"""Return if all bands are selected or not."""
band_indices = interaction.band_indices
num_band = len(interaction.primitive) * 3
if len(band_indices) == num_band:
if (band_indices - np.arange(num_band) == 0).all():
return True
return False
def write_pp(
conductivity: Conductivity,
pp: Interaction,
i,
filename=None,
compression="gzip",
):
"""Write ph-ph interaction strength in hdf5 file."""
grid_point = conductivity.grid_points[i]
sigmas = conductivity.sigmas
sigma_cutoff = conductivity.sigma_cutoff_width
mesh = conductivity.mesh_numbers
triplets, weights, _, _ = pp.get_triplets_at_q()
all_triplets = get_all_triplets(grid_point, pp.bz_grid)
if len(sigmas) > 1:
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,
pp=pp.interaction_strength,
g_zero=pp.zero_value_positions,
grid_point=grid_point,
triplet=triplets,
weight=weights,
triplet_all=all_triplets,
sigma=sigmas[-1],
sigma_cutoff=sigma_cutoff,
filename=filename,
compression=compression,
)

2349
phono3py/phonon3/conductivity_LBTE.py
View File

File diff suppressed because it is too large Load Diff

901
phono3py/phonon3/conductivity_RTA.py
View File

@ -36,6 +36,7 @@
import sys
import numpy as np
from phonopy.structure.tetrahedron_method import TetrahedronMethod
from phonopy.phonon.group_velocity import GroupVelocity
from phono3py.file_IO import (
write_kappa_to_hdf5,
read_gamma_from_hdf5,
@ -44,460 +45,18 @@ from phono3py.file_IO import (
)
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.interaction import Interaction
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
def get_thermal_conductivity_RTA(
interaction,
temperatures=None,
sigmas=None,
sigma_cutoff=None,
mass_variances=None,
grid_points=None,
is_isotope=False,
boundary_mfp=None, # in micrometre
use_ave_pp=False,
is_kappa_star=True,
gv_delta_q=None,
is_full_pp=False,
write_gamma=False,
read_gamma=False,
is_N_U=False,
write_kappa=False,
write_pp=False,
read_pp=False,
write_gamma_detail=False,
compression="gzip",
input_filename=None,
output_filename=None,
log_level=0,
):
"""Run RTA thermal conductivity calculation."""
if temperatures is None:
_temperatures = np.arange(0, 1001, 10, dtype="double")
else:
_temperatures = temperatures
if log_level:
print(
"-------------------- Lattice thermal conducitivity (RTA) "
"--------------------"
)
br = Conductivity_RTA(
interaction,
grid_points=grid_points,
temperatures=_temperatures,
sigmas=sigmas,
sigma_cutoff=sigma_cutoff,
is_isotope=is_isotope,
mass_variances=mass_variances,
boundary_mfp=boundary_mfp,
use_ave_pp=use_ave_pp,
is_kappa_star=is_kappa_star,
gv_delta_q=gv_delta_q,
is_full_pp=is_full_pp,
read_pp=read_pp,
store_pp=write_pp,
pp_filename=input_filename,
is_N_U=is_N_U,
is_gamma_detail=write_gamma_detail,
log_level=log_level,
)
if read_gamma:
if not _set_gamma_from_file(br, filename=input_filename):
print("Reading collisions failed.")
return False
for i in br:
if write_pp:
_write_pp(
br, interaction, i, compression=compression, filename=output_filename
)
if write_gamma:
_write_gamma(
br,
interaction,
i,
compression=compression,
filename=output_filename,
verbose=log_level,
)
if write_gamma_detail:
_write_gamma_detail(
br,
interaction,
i,
compression=compression,
filename=output_filename,
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,
interaction.primitive.volume,
compression=compression,
filename=output_filename,
log_level=log_level,
)
return br
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()
grid_points = br.get_grid_points()
gp = grid_points[i]
sigmas = br.get_sigmas()
sigma_cutoff = br.get_sigma_cutoff_width()
triplets, weights, _, _ = interaction.get_triplets_at_q()
all_triplets = get_all_triplets(gp, interaction.bz_grid)
if all_bands_exist(interaction):
for j, sigma in enumerate(sigmas):
write_gamma_detail_to_hdf5(
temperatures,
mesh,
gamma_detail=gamma_detail,
grid_point=gp,
triplet=triplets,
weight=weights,
triplet_all=all_triplets,
sigma=sigma,
sigma_cutoff=sigma_cutoff,
compression=compression,
filename=filename,
verbose=verbose,
)
else:
for j, sigma in enumerate(sigmas):
for k, bi in enumerate(interaction.get_band_indices()):
write_gamma_detail_to_hdf5(
temperatures,
mesh,
gamma_detail=gamma_detail[:, :, k, :, :],
grid_point=gp,
triplet=triplets,
weight=weights,
band_index=bi,
sigma=sigma,
sigma_cutoff=sigma_cutoff,
compression=compression,
filename=filename,
verbose=verbose,
)
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
gv_by_gv = br.gv_by_gv
mode_heat_capacities = br.mode_heat_capacities
ave_pp = br.averaged_pp_interaction
mesh = br.mesh_numbers
temperatures = br.temperatures
gamma = br.gamma
gamma_isotope = br.gamma_isotope
sigmas = br.sigmas
sigma_cutoff = br.sigma_cutoff_width
volume = interaction.primitive.volume
gamma_N, gamma_U = br.get_gamma_N_U()
gp = grid_points[i]
if all_bands_exist(interaction):
if ave_pp is None:
ave_pp_i = None
else:
ave_pp_i = ave_pp[i]
frequencies = interaction.get_phonons()[0][gp]
for j, sigma in enumerate(sigmas):
if gamma_isotope is not None:
gamma_isotope_at_sigma = gamma_isotope[j, i]
else:
gamma_isotope_at_sigma = None
if gamma_N is None:
gamma_N_at_sigma = None
else:
gamma_N_at_sigma = gamma_N[j, :, i]
if gamma_U is None:
gamma_U_at_sigma = None
else:
gamma_U_at_sigma = gamma_U[j, :, i]
write_kappa_to_hdf5(
temperatures,
mesh,
frequency=frequencies,
group_velocity=group_velocities[i],
gv_by_gv=gv_by_gv[i],
heat_capacity=mode_heat_capacities[:, i],
gamma=gamma[j, :, i],
gamma_isotope=gamma_isotope_at_sigma,
gamma_N=gamma_N_at_sigma,
gamma_U=gamma_U_at_sigma,
averaged_pp_interaction=ave_pp_i,
grid_point=gp,
sigma=sigma,
sigma_cutoff=sigma_cutoff,
kappa_unit_conversion=unit_to_WmK / volume,
compression=compression,
filename=filename,
verbose=verbose,
)
else:
for j, sigma in enumerate(sigmas):
for k, bi in enumerate(interaction.band_indices):
if ave_pp is None:
ave_pp_ik = None
else:
ave_pp_ik = ave_pp[i, k]
frequencies = interaction.get_phonons()[0][gp, bi]
if gamma_isotope is not None:
gamma_isotope_at_sigma = gamma_isotope[j, i, k]
else:
gamma_isotope_at_sigma = None
if gamma_N is None:
gamma_N_at_sigma = None
else:
gamma_N_at_sigma = gamma_N[j, :, i, k]
if gamma_U is None:
gamma_U_at_sigma = None
else:
gamma_U_at_sigma = gamma_U[j, :, i, k]
write_kappa_to_hdf5(
temperatures,
mesh,
frequency=frequencies,
group_velocity=group_velocities[i, k],
gv_by_gv=gv_by_gv[i, k],
heat_capacity=mode_heat_capacities[:, i, k],
gamma=gamma[j, :, i, k],
gamma_isotope=gamma_isotope_at_sigma,
gamma_N=gamma_N_at_sigma,
gamma_U=gamma_U_at_sigma,
averaged_pp_interaction=ave_pp_ik,
grid_point=gp,
band_index=bi,
sigma=sigma,
sigma_cutoff=sigma_cutoff,
kappa_unit_conversion=unit_to_WmK / volume,
compression=compression,
filename=filename,
verbose=verbose,
)
def _show_kappa(br, log_level):
temperatures = br.temperatures
sigmas = br.sigmas
kappa = br.kappa
num_ignored_phonon_modes = br.get_number_of_ignored_phonon_modes()
num_band = br.frequencies.shape[1]
num_phonon_modes = br.get_number_of_sampling_grid_points() * num_band
for i, sigma in enumerate(sigmas):
text = "----------- Thermal conductivity (W/m-k) "
if sigma:
text += "for sigma=%s -----------" % sigma
else:
text += "with tetrahedron method -----------"
print(text)
if log_level > 1:
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)
)
)
else:
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("")
def _write_kappa(br, volume, compression="gzip", filename=None, log_level=0):
temperatures = br.temperatures
sigmas = br.sigmas
sigma_cutoff = br.sigma_cutoff_width
gamma = br.gamma
gamma_isotope = br.gamma_isotope
gamma_N, gamma_U = br.get_gamma_N_U()
mesh = br.mesh_numbers
frequencies = br.frequencies
gv = br.group_velocities
gv_by_gv = br.gv_by_gv
mode_cv = br.mode_heat_capacities
ave_pp = br.averaged_pp_interaction
qpoints = br.qpoints
weights = br.grid_weights
kappa = br.kappa
mode_kappa = br.mode_kappa
for i, sigma in enumerate(sigmas):
kappa_at_sigma = kappa[i]
if gamma_isotope is not None:
gamma_isotope_at_sigma = gamma_isotope[i]
else:
gamma_isotope_at_sigma = None
if gamma_N is None:
gamma_N_at_sigma = None
else:
gamma_N_at_sigma = gamma_N[i]
if gamma_U is None:
gamma_U_at_sigma = None
else:
gamma_U_at_sigma = gamma_U[i]
write_kappa_to_hdf5(
temperatures,
mesh,
frequency=frequencies,
group_velocity=gv,
gv_by_gv=gv_by_gv,
heat_capacity=mode_cv,
kappa=kappa_at_sigma,
mode_kappa=mode_kappa[i],
gamma=gamma[i],
gamma_isotope=gamma_isotope_at_sigma,
gamma_N=gamma_N_at_sigma,
gamma_U=gamma_U_at_sigma,
averaged_pp_interaction=ave_pp,
qpoint=qpoints,
weight=weights,
sigma=sigma,
sigma_cutoff=sigma_cutoff,
kappa_unit_conversion=unit_to_WmK / volume,
compression=compression,
filename=filename,
verbose=log_level,
)
def _set_gamma_from_file(br, filename=None, verbose=True):
"""Read kappa-*.hdf5 files for thermal conductivity calculation."""
sigmas = br.get_sigmas()
sigma_cutoff = br.get_sigma_cutoff_width()
mesh = br.get_mesh_numbers()
grid_points = br.get_grid_points()
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_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")
is_gamma_N_U_in = False
is_ave_pp_in = False
read_succeeded = True
for j, sigma in enumerate(sigmas):
data = read_gamma_from_hdf5(
mesh,
sigma=sigma,
sigma_cutoff=sigma_cutoff,
filename=filename,
verbose=verbose,
)
if 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:
is_ave_pp_in = True
ave_pp[:] = data["ave_pp"]
else:
for i, gp in enumerate(grid_points):
data_gp = read_gamma_from_hdf5(
mesh,
grid_point=gp,
sigma=sigma,
sigma_cutoff=sigma_cutoff,
filename=filename,
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:
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:
is_ave_pp_in = True
ave_pp[i] = data_gp["ave_pp"]
else:
for bi in range(num_band):
data_band = read_gamma_from_hdf5(
mesh,
grid_point=gp,
band_index=bi,
sigma=sigma,
sigma_cutoff=sigma_cutoff,
filename=filename,
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:
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:
is_ave_pp_in = True
ave_pp[i, bi] = data_band["ave_pp"]
else:
read_succeeded = False
if read_succeeded:
br.set_gamma(gamma)
if is_ave_pp_in:
br.set_averaged_pp_interaction(ave_pp)
if is_gamma_N_U_in:
br.set_gamma_N_U(gamma_N, gamma_U)
return True
else:
return False
class Conductivity_RTA(Conductivity):
"""Lattice thermal conductivity calculation with RTA."""
def __init__(
self,
interaction,
interaction: Interaction,
grid_points=None,
temperatures=None,
sigmas=None,
@ -518,8 +77,10 @@ class Conductivity_RTA(Conductivity):
log_level=0,
):
"""Init method."""
self._pp = None
self._gv_obj = None
self._pp: Interaction
self._grid_point_count: int
self._num_sampling_grid_points: int
self._gv_obj: GroupVelocity
self._temperatures = None
self._sigmas = None
self._sigma_cutoff = None
@ -553,17 +114,14 @@ class Conductivity_RTA(Conductivity):
self._use_const_ave_pp = None
self._averaged_pp_interaction = None
self._num_ignored_phonon_modes = None
self._num_sampling_grid_points = None
self._conversion_factor = None
self._is_isotope = None
self._isotope = None
self._mass_variances = None
self._grid_point_count = None
Conductivity.__init__(
self,
super().__init__(
interaction,
grid_points=grid_points,
temperatures=temperatures,
@ -1010,3 +568,446 @@ class Conductivity_RTA(Conductivity):
else:
text += " (dq=%3.1e)" % self._gv_obj.q_length
print(text)
def get_thermal_conductivity_RTA(
interaction: Interaction,
temperatures=None,
sigmas=None,
sigma_cutoff=None,
mass_variances=None,
grid_points=None,
is_isotope=False,
boundary_mfp=None, # in micrometre
use_ave_pp=False,
is_kappa_star=True,
gv_delta_q=None,
is_full_pp=False,
write_gamma=False,
read_gamma=False,
is_N_U=False,
write_kappa=False,
write_pp=False,
read_pp=False,
write_gamma_detail=False,
compression="gzip",
input_filename=None,
output_filename=None,
log_level=0,
):
"""Run RTA thermal conductivity calculation."""
if temperatures is None:
_temperatures = np.arange(0, 1001, 10, dtype="double")
else:
_temperatures = temperatures
if log_level:
print(
"-------------------- Lattice thermal conducitivity (RTA) "
"--------------------"
)
br = Conductivity_RTA(
interaction,
grid_points=grid_points,
temperatures=_temperatures,
sigmas=sigmas,
sigma_cutoff=sigma_cutoff,
is_isotope=is_isotope,
mass_variances=mass_variances,
boundary_mfp=boundary_mfp,
use_ave_pp=use_ave_pp,
is_kappa_star=is_kappa_star,
gv_delta_q=gv_delta_q,
is_full_pp=is_full_pp,
read_pp=read_pp,
store_pp=write_pp,
pp_filename=input_filename,
is_N_U=is_N_U,
is_gamma_detail=write_gamma_detail,
log_level=log_level,
)
if read_gamma:
if not _set_gamma_from_file(br, filename=input_filename):
print("Reading collisions failed.")
return False
for i in br:
if write_pp:
_write_pp(
br, interaction, i, compression=compression, filename=output_filename
)
if write_gamma:
_write_gamma(
br,
interaction,
i,
compression=compression,
filename=output_filename,
verbose=log_level,
)
if write_gamma_detail:
_write_gamma_detail(
br,
interaction,
i,
compression=compression,
filename=output_filename,
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,
interaction.primitive.volume,
compression=compression,
filename=output_filename,
log_level=log_level,
)
return br
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()
grid_points = br.get_grid_points()
gp = grid_points[i]
sigmas = br.get_sigmas()
sigma_cutoff = br.get_sigma_cutoff_width()
triplets, weights, _, _ = interaction.get_triplets_at_q()
all_triplets = get_all_triplets(gp, interaction.bz_grid)
if all_bands_exist(interaction):
for j, sigma in enumerate(sigmas):
write_gamma_detail_to_hdf5(
temperatures,
mesh,
gamma_detail=gamma_detail,
grid_point=gp,
triplet=triplets,
weight=weights,
triplet_all=all_triplets,
sigma=sigma,
sigma_cutoff=sigma_cutoff,
compression=compression,
filename=filename,
verbose=verbose,
)
else:
for j, sigma in enumerate(sigmas):
for k, bi in enumerate(interaction.get_band_indices()):
write_gamma_detail_to_hdf5(
temperatures,
mesh,
gamma_detail=gamma_detail[:, :, k, :, :],
grid_point=gp,
triplet=triplets,
weight=weights,
band_index=bi,
sigma=sigma,
sigma_cutoff=sigma_cutoff,
compression=compression,
filename=filename,
verbose=verbose,
)
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
gv_by_gv = br.gv_by_gv
mode_heat_capacities = br.mode_heat_capacities
ave_pp = br.averaged_pp_interaction
mesh = br.mesh_numbers
temperatures = br.temperatures
gamma = br.gamma
gamma_isotope = br.gamma_isotope
sigmas = br.sigmas
sigma_cutoff = br.sigma_cutoff_width
volume = interaction.primitive.volume
gamma_N, gamma_U = br.get_gamma_N_U()
gp = grid_points[i]
if all_bands_exist(interaction):
if ave_pp is None:
ave_pp_i = None
else:
ave_pp_i = ave_pp[i]
frequencies = interaction.get_phonons()[0][gp]
for j, sigma in enumerate(sigmas):
if gamma_isotope is not None:
gamma_isotope_at_sigma = gamma_isotope[j, i]
else:
gamma_isotope_at_sigma = None
if gamma_N is None:
gamma_N_at_sigma = None
else:
gamma_N_at_sigma = gamma_N[j, :, i]
if gamma_U is None:
gamma_U_at_sigma = None
else:
gamma_U_at_sigma = gamma_U[j, :, i]
write_kappa_to_hdf5(
temperatures,
mesh,
frequency=frequencies,
group_velocity=group_velocities[i],
gv_by_gv=gv_by_gv[i],
heat_capacity=mode_heat_capacities[:, i],
gamma=gamma[j, :, i],
gamma_isotope=gamma_isotope_at_sigma,
gamma_N=gamma_N_at_sigma,
gamma_U=gamma_U_at_sigma,
averaged_pp_interaction=ave_pp_i,
grid_point=gp,
sigma=sigma,
sigma_cutoff=sigma_cutoff,
kappa_unit_conversion=unit_to_WmK / volume,
compression=compression,
filename=filename,
verbose=verbose,
)
else:
for j, sigma in enumerate(sigmas):
for k, bi in enumerate(interaction.band_indices):
if ave_pp is None:
ave_pp_ik = None
else:
ave_pp_ik = ave_pp[i, k]
frequencies = interaction.get_phonons()[0][gp, bi]
if gamma_isotope is not None:
gamma_isotope_at_sigma = gamma_isotope[j, i, k]
else:
gamma_isotope_at_sigma = None
if gamma_N is None:
gamma_N_at_sigma = None
else:
gamma_N_at_sigma = gamma_N[j, :, i, k]
if gamma_U is None:
gamma_U_at_sigma = None
else:
gamma_U_at_sigma = gamma_U[j, :, i, k]
write_kappa_to_hdf5(
temperatures,
mesh,
frequency=frequencies,
group_velocity=group_velocities[i, k],
gv_by_gv=gv_by_gv[i, k],
heat_capacity=mode_heat_capacities[:, i, k],
gamma=gamma[j, :, i, k],
gamma_isotope=gamma_isotope_at_sigma,
gamma_N=gamma_N_at_sigma,
gamma_U=gamma_U_at_sigma,
averaged_pp_interaction=ave_pp_ik,
grid_point=gp,
band_index=bi,
sigma=sigma,
sigma_cutoff=sigma_cutoff,
kappa_unit_conversion=unit_to_WmK / volume,
compression=compression,
filename=filename,
verbose=verbose,
)
def _show_kappa(br, log_level):
temperatures = br.temperatures
sigmas = br.sigmas
kappa = br.kappa
num_ignored_phonon_modes = br.get_number_of_ignored_phonon_modes()
num_band = br.frequencies.shape[1]
num_phonon_modes = br.get_number_of_sampling_grid_points() * num_band
for i, sigma in enumerate(sigmas):
text = "----------- Thermal conductivity (W/m-k) "
if sigma:
text += "for sigma=%s -----------" % sigma
else:
text += "with tetrahedron method -----------"
print(text)
if log_level > 1:
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)
)
)
else:
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("")
def _write_kappa(br, volume, compression="gzip", filename=None, log_level=0):
temperatures = br.temperatures
sigmas = br.sigmas
sigma_cutoff = br.sigma_cutoff_width
gamma = br.gamma
gamma_isotope = br.gamma_isotope
gamma_N, gamma_U = br.get_gamma_N_U()
mesh = br.mesh_numbers
frequencies = br.frequencies
gv = br.group_velocities
gv_by_gv = br.gv_by_gv
mode_cv = br.mode_heat_capacities
ave_pp = br.averaged_pp_interaction
qpoints = br.qpoints
weights = br.grid_weights
kappa = br.kappa
mode_kappa = br.mode_kappa
for i, sigma in enumerate(sigmas):
kappa_at_sigma = kappa[i]
if gamma_isotope is not None:
gamma_isotope_at_sigma = gamma_isotope[i]
else:
gamma_isotope_at_sigma = None
if gamma_N is None:
gamma_N_at_sigma = None
else:
gamma_N_at_sigma = gamma_N[i]
if gamma_U is None:
gamma_U_at_sigma = None
else:
gamma_U_at_sigma = gamma_U[i]
write_kappa_to_hdf5(
temperatures,
mesh,
frequency=frequencies,
group_velocity=gv,
gv_by_gv=gv_by_gv,
heat_capacity=mode_cv,
kappa=kappa_at_sigma,
mode_kappa=mode_kappa[i],
gamma=gamma[i],
gamma_isotope=gamma_isotope_at_sigma,
gamma_N=gamma_N_at_sigma,
gamma_U=gamma_U_at_sigma,
averaged_pp_interaction=ave_pp,
qpoint=qpoints,
weight=weights,
sigma=sigma,
sigma_cutoff=sigma_cutoff,
kappa_unit_conversion=unit_to_WmK / volume,
compression=compression,
filename=filename,
verbose=log_level,
)
def _set_gamma_from_file(br, filename=None, verbose=True):
"""Read kappa-*.hdf5 files for thermal conductivity calculation."""
sigmas = br.get_sigmas()
sigma_cutoff = br.get_sigma_cutoff_width()
mesh = br.get_mesh_numbers()
grid_points = br.get_grid_points()
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_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")
is_gamma_N_U_in = False
is_ave_pp_in = False
read_succeeded = True
for j, sigma in enumerate(sigmas):
data = read_gamma_from_hdf5(
mesh,
sigma=sigma,
sigma_cutoff=sigma_cutoff,
filename=filename,
verbose=verbose,
)
if 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:
is_ave_pp_in = True
ave_pp[:] = data["ave_pp"]
else:
for i, gp in enumerate(grid_points):
data_gp = read_gamma_from_hdf5(
mesh,
grid_point=gp,
sigma=sigma,
sigma_cutoff=sigma_cutoff,
filename=filename,
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:
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:
is_ave_pp_in = True
ave_pp[i] = data_gp["ave_pp"]
else:
for bi in range(num_band):
data_band = read_gamma_from_hdf5(
mesh,
grid_point=gp,
band_index=bi,
sigma=sigma,
sigma_cutoff=sigma_cutoff,
filename=filename,
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:
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:
is_ave_pp_in = True
ave_pp[i, bi] = data_band["ave_pp"]
else:
read_succeeded = False
if read_succeeded:
br.set_gamma(gamma)
if is_ave_pp_in:
br.set_averaged_pp_interaction(ave_pp)
if is_gamma_N_U_in:
br.set_gamma_N_U(gamma_N, gamma_U)
return True
else:
return False

13
phono3py/phonon3/displacement_fc3.py
View File

@ -40,11 +40,16 @@ from phonopy.harmonic.displacement import (
get_displacement,
is_minus_displacement,
)
from phonopy.structure.atoms import PhonopyAtoms
from phonopy.structure.symmetry import Symmetry
from phonopy.structure.cells import get_smallest_vectors
def direction_to_displacement(
direction_dataset, displacement_distance, supercell, cutoff_distance=None
direction_dataset,
displacement_distance,
supercell: PhonopyAtoms,
cutoff_distance=None,
):
"""Convert displacement directions to those in Cartesian coordinates.
@ -125,7 +130,7 @@ def direction_to_displacement(
def get_third_order_displacements(
cell, symmetry, is_plusminus="auto", is_diagonal=False
cell: PhonopyAtoms, symmetry: Symmetry, is_plusminus="auto", is_diagonal=False
):
"""Create dispalcement dataset.
@ -277,11 +282,11 @@ def get_least_orbits(atom_index, cell, site_symmetry, symprec=1e-5):
def get_smallest_vector_of_atom_pair(
atom_number_supercell, atom_number_primitive, supercell, symprec
atom_number_supercell, atom_number_primitive, supercell: PhonopyAtoms, symprec
):
"""Return smallest vectors of an atom pair in supercell."""
s_pos = supercell.scaled_positions
svecs, multi = get_smallest_vectors(
svecs, _ = get_smallest_vectors(
supercell.cell,
[s_pos[atom_number_supercell]],
[s_pos[atom_number_primitive]],

23
phono3py/phonon3/fc3.py
View File

@ -50,13 +50,20 @@ from phono3py.phonon3.displacement_fc3 import (
get_bond_symmetry,
get_smallest_vector_of_atom_pair,
)
from phonopy.structure.cells import compute_all_sg_permutations
from phonopy.structure.atoms import PhonopyAtoms
from phonopy.structure.cells import compute_all_sg_permutations, Primitive
from phonopy.structure.symmetry import Symmetry
logger = logging.getLogger(__name__)
def get_fc3(
supercell, primitive, disp_dataset, symmetry, is_compact_fc=False, verbose=False
supercell: PhonopyAtoms,
primitive: Primitive,
disp_dataset,
symmetry: Symmetry,
is_compact_fc=False,
verbose=False,
):
"""Calculate fc3."""
# fc2 has to be full matrix to compute delta-fc2
@ -75,7 +82,7 @@ def get_fc3(
print("Expanding fc3.")
first_disp_atoms = np.unique([x["number"] for x in disp_dataset["first_atoms"]])
rotations = symmetry.get_symmetry_operations()["rotations"]
rotations = symmetry.symmetry_operations["rotations"]
lattice = supercell.cell.T
permutations = symmetry.atomic_permutations
@ -698,14 +705,16 @@ def _third_rank_tensor_rotation_elem(rot, tensor, ll, m, n):
return sum_elems
def _get_fc3_done(supercell, disp_dataset, symmetry, array_shape):
def _get_fc3_done(
supercell: PhonopyAtoms, disp_dataset, symmetry: Symmetry, array_shape
):
num_atom = len(supercell)
fc3_done = np.zeros(array_shape, dtype="byte")
symprec = symmetry.tolerance
lattice = supercell.cell.T
positions = supercell.scaled_positions
rotations = symmetry.get_symmetry_operations()["rotations"]
translations = symmetry.get_symmetry_operations()["translations"]
rotations = symmetry.symmetry_operations["rotations"]
translations = symmetry.symmetry_operations["translations"]
atom_mapping = []
for rot, trans in zip(rotations, translations):
@ -719,7 +728,7 @@ def _get_fc3_done(supercell, disp_dataset, symmetry, array_shape):
first_atom_num = dataset_first_atom["number"]
site_symmetry = symmetry.get_site_symmetry(first_atom_num)
direction = np.dot(
dataset_first_atom["displacement"], np.linalg.inv(supercell.get_cell())
dataset_first_atom["displacement"], np.linalg.inv(supercell.cell)
)
reduced_site_sym = get_reduced_site_symmetry(site_symmetry, direction, symprec)
least_second_atom_nums = []

2
phono3py/phonon3/interaction.py
View File

@ -51,7 +51,7 @@ from phono3py.phonon3.reciprocal_to_normal import ReciprocalToNormal
from phono3py.phonon3.triplets import get_triplets_at_q, get_nosym_triplets_at_q
class Interaction(object):
class Interaction:
"""Calculate ph-ph interaction and phonons on grid.
This class instance is the heart of phono3py calculation.