In this section, we will demonstrate how to add the quadrupole correction to the e-ph coupling matrices from first-principles.
This tutorial includes one extra step compared to the regular Perturbo workflow. First, the dynamical quadrupole rank-3 tensor (Q) needs to be computed with linear response theory. Though Perturbo is mainly interfaced to QE, this quantity Q can only be calculated via Abinit currently. After the calculation of Abinit finishes, a Fortran program provided in the example folder can be used to extract the Q tensor into a hdf5 file which will be read automatically by qe2pert.x.
Methodology
We express the quadrupole correction to the e-ph coupling interactions as,
\[g^{quad}_{mnv}(k,q) = \frac{e^2}{\Omega \epsilon_0} \sum_\kappa \Big(\frac{\hbar}{2\omega_{vq}M_\kappa}\Big)^{1\over2} \sum_{G\ne -q} {1 \over 2} \frac{(q+G)_\alpha(Q_{\kappa,\alpha\beta\gamma}e^{\kappa}_{vq,\gamma})(q+G)_\beta}{(q+G)_\alpha\epsilon_{\alpha\beta}(q+G)_\beta}\times\langle m,k+q|e^{i(q+G)(r-\tau_\kappa)}\vert nk \rangle\]where \(Q_{\kappa,\alpha\beta\gamma}\) is the calculated dynamical quadrupole tensor. For 3D cases, we adopt the value from first-principles directly. For 2D cases, a physical dynamical quadrupole tensor (\(\hat{Q}\)) will be derived from its 3D analogue (Q) which are computed from linear response theory in a supercell setup (vacuum layers).
For out-of-plane components,
\[\hat{Q}^{zz}_{\kappa\beta} = \frac{Q^{zz}_{\kappa\beta} + 2\tau_{\kappa z} Z^{z}_{\kappa\beta}}{\epsilon_{zz}}\]and mixed components,
\[\hat{Q}^{z\alpha}_{\kappa\beta} = \frac{Q^{z\alpha}_{\kappa\beta} + \tau_{\kappa z} Z^{\alpha}_{\kappa\beta}}{\epsilon_{zz}}\]While for in-plane components, we have,
\[\hat{Q}^{\alpha\gamma}_{\kappa\beta} =\delta_{\alpha\gamma}\hat{Q}^{zz}_{\kappa\beta} Q^{\alpha\gamma}_{\kappa\beta}-\epsilon_{\alpha\gamma}\hat{Q}^{zz}_{\kappa\beta}\]Please refer to this paper for details.
The application of quadrupole correction in GaN is present in this paper:
Phys. Rev. Lett, 125, 136602, (2020).
The theory of quadrupole correction for electron-phonon interactions, along with two examples (silicon and PbTiO3), is presented in this paper:
Phys. Rev. B, 102, 125203, (2020).
Calculate dynamical quadrupole tensor Q
Basically, you need to perform linear response calculations using Abinit. Please visit Abinit webpage for details.
After the calculation, we will get the output file from Abinit. In the example, we provide a demo with the name of tlw_3.abo.
Next, a Fortran program qtensor.f90 is provided in the example folder to extract the dynamical quadrupole tensor from tlw_3.abo file.
We provide a Makefile for gfortran to compile the code. HDF5 and h5fortran library need to be installed before compiling it. Please see this section for more details. After specifying the library path there, then
make
Then an executable qtensor will be produced. Then, copy tlw_3.abo from the Abinit calculations.
~$ ls
qtensor qtensor.h5 tlw_3.abo
~$ ./qtensor
Please give the anaddb output file name: tlw_3.abo
The number of atom in your system is: 2
The Quadrupole tensor has read:
1 x -0.0000000 -0.0000000 0.0000000 -14.1085230 0.0000000 0.0000000
1 y -0.0000000 0.0000000 -0.0000000 0.0000000 -14.1085230 -0.0000000
1 z 0.0000000 0.0000000 0.0000000 0.0000000 -0.0000000 -14.1085230
2 x 0.0000000 -0.0000000 -0.0000000 14.1085230 -0.0000000 -0.0000000
2 y 0.0000000 -0.0000000 0.0000000 -0.0000000 14.1085230 0.0000000
2 z -0.0000000 -0.0000000 -0.0000000 -0.0000000 0.0000000 14.1085230
The Quadrupole tensor has been written to qtensor.h5! Please rename it for your own system!!!
and you maybe asked to fill in the name of the output file (here we have tlw_3.abo) to continue.
After that, you will get an output file like qtensor.h5, which is also the input file of qe2pert.x.
Perturbo workflow
The last step is to copy this file qtensor.h5 to the folder where you run qe2pert.x:
cp qtensor.h5 /path/to/qe2pert
No additional parameters need to be specified in either qe2pert.in or pert.in. The code will detect whether qtensor.h5 file exists and adds the information to prefix_epr.h5 file automatically.
E-ph matrix elements with quadrupole correction
Finally, we are good to run Perturbo calculation with this prefix_epr.h5 file. We will focus on the ephmat calculation of silicon for the lowest valence band.
Run perturbo.x:
mpirun -n 1 perturbo.x -npools 1 -i pert.in > pert.out
We can observe the remarkable correction of the deformation potential for the longitudinal optical phonon around the zone center with the value of around 3 eV/Angstrom. For comparison, without quadrupole correction, it will be about 0.01 eV/Angstrom.