2a_convergence_serial.py — AF Smearing Convergence Study¶
Runs the Allen-Feldman diffusivity calculation over a range of Lorentzian smearing values η and temperatures. For each (η, T) pair it evaluates the AF thermal conductivity using the velocity operators produced by step 1a.
Results are saved as AFC_convergence_results.npz and consumed by
2b_save_convergence.py.
Run after 1a_calc_vel_ops.py.
cd tutorials/diffusivity
python 2a_convergence_serial.py
1# import needed libraries
2import numpy as np
3import os, sys
4import h5py
5from ase.io import read
6from typing import Tuple, List
7import glob
8from tqdm import tqdm
9
10# CONSTANTS
11num_modes_start = 0
12rY_TO_CMM1 = 109737.315701113
13hARTREE_SI = 4.35974394E-18
14rYDBERG_SI = hARTREE_SI / 2.0
15H_PLANCK_SI = 6.62606896E-34
16ryau_sec = H_PLANCK_SI / rYDBERG_SI
17k_BOLTZMANN_SI = 1.3806504E-23 # J K^-1
18aU_SEC = H_PLANCK_SI / (2. * np.pi) / hARTREE_SI
19aU_TERAHERTZ = aU_SEC * 1.0E+12
20rY_TO_THZ = 1.0 / aU_TERAHERTZ / (4 * np.pi)
21k_BOLTZMANN_RY = k_BOLTZMANN_SI / rYDBERG_SI
22BOHR_TO_M = 0.52917720859E-10
23ryvel_si = BOHR_TO_M / ryau_sec
24SpeedOfLight = 299792458 # [m/s]
25THzToCm = 1.0e12 / (SpeedOfLight * 100) # [cm^-1] 33.356410
26Angstrom = 1.0e-10 # [m]
27THz = 1.0e12 # [/s]
28tpi = 2.0 * np.pi
29
30temperature_list = [100, 300, 1500]
31
32
33# -------------------------------------------------------------
34# define function for calculating the conductivity
35def is_acoustic_mode_at_Gamma(iq, im):
36 res = False
37 if (iq == 0 and im < 3):
38 res = True
39 return res
40
41
42def conductivity_q(iq: int, sigma_Gauss: float, temp_list: List = temperature_list):
43 """
44 :param iq: q-point on a mesh for calculation of the conductivity
45 :param sigma_Gauss: sigma which broadens the delta function into a Gaussian [Rydberg]
46 :param temp_list: list of temperatures for which to calculate conductivities
47 :return: tuple (specific_heat_q, tc_Allen_q, tc_Allen_Hardy_q,
48 tc_Allen_Lorentz_q, tc_Allen_Lorentz_Hardy_q, temp_array, iq)
49 """
50 smearing = sigma_Gauss * rY_TO_CMM1 / np.sqrt(np.pi / 2.0)
51
52 filename_vel = './velocity_operators/save_%d.hdf5' % (iq)
53 f_vel = h5py.File(filename_vel, 'r')
54
55 # read velocity operator in phonopy format and convert to QE atomic units
56 array_V_operator = np.asarray(f_vel['velocity_operator'][:]) * (THz * Angstrom / (2.0 * np.pi * ryvel_si))
57 array_frequency = np.asarray(f_vel['frequency'])
58 weight = int(np.copy(f_vel['weight']))
59
60 shape_arr = array_V_operator.shape
61 num_modes = shape_arr[1]
62 n_temp = len(temp_list)
63
64 # setting up arrays for thermal conductivity
65 thermal_conductivity_Allen_q = np.zeros(n_temp)
66 thermal_conductivity_Allen_q_Hardy = np.zeros(n_temp)
67
68 thermal_conductivity_Allen_q_Lorentz = np.zeros(n_temp)
69 thermal_conductivity_Allen_q_Lorentz_Hardy = np.zeros(n_temp)
70
71 specific_heat_q = np.zeros(n_temp)
72
73 lw_Ry_1_approx = sigma_Gauss / (np.sqrt(np.pi / 2.0))
74 lw_Ry_2_approx = lw_Ry_1_approx
75
76 conv_diffusivity = ((tpi * ryvel_si) ** 2) * ryau_sec / tpi
77 conv_spec_heat = (1.0 / k_BOLTZMANN_SI) * (rYDBERG_SI) ** 2
78
79 for id_temp in range(0, n_temp):
80 temperature = temp_list[id_temp]
81 tempm1 = 1. / (temperature * k_BOLTZMANN_RY)
82
83 frequencies_Ry = (array_frequency * THzToCm) / rY_TO_CMM1
84 bose_all = 1.0 / (np.exp(frequencies_Ry * tempm1) - 1.0)
85 f_bose_all = bose_all * (bose_all + 1.0)
86 spec_heat_all = (1.0 / (temperature ** 2)) * np.square(frequencies_Ry) * f_bose_all
87
88 if iq == 0:
89 specific_heat_q[id_temp] = spec_heat_all[3:].sum()
90 else:
91 specific_heat_q[id_temp] = spec_heat_all.sum()
92
93 for im1 in range(0, num_modes): # loop maintained so the memory requirements for large models are not huge
94
95 if is_acoustic_mode_at_Gamma(iq, im1):
96 continue
97
98 omega_1_Ry = frequencies_Ry[im1]
99 spec_heat_s1 = spec_heat_all[im1]
100
101 # define Gaussian and Lorentzian distributions
102 Gaussian = (1.0 / (np.sqrt(2.0 * np.pi) * sigma_Gauss)) * np.exp(
103 -0.5 * ((frequencies_Ry - omega_1_Ry) / sigma_Gauss) ** 2)
104 Lorenztian = (1.0 / np.pi) * ((0.5 * (lw_Ry_1_approx + lw_Ry_2_approx)) / (
105 (frequencies_Ry - omega_1_Ry) ** 2 + 0.25 * (lw_Ry_1_approx + lw_Ry_2_approx) ** 2))
106
107 # calculate square mod average of velocities
108 square_mod_avg_Ry = ((array_V_operator[im1] * np.conjugate(array_V_operator[im1])).sum(axis=1) / 3.0).real
109
110 # used for conversion to Hardy conductivity
111 factor_velocity = (frequencies_Ry + omega_1_Ry) / (2.0 * np.sqrt(frequencies_Ry * omega_1_Ry))
112 factor_velocity_square = np.square(factor_velocity)
113
114 # Allen Feldman conductivity with Gaussian + its Hardy counterpart
115 contr_Allen = (
116 ((frequencies_Ry + omega_1_Ry) / 4.0) *
117 ((spec_heat_all / frequencies_Ry) + (spec_heat_s1 / omega_1_Ry)) *
118 square_mod_avg_Ry *
119 np.pi * Gaussian)
120
121 if iq == 0: # avoid the zero translation modes
122 thermal_conductivity_Allen_q[id_temp] += contr_Allen[3:].sum()
123 thermal_conductivity_Allen_q_Hardy[id_temp] += (contr_Allen[3:] * factor_velocity_square[3:]).sum()
124 else:
125 thermal_conductivity_Allen_q[id_temp] += contr_Allen.sum()
126 thermal_conductivity_Allen_q_Hardy[id_temp] += (contr_Allen * factor_velocity_square).sum()
127
128 # Allen and Feldman conductivity with Lorentzian function + Hardy counterpart
129 contr_Lorentz = (
130 ((frequencies_Ry + omega_1_Ry) / 4.0) *
131 ((spec_heat_all / frequencies_Ry) + (spec_heat_s1 / omega_1_Ry)) *
132 square_mod_avg_Ry *
133 np.pi * Lorenztian)
134
135 if iq == 0:
136 thermal_conductivity_Allen_q_Lorentz[id_temp] += contr_Lorentz[3:].sum()
137 thermal_conductivity_Allen_q_Lorentz_Hardy[id_temp] += (
138 contr_Lorentz[3:] * factor_velocity_square[3:]).sum()
139 else:
140 thermal_conductivity_Allen_q_Lorentz[id_temp] += contr_Lorentz.sum()
141 thermal_conductivity_Allen_q_Lorentz_Hardy[id_temp] += (contr_Lorentz * factor_velocity_square).sum()
142
143 return (
144 specific_heat_q * conv_spec_heat,
145 thermal_conductivity_Allen_q * conv_spec_heat * conv_diffusivity,
146 thermal_conductivity_Allen_q_Hardy * conv_spec_heat * conv_diffusivity,
147 thermal_conductivity_Allen_q_Lorentz * conv_spec_heat * conv_diffusivity,
148 thermal_conductivity_Allen_q_Lorentz_Hardy * conv_spec_heat * conv_diffusivity,
149 np.asarray(temp_list),
150 iq)
151
152
153num_q = len(glob.glob("velocity_operators/save*.hdf5"))
154list_smear = [0.025, 0.05, 0.075, 0.1, 0.2, 0.5, 1.0, 2.0, 4.0, 8.0, 12.0, 16.0, 20.0, 30, 40, 60, 80, 100]
155n_smear = len(list_smear)
156n_temp = len(temperature_list)
157
158# preallocate result arrays indexed [iq, i_smear, i_temp]
159specific_heat = np.zeros((num_q, n_smear, n_temp))
160tc_Allen = np.zeros((num_q, n_smear, n_temp))
161tc_Allen_Hardy = np.zeros((num_q, n_smear, n_temp))
162tc_Allen_Lorentz = np.zeros((num_q, n_smear, n_temp))
163tc_Allen_Lorentz_Hardy = np.zeros((num_q, n_smear, n_temp))
164
165for iq in tqdm(range(num_q), desc='q-points'):
166 for i_smear, smear in enumerate(list_smear):
167 sigma_Gauss_Ry = smear * np.sqrt(np.pi / 2.0) / rY_TO_CMM1
168 result = conductivity_q(iq=iq, sigma_Gauss=sigma_Gauss_Ry)
169 specific_heat[iq, i_smear] = result[0]
170 tc_Allen[iq, i_smear] = result[1]
171 tc_Allen_Hardy[iq, i_smear] = result[2]
172 tc_Allen_Lorentz[iq, i_smear] = result[3]
173 tc_Allen_Lorentz_Hardy[iq, i_smear] = result[4]
174
175np.savez('AFC_convergence_results.npz',
176 specific_heat=specific_heat,
177 tc_Allen=tc_Allen,
178 tc_Allen_Hardy=tc_Allen_Hardy,
179 tc_Allen_Lorentz=tc_Allen_Lorentz,
180 tc_Allen_Lorentz_Hardy=tc_Allen_Lorentz_Hardy,
181 temperatures=np.array(temperature_list),
182 smearings=np.array(list_smear))