mirror of https://github.com/phonopy/phonopy.git
Fix ave_pp write out problem and add kappa_conversion_factor to kappa-xxx.hdf5
This commit is contained in:
Atsushi Togo
parent
68fbd2b5c8
commit
893ec5470d
|
|
@ -637,6 +637,7 @@ def write_kappa_to_hdf5(temperature,
|
|||
grid_point=None,
|
||||
band_index=None,
|
||||
sigma=None,
|
||||
kappa_unit_conversion=None,
|
||||
filename=None,
|
||||
verbose=True):
|
||||
if band_index is None:
|
||||
|
|
@ -671,6 +672,9 @@ def write_kappa_to_hdf5(temperature,
|
|||
w.create_dataset('qpoint', data=qpoint)
|
||||
if weight is not None:
|
||||
w.create_dataset('weight', data=weight)
|
||||
if kappa_unit_conversion is not None:
|
||||
w.create_dataset('kappa_unit_conversion',
|
||||
data=kappa_unit_conversion)
|
||||
|
||||
if verbose:
|
||||
text = ""
|
||||
|
|
|
|||
|
|
@ -4,7 +4,7 @@ from phonopy.units import THzToEv, THz, Angstrom
|
|||
from phonopy.phonon.thermal_properties import mode_cv as get_mode_cv
|
||||
from anharmonic.file_IO import (write_kappa_to_hdf5, write_triplets,
|
||||
read_gamma_from_hdf5, write_grid_address)
|
||||
from anharmonic.phonon3.conductivity import Conductivity
|
||||
from anharmonic.phonon3.conductivity import Conductivity, unit_to_WmK
|
||||
from anharmonic.phonon3.imag_self_energy import ImagSelfEnergy
|
||||
from anharmonic.phonon3.triplets import get_grid_points_by_rotations
|
||||
|
||||
|
|
@ -58,7 +58,7 @@ def get_thermal_conductivity_RTA(
|
|||
if not _set_gamma_from_file(br, filename=input_filename):
|
||||
print("Reading collisions failed.")
|
||||
return False
|
||||
|
||||
|
||||
for i in br:
|
||||
if write_gamma:
|
||||
_write_gamma(br,
|
||||
|
|
@ -72,7 +72,10 @@ def get_thermal_conductivity_RTA(
|
|||
if write_kappa:
|
||||
if (grid_points is None and _all_bands_exist(interaction)):
|
||||
br.set_kappa_at_sigmas()
|
||||
_write_kappa(br, filename=output_filename, log_level=log_level)
|
||||
_write_kappa(br,
|
||||
interaction,
|
||||
filename=output_filename,
|
||||
log_level=log_level)
|
||||
|
||||
return br
|
||||
|
||||
|
|
@ -87,9 +90,14 @@ def _write_gamma(br, interaction, i, filename=None, verbose=True):
|
|||
gamma = br.get_gamma()
|
||||
gamma_isotope = br.get_gamma_isotope()
|
||||
sigmas = br.get_sigmas()
|
||||
volume = interaction.get_primitive().get_volume()
|
||||
|
||||
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:
|
||||
|
|
@ -101,18 +109,22 @@ def _write_gamma(br, interaction, i, filename=None, verbose=True):
|
|||
frequency=frequencies,
|
||||
group_velocity=group_velocities[i],
|
||||
heat_capacity=mode_heat_capacities[:, i],
|
||||
kappa=None,
|
||||
gamma=gamma[j, :, i],
|
||||
gamma_isotope=gamma_isotope_at_sigma,
|
||||
averaged_pp_interaction=ave_pp[i],
|
||||
averaged_pp_interaction=ave_pp_i,
|
||||
mesh_divisors=mesh_divisors,
|
||||
grid_point=gp,
|
||||
sigma=sigma,
|
||||
kappa_unit_conversion=unit_to_WmK / volume,
|
||||
filename=filename,
|
||||
verbose=verbose)
|
||||
else:
|
||||
for j, sigma in enumerate(sigmas):
|
||||
for k, bi in enumerate(interaction.get_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, k]
|
||||
if gamma_isotope is not None:
|
||||
gamma_isotope_at_sigma = gamma_isotope[j, i, k]
|
||||
|
|
@ -124,19 +136,17 @@ def _write_gamma(br, interaction, i, filename=None, verbose=True):
|
|||
frequency=frequencies,
|
||||
group_velocity=group_velocities[i, k],
|
||||
heat_capacity=mode_heat_capacities[:, i, k],
|
||||
kappa=None,
|
||||
gamma=gamma[j, :, i, k],
|
||||
gamma_isotope=gamma_isotope_at_sigma,
|
||||
averaged_pp_interaction=ave_pp[i, k],
|
||||
averaged_pp_interaction=ave_pp_ik,
|
||||
mesh_divisors=mesh_divisors,
|
||||
grid_point=gp,
|
||||
band_index=bi,
|
||||
sigma=sigma,
|
||||
kappa_unit_conversion=unit_to_WmK / volume,
|
||||
filename=filename,
|
||||
verbose=verbose)
|
||||
|
||||
|
||||
|
||||
def _all_bands_exist(interaction):
|
||||
band_indices = interaction.get_band_indices()
|
||||
num_band = interaction.get_primitive().get_number_of_atoms() * 3
|
||||
|
|
@ -157,7 +167,7 @@ def _write_triplets(interaction, filename=None):
|
|||
filename=filename)
|
||||
write_grid_address(grid_address, mesh, filename=filename)
|
||||
|
||||
def _write_kappa(br, filename=None, log_level=0):
|
||||
def _write_kappa(br, interaction, filename=None, log_level=0):
|
||||
temperatures = br.get_temperatures()
|
||||
sigmas = br.get_sigmas()
|
||||
gamma = br.get_gamma()
|
||||
|
|
@ -175,7 +185,8 @@ def _write_kappa(br, filename=None, log_level=0):
|
|||
num_ignored_phonon_modes = br.get_number_of_ignored_phonon_modes()
|
||||
num_band = br.get_frequencies().shape[1]
|
||||
num_phonon_modes = br.get_number_of_sampling_grid_points() * num_band
|
||||
|
||||
volume = interaction.get_primitive().get_volume()
|
||||
|
||||
for i, sigma in enumerate(sigmas):
|
||||
kappa_at_sigma = kappa[i]
|
||||
if gamma_isotope is not None:
|
||||
|
|
@ -194,7 +205,7 @@ def _write_kappa(br, filename=None, log_level=0):
|
|||
("T(K)", "xx", "yy", "zz", "yz", "xz", "xy"))
|
||||
for j, (t, k) in enumerate(zip(temperatures, kappa_at_sigma)):
|
||||
print(("%7.1f" + " %10.3f" * 6 + " %d/%d") %
|
||||
((t,) + tuple(k) +
|
||||
((t,) + tuple(k) +
|
||||
(num_ignored_phonon_modes[i, j], num_phonon_modes)))
|
||||
else:
|
||||
print(("#%6s " + " %-10s" * 6) %
|
||||
|
|
@ -217,9 +228,10 @@ def _write_kappa(br, filename=None, log_level=0):
|
|||
weight=weights,
|
||||
mesh_divisors=mesh_divisors,
|
||||
sigma=sigma,
|
||||
kappa_unit_conversion=unit_to_WmK / volume,
|
||||
filename=filename,
|
||||
verbose=log_level)
|
||||
|
||||
|
||||
def _set_gamma_from_file(br, filename=None, verbose=True):
|
||||
sigmas = br.get_sigmas()
|
||||
mesh = br.get_mesh_numbers()
|
||||
|
|
@ -335,7 +347,7 @@ class Conductivity_RTA(Conductivity):
|
|||
self._symmetry = None
|
||||
self._point_operations = None
|
||||
self._rotations_cartesian = None
|
||||
|
||||
|
||||
self._grid_points = None
|
||||
self._grid_weights = None
|
||||
self._grid_address = None
|
||||
|
|
@ -352,7 +364,7 @@ class Conductivity_RTA(Conductivity):
|
|||
self._averaged_pp_interaction = None
|
||||
self._num_ignored_phonon_modes = None
|
||||
self._num_sampling_grid_points = None
|
||||
|
||||
|
||||
self._mesh = None
|
||||
self._mesh_divisors = None
|
||||
self._coarse_mesh = None
|
||||
|
|
@ -387,16 +399,16 @@ class Conductivity_RTA(Conductivity):
|
|||
def set_kappa_at_sigmas(self):
|
||||
num_band = self._primitive.get_number_of_atoms() * 3
|
||||
self._num_sampling_grid_points = 0
|
||||
|
||||
|
||||
for i, grid_point in enumerate(self._grid_points):
|
||||
cv = self._cv[:, i, :]
|
||||
gp = self._grid_points[i]
|
||||
frequencies = self._frequencies[gp]
|
||||
|
||||
|
||||
# Outer product of group velocities (v x v) [num_k*, num_freqs, 3, 3]
|
||||
gv_by_gv_tensor, order_kstar = self._get_gv_by_gv(i)
|
||||
self._num_sampling_grid_points += order_kstar
|
||||
|
||||
|
||||
# Sum all vxv at k*
|
||||
gv_sum2 = np.zeros((6, num_band), dtype='double')
|
||||
for j, vxv in enumerate(
|
||||
|
|
@ -433,7 +445,7 @@ class Conductivity_RTA(Conductivity):
|
|||
|
||||
def set_averaged_pp_interaction(self, ave_pp):
|
||||
self._averaged_pp_interaction = ave_pp
|
||||
|
||||
|
||||
def _run_at_grid_point(self):
|
||||
i = self._grid_point_count
|
||||
self._show_log_header(i)
|
||||
|
|
@ -453,14 +465,14 @@ class Conductivity_RTA(Conductivity):
|
|||
print("Calculating interaction...")
|
||||
|
||||
self._set_gamma_at_sigmas(i)
|
||||
|
||||
|
||||
if self._isotope is not None and not self._read_gamma_iso:
|
||||
self._set_gamma_isotope_at_sigmas(i)
|
||||
|
||||
freqs = self._frequencies[grid_point][self._pp.get_band_indices()]
|
||||
self._cv[:, i, :] = self._get_cv(freqs)
|
||||
self._set_gv(i)
|
||||
|
||||
|
||||
if self._log_level:
|
||||
self._show_log(self._qpoints[i], i)
|
||||
|
||||
|
|
@ -495,7 +507,7 @@ class Conductivity_RTA(Conductivity):
|
|||
self._collision = ImagSelfEnergy(
|
||||
self._pp,
|
||||
unit_conversion=self._gamma_unit_conversion)
|
||||
|
||||
|
||||
def _set_gamma_at_sigmas(self, i):
|
||||
for j, sigma in enumerate(self._sigmas):
|
||||
if self._log_level:
|
||||
|
|
@ -522,14 +534,14 @@ class Conductivity_RTA(Conductivity):
|
|||
self._collision.set_temperature(t)
|
||||
self._collision.run()
|
||||
self._gamma[j, k, i] = self._collision.get_imag_self_energy()
|
||||
|
||||
|
||||
def _get_gv_by_gv(self, i):
|
||||
rotation_map = get_grid_points_by_rotations(
|
||||
self._grid_address[self._grid_points[i]],
|
||||
self._point_operations,
|
||||
self._mesh)
|
||||
gv_by_gv = np.zeros((len(self._gv[i]), 3, 3), dtype='double')
|
||||
|
||||
|
||||
for r in self._rotations_cartesian:
|
||||
gvs_rot = np.dot(self._gv[i], r.T)
|
||||
gv_by_gv += [np.outer(r_gv, r_gv) for r_gv in gvs_rot]
|
||||
|
|
@ -586,7 +598,7 @@ class Conductivity_RTA(Conductivity):
|
|||
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:
|
||||
|
|
|
|||
Loading…
Reference in New Issue