This notebook can be run on Google Colab.

Open in Colab

In Colab, you can enable the GPU acceleration from Edit > Notebook Settings > Accelerator > GPU.

In Colab, you can enable the TPU acceleration from Edit > Notebook Settings > Accelerator > TPU.

Silicon diamond tutorial

Complie LAMMPS as shared-library of python (\sim 8 min)

[ ]:

%%capture
# 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
! wget https://download.lammps.org/tars/lammps-4May2022.tar.gz
! tar xzvf lammps-4May2022.tar.gz
! mv lammps-4May2022 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

[5]:

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 \kappa_{per \ mode} and \kappa_{cum}

[9]:

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()