(optics-ipa-tutorial-target)= # Optics at the independent particle level ```{todo} Add links to input variables page. ``` In this tutorial you will learn the first steps for calculating optical properties with Yambo. In particular we will perform optics calculations at the independent particle (IPA) level for bulk hexagonal boron nitride (hBN). ```{admonition} Requirements - `hBN` folder: download and extract `hBN.tar.gz` from the [Tutorial files](#tutorial-files-target) page - the `SAVE` folder (present in `hBN/YAMBO`) - `yambo` executable - `gnuplot` for plotting spectra ``` The dielectric function in the long-wavelength limit at the independent particle level (RPA without local fields) is essentially given by the following: ```{math} :label: eps00 \epsilon_{\alpha, \alpha}(\omega)= 1+\frac{16 \pi}{\Omega} \sum_{c, v} \sum_{\mathbf{k}}\frac{1}{E_{c \mathbf{k}}-E_{v \mathbf{k}}} \frac{\left|\left\langle v \mathbf{k}\left|\mathbf{p}_{\alpha}+\mathrm i\left[V^{\mathrm{NL}}, \mathbf{r}_{\alpha}\right]\right| c \mathbf{k}\right\rangle\right|^{2} }{\left(E_{c \mathbf{k}}-E_{v \mathbf{k}}\right)^{2}-(\omega + i \gamma)^{2}} ``` In practice, Yambo does not use this expression directly but solves the Dyson equation for the susceptibility $\chi$ (see tutorial on [local fields](#local-fields-tutorial-target)). We will use bulk hBN as an example, so download and unpack the files and then go to the `YAMBO` directory: ```{code-block} none $ cd hBN/YAMBO ``` Remember that the first thing to do is to initialize the `SAVE` folder (or check if it has been already initialized): ```{code-block} none $ ls SAVE/ndb.* No such file or directory $ yambo [...] $ ls SAVE l_setup r_setup $ ls SAVE SAVE/ndb.gops SAVE/ndb.kindx ``` ## Input file and parameters We want to perform a IPA optical spectrum calculation: the option to generate the input file is `-optics c` (or `-o c`). Let's also add the option to choose the input file name: ```{code-block} none $ yambo -F yambo.in_IP -o c ``` ```{tip} :class: dropdown Remember the `yambo -h` command to print all the available options. ``` ````{important} In the input file the independent particle approximation for the optics is indicated by ```{code-block} none Chimod= "IP" # [X] IP/Hartree/ALDA/LRC/PF/BSfxc ``` ```` ## Optics runlevel For optical properties we are interested just in the long-wavelength limit $\mathbf{q} \to 0$. This always corresponds to the first $\mathbf{q}$-point in the set of possible $\mathbf{q} =\mathbf{k} - \mathbf{k}'$ points. Set the following input variables in `yambo.in_IP` to: ```{raw} html 
 % QpntsRXd   
  1 | 1 |                    # [Xd] Transferred momenta
 %
 ETStpsXd= 1001              # [Xd] Total Energy steps
``` Setting the range of the variable `QpntsRXd` to `1 | 1` selects just the first $\mathbf{q}$. `ETStpsXd` instead selects the number of calculated frequencies: a large value ensures a smooth spectrum. We can now run the calculation, specifying a job label using the `-J` option: ```{code-block} none $ yambo -F yambo.in_IP -J Full ``` ````{dropdown} View l-Full_optics_dipoles_chi ```{code-block} none <---> [01] MPI/OPENMP structure, Files & I/O Directories <---> MPI assigned to GPU : 0 <---> [02] CORE Variables Setup <---> [02.01] Unit cells <---> [02.02] Symmetries <---> [02.03] Reciprocal space <---> [02.04] K-grid lattice <---> Grid dimensions : 6 6 2 <---> [02.05] Energies & Occupations <---> [03] Transferred momenta grid and indexing <---> [04] Dipoles <---> [04.01] Setup: observables and procedures <---> [DIP] Database not correct or missing. To be computed <---> [x,Vnl] computed using 4 projectors <02s> R V P [g-space] |########################################| [100%] --(E) --(X) <02s> [WARNING] [r,Vnl^pseudo] included in position and velocity dipoles. <02s> [WARNING] In case H contains other non local terms, these are neglected <02s> [05] Optics <02s> [LA] SERIAL linear algebra <02s> [X-CG] R(p) Tot o/o(of R): 5499 52992 100 <02s> Xo@q[1] |########################################| [100%] --(E) --(X) <02s> [06] Timing Overview <02s> [07] Memory Overview <02s> [08] Game Over & Game summary ``` ```` ```{note} :class: dropdown Yambo runs the calculation when there are no **lowercase** options after `yambo`. ``` When the calculation has finished you should see these files: ```{code-block} none $ ls Full SAVE yambo.in_IP r_setup l_setup l-Full_optics_dipoles_chi o-Full.eel_q1_ip o-Full.eps_q1_ip r-Full_optics_dipoles_chi ``` Let's take a moment to understand what Yambo has done inside the Optics runlevel: 1. Compute the $[\mathbf{r}, V^\mathrm{NL}]$ term 2. Read the wavefunctions from disc `[WF]` 3. Compute the _dipoles_, i.e. matrix elements of $\mathbf{p}$ 4. Write the dipoles to disk as `Full/ndb.dip*` databases. You can check this from the report file: ``````{tab-set} `````{tab-item} v5 ````{code-block} none $ grep -A16 "WR" r-Full_optics_* ```` ````{dropdown} grep output ```{literalinclude} snippets/grep-WR-r-Full_optics_v5.txt :language: none ``` ```` ````` `````{tab-item} v4.5 (or older) ````{code-block} none $ grep -A16 "WR" r-Full_optics_* ```` ````{dropdown} grep output ```{literalinclude} snippets/grep-WR-r-Full_optics_v4.txt :language: none ``` ```` ````` `````` 5. Finally, Yambo computes the non-interacting susceptibility {math}`\chi` for this {math}`\mathbf{q}`, and writes the dielectric function inside the `o-Full.eps_q1_ip` file for plotting ```{important} If `yambo` is computing new databases you will see `WR.` at the corresponding runlevel in the report file. If `yambo` finds a compatible database, this will be read and you will see `RD.` in place of `WR.`. This is useful to check wheter `yambo` is actually computing the database, especially after you have changed some input variables. ``` ````{dropdown} Note on the energy cutoff The previous calculation used **all the G-vectors** in expanding the wavefunctions, up to 1491 (~1016 components). This corresponds roughly to the cut off energy of 40 Ry we used in the DFT calculation. Generally, however, we can use a smaller value. We use the verbosity to switch on this variable: ```{raw} html 
$ yambo -F yambo.in_IP  -V RL -o c
``` Change the **value** of `FFTGvecs` and also its **unit** from `RL` (number of G-vectors) to `Ry` (energy in Rydberg): ```{raw} html 
FFTGvecs= 6 Ry  # [FFT] Plane-waves
``` Save the input file and launch the code again with a new `-J` flag to avoid reading the database from before: ```{raw} html 
$ yambo -F yambo.in_IP -J 6Ry
``` The result will be printed in the file `o-6Ry.eps_q1_ip`. ```{caution} In general the best practice is to keep the value of `FFTGvecs` at its default (maximum) value. ``` ```` ## Plotting the output To plot the result in `o-Full.eps_q1_ip` file type: ```{code-block} none $ gnuplot gnuplot> plot "o-Full.eps_q1_ip" w l ``` The figure below also shows the results with lower energy cutoff. If you also did this calculation you can add the plot on top of the previous one with `replot "o-6Ry.eps_q1_ip" w p`. ```{figure} ./images/hBN-IPA-eps_Full-vs-6Ry.png :width: 85% Comparison between the IPA dielectric functions in the long-wavelenght limit, i.e. $\Im(\varepsilon_{\alpha,\alpha}(\omega))$ (Eq.{eq}`eps00`), computed with the total number available G-vectors (Full) and a cutoff value of 6 Ry. ``` ````{note} In this case the direction of the electric field (i.e. $\alpha,\alpha$) is the one of the lattice vector $\boldsymbol{\mathrm{a}}$ of the unit cell -- which is also parallel to the $x$ axis -- because in the input file we have ```{code-block} none % LongDrXd 1.000000 | 0.000000 | 0.000000 | # [Xd] [cc] Electric Field % ``` ```` %```{attention} %:class: dropdown % %If you cannot reproduce this figure, you should check the value of `FFTGvecs` %in the input file, delete the `6Ry` folder, and try again - taking care to plot %the right file! %Indeed keep in mind that if Yambo finds another existing file named `o-6Ry.eps_q1_ip`, %then it will append a `01` suffix to the second output, e.g. `o-6Ry.eps_q1_ip_01`. %``` ```{important} We can see that there is very little difference between the two spectra. This highlights an important point in calculating excited state properties: generally, fewer G-vectors are needed than what are needed in DFT calculations. **However we suggest to never change the value of `FFTGvecs`** (unless you have very limited memory, in which case carefully check how changing this value affects the results). ``` Regarding the spectrum itself, the first peak occurs at about 4.4 eV. This is consistent with the minimum direct gap reported by Yambo: 4.28 eV. The comparison with experiment (not shown) is very poor however. This is because in most materials more advanced techniques than IPA/RPA are needed. ## q-direction ```{figure} ../basic/images/hBN-bulk-3x3.png :width: 50% Atomic structure of bulk hBN ($z \parallel c$ is the stacking direction). ``` Now let's do the calculation for the {math}`\epsilon_{z,z}(\omega)` component of the dielectric tensor. To do so we modify `LongDrXd` in the input file ```{raw} html 
$ yambo -F yambo.in_IP -V RL -o c
...
% LongDrXd
0.000000| 0.000000 | 1.000000 |        # [Xd] [cc] Electric Field
%
...
``` and run the calculation again ```{code-block} none $ yambo -F yambo.in_IP -J Full ``` This time yambo reads from the `Full` folder, so it does not need to compute the dipole matrix elements again and the calculation is faster. Let's plot the result: ```{code-block} none $ gnuplot gnuplot> plot "o-Full.eps_q1_ip" t "q || x-axis" w l,"o-Full.eps_q1_ip_01" t "q || z-axis" w l ``` ```{figure} ./images/hBN-IPA-eps_x-vs-z.png :width: 85% Comparison between $\Im(\varepsilon_{x,x}(\omega))$ and $\Im(\varepsilon_{z,z}(\omega))$ (Eq.{eq}`eps00`). ``` The absorption is suppressed in the stacking direction. As the interplanar spacing is increased, we would eventually arrive at the absorption of the BN sheet (see [local fields](#local-fields-tutorial-target) tutorial). ## Non-local commutator Last, we show the effect of switching off the non-local commutator term (the term with {math}`V^\mathrm{NL}` in Eq. {eq}`eps00`) due to the pseudopotential. As there is no option to do this inside yambo, you need to hide the database file, for example renaming it like this: ```{code-block} none $ mv SAVE/ns.kb_pp_pwscf SAVE/ns.kb_pp_pwscf_OFF ``` Then open the input file and modify it back to do the calculation along the {math}`q \parallel (1 0 0)` direction ```{raw} html 
 % LongDrXd
 1.000000  | 0.000000 |  0.000000  | 
 %
``` and launch yambo with a different `-J` label: ```{code-block} none $ yambo -F yambo.in_IP -J Full_NoVnl ``` ````{important} Note the warning in the log: ```{code-block} none <---> [WARNING] [r,Vnl^pseudo] not included in position and velocity dipoles ``` which also appears in the report file `r-Full_NoVnl_optics_dipoles_chi` as ```{code-block} none [WR./Full_NoVnl//ndb.dipoles]---------------------------------------------------- ... [r,Vnl] included : no ... ``` ```` Plotting the result we see that the difference is small: ```{code-block} none $ gnuplot gnuplot> plot "o-Full.eps_q1_ip" w l t "w/ Vnl", "o-Full_NoVnl.eps_q1_ip" w p t "no Vnl" ``` ```{figure} ./images/hBN-IPA-eps_w-vs-noVnl.png :width: 85% Comparison between $\Im(\varepsilon_{x,x}(\omega))$ computed with and without the {math}`V^{\mathrm{NL}}` term (Eq.{eq}`eps00`). ``` ```{tip} :class: dropdown When you have a large sistem, either with more projectors in the pseudopotential or more k-points, the inclusion of {math}`V^\mathrm{NL}` can make a huge difference in the computational load. So it is always worth checking to see if these terms are important for your system. ```