2b_save_convergence.py — Save Convergence Summary

Loads the raw convergence data from AFC_convergence_results.npz, applies BZ-weight averaging over all irreducible q-points, and saves summary arrays (conductivity vs smearing for each temperature).

Run after 2a_convergence_serial.py.

cd tutorials/diffusivity
python 2b_save_convergence.py

 1import numpy as np
 2import os
 3from ase.io import read
 4from scipy.constants import physical_constants
 5
 6cell = read('POSCAR')
 7volume_A3 = cell.get_volume()
 8total_mass = cell.get_masses().sum()  # in atomic mass units
 9
10amu = physical_constants['atomic mass constant'][0]
11conversion_factor = amu * 1e27
12
13# converting from atomic mass units per angstrom^3 to kg/m^3
14density_kg_m3 = total_mass / volume_A3 * conversion_factor * 1000
15
16grid = [5, 5, 5]
17num_q = int((grid[0] * grid[1] * grid[2] + 1) / 2)
18
19weight = np.ones(num_q) * 2  # for the 5x5x5 grid, all q-points different from Gamma have weight 2
20weight[0] = 1.0              # Gamma has weight 1
21
22# load single results file written by 2a_convergence_serial.py
23data = np.load('AFC_convergence_results.npz')
24specific_heat_all        = data['specific_heat']          # (n_q, n_smear, n_temp)
25tc_Allen_all             = data['tc_Allen']
26tc_Allen_Hardy_all       = data['tc_Allen_Hardy']
27tc_Allen_Lorentz_all     = data['tc_Allen_Lorentz']
28tc_Allen_Lorentz_Hardy_all = data['tc_Allen_Lorentz_Hardy']
29temperature_array        = data['temperatures']
30smearings                = data['smearings']
31
32n_smear = len(smearings)
33n_temp  = len(temperature_array)
34
35kappas_G       = np.zeros((n_smear, n_temp))
36kappas_L       = np.zeros((n_smear, n_temp))
37kappas_G_gamma = np.zeros((n_smear, n_temp))
38kappas_L_gamma = np.zeros((n_smear, n_temp))
39
40normalization_NV = grid[0] * grid[1] * grid[2] * volume_A3 * (1E-10 ** 3)
41
42for i_smear in range(n_smear):
43    specific_heat_T        = (specific_heat_all[:, i_smear, :]        * weight[:, None]).sum(axis=0)
44    conductivity_Allen     = (tc_Allen_all[:, i_smear, :]             * weight[:, None]).sum(axis=0)
45    conductivity_Allen_Hardy     = (tc_Allen_Hardy_all[:, i_smear, :] * weight[:, None]).sum(axis=0)
46    conductivity_Allen_Lorentz   = (tc_Allen_Lorentz_all[:, i_smear, :] * weight[:, None]).sum(axis=0)
47    conductivity_Allen_Lorentz_Hardy = (tc_Allen_Lorentz_Hardy_all[:, i_smear, :] * weight[:, None]).sum(axis=0)
48
49    specific_heat_T              /= normalization_NV * density_kg_m3
50    conductivity_Allen           /= normalization_NV
51    conductivity_Allen_Hardy     /= normalization_NV
52    conductivity_Allen_Lorentz   /= normalization_NV
53    conductivity_Allen_Lorentz_Hardy /= normalization_NV
54
55    kappas_G[i_smear]       = conductivity_Allen
56    kappas_L[i_smear]       = conductivity_Allen_Lorentz
57    kappas_G_gamma[i_smear] = tc_Allen_all[0, i_smear, :]       * weight[0] / (volume_A3 * (1E-10 ** 3))
58    kappas_L_gamma[i_smear] = tc_Allen_Lorentz_all[0, i_smear, :] * weight[0] / (volume_A3 * (1E-10 ** 3))
59
60os.makedirs('results', exist_ok=True)
61np.save('./results/convergence_test.npy',
62        np.array([kappas_G, kappas_L, kappas_G_gamma, kappas_L_gamma], dtype=object))