diff --git a/anharmonic/file_IO.py b/anharmonic/file_IO.py index d3f99211..668022f1 100644 --- a/anharmonic/file_IO.py +++ b/anharmonic/file_IO.py @@ -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 = "" diff --git a/anharmonic/phonon3/conductivity_RTA.py b/anharmonic/phonon3/conductivity_RTA.py index 1ffc21d7..21177873 100644 --- a/anharmonic/phonon3/conductivity_RTA.py +++ b/anharmonic/phonon3/conductivity_RTA.py @@ -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: