(local-fields-tutorial-target)= # RPA in Yambo: local field effects In this tutorial you will learn how to calculate the macroscopic dielectric function with Yambo using two-dimensional (2D) hBN as an example. ```{admonition} Requirements - `SAVE` folder of 2D hBN: download and extract `hBN-2D.tar.gz` from the [Tutorial files](#tutorial-files-target) page - `yambo` executable - `gnuplot` for plotting spectra ``` The macroscopic dielectric function is obtained by including the so called local field effects (LFE) in the calculation of the response function. Within the time-dependent DFT formalism this is achieved by solving the Dyson equation for the susceptibility {math}`\chi`. In reciprocal space this is given by: ```{math} :label: chi_dyson \chi_{\mathbf{G}, \mathbf{G}^{\prime}}(\mathbf{q}, \omega) = \chi_{\mathbf{G}, \mathbf{G}^{\prime}}^{0}(\mathbf{q}, \omega)+\sum_{\mathbf{G}_{1}, \mathbf{G}_{2}} \chi_{\mathbf{G}, \mathbf{G}_{1}}^{0}(\mathbf{q}, \omega) \left[ v_{\mathbf{G}_{1}}(\mathbf{q}) \delta_{\mathbf{G}_{1}, \mathbf{G}_{2}}+f_{\mathbf{G}_{1}, \mathbf{G}_{2}}^{x c} \right] \chi_{\mathbf{G}_{2}, \mathbf{G}^{\prime}}(\mathbf{q}, \omega) \, . ``` The microscopic dielectric function is related to {math}`\chi` by ```{math} :label: eps_micro \epsilon_{\mathbf{G}, \mathbf{G}^{\prime}}^{-1}(\mathbf{q}, \omega) = \delta_{\mathbf{G}, \mathbf{G}^{\prime}} + v_{\mathbf{G}}(\mathbf{q}) \chi_{\mathbf{G}, \mathbf{G}^{\prime}}(\mathbf{q}, \omega) ``` and the macroscopic dielectric function is obtained by taking the {math}`\mathbf{G}=\mathbf{G^{\prime}}=0` component of the inverse microscopic one ```{math} :label: eps_macro \epsilon_{M}(\omega) = \lim_{\mathrm{q} \rightarrow 0} \frac{1}{\epsilon_{\mathrm{G}=0, \mathrm{G}^{\prime}=0}^{-1}(\mathbf{q}, \omega)} \, . ``` Experimental observables like the optical absorption and the electron energy loss are obtained from the macroscopic dielectric function: ```{math} :label: abs_and_eels \operatorname{Abs}(\omega) = \operatorname{Im} \epsilon_{M}(\omega) \; \; , \; \; \operatorname{EELS}(\omega) = -\operatorname{Im} \frac{1}{\epsilon_{M}(\omega)} \, . ``` In the following we will neglect the {math}`f_{\mathbf{G}_1,\mathbf{G}_2}^{xc}` term in the kernel: we perform the calculation at the RPA level including just the Hartree term (from {math}`v_G`). If we also neglect the Hartree term we go back to the independent particle approximation since we effectively remove the kernel and set {math}`\chi = \chi_0`. Before starting go to the `hBN-2D/YAMBO` directory and initialize the `SAVE` folder: ```{code-block} none $ cd hBN-2D/YAMBO $ ls SAVE $ yambo ``` ## Input file and parameters After the initialization we can generate the input file. To include the local fields variables use `-o c -k hartree`. Check the meaning of the options with the command-line interface if you are not familiar with them. Let's start doing a calculation with the light polarization {math}`\mathbf{q}` in the plane of the BN sheet. Generate the input like this ```{code-block} none $ yambo -F yambo.in_RPA -V RL -o c -k hartree ``` (thus naming it `yambo.in_RPA`) and then set the following variables to these values: ```{raw} html 
 FFTGvecs=  3951            RL       # [FFT] Plane-waves
 Chimod= "Hartree"                   # [X] IP/Hartree/ALDA/LRC/BSfxc
 NGsBlkXd=  3        Ry              # [Xd] Response block size
 % QpntsRXd
  1 | 1 |                            # [Xd] Transferred momenta
 %
 % EnRngeXd
  0.00000 | 20.00000 | eV            # [Xd] Energy range
 %
 % DmRngeXd
  0.20000 | 0.200000 | eV            # [Xd] Damping range
 %
 ETStpsXd= 2001                      # [Xd] Total Energy steps
 % LongDrXd
 1.000000 | 0.000000 | 0.000000 |    # [Xd] [cc] Electric Field
 %
``` In this input file: - We evaluate the {math}`\mathbf q \to 0` response function choosing the direction in the limit parallel to the plane of the hBN sheet; - We set the energy range for the spectrum to $0-20$ eV; - We set the broadening through the variable `DmRngeXd`; - We select the Hartree kernel; - We included G-vectors in the screening up to 3 Ry (about 85 G-vectors for this system[^85Gvecs]). [^85Gvecs]: This value does not appear in the report file of an "optics,dipoles,kernel,chi" calculation, but you can find it as `Hartree size` in the header of the `o-` files. Otherwise you can get it from the report of a "screen,dipoles,em1s" calulation (input file from `yambo -X s`) grepping `'X matrix size'`. ```{note} While `FFTGvecs` controls the cutoff for the plane wave expansion of the wavefuntions (hidden in dipole matrix elements appearing in the expression of $\chi_0$), `NGsBlkXd` determines the size of the dielectric matrix in G-space (or equivalently the cutoff on the sum over $\mathbf{G}_1, \mathbf{G}_2$ in Eq. {eq}`chi_dyson`). Remember that `NGsBlkXd` is a typical parameter to be converged while it is better to **not** modify `FFTGvecs` and leave it at its default (maximum) value. ``` ## LFE in the periodic direction ```{figure} ../basic/images/hBN-2D-3x3.png :width: 25% Atomic structure of 2D hBN ($z \parallel c$ is the non-periodic direction). ``` Let's now run the calculation specifying the electric field direction in the job string label: ```{code} $ yambo -F yambo.in_RPA -J q100 ``` This calculation should take about 4 minutes in serial. ````{dropdown} yambo -F yambo.in_RPA -J q100 (v5.3 log) ```{literalinclude} snippets/l-q100_optics_dipoles_kernel_chi_v5.3 :language: none ``` ```` Let's now look at the results. We want to compare the absorption with and without the local fields: by inspecting the output file `o-q100.eps_q1_inv_rpa_dyson` we find that this information is given in the 2$^{nd}$ and 4$^{th}$ columns, respectively: ```{code} $ head -n 50 o-q100.eps_q1_inv_rpa_dyson | grep '#' | tail -n 3 # # E[1] [eV] Im(eps) Re(eps) Im(eps_o) Re(eps_o) # ``` Plot the RPA and the IPA response: ```{code-block} $ gnuplot $ gnuplot> plot "o-q100.eps_q1_inv_rpa_dyson" u 1:2 w l t "RPA-LFE","o-q100.eps_q1_inv_rpa_dyson" u 1:4 w l t "noLFE", "o-q100.eel_q1_inv_rpa_dyson" u 1:4 w l ls 7 dt 2 t "EELS" ``` ```{figure} ./images/hBN-2D-LFE_qx.png :width: 85% Absorption and EELS spectra from the macroscopic dielectric function (Eq. {eq}`abs_and_eels`) with the electric field in the in-plane direction. ``` There is very little influence of local fields for 2D hBN. This is generally the case for semiconductors or materials with a smoothly varying electronic density. In the same figure, we have also shown the EELS spectrum (`o-q100.eel_q1_inv_rpa_dyson`) for comparison. ## LFE in the non-periodic direction Now let's orientate {math}`\mathbf{q}` in the direction perpendicular to the BN plane. Open the input file `yambo.in_RPA` and change `LongDrXd` to: ```{code-block} none % LongDrXd 0.000000 | 0.000000 | 1.000000 | # [Xd] [cc] Electric Field % ``` Run the calculation (don't forget to change the `-J` label): ```{code} $ yambo -F yambo.in_RPA -J q001 >& /dev/null & ``` The `>& /dev/null &` makes Yambo run in the background and print the log into a file. You can follow the progress of the calculation with ```{code-block} none $ tail -f l-q001_optics_dipoles_kernel_chi ``` When the calculation is done we can plot the new output: ```{code-block} none $ gnuplot gnuplot> plot "o-q001.eps_q1_inv_rpa_dyson" u 1:2 w l,"o-q001.eps_q1_inv_rpa_dyson" u 1:4 w l ``` ```{figure} ./images/hBN-2D-LFE_qz.png :width: 85% Absorption spectra from the macroscopic dielectric function (Eq. {eq}`abs_and_eels`) with the electric field in the out-of-plane direction. ``` In the out of plane direction, the absorption is strongly blueshifted with respect to the in-plane absorption. Furthermore, the influence of local fields is striking and quenches the spectrum strongly. This is the well known depolarization effect. Local field effects are much stronger in the perpendicular direction because the charge inhomogeneity is dramatic. Also keep in mind that many G-vectors are needed to account for the sharp change in the potential across the BN-vacuum interface. ### Absorption vs EELS In order to understand this further, we plot the electron energy loss spectrum for this component and compare with the absorption: ```{code-block} none $ gnuplot $ gnuplot > plot "o-q001.eps_q1_inv_rpa_dyson" w l,"o-q001.eel_q1_inv_rpa_dyson" w l ``` ```{figure} ./images/hBN-2D-LFE_qz_abs-vs-eels.png :width: 85% Absorption and EELS spectra from the macroscopic dielectric function (Eq. {eq}`abs_and_eels`) with the electric field in the out-of-plane direction. ``` The conclusion is that the absorption and EELS coincide for isolated systems. To understand this result, you need to consider the role of the _macroscopic_ screening in the response function and the long-range part of the Coulomb potential. See e.g. [this article](https://doi.org/10.1002/qua.20486)[^SottileTDDFT2005] [^SottileTDDFT2005]: F. Sottile, et al., "TDDFT from molecules to solids: The role of long-range interactions", International journal of quantum chemistry 102.5 (2005): 684-701.