Polarons¶
Note
Hands-on based on Quantum ESPRESSO v7.6 and EPW v6.1
Introduction¶
In this session, you will learn how to perform calculations of self-trapped polarons using the EPW code.
For the description for all input flags please follow the link:
https://docs.epw-code.org/Inputs/Inputs.html
Theoretical background¶
In this short introduction, we briefly present the main concepts and equations relevant to this tutorial. A complete description of the theory and computational method can be found in Phys. Rev. Lett. **122**, 246403 (2019) and Phys. Rev. B **99**, 235139 (2019).
The ground state wave function \(\psi(\mathbf{r})\) and atomic displacements \(\Delta \tau_{\kappa\alpha p}\) forming a polaron can be found by minimizing the total DFT energy functional of an excess electron added to a crystal, which translates into the solution of the following coupled system of equations:
In these expressions, \(\tau_{\kappa\alpha p}\) represents the Cartesian coordinate of the atom \(\kappa\) in the unit cell \(p\) along the direction \(\alpha\), \(C^{0}_{\kappa\alpha p,\kappa'\alpha'p'}\) is the matrix of interatomic force constants, and \(\hat{H}_{\mathrm{KS}}^{0}\) and \(V_{\mathrm{KS}}^{0}\) represent the Kohn-Sham Hamiltonian and the self-consistent potential, respectively. The superscript \(^{0}\) indicates that the quantities are evaluated in the ground state of the pristine system without extra electron. The real space integral is performed over a Born-von Karman (BvK) supercell of the crystal containing \(N_p\) unit cells. We will refer to \(\varepsilon\) as the polaron eigenvalue.
By expanding the polaron wave function in terms of the single-particle Kohn-Sham states \(\psi_{n\mathbf{k}}\) with eigenvalues \(\varepsilon_{n\mathbf{k}}\):
and the atomic displacements in terms of the phonon eigenmodes \(e_{\kappa\alpha,\nu}(\mathbf{q})\) with frequencies \(\omega_{\mathbf{q}\nu}\):
where \(M_\kappa\) is the mass of the atom \(\kappa\) and \(\mathbf{R}_p\) is the lattice vector of the unit cell \(p\), we can transform Eq. (1) and (2) into a coupled set of equations for the expansion coefficients in reciprocal space:
In these expressions, \(\varepsilon_{n\mathbf{k}}\) are the Kohn-Sham eigenvalues, \(\omega_{\mathbf{q}\nu}\) are the phonon frequencies and \(g_{mn\nu}(\mathbf{k},\mathbf{q})\) are the electron-phonon matrix elements. All these quantities can be obtained by standard DFT and DFPT calculations.
Eq. (5) and (6) can be rewritten more conveniently as follows:
where
In practice, Eqs. (eq:eq4)-((8)) are solved iteratively until self-consistency is reached. The final atomic displacements associated with the polaron will be obtained from Eq. (4), and the wave function in real space can be conveniently obtained from the maximally localized wannier functions \(w\) as:
having defined
and
where \(U^{\dagger}_{mn\mathbf{k}}\) is the unitary matrix that generates the smooth Bloch gauge (see e.g. Rev. Mod. Phys. **84**, 1419 (2012)).
The polaron formation energy \(\Delta E_{f}\), defined as the energy required to trap a conduction band state with eigenvalue \(\varepsilon_{\mathrm{CBM}}\) into a localized polaron state, can be obtained from the expansion coefficients that solve Eq. (5) and (6) by:
We will refer to the first and second terms on the right hand side as the electron and phonon parts of the formation energy, respectively.
From Eq. (5) and (6), we observe that the quantities required to solve the polaron equations are the Kohn-Sham eigenvalues, phonon frequencies, and electron-phonon matrix elements on the \(\mathbf{k}\)- and \(\mathbf{q}\)-point grids corresponding to the equivalent BvK supercell in which Eq. (1) and (2) are defined. To obtain the formation energy of an isolated polaron, we solve Eq. (5) and (6) on progressively denser grids and extrapolate the results to the infinite-supercell limit.
Exercise 1¶
In this exercise you will calculate and analyze the formation of hole polarons in LiF.
You are advised to run the calculations in your scratch folder ($SCRATCH). To do so, navigate to your scratch directory, download the tar file which contains all necessary files for this tutorial, and extract it:
$ cd $SCRATCH
$ wget --content-disposition \
"https://drive.google.com/uc?export=download&id=1ZrlYL8Yh6Kib9dyS8q--hCurBnfFba-L" .
$ tar -xvf Wed.4.LafuenteBartolome.tar; cd Wed.4.LafuenteBartolome
A quick inspection of the directory shows the four subdirectories containing the files required for the four exercises in this hands-on.
$ ls
exercise1 exercise2 exercise3 exercise4
Begin by navigating to the corresponding directory:
$ cd exercise1
The initial steps of the calculation are similar to the ones followed in the previous tutorials. First, we need to perform an electronic ground state calculation using Quantum ESPRESSO, and a phonon calculation using the PHonon code. The corresponding input files (lif.scf.in and lif.ph.in) are the following:
lif.scf.in
&control
calculation = 'scf'
prefix = 'lif'
pseudo_dir = './'
outdir = './'
tprnfor = .true.
tstress = .true.
/
&system
ibrav = 2
celldm(1) = 7.19084
nat = 2
ntyp = 2
ecutwfc = 60
/
&electrons
conv_thr = 1.0d-16
/
ATOMIC_SPECIES
Li 6.941 Li.upf
F 18.9984 F.upf
ATOMIC_POSITIONS crystal
Li 0.0000 0.0000 0.0000
F 0.5000 0.5000 0.5000
K_POINTS automatic
8 8 8 0 0 0
lif.ph.in
&inputph
prefix = 'lif'
outdir = './'
epsil = .true.
zeu = .true.
ldisp = .true.
fildyn = 'lif.dyn'
fildvscf = 'dvscf'
tr2_ph = 1.0d-16,
nq1 = 4,
nq2 = 4,
nq3 = 4
/
You can use the following submission job (submit1.sh) to perform the calculation in Frontera: TO BE UPDATED ONCE WE HAVE THE RESERVATION CODE ETC.
submit1.sh
#!/bin/bash
#SBATCH -J myjob # Job name
#SBATCH -N 1 # Total # of nodes
#SBATCH --ntasks-per-node 32
#SBATCH -t 01:00:00 # Run time (hh:mm:ss)
#SBATCH -p small
#SBATCH --reservation=XXX
# echo loaded modules, current directory, and starting time
module list
pwd
date
# Export the path which contains executable file
export PATHQE=/work2/06868/giustino/SCHOOL/q-e-qe-7.6
# Launch MPI code...
# Total # of parallel tasks
MPIOPT="-np "$SLURM_NTASKS
# kpoint parallel groups
KPTPRL="-npool 4"
# run jobs
# scf calculation
ibrun ${MPIOPT} ${PATHQE}/bin/pw.x ${KPTPRL} -input lif.scf.in > lif.scf.out
# ph calculation
ibrun ${MPIOPT} ${PATHQE}/bin/ph.x ${KPTPRL} -input lif.ph.in > lif.ph.out
# echo finishing time
date
Action. Run a self-consistent calculation and a phonon calculation on a homogeneous \(8\times8\times8\) \(\mathbf{k}\) and \(4\times4\times4\) \(\mathbf{q}\)-point grid.
$ sbatch submit1.sh
Note1: These coarse \(\mathbf{k}\)- and \(\mathbf{q}\)-point grids and the plane-wave expansion cutoff have been chosen to provide a reasonable balance between accuracy and computational cost. Fully converged calculations require denser coarse \(\mathbf{k}\) and \(\mathbf{q}\) grids, as well as a larger cutoff.
The keyword fildvscf instructs the code to write the change of the self-consistent potential due to phonon perturbations, \(\partial_{\mathbf{q}\nu}V^{\rm scf}\), which is needed to compute the electron-phonon matrix elements. In the output file lif.dyn0, you can find the list of 16 irreducible \(\mathbf{q}\) points in the Brillouin Zone (IBZ). If you type ls, you can see the lif.dynX files containing the dynamical matrix for each irreducible \(\mathbf{q}\) point. The dvscf files are all named lif.dvscf1 and are located inside the _ph0/lif.q_X/ folders, except for the one corresponding to the first \(\mathbf{q}\) point (\(\Gamma\)) that is located in _ph0/.
Action. Gather the .dyn and .dvscf files into a new save/ directory which EPW will read.
The files in _ph0/lif.phsave/ containing the displacement patterns are also needed. This can easily be done using the pp.py python script which is included in the EPW distribution:
pp.py
$ python3 /work2/06868/giustino/SCHOOL/q-e-qe-7.6/EPW/bin/pp.py
The script will ask you to prompt the prefix of your calculation (in this case lif): {.. code-block:: text
Enter the prefix used for PH calculations (e.g. diam) lif
}
Action. Run a non self-consistent calculation on a full homogeneous \(\mathbf{k}\)-point grid, and an EPW calculation to obtain the electron-phonon matrix elements in the Wannier representation.
To proceed with the calculation, we first construct the Wannier functions associated with the valence-band manifold using wannier90 as an internal library. We then compute the electron-phonon matrix elements and transform them to the Wannier representation using EPW.
Begin by performing a non self-consistent calculation on a homogeneous \(4\times4\times4\) \(\mathbf{k}\)-point grid (lif.nscf.in):
lif.nscf.in
&control
calculation = 'bands'
prefix = 'lif'
pseudo_dir = './'
outdir = './'
/
&system
ibrav = 2
celldm(1) = 7.19084
nat = 2
ntyp = 2
ecutwfc = 60
nbnd = 15
/
&electrons
conv_thr = 1.0d-16
/
ATOMIC_SPECIES
Li 6.941 Li.upf
F 18.9984 F.upf
ATOMIC_POSITIONS crystal
Li 0.0000 0.0000 0.0000
F 0.5000 0.5000 0.5000
K_POINTS crystal
64
0.00000000 0.00000000 0.00000000 1.562500e-02
0.00000000 0.00000000 0.25000000 1.562500e-02
0.00000000 0.00000000 0.50000000 1.562500e-02
0.00000000 0.00000000 0.75000000 1.562500e-02
0.00000000 0.25000000 0.00000000 1.562500e-02
...
Note: The \(\mathbf{k}\)-point list under K_POINTS crystal in the file above is for illustrative purposes and is not complete. The complete positive-definite homogeneous \(4\times4\times4\) \(\mathbf{k}\)-point grid between 0 and 1 can be generated by using the script kmesh.pl included in the wannier90 package:
kmesh.pl
$ /work2/06868/giustino/SCHOOL/q-e-qe-7.6/external/wannier90-3.1.0/utility/kmesh.pl 4 4 4
You can find the complete list in the input file lif.nscf.in.
Next, prepare the EPW input file to calculate the electron-phonon matrix elements in the Wannier representation (lif.epw1.in):
lif.epw1.in
&inputepw
prefix = 'lif'
outdir = './'
elph = .true.
epwwrite = .true.
lpolar = .true.
nbndsub = 3
dvscf_dir = './save/'
etf_mem = 0
bands_skipped = 'exclude_bands = 1:2, 6:15'
wannierize = .true.
num_iter = 500
iprint = 2
proj(1) = 'F:p'
wannier_plot = .true.
wannier_plot_supercell = 4 4 4
nk1 = 4
nk2 = 4
nk3 = 4
nq1 = 4
nq2 = 4
nq3 = 4
band_plot = .true.
filkf = './path.kpt'
filqf = './path.kpt'
/
As shown in the input file, we consider three Wannier functions associated with the three isolated valence bands below the Fermi level. Consistent with their orbital character, the initial projections are chosen as three \(p\) orbitals centered on the F atom.
You can now submit the calculation:
$ sbatch submit2.sh
Note: The submission script submit2.sh can be obtained from the previous submit1.sh by replacing the two lines containing ibrun by the following ones:
submit1.sh
# nscf calculation
ibrun ${MPIOPT} ${PATHQE}/bin/pw.x ${KPTPRL} -input lif.nscf.in > lif.nscf.out
# epw calculation
ibrun ${MPIOPT} ${PATHQE}/bin/epw.x -npool ${SLURM_NTASKS} -input lif.epw1.in > lif.epw1.out
Action. Interpolate the electron-phonon matrix elements to a fine \(\mathbf{k}\)- and \(\mathbf{q}\)-point grid and solve the self-consistent polaron equations.
Begin by preparing the EPW input file to perform the interpolation and solve the polaron equations (lif.epw2.in):
lif.epw2.in
&inputepw
prefix = 'lif'
outdir = './'
elph = .true.
epwread = .true.
lpolar = .true.
nbndsub = 3
dvscf_dir = './save/'
etf_mem = 0
wannierize = .false.
plrn = .true.
restart_plrn = .false.
type_plrn = 1
init_plrn = 1
init_sigma_plrn = 1.0
niter_plrn = 500
conv_thr_plrn = 1E-4
ethrdg_plrn = 1E-5
nk1 = 4
nk2 = 4
nk3 = 4
nq1 = 4
nq2 = 4
nq3 = 4
nkf1 = 6
nkf2 = 6
nkf3 = 6
nqf1 = 6
nqf2 = 6
nqf3 = 6
/
You can observe that we use a \(6\times 6\times 6\) fine \(\mathbf{k}\)- and \(\mathbf{q}\)-point grid in this example. The tag plrn = .true. activates the polaron calculations, and the tag type_plrn = 1 indicates that a hole polaron is to be calculated. The self-consistent polaron equations will be initialized with a Gaussian of width \(\sigma=10.0 \mathrm{bohr}\) (init_sigma_plrn), and the solution will be considered converged when the greatest difference in the magnitude of the atomic displacements between two subsequent iterations is lower than \(10^{-4} \mathrm{bohr}\) (conv_thr_plrn). The input variable ethrdg_plrn controls the convergence threshold (Ry) used in the iterative diagonalization at each step in the self-consistent solution of the polaron equations.
You can now submit the calculation:
$ sbatch submit3.sh
Note: The submission script submit3.sh can be obtained from the previous submit2.sh by replacing the two lines containing ibrun by the following one:
submit2.sh
ibrun ${MPIOPT} ${PATHQE}/bin/epw.x -npool ${SLURM_NTASKS} -input lif.epw2.in > lif.epw2.out
You can analyze the iterative self-consistent process in the output:
{.. code-block:: text
- iter Eigval/eV Phonon/eV Electron/eV Formation/eV Error/eV
1 0.4145E+00 -0.1735E+00 -0.3746E+00 0.2011E+00 0.1467E+00 2 0.1747E+01 -0.9373E+00 -0.7438E+00 -0.1935E+00 0.2134E+00 3 0.2662E+01 -0.1687E+01 -0.8773E+00 -0.8101E+00 0.1456E+00 4 0.2839E+01 -0.1857E+01 -0.9008E+00 -0.9566E+00 0.2514E-01 5 0.2868E+01 -0.1883E+01 -0.9057E+00 -0.9776E+00 0.7096E-02
… 100 0.2938E+01 -0.1919E+01 -0.8998E+00 -0.1019E+01 0.1103E-03 101 0.2938E+01 -0.1919E+01 -0.8998E+00 -0.1019E+01 0.1023E-03 102 0.2938E+01 -0.1919E+01 -0.8998E+00 -0.1019E+01 0.1086E-03
End of self-consistent cycle
}
All the information related to the energetics of the converged polaron solution is printed at the end of the self-consistent process:
{.. code-block:: text
Eigenvalue (eV): 2.9379456
Phonon part (eV): -1.9189683
Electron part (eV): 0.8998334
Formation Energy (eV): -1.0191349
}
The breakdown of the formation energy in the different terms is briefly explained in the hyperref[sec:intro]{introduction}. More details and further discussion can be found in Phys. Rev. B **99**, 235139 (2019).
A quick inspection of the output files will show that several files have been written with the .plrn extension. The most relevant output files for this exercise are Amp.plrn and dtau.plrn, which contain the polaron wave function coefficients in the Wannier basis and the atomic displacements, respectively. These files are needed to postprocess and analyze the polaron solution in more detail.
Action. Plot and analyze the polaron wave function.
The post-processing of the polaron solution can be activated by the following input file (lif.epw3.in):
lif.epw3.in
&inputepw
prefix = 'lif'
outdir = './'
elph = .true.
epwread = .true.
lpolar = .true.
nbndsub = 3
dvscf_dir = './save/'
etf_mem = 0
wannierize = .false.
plrn = .true.
restart_plrn = .true.
type_plrn = 1
cal_psir_plrn = .true.
interp_Ank_plrn = .true.
interp_Bqu_plrn = .true.
filkf = './path.kpt'
filqf = './path.kpt'
nk1 = 4
nk2 = 4
nk3 = 4
nq1 = 4
nq2 = 4
nq3 = 4
nkf1 = 6
nkf2 = 6
nkf3 = 6
nqf1 = 6
nqf2 = 6
nqf3 = 6
/
The input tag restart_plrn = .true. skips the self-consistent solution of the polaron equations and triggers the post-processing of the previously converged solution by reading the Amp.plrn and dtau.plrn files.
The input tag cal_psir_plrn = .true. calls the subroutine to plot the polaron wave function in real space, which will be printed in the psir_plrn.xsf file. This file can be plotted with VESTA. First submit the calculation:
$ sbatch submit4.sh
and then copy the file containing the wave function to your own computer:
$ scp username@frontera.tacc.utexas.edu:path_to_exercise/psir_plrn.xsf .
As when logging in, enter your TACC password at the password prompt. At the TACC Token prompt, enter your 6-digit code. Then you can visualize the polaron wave function by running VESTA (https://jp-minerals.org/vesta/en/download.html) locally on your computer:
(local)$ VESTA psir_plrn.xsf
You should obtain something similar to this:
The polaron wave function is plotted as an isosurface, and the red arrows indicate the magnitude and direction of the atomic displacements in the polaronic configuration. You can observe that the \(6\times 6\times 6\) \(\mathbf{k}\)-point grid transforms to a \(6\times 6\times 6\) supercell in real space, with periodic boundary conditions. From this figure it is also appreciated that holes form small polarons in LiF, localized in one unit cell, and have the shape of a \(p\)-orbital centered on a F-atom.
Next we analyze the polaron wave function coefficients in the single-particle basis of the DFT-Kohn Sham states. This calculation is activated by the interp_Ank_plrn = .true. tag. Similarly, the calculation of the expansion coefficients of the atomic displacements in the phonon basis is activated by interp_Bqu_plrn = .true. . The files containing the \(\mathbf{k}\)- and \(\mathbf{q}\)-points in which the polaron coefficients will be interpolated have to be indicated in the filkf and filqf input tags, respectively. In this example we use a path along the high-symmetry W-L-\(\Gamma\)-X-W-K points, as contained in the path.kpt file.
The output files containing the interpolated coefficients are Ank.band.plrn and Bmat.band.plrn. You can visualize your results by using, for example, the following script (plot_Ank.py):
plot_Ank.py
## Script to plot Ank coefficients
import matplotlib as mpl
mpl.use('pdf')
import matplotlib.pyplot as plt
import numpy as np
# change font on mathematical expressions on plots
mpl.rcParams['mathtext.fontset'] = 'cm'
# Read kpath file
kpath = np.loadtxt('path.kpt', skiprows=1)
# Read Ank file
ik, ibnd, ek0, ReAnk, ImAnk, AbsAnk = np.loadtxt('Ank.band.plrn',
unpack=True, skiprows=1)
# Separate data in different bands
nbnd=int(max(ibnd))
maxik=int(max(ik))
Ak=np.zeros((maxik,nbnd))
ek=np.zeros((maxik,nbnd))
iklist=np.zeros((maxik,nbnd))
for i in range(maxik):
for ibnd in range(nbnd):
Ak[i][ibnd]=AbsAnk[i*nbnd+ibnd]
ek[i][ibnd]=ek0[i*nbnd+ibnd]
iklist[i][ibnd]=i
# Get vbm and bandwidth
vbm=max(ek[:,nbnd-1])
bandwidth=vbm-min(ek[:,0])
## Plot bands
f, ax = plt.subplots(figsize=(10,6))
ax.plot(iklist, ek, color='blue')
# Plot Ank
ax.scatter(iklist, ek, 100*Ak, color='gold', edgecolors='gray', alpha=0.8,
label=r'$|A_{n\mathbf{k}}|$')
# Set tick params etc.
ax.set_ylim(-bandwidth-0.1*bandwidth,0.1*bandwidth)
ax.set_xlim((min(iklist[:,0]),max(iklist[:,0])))
ax.set_xticklabels([])
ax.tick_params(axis='x', color='black', labelsize='0', pad=0, length=0, width=0)
ax.tick_params(axis='y', color='black', labelsize='18', pad=5, length=5, width=1)
ax.set_ylabel(r'$E-E_{\mathrm{VBM}} ~ (\mathrm{eV})$', fontsize=25, labelpad=10)
ax.legend(loc='upper right', fontsize=25)
plt.savefig('Ank.pdf')
#plt.show()
Note: You can find a similar script to plot the \(B_{\mathbf{q}\nu}\) coefficients in the example directory (plot_Bqv.py).
Executing these scripts:
$ python3 plot_Ank.py
$ python3 plot_Bqv.py
should generate two .pdf output files (Ank.pdf and Bqv.pdf) containing the corresponding figures. You can visualize the plots in Stampede3 if you have accessed with X11 forwarding:
$ evince Ank.pdf
$ evince Bqv.pdf
Note: Alternatively, you can copy the output files to your own computer (enter your TACC password and token as before):
$ scp username@frontera.tacc.utexas.edu:path_to_exercise/*.pdf .
and open them with your preferred .pdf file reader.
Your figures should be similar to:
As you can observe in these results, due to the strongly localized nature of the hole polaron in LiF, its wave function coefficients in the Bloch basis are spread across the entire Brillouin zone (see Eq. (3) and (4) in the hyperref[sec:intro]{introduction}).
Action. Investigate the dependence of the polaron formation energy on the supercell size, and extrapolate the results to the infinite-supercell limit corresponding to an isolated polaron.
In the previous calculation we have seen that the \(6\times 6\times 6\) \(\mathbf{k}\)- and \(\mathbf{q}\)-point grid in momentum space corresponds to a \(6\times 6\times 6\) BvK supercell in real space. Due to the interaction between replicated polarons in the neighbouring supercells, the formation energy of the polaron will depend on the momentum-grid (supercell) size.
You can analyze the dependence of the polaron formation energy on the supercell size by using the following script (submit5.sh):
submit5.sh
#!/bin/bash
#SBATCH -J myjob # Job name
#SBATCH -N 1 # Total # of nodes
#SBATCH --ntasks-per-node 32
#SBATCH -t 01:00:00 # Run time (hh:mm:ss)
#SBATCH -p small
#SBATCH --reservation=XXX
# echo loaded modules, current directory, and starting time
module list
pwd
date
# export the path which contains executable file
export PATHQE=/work2/05193/sabyadk/stampede3/EPWSchool2024/q-e
# Launch MPI code...
echo "# nk Eform" > "E_vs_nk.dat"
for i in 8 10 12 14 16
do
sed -e "32s/.*/ nkf1 = $i/" \
-e "33s/.*/ nkf2 = $i/" \
-e "34s/.*/ nkf3 = $i/" \
-e "35s/.*/ nqf1 = $i/" \
-e "36s/.*/ nqf2 = $i/" \
-e "37s/.*/ nqf3 = $i/" \
"lif.epw4.in" > "lif.epw4.$i.in"
# run jobs
j=$(( 4*$i ))
ibrun -np $j ${PATHQE}/bin/epw.x -npool $j -input lif.epw4.$i.in > lif.epw4.$i.out
grep 'Formation Energy (eV):' "lif.epw4.$i.out" >> "E_vs_nk.dat"
sed -i "s/.Formation Energy (eV):/ $i/" "E_vs_nk.dat"
done
# echo finishing time
date
You can submit this job by:
$ sbatch submit5.sh
For each of the desired \(\mathbf{k}\)-point grids (in this case \(8^{3}\), \(10^{3}\), \(12^{3}\), \(14^{3}\), and \(16^{3}\)), this script creates a copy of the input file lif.epw4.in with modified nkf and nqf input variables into lif.epw4.*.in, launches the calculation, and copies the resulting formation energy to the E_vs_nk.dat file:
{.. code-block:: text
- # nk Eform
8 -1.2050655 10 -1.3023245 12 -1.3542913 14 -1.3927680 16 -1.4178908
}
You can plot the formation energy as a function of the inverse supercell size (here defined as \(L^{-1}\) where \(L^{3}\) is the supercell volume) and extrapolate to the infinite supercell by using, for example, the following script (plot_E_s_nk.py):
plot_E_s_nk.py
import matplotlib as mpl
mpl.use('pdf')
import matplotlib.pyplot as plt
from matplotlib.ticker import MultipleLocator,FormatStrFormatter
import numpy as np
# Read E vs nk
Nk, Eform = np.genfromtxt('E_vs_nk.dat', unpack=True)
# Set unit cell volume to convert Nk to inverse supercell size
ucell_volume = 92.9563
# Start figure
fig, ax = plt.subplots(figsize=(12,8))
# Plot Eform
ax.scatter(1/(Nk*ucell_volume**(1/3)), Eform, s=50, marker='o',
color='darkred', edgecolors='black')
# Perform linear fit and plot
mf, bf = np.polyfit(1/(Nk*ucell_volume**(1/3)), Eform, 1)
xlist = np.linspace(0.0, np.max(1/(Nk*ucell_volume**(1/3))), 100)
print("Extrapolation to isolated polaron formation energy = ", "%.3f" % bf, "eV")
ax.plot(xlist, mf*xlist+bf, '--', color='gray')
# Set tick params etc.
ax.set_xlabel(r'Inverse supercell size ($\mathrm{bohr}^{-1}$)',fontsize=20)
ax.set_ylabel('Formation energy (eV)',fontsize=20, labelpad=5)
ax.tick_params(axis='x', color='black', labelsize='20', pad=5, length=5, width=2)
ax.tick_params(axis='y', color='black', labelsize='20', pad=5, length=5, width=2,
right=True)
ax.yaxis.set_minor_locator(MultipleLocator(0.1))
ax.tick_params(axis='y', which='minor', color='black', labelsize='20', pad=5,
length=3, width=1.2, right=True)
ax.set_xlim(0.0, np.max(1/(Nk*ucell_volume**(1/3)))+0.01)
ax.set_title('LiF hole polaron', fontsize=20)
plt.savefig("E_vs_nk.pdf")
#plt.show()
Executing this script :
$ python3 plot_E_vs_nk.py
and opening the .pdf file:
$ evince E_vs_nk.pdf
You should obtain something similar to this:
In this example, the extrapolation to infinite supercell gives a formation energy of \(\Delta E_{f}=-1.635 \mathrm{eV}\) for the isolated hole polaron in LiF. Note that the present results are not fully converged (see Phys. Rev. B **99**, 235139 (2019) for fully converged calculation parameters).
Exercise 2¶
In this exercise you will repeat a similar procedure as before, but this time to calculate and analyze the formation of electron polarons in LiF.
Start by navigating to the directory for this exercise, and creating symbolic links to the lif.save and save directories generated in the previous exercise, so that we do not need to run the scf, ph, and nscf calculations again:
$ cd ../exercise2/
$ ln -s ../exercise1/lif.save/ .
$ ln -s ../exercise1/save/ .
You can follow the calculation steps one by one as in Exercise 1, and analyze the intermediate results. Alternatively, you can submit the following script that will perform all the calculations automatically (submit1.sh):
submit1.sh
#!/bin/bash
#SBATCH -J myjob # Job name
#SBATCH -N 1 # Total # of nodes
#SBATCH --ntasks-per-node 32
#SBATCH -t 01:00:00 # Run time (hh:mm:ss)
#SBATCH -p small
#SBATCH --reservation=XXX
# echo loaded modules, current directory, and starting time
module list
pwd
date
# Export the path which contains executable file
export PATHQE=/work2/06868/giustino/SCHOOL/q-e-qe-7.6
# run jobs
# epw1: g(Re,Rp)
ibrun -np ${SLURM_NTASKS} ${PATHQE}/bin/epw.x -npool ${SLURM_NTASKS} -input lif.epw1.in > lif.epw1.out
# epw2: polaron for different nk
echo "# nk Eform" > "E_vs_nk.dat"
for i in 6 8 10 12 14 16
do
sed -e "33s/.*/ nkf1 = $i/" \
-e "34s/.*/ nkf2 = $i/" \
-e "35s/.*/ nkf3 = $i/" \
-e "36s/.*/ nqf1 = $i/" \
-e "37s/.*/ nqf2 = $i/" \
-e "38s/.*/ nqf3 = $i/" \
"lif.epw2.in" > "lif.epw2.$i.in"
# run jobs
j=$(( 2*$i ))
ibrun -np $j ${PATHQE}/bin/epw.x -npool $j -input lif.epw2.$i.in > lif.epw2.$i.out
grep 'Formation Energy (eV):' "lif.epw2.$i.out" >> "E_vs_nk.dat"
sed -i "s/.Formation Energy (eV):/ $i/" "E_vs_nk.dat"
done
Note: Analyze the input file lif.epw2.in. You will see that in this case we consider one conduction band for the Wannierization, and that the electron polaron calculation is set by the type_plrn = -1 input tag.
Action. Study the change of the electron polaron formation energy as a function of the supercell size.
You can submit the calculation with:
$ sbatch submit1.sh
Note: This calculation is computationally more demanding and should take around 10 minutes to complete.
The dependence of the polaron formation energy with respect to the momentum grid (BvK supercell) size in this case is (E_vs_nk.dat):
{.. code-block:: text
# nk Eform 6 0.0000046 8 0.0000070 10 0.0000667 12 -0.0078670 14 -0.0360222 16 -0.0580386
}
Exercise: Why do we obtain negative formation energies only for supercell sizes larger than \(10\times 10 \times 10\)?
As in the previous exercise, we can postprocess and analyze the converged polaron solution by studying its wave function in real space and the distribution of the expansion coefficients in reciprocal space. We will analyze the polaron solution obtained in the last \(16\times 16\times 16\) \(\mathbf{k}\)-point grid calculation:
$ sbatch submit2.sh
Note: The calculation of the polaron wave function in the large 16:math:times`16:math:times`16 real space supercell should take around 5 min to complete.
Note: You can obtain the lif.epw3.in input file by modifying lif.epw2.16.in in a similar way as we have done before for lif.epw3.in in Exercise 1. The input tag step_wf_grid_plrn = 2 reduces the amount of real-space points in which the polaron wave function is computed, and it is only used here to reduce the computational time. You can obtain the submit2.sh submission file by copying submit4.sh from ../exercise1/
The distribution of the polaron wave function coefficients along the high-symmetry W-L-\(\Gamma\)-X-W-K points should be similar to:
Note: You can obtain plot the polaron wave function coefficients by copying plot_Ank.py and plot_Bqv.py from ../exercise1/, and modifying the following lines in plot_Ank.py: ax.set_ylim(-0.1*bandwidth,bandwidth+0.1*bandwidth) ax.set_ylabel(r':math:`E-E_{\mathrm{CBM}} (\mathrm{eV})`', fontsize=25, labelpad=10)
And the real space plot of the polaron wave function should look similar to:
As you can observe in these results, electrons in LiF form large polarons in real space, which correspond to wave function coefficients localized around the conduction band bottom.
Exercise 3¶
In this exercise, we illustrate how to perform polaron calculations in non-diagonal supercells with any shape.
We will take the hole polaron in rocksalt MgO as an example.
Background and notation¶
We denote by \(\mathbf{a}_{s1}\), \(\mathbf{a}_{s2}\), \(\mathbf{a}_{s3}\) the primitive lattice vectors of a BvK supercell, and by \(\mathbf{a}_{p1}\), \(\mathbf{a}_{p2}\), \(\mathbf{a}_{3}\) the primitive lattice vectors of the crystal unit cell. These vectors must be related as follows:
where the matrix elements \(S_{ij}\) are integers. The transformation between the corresponding reciprocal lattice vectors is:
where \(\mathbf{b}_{s1}\), \(\mathbf{b}_{s2}\), \(\mathbf{b}_{s3}\) are the primitive reciprocal lattice vectors of the BvK supercell, and \(\mathbf{b}_{p1}\), \(\mathbf{b}_{p2}\), \(\mathbf{b}_{3}\) are the primitive vectors of the reciprocal lattice of the crystal unit cell. The matrix elements of these two transformations are related by \(\bar{S}=(S^{-1})^{T}\).
When the off-diagonal elements of the matrix \(S\) in Eq. (10) are non-zero, one obtains a so-called ``non-diagonal’’ BvK supercell. A non-diagonal supercell differs from a standard supercell in that it does not have the same shape as the primitive unit cell. This observation was exploited in Phys. Rev. B **92**, 184301 (2015) to generate supercells for computing phonons at desired \(q\)-vectors from finite-difference calculations. Here, we employ the same strategy but in reverse: first we choose the shape of the BvK supercell, then we use Eq. (11) to determine the Brillouin zone grid that we need in Eq. (5)-(6) to obtain polarons in such a supercell. A more detailed discussion on polaron calculations in non-diagonal supercells can be found at Proc. Natl. Acad. Sci. USA 121, e2318151121 (2024).
Example: FCC primitive to conventional cell¶
The lattice vectors of the MgO primitive cell are:
The lattice vectors corresponding to a conventional cubic supercell of MgO are given by:
where \(N\) is a positive integer. Such a supercell is obtained by using Eq. (10) with the transformation matrix:
which contains \(|\mathrm{det} S| = 4N\) crystal unit cells. The inverse transformation \(\bar S\) is given by:
The \(k\)-point (and \(q\)-point) grids to be used in the solution of (5) and (6) are given by the points \(\mathbf{k} = i\mathbf{b}_{s1}+j\mathbf{b}_{s2}+k\mathbf{b}_{s3}\) that fall within the first Brillouin zone of the crystal unit cell, being \(i\), \(j\), \(k\) integer numbers and \(\mathbf{b}_{s1}\), \(\mathbf{b}_{s2}\), and \(\mathbf{b}_{s3}\) the BvK supercell reciprocal lattice vectors obtained from (11).
Action. Calculate the polaron ground state.
We begin by calculating the formation energy and atomic displacements of the polaron ground state. To do so, we follow a workflow similar to that used in Exercise 1: ground-state DFT electronic-structure calculation with PW; phonon calculations within DFPT using PH; calculation of the electron-phonon matrix elements on the coarse grid and transformation to the Wannier representation using EPW; interpolation onto a fine grid and self-consistent solution of the polaron equations with EPW; post-processing and visualization of the wave function.
All these calculations can be run using the following script (submit1.sh):
submit1.sh
#!/bin/bash
#SBATCH -J myjob # Job name
#SBATCH -N 1 # Total # of nodes
#SBATCH --ntasks-per-node 48
#SBATCH -t 01:00:00 # Run time (hh:mm:ss)
##SBATCH -A EPSchool2026
#SBATCH -A DMR21002
#SBATCH -p development
##SBATCH -p small
##SBATCH --reservation=XXX
# echo loaded modules, current directory, and starting time
module list
pwd
date
# Export the path which contains executable file
##export PATHQE=/work2/06868/giustino/SCHOOL/q-e-qe-7.6
export PATHQE=/work2/11514/acarrasc/frontera/q-e
# Total # of parallel tasks
MPIOPT="-np "$SLURM_NTASKS
# kpoint parallel groups
KPTPRL="-npool 4"
# Launch MPI jobs
# scf calculation
ibrun ${MPIOPT} $PATHQE/bin/pw.x ${KPTPRL} -input mgo.scf.in > mgo.scf.out
# ph calculation
ibrun ${MPIOPT} $PATHQE/bin/ph.x ${KPTPRL} -input mgo.ph.in > mgo.ph.out
# Convert ph output to EPW-readable format
printf "mgo\n" | python3 $PATHQE/EPW/bin/pp.py > pp.out
# nscf calculation
ibrun ${MPIOPT} $PATHQE/bin/pw.x ${KPTPRL} -input mgo.nscf.in > mgo.nscf.out
# epw calculation
ibrun ${MPIOPT} $PATHQE/bin/epw.x -npool ${SLURM_NTASKS} -input mgo.epw1.in > mgo.epw1.out
ibrun ${MPIOPT} $PATHQE/bin/epw.x -npool ${SLURM_NTASKS} -input mgo.epw2.in > mgo.epw2.out
ibrun ${MPIOPT} $PATHQE/bin/epw.x -npool ${SLURM_NTASKS} -input mgo.epw3.in > mgo.epw3.out
# echo finishing time
date
Navigate to the exercise directory and submit the calculation:
$ cd ../exercise3/
$ sbatch submit1.sh
Inspecting the fine-grid EPW input file (mgo.epw2.in), you can observe that this calculation employs a homogeneous \(6\times6\times6\) supercell for the polaron calculation. Once the calculation has finished, plot the displacements and wave function corresponding to this supercell, similar to what we did in Exercise 1:
mgo.epw2.in
$ scp username@frontera.tacc.utexas.edu:path_to_exercise/psir_plrn.xsf .
(local)$ VESTA psir_plrn.xsf
You should obtain something similar to this:
Action. Solve the polaron equations for an electron polaron in a conventional \(4\times 4\times 4\) cubic supercell of MgO.
Inspect the mgo.epw4.in file:
mgo.epw4.in
&inputepw
prefix = 'mgo'
outdir = './'
elph = .true.
epwread = .true.
lpolar = .true.
nbndsub = 3
dvscf_dir = './save/'
etf_mem = 0
wannierize = .false.
plrn = .true.
restart_plrn = .false.
type_plrn = 1
init_plrn = 1
init_sigma_plrn = 0.1
niter_plrn = 500
conv_thr_plrn = 1E-4
ethrdg_plrn = 1E-6
nk1 = 4
nk2 = 4
nk3 = 4
nq1 = 4
nq2 = 4
nq3 = 4
scell_mat_plrn = .true.
scell_mat(1, 1:3) = -4, 4, -4
scell_mat(2, 1:3) = -4, 4, 4
scell_mat(3, 1:3) = 4, 4, -4
/
The calculation of the polaron equations using general non-diagonal supercells is activated by the input flag scell_mat_plrn. The size and shape of transformation matrix in Eq. (10) is determined by the input flags scell_mat(:,:). You can notice that no nkf1, nkf2, nkf3 are specified in the input file. The \(k\)-point (and \(q\)-point) grids to be used in the solution of (5) and (6) corresponding to the supercell specified in the input will be automatically generated. You can inspect such grid points by looking at the output file kgrid.scell.plrn. Conversely, the list of lattice vectors forming the supercell in real space are written to the output file Rp.scell.plrn.
The job can be submitted by following a similar script as in the previous exercises:
$ sbatch submit2.sh
All the relevant information regarding the supercell transformation matrices is given in the output mgo.epw4.out:
{.. code-block:: text
Supercell transformation activated (k), as=S*at S(1, 1:3): -4 4 -4 S(2, 1:3): -4 4 4 S(3, 1:3): 4 4 -4 Transformed lattice vectors (alat units): as(1, 1:3): 4.000000 0.000000 0.000000 as(2, 1:3): 0.000000 4.000000 0.000000 as(3, 1:3): 0.000000 0.000000 4.000000 Reciprocal lattice transformation matrix, Sbar = (S^{-1})^{t}: Sbar(1, 1:3): -0.125000 0.000000 -0.125000 Sbar(2, 1:3): 0.000000 0.125000 0.125000 Sbar(3, 1:3): 0.125000 0.125000 0.000000 Transformed reciprocal lattice vectors (2pi/alat units): bs(1, 1:3): 0.250000 0.000000 0.000000 bs(2, 1:3): 0.000000 0.250000 0.000000 bs(3, 1:3): 0.000000 0.000000 0.250000
Number of unit cells within supercell: 256 Number of k-points needed: 256
}
After setting the \(\mathbf{k}\) and \(\mathbf{q}\) point grids the polaron calculation proceeds as before. You can visualize the conventional cubic cell by plotting the final atomic displacements contained in the dtau.plrn.xsf file:
dtau.pl
$ scp username@frontera.tacc.utexas.edu:path_to_exercise/dtau.plrn.xsf .
Then you can visualize the displacements by running VESTA locally on your computer:
(local)$ VESTA dtau.plrn.xsf
You should obtain something similar to this:
Exercise 4¶
In this exercise, we illustrate how to perform calculations of polaron energy landscapes.
The ground state energy of the polaron for a given configuration of the atomic displacements, \(\{\Delta \tau_{\kappa\alpha p}\}\), is obtained by performing the variational minimization of the energy functional with respect to the polaron wave function, leading to (1). The dependence on the atomic displacements allows us to explore the adiabatic potential energy surface of the polaron. An explicit relation between the eigenvalue in (1) and the formation energy can be obtained by:
where the \(B_{\mathbf{q}\nu}\) coefficients and the atomic displacements are related by Eq. (4). In Eq. (12), we use the tilde symbol in \(\Delta \tilde E_{\mathrm{f}}\) to emphasize that this energy is a minimum with respect to the electronic degrees of freedom, but the forces on the atoms are generally nonvanishing in this configuration. If we minimize this quantity also with respect to the atomic displacements, then we obtain the formation energy \(\Delta E_{\mathrm{f}}\) of (9).
In order to explore polaron energy surfaces we proceed as follows: First, we define a given displacement configuration \(\{\Delta \tau_{\kappa\alpha p}\}\); then, we perform the transformation from \(\{\Delta \tau_{\kappa\alpha p}\}\) to \(\{ B_{\mathbf{q}\nu}\}\) using (4); we diagonalize the effective Hamiltonian in (7); and we obtain the associated formation energy using (12).
Action. Calculate the polaron ground state.
We begin by calculating the formation energy and atomic displacements of the polaron ground state. The workflow is identical to that used in Exercise 3, except that here the polaron ground state is calculated directly in a non-diagonal conventional \(4\times4\times4\) supercell.
Navigate to the directory for this exercise and launch the calculation:
$ cd ../exercise3/
$ sbatch submit1.sh
You can inspect the value of the polaron formation energy in this supercell:
{.. code-block:: text
- End of self-consistent cycle
Eigenvalue (eV): 1.3341956
Phonon part (eV): -1.1248191
Electron part (eV): 0.9157267
Formation Energy (eV): -0.2090924
}
Action. Generate a path along the configuration space.
We will consider a path which corresponds to a hopping of hole polarons in MgO across nearest-neighbor sites.
As the initial configuration-coordinate point, we use the polaron displacement pattern obtained in the previous step, contained in the dtau.plrn output file. The final configuration is constructed by shifting the original small-polaron distortion by one unit cell along the \([1,\bar{1},0]\) direction. The hopping path is then defined through a linear interpolation between the initial and final configurations. An automated script to generate the displacement configurations along the hopping path is provided in the directory:
dtau.pl
$ python3 hopping_path.py
This script will output several dtau_disp.plrn_* files, which contain the displacement configurations for each of the points along the path. It also outputs a initial_final_displacements.pdf file, which contains figures of the initial and final displacement configurations.
Action. Compute the polaron energy landscape along the configuration space.
The input file corresponding to the polaron energy landscape calculation is mgo.epw3.in:
mgo.epw3.in
&inputepw
prefix = 'mgo'
outdir = './'
elph = .true.
epwread = .true.
lpolar = .true.
nbndsub = 3
dvscf_dir = './save/'
etf_mem = 0
wannierize = .false.
plrn = .true.
restart_plrn = .false.
type_plrn = 1
init_plrn = 6
init_ntau_plrn = 21
niter_plrn = 1
full_diagon_plrn = .true.
nk1 = 4
nk2 = 4
nk3 = 4
nq1 = 4
nq2 = 4
nq3 = 4
scell_mat_plrn = .true.
scell_mat(1, 1:3) = -4, 4, -4
scell_mat(2, 1:3) = -4, 4, 4
scell_mat(3, 1:3) = 4, 4, -4
/
You can notice that, in contrast to the input files used in the previous calculations, the flag init_plrn is here set to 6. This activates the initialization of the polaron equations from displacement configurations \(\{\Delta \tau_{\kappa\alpha p}\}\), which must be provided as input files. The flag init_ntau_plrn specifies the number of displacement configurations to be considered.
As discussed at the beginning of the exercise, computing the polaron energy landscape only requires a single diagonalization of the effective polaron Hamiltonian. Therefore, the flag niter_plrn is set to 1 and no self-consistency cycle is performed.
Additionally, the flag full_diagon_plrn=.true. activates the serial LAPACK diagonalization of the effective polaron Hamiltonian. For small supercells, this improves numerical aaccuracy, although at a significantly higher computational cost.
This calculation can be launched by:
$ sbatch submit2.sh
A quick inspection of the output file mgo.epw3.out shows all the relevant information about the polaron energetics for each displacement configuration along the configuration coordinate path:
{.. code-block:: text
- iter Eigval/eV Phonon/eV Electron/eV Formation/eV Error/eV
- 1 0.1334E+01 -0.1125E+01 -0.9157E+00 -0.2091E+00 0.2065E+00
Eigenvalue (eV): 1.3341949
Phonon part (eV): -1.1248183
Electron part (eV): 0.9157264
- Formation Energy at this dtau (eV): -0.2093765
- 1 0.1262E+01 -0.1055E+01 -0.8932E+00 -0.1622E+00 0.1149E-01
Eigenvalue (eV): 1.2619775
Phonon part (eV): -1.0553876
Electron part (eV): 0.8931892
- Formation Energy at this dtau (eV): -0.2065898
- 1 0.1192E+01 -0.9933E+00 -0.8688E+00 -0.1245E+00 0.1149E-01
Eigenvalue (eV): 1.1917985
Phonon part (eV): -0.9932654
Electron part (eV): 0.8688096
Formation Energy at this dtau (eV): -0.1985331 …
}
As you can observe, the polaron energy increases as we move away from the original ground state configuration, and then decreases again as we get closer to the final configuration after the hopping. This information about the energy barrier can be extracted from the output file and plotted with the script plot_barrier.py:
plot_barrier.py
$ python3 plot_hopping_barrier.py
The resulting hopping barrier height is given as an output: {.. code-block:: text
Hopping barrier relative to initial state: 132.12 meV
} and a figure with the energy landscape is written to the hopping_barrier.pdf file. Your figures should be similar to: