Phonon-assisted luminescence by finite atomic displacements

Finite atomic displacements

In recent years, several articles have been published with Yambo, who calculated the phonon-assisted luminescence by means of finite atomic displacements. In this tutorial we will show how to use YamboPy in combination with Yambo to perform these kind of calculations. For a general review of finite displacements calculations of vibrational properties we recommend this article.
Notice that calculation of phonon-assisted luminescence by finite atomic displacements are quite involved and expensive from a computational point of view, for this reason they are limited to relatively simple systems.
In this tutorial we will follow the approach of Ref., and calculate luminescence using the Roosbroeck–Shockley (RS) relation applied to the excitonic case.
Important: This tutorial requires that you are familiar with Bethe-Salpeter calculations with Yambo, and phonons with Quantum-Espresso, and YamboPy. If it is not the case please study the corresponding tutorial before continuing.

Locate the lowest indirect exciton (or indirect band gap)

In this tutorial as example we will consider the hexagonal-BN, that is an indirect insulator. The first step to calculate the phonon-assisted luminescence is to locate the lowest indirect exciton. The cleanest strategy to identify the lowest excitation would be to calculate the excitons along the whole Irreducible Brillouin zone(IBZ), then interpolate them as explained in Interpolate exciton dispersion and find the q point that correspond to the minimum.
A simpler approach is to use the GW (or Kohn-Sham) band structure to find the momentum responsible for the indirect band gap, in the major part of cases this correspond to the lowest exciton.
In the h-BN case we have:

The indirect gap is between a point close to K and the M point, see panel (a) of the figure above. We will approximate the momentum responsible for the indirect emission with q=K-M. Notice that then we should build supercells that contain this moment, so it is good to approximate it with fractions of integers that are not too large so as not to have giant supercells. In this case we get q=(1/3,-1/6, 0).
Notice that if you have a direct band gap material, the momentum q corresponding to the lowest transition is q=0 and you do not need to construct special supercells to sample the corresponding vibration, it will be sufficient to displace atoms in the primitive cell.

Calculate phonons at momentum q of the indirect transition

Now we need to calculate all atomic vibrations compatible with the momentum of the lowest indirect transition. In order to do so we first calculate the phonons with QE, here the inputs hBN_phonons.tgz, then we calculate the force constants with q3r.x and now we can interpolate phonons along all the BZ as shown in panel (b) of Figure 1. Then we use matdyn.x to calculate phonons at momentum q=K-M = (1/3, -1/6, 0):

&input
   asr='crystal',
   flfrc='bn.fc',
   flfrq='bn.freq',
   deltaE=1.d0,
   q_in_cryst_coord=.true.
   fleig='bn.eig'
/
2
0 0 0
0.3333333333333333 -0.16666666666666666666 0.0

The obtained phonon modes correspond phonons at the vertical red line in panel(b) of Fig. 1. They are written in the matdyn.modes file with the corresponding eigenvectors and will be later used to displace atoms along the phonon modes. If you open the 'matdyn.modes' file you will find phonon modes for the q=0 and q=(1/3, -1/6, 0) that in cartesian coordinates is (1/3,0,0). We copy the modes corresponding to the second q=(1/3, -1/6, 0) in a separate file, called Q2.modes that will look like:

     diagonalizing the dynamical matrix ...

q =       0.3333      0.0000      0.0000
**************************************************************************
    freq (    1) =       5.194139 [THz] =     173.257815 [cm-1]
(  0.000000  -0.000000    -0.000000  -0.000000    -0.274372  -0.475224   )
( -0.000000   0.000000    -0.000000  -0.000000    -0.445962   0.000001   )
( -0.000000  -0.000000    -0.000000  -0.000000    -0.222981  -0.386214   )
(  0.000000   0.000000    -0.000000   0.000000    -0.548742   0.000000   )
...........................

Generate the supercell

Now we use YamboPy to generate the supercell corresponding to the momentum q=(1/3, -1/6, 0). A tutorial is already present in the file tutorial/supercell/generate_supercells.py. You can edit the file to specify the momentum q, the name of the SCF initial file (uc_filnm) and the name of the phonon mode file (mode_file).

   .....
   uc_filnm = 'hBN.scf.in' #pw input
   Q = [[1,-1,0],[3,6,0]] #qpoint to fold in nondiagonal supercell: Q=(1/3,-1/6,0)
   kpoints = [12,12,4]    #NB: these fractional crystal coordinates must exactly divide the kpoint mesh!
   modes_file = 'Q2.modes' #file with phonon eigenmodes
   .....

run the script to generate the supercell without displace atoms python generate_supercells.py -N. The script will generate a non-diagonal supercell in the file "sc_nondiagonal.bn.scf' that is shown below:

The YamboPy script automatically reduce the number of k-point in such a way to have the same sampling of the corresponding primitive cell, however we advice you to check if the new k-point grid is what you expect, and modify it in case of incorrect grid.

Calculate the optical spectra in the supercell

Now you have to calculate the optical spectra and exciton in the new supercell. Do not forget to increase the number of conduction bands in both the screening the Bethe-Salpeter equation. For example in the supercell shown above there are 24 atoms, 6 times more than in the original cell, so you will have to multiply the number of bands by 6, if in the base cell you have used 40 bands here you will have to put 240 to have the same convergence parameters.
N. B. : in order to not recalculate W for each atomic displacement you can perform calculations without symmetries in the pristine supercell and then copy the corresponding dielectric constant in the calculations with displaced atoms, this is usually a safe approximation described in Ref..
After solving the BSE in the supecell you will find the excitons that were already present in the primary cell (within the numerical noise) and new lower energy excitons, with zero dipole matrix elements, that correspond to the indirect excitons that havenow been mapped at q=0, something like:

#    Maximum Residual Value =  0.38659E+01
#   
#    E [ev]             Strength           Index
#
    5.7165937         0.47290171E-11      1.0000000         <-- Indirect exciton mapped at q=0
    5.7166047         0.58239158E-09      2.0000000         <-- Indirect exciton mapped at q=0
    5.7367945         0.31830746E-09      3.0000000         <-- Indirect exciton mapped at q=0
    5.7368207         0.33674005E-10      4.0000000         <-- Indirect exciton mapped at q=0
    5.80542803        0.100763231E-6      5.00000000        <-- Direct exciton
    5.80598068        0.384604526E-8      6.00000000        <-- Direct exciton
    5.88944197         1.00000000         7.00000000        <-- Direct exciton
    5.88967037        0.102969212E-1      8.00000000        <-- Direct exciton
    ................

This result shows that indirect exciton without the coupling with phonon do not contribute the optical response.