Phonon-limited carrier mobility in semiconductors#
Author: S. Tiwari (v1, 06/01/2024) Revision: F. Giustino (v1.2, 07/03/2024)
In this Noteboook, we compute the phonon-limited carrier mobility of silicon using the ab initio Boltzmann Transport Equation (BTE). Electrons and phonons are computed with Quantum ESPRESSO (QE), maximally-localized Wannier functions are computed with Wannier90, and transport properties are computed with EPW.
Below we define constants that will remail unchanged throughout the Notebook. The object silicon is created as an instance of the EPWpy class. This object will contain everything that we need to execute and analyze the calculations.
[4]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rc
import time, sys, os
import EPWpy
from EPWpy import EPWpy
from EPWpy.plotting import plot_bands
import plotly.io
plotly.io.renderers.default = "sphinx_gallery"
silicon=EPWpy.EPWpy({'prefix':'si',
'calculation':'\'scf\'',
'structure_mp':"mp-149",
'ecutwfc':'40',
'celldm(1)':'10.262',
'pseudo_auto':True,
},
env='mpirun')
silicon.run_serial=True
# Summary
pseudopot=silicon.__dict__['pw_atomic_species']['pseudo'][0]
print('Pseudopotential:', silicon.__dict__['pw_atomic_species']['pseudo'][0])
print('Pseudopotential directory:', silicon.__dict__['pw_control']['pseudo_dir'])
print('Prefix:',silicon.__dict__['prefix'])
app = silicon.display_lattice(supercell=[2,2,1])
app.run()
-*#*- ...............-
.+*= .+%*-=%%: .=#*- -===============-:.
:*%=*%%- *%* #%* :+%+-%%+ .:. -=. :==-.
-%S -%%*: :#%. -%%-. -##: #%* -=. :==-
.. .%S: +%%%%*. :*%%%#= %%= -=. :==-
:=#%%*- .#S- .. .%%= :*#*: -=. :==-
-%S:.=#%%*==%# *%%=::=##-+%%. . .=-. :==- .=
:%%- .-+++: -+##*=. =%S :-::==: .==- --.
*%# #%+ -=--:. .----:.
:%%- -%S.
.-=*####SS%#########: -###########*+: +#######= =%SS####+::.
=+#**%SSSSSSSSSSSSSSSS= =SSSSSSSSSSSSSSS= #SSSSSSS*. +SSSSSSSS%%%%-
*%% =SSS.. .SSS= SSS=. .:+SSS+ -SSS:. .-::. .SSS+. #%*
#%#. =SSS. *S# *S%: SSS= %SS%. .SSS= #SSS%. :SSS: =%#.
*%%: =SSS#*#SSS- SSS= .+SSS+. %SS*.=SSSSS+ +SS%.:+%+
+%%:=SSSSSSSSS- SSSSSSSSSSSSS=. +SSS:SSS%SSS:%SS*-#%=
....#S==SSS:..SSS: =+- SSS%######+=. -SSS%SS%.#SS%SSS==%%=.
.:+##%%#*- =SSS. ::. :SSS. SSS=. SSSSSS:..SSSSSS: :*%%%*=-
#%+. -+#SSS*+++++++*SSS: .=+SSS#++++: %SSSS+. =SSSS%. :-+#%%+
#%* .SSSSSSSSSSSSSSSSSS: *SSSSSSSSSSS. +SSS%. %SSS*. . .%%=
=%S :::::::::::::::::. .:::::::::. :::. :::. =+=-.#%+
-%S: =+===#S:
==*------------------------------=========+++++++++++++++++++++++========++-+##
=+++++++++++*******++++++++++++++++========------------------==========+++++++-
-- -- -- -- -- -- -- -- -- -- -- -- Structure Info -- -- -- -- -- -- -- -- -- -- -- -- --
lattice vector(1): [3.333573 0. 1.924639]
lattice vector(2): [1.111191 3.142924 1.924639]
lattice vector(3): [0. 0. 3.849278]
atom(1): Si [0.875 0.875 0.875]
atom(2): Si [0.125 0.125 0.125]
1
https://raw.githubusercontent.com/PseudoDojo/ONCVPSP-PBE-FR-PDv0.4/master/Si/Si.upf
https://raw.githubusercontent.com/PseudoDojo/ONCVPSP-PBE-FR-PDv0.4/master/Si/Si_r.upf
pseudo found at pseudodojo : ONCVPSP-PBE-FR-PDv0.4/Si_r.upf
-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
Pseudopotential: Si_r.upf
Pseudopotential directory: '/mnt/storage/sabya/For_video/epwpy/notebooks_basic/pseudo/'
Prefix: si
Self-consistent field (SCF) calculations#
Here we perform the self-consistent field calculation to obtain the electron charge density of silicon in the ground state. The calculation consists of three separate steps: 1. Apply the method scf to the object silicon. This step specifies runtime parameters for an SCF calculation on siicon 2. Based on the properties defined at step 1 as well as other properties that are set by default within EPWpy, the method prepare creates the input file needed by QE 3. The method run applied
to the object silicon instructs QE to perform the SCF calculation
[2]:
silicon.scf(name='scf',kpoints={'kpoints':[[6,6,6]]})
silicon.prepare(1,type_run='scf')
silicon.run(4)
silicon.pw_util = silicon.PW_utilities()
-- -- -- -- -- -- -- -- -- -- -- Calculation: scf -- -- -- -- -- -- -- -- -- -- --
Running scf |████████████████████████████████████████| in 3.9s (0.42/s)
-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
Band structure calculation#
In this step, we compute the band structure of silicon along some high-symmetry lines in the Brillouin zone.
This calculation is not strictly necessary to compute the mobility, but it is useful to understand the electronic structure of the system under consideration.
Also in this case, we use three instructions to specify runtime parameters, prepare the input file, and execute the QE calculation.
[3]:
silicon.scf(control={'calculation':'\'bands\''},
system={'nbnd':12},
kpoints={'kpoints':[
['0.5','0.50','0.50','20'],
['0.0','0.00','0.00','20'],
['0.5','0.25','0.75','20']
],
'kpoints_type':'{crystal_b}'
},
name='bs')
silicon.prepare(type_run='bs')
silicon.run(4,type_run='bs')
-- -- -- -- -- -- -- -- -- -- -- Calculation: bs -- -- -- -- -- -- -- -- -- -- --
Running bs |████████████████████████████████████████| in 4.7s (0.34/s)
-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
Band structure plot#
We now plot the electronic band structure computed at the previous step. The zero of the energy axis is set to the value specified manually via ef0.
[4]:
ef_from_file = silicon.pw_util.efermi
Band=plot_bands.plot_band_scf(f'./{silicon.prefix}/bs/bs.out')
plot_bands.plot_band_prod(Band,
ef0=ef_from_file,
xticks=['X','$\Gamma$','L'],
xlabel = 'Wavevector',
ylabel = 'Energy (eV)'
)
Phonon dispersion relations#
To compute phonon-limited mobilities, we need to determine vibrational frequencies and eigenmodes. This operation consists of three steps 1. We compute these properties on a uniform and centered Brillouin zone grid 2. We perform a Fourier transform of the results in order to obtain the interatomic force constants (IFCs) 3. We interpolate the IFCs along specified Brillouin zone paths to obtain phonon dispersions.
This plot of phonon dispersions is only meant for us to develop a qualitative understanding of the vibrational properties of the system under consideration. The phonon interpolation needed for transport calculations is performed once again later by EPW.
Step 1: Calculations of phonons on uniform Brillouin zone grid#
[5]:
silicon.ph(phonons={'fildyn':'\'si.dyn\'',
'nq1':3,
'nq2':3,
'nq3':3,
'fildvscf':'\'dvscf\''}
)
silicon.prepare(type_run='ph')
silicon.run(16,type_run='ph')
-- -- -- -- -- -- -- -- -- -- -- Calculation: ph -- -- -- -- -- -- -- -- -- -- --
Running ph |████████████████████████████████████████| in 5:28.7 (0.00/s)
-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
Step 2: Generation of IFCs#
[6]:
silicon.q2r(name='q2r')
silicon.prepare(type_run='q2r')
silicon.run(1,type_run='q2r')
-- -- -- -- -- -- -- -- -- -- -- Calculation: q2r -- -- -- -- -- -- -- -- -- -- --
Running q2r |████████████████████████████████████████| in 2.1s (1.08/s)
-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
Step 3: Interpolation of IFCs and generation of phonon dispersions plot#
The logic and syntax of this operation are the same as for the band structure plot above: three instructions to execute matdyn.x and then plotting.
[7]:
silicon.matdyn(name='matdyn',
kpoints={'kpoints':[
['0.5','0.50','0.50','20'],
['0.0','0.00','0.00','20'],
['0.5','0.25','0.75','20']
],
'kpoints_type':'{crystal_b}'
},
)
silicon.prepare(type_run='matdyn')
silicon.run(1,type_run='matdyn')
-- -- -- -- -- -- -- -- -- -- -- Calculation: matdyn -- -- -- -- -- -- -- -- -- -- --
Running matdyn |████████████████████████████████████████| in 1.3s (3.28/s)
-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
[8]:
Band=plot_bands.plot_band_eig(f'./{silicon.prefix}/ph/si.freq')
plot_bands.plot_band_freq(Band,ylabel='E (meV)',ef0=0,xticks=['L','$\Gamma$','X'],color='royalblue')
Transformation of electrons and phonons to Wannier basis with EPW#
Now we have Kohn-Sham wavefunctions and variations of the self-consistent Kohn-Sham potential on coarse Brillouin zone grid. We will generate the electron Hamiltonian, the IFCs, and the electron-phonon matrix elements in the Wannier representation using EPW. Details on the underlying formalism can be found here (free version) or here (journal version).
This operation involves four logicals steps: 1. Compute Kohn-Sham states on a uniform centered Brillouin zone grid (QE) 2. Use EPW to load these states and call Wannier90 to generate Wannier functions 3. Use EPW to load IFCs and potential variations, and combine with 2. to compute electron-phonon matrix elements in the Bloch representation 4. Use EPW to perform the transformation to the Wannier basis and write to file
Step 1: Calculations of Kohn-Sham states on uniform Brillouin zone grid#
[9]:
silicon.nscf(system={'nbnd':8},
kpoints={'grid':[6,6,6],'kpoints_type': 'crystal'})
silicon.prepare(type_run='nscf')
silicon.run(4,type_run='nscf')
-- -- -- -- -- -- -- -- -- -- -- Calculation: nscf -- -- -- -- -- -- -- -- -- -- --
Running nscf |████████████████████████████████████████| in 11.0s (0.13/s)
-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
Steps 2-4: Load Bloch representation, Wannierize, write to file quantities in Wannier representation#
[18]:
# File with k-path for sanity checks
silicon.epw(epwin={'wdata':['guiding_centres = .true.',
'dis_num_iter = 500',
'num_print_cycles = 10',
'dis_mix_ratio = 1',
'use_ws_distance = T'],
'proj':['\'Si : sp3\''],
'band_plot':'.true.',
'filkf':silicon.filkf_file,
'filqf':silicon.filkf_file,
'etf_mem':0,
'fsthick':12.0,
'wannierize':'.true.',
'elph':'.true.',
'num_iter':500
},
name='epw1')
silicon.filkf_file = 'LGX.txt'
silicon.filkf(path=[[0.5,0.5,0.5],
[0.0,0.0,0.0],
[0.5,0.25,0.75]],
length=[51,51],
)
silicon.prepare(type_run='epw1')
silicon.run(8,type_run='epw1')
(4, 3)
[51, 51, 51]
-- -- -- -- -- -- -- -- -- -- -- Calculation: epw1 -- -- -- -- -- -- -- -- -- -- --
Running epw1 |████████████████████████████████████████| in 1:00.7 (0.02/s)
-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
Sanity check: Interpolated bands and phonons from EPW#
At this point we have all necessary quantities in the Wannier representation on file. As a sanity check, we perform a simple interpolation of bands and phonons to make sure that we reproduce the results found above without Wannier interpolation.
[30]:
# Electrons
Band_EPW=plot_bands.plot_band_eig(f'./{silicon.prefix}/epw/band.eig')
Band_QE=plot_bands.plot_band_scf(f'./{silicon.prefix}/bs/bs.out')
plot_bands.plot_band_prod(Band_EPW,
ef0=ef_from_file,
xlabel='Wavevector',
ylabel='Electron energy (eV)',
xticks=['L','$\Gamma$','X'],linestyle='--',color_c='b',color_v='b',first = True)
plot_bands.plot_band_prod(Band_QE,
ef0=ef_from_file,
xlabel='Wavevector',
ylabel='Electron energy (eV)',
xticks=['L','$\Gamma$','X'],first = False) # False controls if this is the first set of plots
# Phonons
PH_epw=plot_bands.plot_band_eig(f'./{silicon.prefix}/epw/phband.freq')
PH_matdyn=plot_bands.plot_band_eig(f'./{silicon.prefix}/ph/si.freq')
PH_matdyn=PH_matdyn*0.124
plot_bands.plot_band_freq(PH_epw,
xlabel='Wavevector',
ylabel='Phonon energy (meV)',
ef0=0,
xticks=['L','$\Gamma$','X'],linestyle='--',first = True,color='blue')
plot_bands.plot_band_freq(PH_matdyn,
xlabel='Wavevector',
ylabel='Phonon energy (meV)',
ef0=0,
xticks=['L','$\Gamma$','X'],first = False)
Calculation of carrier mobility#
In order to compute the carrier mobility, we perform the following operations: 1. We interpolate the electrons, phonons, and electron-phonon couplings onto a fine Brillouin zone grid (60 x 60 x 60 for electrons and phonons in this example) 2. We use these data to solve the BTE within EPW
Both steps are performed within a single call of EPW. Note the keyword temps which is used to perform calculations for multiple temperatures (300 K, …, 100 K in this example).
[33]:
#silicon.reset()
silicon.epw(epwin={'elph':'.true.',
'etf_mem': '3',
'nkf1':40,
'nkf2':40,
'nkf3':40,
'nqf1':40,
'nqf2':40,
'nqf3':40,
'mp_mesh_k':'.true.',
'efermi_read':'.true.',
'fsthick': 0.3,
'fermi_energy':6.5,
'temps':'300 250 200 150 100',
'degaussw':0.0,
'scattering':'.true.',
'int_mob':'.false.',
'carrier':'.true.',
'ncarrier' :'1E13',
'iterative_bte':'.true.',
'nstemp': 5,
'clean_transport':None},
name='epw2')
silicon.prepare(type_run='epw2')
silicon.run(16,type_run='epw2')
rm: cannot remove './epw/si.epmatkq1': No such file or directory
-- -- -- -- -- -- -- -- -- -- -- -- -- Warning -- -- -- -- -- -- -- -- -- -- -- -- -- --
Refreshing EPW input (remove refresh from epw_save.json if not needed)
-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
-- -- -- -- -- -- -- -- -- -- -- Calculation: epw2 -- -- -- -- -- -- -- -- -- -- --
Running epw2 |████████████████████████████████████████| in 6:26.9 (0.00/s)
-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
[53]:
print(os.getcwd())
silicon.reset()
silicon.epw_fold = 'epw'
silicon.epw_file = 'epw2'
temps=[300, 250, 200, 150, 100]
mob=[]
font=16
for T in temps:
silicon.epw_params['temps']=T
print(silicon.ibte_mobility_e[0,0],T)
mob.append(silicon.ibte_mobility_e[1,1])
plt.scatter(T,silicon.ibte_mobility_e[1,1], color = 'crimson')
plt.scatter(T,silicon.serta_mobility_e[1,1], color = 'royalblue')
print(mob)
plt.plot(temps[::-1],mob[::-1], color = 'crimson')
plt.xticks([100,150,200,250,300],fontsize=font)
plt.yticks([1000,3000,5000,8000],fontsize=font)
plt.xlabel('Temperature (K)',fontsize=16)
plt.ylabel('$\mu$ (cm$^2$V$^{-1}$s$^{-1}$)',fontsize=16)
plt.xticks(fontsize=16)
plt.yticks(fontsize=16)
plt.show()
/mnt/storage/sabya/For_video/epwpy/notebooks_basic/si
obtaining nscf and ph attributes
916.544 300
1411.82 250
2516.74 200
5680.29 150
18524.0 100
[1106.83, 1612.69, 2532.93, 4305.52, 7503.36]
\(\underline{\large \rm We\, now\, calculate\, the\, drift\, current\, in\, linear\, regime}\)
\(J_{\rm drift} = e n \mu E\)
[37]:
E=np.linspace(0,1,50)*1e6 # V/cm
n=1e15 # cm^-3
e = 1.6*10**-19 # Coulomb
J = []
for mu in mob:
J.append(e*n*mu*E) # Ampere/cm^2
for i,T in enumerate(temps):
print(T)
plt.plot(E,J[i])
plt.text(E[-15],J[i][-15],str(T)+' K')
plt.xlabel('$F \\rm (V/cm)$',fontsize=20)
plt.ylabel('$J \\rm (A/cm^2)$',fontsize=20)
plt.show()
300
250
200
150
100
\(\newcommand{\mrm}[1]{\mathrm{#1}}\) \(\newcommand{\vec}[1]{{\mathbf{\mathrm{#1}}}}\) \(\newcommand{\Ver}[1]{\)V_:nbsphinx-math:mathrm{er}\(}\) $:nbsphinx-math:def>{rangle} $ $:nbsphinx-math:def<{langle} $ \(\def\bk{{\bf k}}\) \(\def\bq{{\bf q}}\) \(\def\be{{\bf e}}\) \(\def\bv{{\bf v}}\) \(\def\a{\alpha}\) \(\def\d{\delta}\) \(\def\w{\omega}\) \(\def\ve{\varepsilon}\) \(\def\hc{{\hat{c}^{\phantom{\dagger}}}}\) \(\def\hcd{{\hat{c}^\dagger}}\) \(\def\hH{{\hat{H}}}\) \(\def\hp{{\hat{\bf p}}}\) \(\def\ha{{\hat{a}^{\phantom{\dagger}}}}\) \(\def\had{{\hat{a}^\dagger}}\) \(\def\Vep{\hat{V}_{\rm ep}}\) \(\def\Ver{\hat{V}_{\rm er}}\)
\(\underline{\large{\rm Plotting\,inverse\,relaxation\,time}}\)
We next obtain the relaxation time as a function of energy which is defined as,
\(\frac{1}{\tau_{k}}\)=\(\big|\frac{1}{N_\bq}\sum_{\bq}\langle \psi_{\bk+\bq}|\hat{V}_{\rm ep}|\psi_{bk}\rangle\big|^2\)
[44]:
silicon.temp=400
tau_inv=silicon.inv_taucb
plt.scatter(tau_inv[:,3]-6.6,tau_inv[:,4])
plt.yscale('log')
plt.xlabel('Energy eV')
plt.ylabel('$\\tau$ (ps$^{-1}$)')
plt.axis([0,0.2,1e-2,1e3])
plt.show()
[ ]: