3b_tensor_conductivity_save_process.py — Assemble Conductivity Tensors

Collects all per-q-point, per-mode-batch HDF5 files written by 3a_tensor_conductivity_save.py, applies BZ weighting, and assembles the full Allen-Feldman conductivity tensor dataset into data_save/IC_dataset_tensor.npz.

Run after 3c_launch_serial.py (which calls 3a for every q-point).

cd tutorials/diffusivity
python 3b_tensor_conductivity_save_process.py

  1import os, sys
  2import numpy as np
  3from ase.io import read
  4from scipy.constants import physical_constants
  5import h5py
  6import glob
  7
  8
  9# CONSTANTS
 10rY_TO_CMM1 = 109737.315701113
 11hARTREE_SI = 4.35974394E-18
 12rYDBERG_SI = hARTREE_SI / 2.0
 13H_PLANCK_SI = 6.62606896E-34
 14ryau_sec = H_PLANCK_SI / rYDBERG_SI
 15k_BOLTZMANN_SI = 1.3806504E-23  # J K^-1
 16aU_SEC = H_PLANCK_SI / (2. * np.pi) / hARTREE_SI
 17aU_TERAHERTZ = aU_SEC * 1.0E+12
 18rY_TO_THZ = 1.0 / aU_TERAHERTZ / (4 * np.pi)
 19k_BOLTZMANN_RY = k_BOLTZMANN_SI / rYDBERG_SI
 20BOHR_TO_M = 0.52917720859E-10
 21ryvel_si = BOHR_TO_M / ryau_sec
 22SpeedOfLight = 299792458  # [m/s]
 23THzToCm = 1.0e12 / (SpeedOfLight * 100)  # [cm^-1] 33.356410
 24Angstrom = 1.0e-10  # [m]
 25THz = 1.0e12  # [/s]
 26tpi = 2.0 * np.pi
 27eV_TO_JOULE = 1.602177e-19
 28amu_TO_kg = 1.660539e-27
 29
 30temp_list = np.arange(50, 2001, 50).tolist()
 31n_temp = len(temp_list)
 32
 33
 34# change parameters HERE
 35grid = [5, 5, 5]
 36num_q = int((grid[0] * grid[1] * grid[2] + 1) / 2)
 37gamma_min_plateau_cmm1 = 8.0
 38
 39
 40
 41# get density
 42cell = read('POSCAR')
 43n_atoms_structure = len(cell.get_positions())
 44num_modes = n_atoms_structure * 3
 45amount_of_modes = 1000
 46num_starts = np.arange(0, num_modes + amount_of_modes, amount_of_modes)
 47num_starts[-1] = num_modes
 48
 49
 50volume_A3 = cell.get_volume()
 51total_mass = cell.get_masses().sum()  # in atomic mass units
 52
 53amu = physical_constants['atomic mass constant'][0]
 54conversion_factor = amu * 1e27
 55
 56# converting from atomic mass units per angstrom^3 to kg/m^3
 57density_kg_m3 = total_mass / volume_A3 * conversion_factor * 1000
 58
 59dataset = {'name': f'density {density_kg_m3:.3f}, {n_atoms_structure} atoms', 'weights': [], 'qpoints': [],
 60           'frequencies': [], 'temperatures': temp_list, 'specific_heats': [],
 61           'tc_Allen': [], 'diff_Allen': [],
 62           'gamma_min_plateau': gamma_min_plateau_cmm1, 'density': density_kg_m3}
 63
 64for iq in range(num_q):
 65    file = f"velocity_operators/save_{iq}.hdf5"
 66
 67    with h5py.File(file, 'r') as f:
 68        keys = list(f.keys())
 69        print(f"file: {file}, keys: {keys}")
 70        dataset["frequencies"].append(np.copy(f['frequency'])[:] * THzToCm)
 71
 72        dataset["weights"].append(int(np.copy(f['weight'])))
 73        dataset['qpoints'].append(np.copy(f['qpoint']))
 74
 75    temp = {'specific_heats': [],
 76            'tc_Allen': [],
 77            'diff_Allen': []}
 78
 79    for idx in range(len(num_starts) - 1):
 80        num_start = num_starts[idx]
 81        num_stop = num_starts[idx + 1]
 82        infile = f'./data_save/interpolation_conductivity_dataset/results_iq_{iq}_start_{num_start}_stop_{num_stop}_gamma_{int(1000*gamma_min_plateau_cmm1)}.npy'
 83        results = np.load(infile, allow_pickle=True)
 84
 85        normalization_NV = grid[0] * grid[1] * grid[2] * volume_A3 * (1E-10 ** 3)
 86
 87        temp['specific_heats'].append(results[0] / (normalization_NV * density_kg_m3))
 88        temp['tc_Allen'].append(results[1] / normalization_NV)
 89        temp['diff_Allen'].append(results[2])
 90
 91
 92    for key in ['specific_heats', 'tc_Allen', 'diff_Allen']:
 93        # perform vstack
 94        curr_array = temp[key]
 95        n_stacks = len(curr_array)
 96        temp[key] = np.vstack([curr_array[local_idx] for local_idx in range(n_stacks)])
 97
 98        dataset[key].append(temp[key])
 99
100list_keys = ['weights', 'qpoints', 'frequencies', 'specific_heats', 'temperatures',
101             'tc_Allen', 'diff_Allen']
102
103for key in list_keys:
104    dataset[key] = np.array(dataset[key])
105
106outfile = f'data_save/IC_dataset_tensor.npz'
107np.savez(outfile, **dataset)