Sodium Chloride Upetcalculator¶

# Denote data path
data_folder = "kaldo_runs/"

# Extract supercell structure for later display
atoms = read(data_folder  + 'fd_upet/replicated_atoms.xyz')

# Retrieve elastic constants and compute bulk modulus
# note that the equation works for cubic system only
Cijs = np.load(data_folder + 'Cij.npy')
C11 = Cijs[0, 0, 0, 0]
C12 = Cijs[0, 0, 1, 1]
C44 = Cijs[1, 2, 1, 2]
Bulk_modulus = (C11 + 2 * C12)/3

# Display supercell and bulk modulus side by side
figure(figsize=(4, 3))
set_fig_properties([gca()])
plot_atoms(atoms)
gca().axis('off')
gca().set_title('NaCl supercell')
show()

print('\n')

# Prepare data for pandas
df = pd.DataFrame({
    "Properties": ["C11", "C12", "C44", "Bulk Modulus"],
    "Predictions (GPa)": [C11, C12, C44, Bulk_modulus],
    "DFT references (GPa)": [47.0, 12.0, 12.0, 24.0]
})

# Format the floats to one decimal place and print
pd.set_option('display.float_format', '{:.1f}'.format)
print(df.to_string(index=False))

# Derive symbols for high symmetry directions in Brillouin zone
dispersion = np.loadtxt(data_folder  + 'plots/12_12_12/dispersion')
point_names = np.loadtxt(data_folder  + 'plots/12_12_12/point_names', dtype=str)

point_names_list = []
for point_name in point_names:
    if point_name == 'G':
        point_name = r'$\Gamma$'
    elif point_name == 'U':
        point_name = 'U=K'
    point_names_list.append(point_name)

# Load in the "tick-mark" values for these symbols
q = np.loadtxt(data_folder  + 'plots/12_12_12//q')
Q = np.loadtxt(data_folder  + 'plots/12_12_12/Q_val')

# Load in grid and dos individually
dos_NaCl = np.load(data_folder + 'plots/12_12_12/dos.npy')
pdos_Na = np.load(data_folder + 'plots/12_12_12/pdos_Na.npy')
pdos_Cl = np.load(data_folder + 'plots/12_12_12/pdos_Cl.npy')

# Plot dispersion
fig = figure(figsize=(4,3))
set_fig_properties([gca()])
plot(q[0], dispersion[0, 0], 'r-', ms=1)
plot(q, dispersion, 'r-', ms=1)
for i in range(1, 4):
    axvline(x=Q[i], ymin=0, ymax=2, ls='--',  lw=2, c="k")
ylabel('Frequency (THz)', fontsize=14)
xlabel(r'Wave vector ($\frac{2\pi}{a}$)', fontsize=14)
gca().set_yticks(np.arange(0, 10, 2))
xticks(Q, point_names_list)
ylim([0, 8])
xlim([Q[0], Q[4]])
dosax = fig.add_axes([0.91, .11, .17, .77])
set_fig_properties([gca()])

# Plot per projection
for p_Na in np.expand_dims(pdos_Na[1],0):
        dosax.plot(p_Na, pdos_Na[0], c='m', label='Na')
for p_Cl in np.expand_dims(pdos_Cl[1],0) :
        dosax.plot(p_Cl, pdos_Cl[0], c='orange', label='Cl')
for d in np.expand_dims(dos_NaCl[1],0):
    dosax.plot(d, dos_NaCl[0],c='r', label='NaCl')
dosax.set_yticks([])
dosax.set_xticks([])
dosax.set_xlabel("DOS")
dosax.legend(fontsize=6.5)
ylim([0, 8])
show()

# Load in group velocity, heat_capacity (cv) and frequency data
frequency =  np.load(
    data_folder + 'ALD_NaCl_upet/12_12_12/frequency.npy',
    allow_pickle=True)

cv =  np.load(
    data_folder + 'ALD_NaCl_upet/12_12_12/300/quantum/heat_capacity.npy',
    allow_pickle=True)


group_velocity = np.load(
    data_folder + 'ALD_NaCl_upet/12_12_12/velocity.npy')

phase_space = np.load(data_folder +
 'ALD_NaCl_upet/12_12_12/300/quantum/_ps_and_gamma.npy',
                      allow_pickle=True)[:,0]

# Compute norm of group velocity
# Convert the unit from angstrom / picosecond to kilometer/ second
group_velcotiy_norm = np.linalg.norm(
    group_velocity.reshape(-1, 3), axis=1) / 10.0

# Plot observables in subplot
figure(figsize=(12, 3))
subplot(1,3, 1)
set_fig_properties([gca()])
scatter(frequency.flatten(order='C')[3:], 1e23*cv.flatten(order='C')[3:],
        facecolor='w', edgecolor='r', s=10, marker='8')
ylabel (r"$C_{v}$ ($10^{23}$ J/K)")
xlabel('Frequency (THz)', fontsize=14)
gca().set_xticks(np.arange(0, 10, 2))
ylim(0.9*1e23*cv.flatten(order='C')[3:].min(), 1.05*1e23*cv.flatten(order='C')[3:].max())

subplot(1 ,3, 2)
set_fig_properties([gca()])
scatter(frequency.flatten(order='C'),
        group_velcotiy_norm, facecolor='w', edgecolor='r', s=10, marker='^')
gca().set_xticks(np.arange(0, 10, 2))
xlabel('Frequency (THz)', fontsize=14)
ylabel(r'$|v| \ (\frac{km}{s})$', fontsize=14)

subplot(1 ,3, 3)
set_fig_properties([gca()])
scatter(frequency.flatten(order='C'),
        phase_space, facecolor='w', edgecolor='r', s=10, marker='o')
gca().set_xticks(np.arange(0, 10, 2))
xlabel('Frequency (THz)', fontsize=14)
ylabel('Phase space', fontsize=14)
subplots_adjust(wspace=0.33)
show()

# Load in scattering rate
scattering_rate = np.load(
    data_folder +
    'ALD_NaCl_upet/12_12_12/300/quantum/bandwidth.npy', allow_pickle=True)

# Derive lifetime, which is inverse of scattering rate
life_time = scattering_rate ** (-1)

# Denote lists to intake mean free path in each direction
mean_free_path = []
for i in range(3):
    mean_free_path.append(np.loadtxt(
    data_folder +
    'ALD_NaCl_upet/12_12_12/inverse/300/quantum/mean_free_path_' +
    str(i) + '.dat'))

# Convert list to numpy array and compute norm for mean free path
# Convert the unit from angstrom  to nanometer
mean_free_path = np.array(mean_free_path).T
mean_free_path_norm = np.linalg.norm(
    mean_free_path.reshape(-1, 3), axis=1) / 10.0

# Plot observables in subplot
figure(figsize=(12, 3))
subplot(1,3, 1)
set_fig_properties([gca()])
scatter(frequency.flatten(order='C'),
        life_time, facecolor='w', edgecolor='r', s=10, marker='s')
gca().set_xticks(np.arange(0, 10, 2))
yscale('log')
ylabel(r'$\tau$  (ps)', fontsize=14)
xlabel('Frequency (THz)', fontsize=14)

subplot(1,3, 2)
set_fig_properties([gca()])
scatter(frequency.flatten(order='C'),
        scattering_rate, facecolor='w', edgecolor='r', s=10, marker='d')
gca().set_xticks(np.arange(0, 10, 2))
ylabel(r'$\Gamma$ (THz)', fontsize=14)
xlabel('Frequency (THz)', fontsize=14)

subplot(1,3, 3)
set_fig_properties([gca()])
scatter(frequency.flatten(order='C'),
        mean_free_path_norm, facecolor='w', edgecolor='r', s=10, marker='8')
gca().set_xticks(np.arange(0, 10, 2))
ylabel(r'$\lambda$ (nm)', fontsize=14)
xlabel('Frequency (THz)', fontsize=14)
yscale('log')
ylim([1e-2, 1e4])
subplots_adjust(wspace=0.33)
show()

# Denote zeros to intake kappa tensor
kappa_tensor = np.zeros([mean_free_path.shape[0], 3, 3])
kappa_tensor_iso = np.zeros([mean_free_path.shape[0], 3, 3])

for i in range(3):
    for j in range(3):
        kappa_tensor[:, i, j] = np.loadtxt(data_folder +
    'ALD_NaCl_upet/12_12_12/inverse/300/quantum/conductivity_' +
    str(i) + '_' + str(j) + '.dat')
        kappa_tensor_iso[:, i, j] = np.loadtxt(data_folder +
    'ALD_NaCl_upet/12_12_12/inverse/300/quantum/isotopes/conductivity_' +
    str(i) + '_' + str(j) + '.dat')

# Sum over the 0th dimension to recover 3-by-3 kappa matrix
kappa_matrix = kappa_tensor.sum(axis=0)
print("Bulk thermal conductivity: %.1f W m^-1 K^-1\n"
      %np.mean(np.diag(kappa_matrix)))
print("kappa matrix: ")
print(kappa_matrix)
print('\n')

kappa_matrix_iso = kappa_tensor_iso.sum(axis=0)
print("Bulk thermal conductivity with isotopic scattering: %.1f W m^-1 K^-1\n"
      %np.mean(np.diag(kappa_matrix_iso)))
print("kappa matrix with isotopic scattering: ")
print(kappa_matrix_iso)
print('\n')

# Compute kappa in per mode and cumulative representations
kappa_per_mode = kappa_tensor.sum(axis=-1).sum(axis=1)
freq_sorted, kappa_cum_wrt_freq = cumulative_cond_cal(
    frequency.flatten(order='C'), kappa_tensor)
lambda_sorted, kappa_cum_wrt_lambda = cumulative_cond_cal(
    mean_free_path_norm, kappa_tensor)

kappa_per_mode_iso = kappa_tensor_iso.sum(axis=-1).sum(axis=1)
freq_sorted, kappa_cum_wrt_freq_iso = cumulative_cond_cal(
    frequency.flatten(order='C'), kappa_tensor_iso)
lambda_sorted, kappa_cum_wrt_lambda_iso = cumulative_cond_cal(
    mean_free_path_norm, kappa_tensor_iso)

# Plot observables in subplot
figure(figsize=(12, 3))
subplot(1,3, 1)
set_fig_properties([gca()])
scatter(frequency.flatten(order='C'),
        kappa_per_mode, facecolor='w', edgecolor='r', s=10, marker='>', label='$\kappa_{per \ mode, pure }$')
scatter(frequency.flatten(order='C'),
        kappa_per_mode_iso, facecolor='w', edgecolor='m', s=10, marker='o', label='$\kappa_{per \ mode, iso }$')
gca().axhline(y = 0, color='k', ls='--', lw=1)
gca().set_xticks(np.arange(0, 10, 2))
ylabel(r'$\kappa_{per \ mode}\;\left(\frac{\rm{W}}{\rm{m}\cdot\rm{K}}\right)$',fontsize=14)
xlabel('Frequency (THz)', fontsize=14)
legend(loc=1, fontsize=10)
ylim([-0.02, 0.4])

subplot(1,3, 2)
set_fig_properties([gca()])
plot(freq_sorted, kappa_cum_wrt_freq, 'r',
     label=r'$\kappa_{pure, PET-MAD} \approx 8.0\;\frac{\rm{W}}{\rm{m}\cdot\rm{K}}$')
plot(freq_sorted, kappa_cum_wrt_freq_iso, c='m',
     ls='-.', label=r"$\kappa_{iso,PET-MAD} \approx 7.5\;\frac{\rm{W}}{\rm{m}\cdot\rm{K}}$")
gca().set_yticks(np.arange(0, 10, 2))
gca().set_xticks(np.arange(0, 10, 2))
ylabel(r'$\kappa_{cumulative, \omega}\;\left(\frac{\rm{W}}{\rm{m}\cdot\rm{K}}\right)$',fontsize=14)
xlabel('Frequency (THz)', fontsize=14)
legend(loc=4, fontsize=6)

subplot(1,3, 3)
set_fig_properties([gca()])
plot(lambda_sorted, kappa_cum_wrt_lambda, 'r')
plot(lambda_sorted, kappa_cum_wrt_lambda_iso, c='m')
gca().set_yticks(np.arange(0, 10, 2))
xlabel(r'$\lambda$ (nm)', fontsize=14)
ylabel(r'$\kappa_{cumulative, \lambda}\;\left(\frac{\rm{W}}{\rm{m}\cdot\rm{K}}\right)$',fontsize=14)
xscale('log')
xlim([1e0, 1e3])
subplots_adjust(wspace=0.33)
show()