3a_tensor_conductivity_save.py — Per-q Conductivity Tensor

Computes the Allen-Feldman conductivity tensor contribution for a single irreducible q-point iq and a block of mode indices [num_start, num_stop). Saves the result as a partial HDF5 file postprocessed by 3b_tensor_conductivity_save_process.py.

This script is not run directly; it is dispatched in a loop by 3c_launch_serial.py, which iterates over all q-points and mode batches.

# Called internally by 3c_launch_serial.py:
python 3a_tensor_conductivity_save.py <iq> <num_start> <num_stop> <gamma_min_plateau_cmm1>

  1# import needed libraries
  2import numpy as np
  3import os, sys
  4import h5py
  5from typing import List
  6
  7# CONSTANTS
  8rY_TO_CMM1 = 109737.315701113
  9hARTREE_SI = 4.35974394E-18
 10rYDBERG_SI = hARTREE_SI / 2.0
 11H_PLANCK_SI = 6.62606896E-34
 12ryau_sec = H_PLANCK_SI / rYDBERG_SI
 13k_BOLTZMANN_SI = 1.3806504E-23  # J K^-1
 14aU_SEC = H_PLANCK_SI / (2. * np.pi) / hARTREE_SI
 15aU_TERAHERTZ = aU_SEC * 1.0E+12
 16rY_TO_THZ = 1.0 / aU_TERAHERTZ / (4 * np.pi)
 17k_BOLTZMANN_RY = k_BOLTZMANN_SI / rYDBERG_SI
 18BOHR_TO_M = 0.52917720859E-10
 19ryvel_si = BOHR_TO_M / ryau_sec
 20SpeedOfLight = 299792458  # [m/s]
 21THzToCm = 1.0e12 / (SpeedOfLight * 100)  # [cm^-1] 33.356410
 22Angstrom = 1.0e-10  # [m]
 23THz = 1.0e12  # [/s]
 24tpi = 2.0 * np.pi
 25
 26if len(sys.argv) > 1:
 27    iq, num_start, num_stop, gamma_min_plateau_cmm1 = (int(sys.argv[1]), int(sys.argv[2]),
 28                                                        int(sys.argv[3]), float(sys.argv[4]))
 29else:
 30    print('error: iq, num_start, num_stop, gamma_min_plateau_cmm1 not given')
 31    exit()
 32
 33temperature_list = np.arange(50, 2001, 50).tolist()
 34n_temp = len(temperature_list)
 35
 36gamma_min_plateau_Ry = gamma_min_plateau_cmm1 / rY_TO_CMM1
 37sigma_Gauss_Ry = np.sqrt(np.pi / 2.0) * gamma_min_plateau_Ry
 38
 39
 40# -------------------------------------------------------------
 41# define function for calculating the conductivity
 42def is_acoustic_mode_at_Gamma(iq, im):
 43    res = False
 44    if (iq == 0 and im < 3):
 45        res = True
 46    return res
 47
 48
 49def conductivity_q(iq: int = iq,
 50                   num_start: int = num_start,
 51                   num_stop: int = num_stop,
 52                   sigma_Gauss: float = sigma_Gauss_Ry,
 53                   temp_list: List = temperature_list):
 54    """
 55    Saves and returns specific heat, Allen conductivity and diffusivity, weight, temperatures and q-point.
 56    Quantities are resolved w.r.t. temperature and mode.
 57
 58    :param iq: q-point for calculation
 59    :param num_start: mode to start - used for parallelisation
 60    :param num_stop: mode to stop (not inclusive) - used for parallelisation
 61    :param sigma_Gauss: sigma which broadens the delta function into a Gaussian [Rydberg]
 62    :param temp_list: list of temperatures for which to calculate conductivities
 63
 64    :return: specific_heat_q,             shape: (num_stop-num_start, n_temp)
 65             thermal_conductivity_Allen_q, shape: (num_stop-num_start, n_temp, 3, 3)
 66             diffusivity_Allen_q,          shape: (num_stop-num_start, 3, 3)
 67             weight,
 68             np.asarray(temp_list),        shape: (n_temp,)
 69             iq
 70    """
 71    smearing = sigma_Gauss * rY_TO_CMM1 / np.sqrt(np.pi / 2.0)  # value of the eta parameter used here
 72
 73    # read frequencies and velocities
 74    filename_vel = './velocity_operators/save_%d.hdf5' % (iq)
 75    f_vel = h5py.File(filename_vel, 'r')
 76
 77    # read velocity operator in phonopy format and convert to QE atomic units
 78    array_V_operator = np.asarray(f_vel['velocity_operator'][:]) * (THz * Angstrom / (2.0 * np.pi * ryvel_si))
 79    array_frequency = np.asarray(f_vel['frequency'])
 80    weight = int(np.copy(f_vel['weight']))
 81
 82    # read frequencies
 83    num_modes = len(array_frequency)
 84    n_temp = len(temp_list)
 85
 86    # setting up arrays for thermal conductivity dataset
 87    array_modes = num_stop - num_start
 88    specific_heat_q = np.zeros((array_modes, n_temp))  # after conversion should be in SI units
 89
 90    thermal_conductivity_Allen_q = np.zeros((array_modes, n_temp, 3, 3))
 91    diffusivity_Allen_q = np.zeros([array_modes, 3, 3])
 92
 93    conv_diffusivity = ((tpi * ryvel_si) ** 2) * ryau_sec / tpi
 94    conv_spec_heat = (1.0 / k_BOLTZMANN_SI) * (rYDBERG_SI) ** 2
 95
 96    # -----------------
 97    # MAIN FOR-LOOP
 98    # -----------------
 99
100    for id_temp in range(0, n_temp):
101        temperature = temp_list[id_temp]
102        tempm1 = 1. / (temperature * k_BOLTZMANN_RY)
103
104        frequencies_Ry = (array_frequency * THzToCm) / rY_TO_CMM1
105        bose_all = 1.0 / (np.exp(frequencies_Ry * tempm1) - 1.0)
106        f_bose_all = bose_all * (bose_all + 1.0)
107        spec_heat_all = (1.0 / (temperature ** 2)) * np.square(frequencies_Ry) * f_bose_all
108
109        specific_heat_q[:, id_temp] = spec_heat_all[num_start:num_stop]
110
111        for im1 in range(num_start, num_stop):  # loop maintained so the memory requirements for large models are not huge
112
113            if is_acoustic_mode_at_Gamma(iq, im1):
114                continue
115
116            omega_1_Ry = frequencies_Ry[im1]
117            spec_heat_s1 = spec_heat_all[im1]
118
119            # vectorisation over second mode
120            local_frequencies = frequencies_Ry
121            local_spec_heat = spec_heat_all
122            n_skip = 3 if iq == 0 else 0  # used for skipping zero modes
123            local_velocity_operator = array_V_operator[im1]
124
125            # distributions
126            Gaussian = (1.0 / (np.sqrt(2.0 * np.pi) * sigma_Gauss)) * np.exp(
127                -0.5 * ((local_frequencies - omega_1_Ry) / sigma_Gauss) ** 2)
128
129            # calculate v_alpha v_beta tensor
130            local_velocity_tensor = np.zeros((3, 3, len(local_velocity_operator)))
131            for alpha in range(3):
132                for beta in range(3):
133                    local_velocity_tensor[alpha, beta, :] = (local_velocity_operator[:,
134                                                            alpha] * local_velocity_operator.conj()[:, beta]).real
135
136            # prefactors
137            prefactor_conductivity_tensor = (((omega_1_Ry + local_frequencies) / 4.0) *
138                                      ((spec_heat_s1 / omega_1_Ry) + (local_spec_heat / local_frequencies)) *
139                                      np.pi * local_velocity_tensor)
140
141            prefactor_diffusivity_tensor = (
142                    ((omega_1_Ry + local_frequencies) / (2.0 * (spec_heat_s1 + local_spec_heat))) *
143                    ((spec_heat_s1 / omega_1_Ry) + (local_spec_heat / local_frequencies)) *
144                    local_velocity_tensor *
145                    np.pi)
146
147            thermal_conductivity_Allen_q[im1 - num_start, id_temp, :, :] = (prefactor_conductivity_tensor * Gaussian)[:, :, n_skip:].sum(axis=2)
148
149            # diffusivity is temperature-independent; compute once at the last temperature step
150            if id_temp == n_temp - 1:
151                diffusivity_Allen_q[im1 - num_start] = (prefactor_diffusivity_tensor * Gaussian)[:, :, n_skip:].sum(axis=2)
152
153    result = (specific_heat_q * conv_spec_heat,
154             thermal_conductivity_Allen_q * conv_spec_heat * conv_diffusivity,
155             diffusivity_Allen_q * conv_diffusivity,
156             weight,
157             np.asarray(temp_list),
158             iq)
159
160    os.makedirs('data_save/interpolation_conductivity_dataset', exist_ok=True)
161    np.save(
162        f'./data_save/interpolation_conductivity_dataset/results_iq_{iq}_start_{num_start}_stop_{num_stop}_gamma_{int(np.round(smearing * 1000))}.npy',
163        np.array(result, dtype=object))
164
165    return result
166
167result = conductivity_q()