This notebook can be run on Google Colab.
In Colab, you can enable the GPU acceleration from Edit > Notebook Settings > Accelerator > GPU.
Silicon diamond tutorial¶
Complie LAMMPS as shared-library of python (
8 min)¶
[ ]:
# Sudo apt-get softwares
! apt-get update
! apt install -y cmake build-essential git ccache openmpi-bin libopenmpi-dev python3.10-venv
! echo "Sudo apt-get finishes!"
# Upgrade python packages
! pip install --upgrade pip
! pip install numpy torch scipy virtualenv psutil pandas tabulate mpi4py Cython
! echo "Python packages installation finishes!"
# Build lammps with cmake
%cd /content
!rm -rf lammps
! git clone https://github.com/lammps/lammps.git lammps
%cd /content/lammps
! rm -rf build
! mkdir build
%cd build
! cmake ../cmake -DLAMMPS_EXCEPTIONS=yes \
-DBUILD_SHARED_LIBS=yes \
-DMLIAP_ENABLE_PYTHON=yes \
-DPKG_PYTHON=yes \
-DPKG_MANYBODY=yes \
-DPKG_KSPACE=yes \
-DPKG_PHONON=yes \
-DPYTHON_EXECUTABLE:FILEPATH=`which python`
# Complie lammps as share-libary of python
! make -j 2
! make install-python
! echo "LAMMPS compilation done!"
# Redirect back to main folder
%cd /content/
Remote fetch and install source code from Github¶
[ ]:
! git lfs install --skip-repo --skip-smudge
! pip install git+https://github.com/nanotheorygroup/kaldo
Remote fetch supplyment files¶
[ ]:
# Remote fetch kaldo resources from drop box
! wget https://www.dropbox.com/s/bvw0qcxy397g25q/kaldo_resources.zip?dl=0
! mv kaldo_resources.zip?dl=0 kaldo_resources.zip
! unzip kaldo_resources.zip
# Unzip files
!unzip forcefields.zip
# Clean workspace
! rm -r forcefields.zip
! rm -r structure_a_si_512.zip
! rm -r kaldo_resources.zip
! rm -r sample_data
Thermal transport simulation for silicon-bulk¶
[ ]:
from ase.build import bulk
from ase.calculators.lammpslib import LAMMPSlib
from kaldo.forceconstants import ForceConstants
import numpy as np
# We start from the atoms object
atoms = bulk('Si', 'diamond', a=5.432)
# Config super cell and calculator input
supercell = np.array([3, 3, 3])
lammps_inputs = {
'lmpcmds': [
'pair_style tersoff',
'pair_coeff * * forcefields/Si.tersoff Si'],
'log_file': 'lammps-si-bulk.log',
'keep_alive':True}
# Create a finite difference object
forceconstants_config = {'atoms':atoms,'supercell': supercell,'folder':'fd'}
forceconstants = ForceConstants(**forceconstants_config)
# Compute 2nd and 3rd IFCs with the defined calculators
forceconstants.second.calculate(LAMMPSlib(**lammps_inputs), delta_shift=1e-3)
forceconstants.third.calculate(LAMMPSlib(**lammps_inputs), delta_shift=1e-3)
Create phonons object¶
[ ]:
from kaldo.phonons import Phonons
# Define k-point grids, temperature
# and the assumption for the
# phonon poluation (i.e classical vs. quantum)
k = 7
kpts = [k, k, k]
temperature = 300
is_classic = False
k_label = str(k) + '_' + str(k) + '_' + str(k)
# Create a phonon object
phonons = Phonons(forceconstants=forceconstants,
kpts=kpts,
is_classic=is_classic,
temperature=temperature,
folder='si-bulk-ald-' + k_label,
storage='numpy')
Calculate conductivities for infinite-size sample¶
[ ]:
from kaldo.conductivity import Conductivity
# Calculate conductivity with direct inversion approach (inverse)
print('\n')
inv_cond_matrix = (Conductivity(phonons=phonons, method='inverse').conductivity.sum(axis=0))
print('Inverse conductivity (W/mK): %.3f'%(np.mean(np.diag(inv_cond_matrix))))
print(inv_cond_matrix)
# Calculate conductivity with relaxation time approximation (rta)
print('\n')
rta_cond_matrix = Conductivity(phonons=phonons, method='rta').conductivity.sum(axis=0)
print('Rta conductivity (W/mK): %.3f'%(np.mean(np.diag(rta_cond_matrix))))
print(rta_cond_matrix)
# Calculate conductivity with self-consistent approach (sc)
print('\n')
sc_cond_matrix = Conductivity(phonons=phonons, method='sc',n_iterations=20).conductivity.sum(axis=0)
print('Self-consistent conductivity (W/mK): %.3f'%(np.mean(np.diag(sc_cond_matrix))))
print(sc_cond_matrix)
Visualize harmonic properties using built-in plotter¶
[ ]:
import kaldo.controllers.plotter as plotter
import matplotlib.pyplot as plt
plt.style.use('seaborn-poster')
# Plot dispersion relation and group velocity in each direction
plotter.plot_dispersion(phonons,n_k_points=int(k_label))
print('\n')
Access and visualize properties calculated during simulations¶
[ ]:
# Direct access to properties
# calculated during the simulation
# Plot heat capacity vs frequency
freq_full = phonons.frequency.flatten(order='C')
cv_1d = phonons.heat_capacity.flatten(order='C')[3:]
print('\n')
plt.figure()
plt.scatter(freq_full[3:],1e23*cv_1d,s=15)
plt.ylabel (r"$C_{v}$ ($10^{23}$ J/K)", fontsize=25, fontweight='bold')
plt.xlabel ("$\\nu$ (THz)", fontsize=25, fontweight='bold')
plt.ylim(0.9*1e23*cv_1d[cv_1d>0].min(),
1.05*1e23*cv_1d.max())
plt.show()
# Plot phonon bandwidth vs frequency
band_width_flatten = phonons.bandwidth.flatten(order='C')
freq = freq_full[band_width_flatten!=0]
print('\n')
plt.figure()
plt.scatter(freq,band_width_flatten[band_width_flatten!=0] ,s=15)
plt.ylabel (r"$\Gamma$ (THz)", fontsize=25, fontweight='bold')
plt.xlabel ("$\\nu$ (THz)", fontsize=25, fontweight='bold')
plt.ylim(0.95*band_width_flatten .min(), 1.05*band_width_flatten .max())
plt.show()
# Plot phase space vs frequency
print('\n')
plt.figure()
plt.scatter(freq_full[3:],phonons.phase_space.flatten(order='C')[3:],s=15)
plt.ylabel ("Phase space", fontsize=25, fontweight='bold')
plt.xlabel ("$\\nu$ (THz)", fontsize=25, fontweight='bold')
plt.ylim(phonons.phase_space.min(), phonons.phase_space.max())
plt.show()
Calculate and visualize 
and 
¶
[ ]:
def cumulative_cond_cal(freq,full_cond,n_phonons):
conductivity = np.einsum('maa->m', 1/3 * full_cond)
conductivity = conductivity.reshape(n_phonons)
cumulative_cond = np.zeros_like(conductivity)
freq_reshaped = freq.reshape(n_phonons)
for mu in range(cumulative_cond.size):
single_cumulative_cond = conductivity[(freq_reshaped < freq_reshaped[mu])].sum()
cumulative_cond[mu] = single_cumulative_cond
return cumulative_cond
[ ]:
# Compute conductivity with per phonon mode basis using different methods
kappa_rta_per_mode = np.einsum('maa->m',1/3*Conductivity(phonons=phonons, method='rta').conductivity)
kappa_inv_per_mode = np.einsum('maa->m',1/3*Conductivity(phonons=phonons, method='inverse').conductivity)
kappa_sc_per_mode = np.einsum('maa->m',1/3*Conductivity(phonons=phonons, method='sc',n_iterations=20).conductivity)
# Compute cumulative conductivity by frequency using different methods
kappa_rta_cum_freq = cumulative_cond_cal(phonons.frequency,Conductivity(phonons=phonons, method='rta').conductivity,phonons.n_phonons)
kappa_sc_cum_freq = cumulative_cond_cal(phonons.frequency,Conductivity(phonons=phonons, method='sc',n_iterations=20).conductivity,phonons.n_phonons)
kappa_inv_cum_freq = cumulative_cond_cal(phonons.frequency,Conductivity(phonons=phonons, method='inverse').conductivity,phonons.n_phonons)
kappa_qhgk_cum_freq = cumulative_cond_cal(phonons.frequency,Conductivity(phonons=phonons, method='qhgk').conductivity,phonons.n_phonons)
print('\n')
# Visualize the cumulative conductivity vs frequency
plt.figure()
plt.plot(freq_full,kappa_rta_per_mode,'r.',label='RTA')
plt.plot(freq_full,kappa_sc_per_mode,'mo',label='Self Consistent',ms=8)
plt.plot(freq_full,kappa_inv_per_mode,'k.',label='Direct Inversion')
plt.xlabel ("$\\nu$ (THz)", fontsize=25, fontweight='bold')
plt.ylabel(r'$\kappa(W/m/K)$', fontsize=25, fontweight='bold')
plt.legend(loc=1,frameon=False)
#plt.grid()
plt.show()
print('\n')
# Visualize the cumulative conductivity vs frequency
plt.figure()
plt.plot(freq_full,kappa_rta_cum_freq,'r.',label='RTA')
plt.plot(freq_full,kappa_sc_cum_freq,'mo',label='Self Consistent',ms=8)
plt.plot(freq_full,kappa_inv_cum_freq,'k.',label='Direct Inversion')
plt.xlabel ("$\\nu$ (THz)", fontsize=25, fontweight='bold')
plt.ylabel(r'$\kappa_{cum}(W/m/K)$', fontsize=25, fontweight='bold')
plt.legend(loc=4,frameon=False)
plt.grid()
plt.show()