Many-Body Kernels for TDDFT Calculations

by Caterina Cocchi & Santiago Rigamonti for exciting nitrogen

Purpose: In this tutorial we will learn how to perform a time-dependent density-functional (TDDFT) calculation using different xc-kernels. Three examples will be proposed: One for the BSE-derived xc kernel, one for the LRC kernel, and another for the RBO (RPA-bootstrap) kernel. As a test case, the optical spectrum of LiF will be studied.

This tutorial requires the use of the following script:

• PLOT-spectra.py: Python script for generating plots of the dielectric function.

Click above on the script name for downloading it. Assuming that the script PLOT-spectra.py has been downloaded in the directory $HOME/Downloads, move it to the directory$EXCITINGTOOLS and make it executable:

$mv$HOME/Downloads/PLOT-spectra.py $EXCITINGTOOLS/PLOT-spectra.py$ chmod u+rwx $EXCITINGTOOLS/PLOT-spectra.py  ### 1. Calculation setup Before starting, be sure that the relevant environment variables are already defined as specified in How to set environment variables for tutorials scripts. Important note: All input parameters that will appear will be given in atomic units! As a preliminary step for this excited-state calculation, a ground-state calculation will be performed. In this tutorial we consider as an example LiF. Create a directory named LiF-TDDFT-kernels and move into it. $ mkdir LiF-TDDFT-kernels
$cd LiF-TDDFT-kernels  Inside the directory LiF-TDDFT-kernels we create a sub-directory GS where we perform the preliminary ground-state calculation: $ mkdir GS
$cd GS  Inside the GS sub-directory we create the input file for LiF. In the structure element we include the lattice parameter and basis vectors of LiF, which has a rock-salt cubic lattice, as well as the positions of the Li and F atoms. In the groundstate element, we include a 10$\times$10$\times$10 k-point mesh (ngridk) and a value of 14.0 for gmaxvr. This value, which is larger than the default, is needed in view of the excited-state calculation (for details on this we refer to Excited States from BSE). The resulting input file is the following: <input> <title>LiF-BSE</title> <structure speciespath="$EXCITINGROOT/species/">
<crystal scale="7.608">
<basevect>0.5 0.5 0.0</basevect>
<basevect>0.5 0.0 0.5</basevect>
<basevect>0.0 0.5 0.5</basevect>
</crystal>
<species speciesfile="Li.xml">
<atom coord="0.0000  0.0000  0.0000" />
</species>
<species speciesfile="F.xml">
<atom coord="0.5000  0.5000  0.5000" />
</species>
</structure>

<groundstate
xctype="GGA_PBE_SOL"
ngridk="10  10  10"
gmaxvr="14.0"/>

</input>


N. B.: Do not forget to replace in the input.xml the string "$EXCITINGROOT" by the actual value of the environment variable$EXCITINGROOT.You can do this by directly editing the input.xml file or by using the following command:

$SETUP-excitingroot.sh  Start now the ground-state SCF calculation and check if it finishes gracefully. $ time excitingser &


If the run ends successfully, the files STATE.OUT and EFERMI.OUT will be present in the directory. These two files are needed as a starting point for the excited-state calculation.

### 2. Theoretical background

There is a large literature dealing with the inclusion of many-body effects into TDDFT kernels, in order to correctly reproduce excitonic features (for further details we refer the seminal review ORR-2002). In the following, we will present examples related to two different approaches. In the first one we will present a TDDFT calculation of the optical spectrum of LiF performed with a xc kernel derived from the solution of the Bethe-Sapeter equation (BSE). In the second part of the tutorial, we will deal with kernels including the so-called long-range correction (LRC).

In the example treating the BSE-derived xc kernel, the scheme proposed in MDR-2003 is adopted. In this approach, a nonlocal exchange-correlation functional is derived by imposing TDDFT to reproduce the many-body diagrammatic expansion of the Bethe-Salpeter polarization function. In this way, it is shown that the TDDFT kernel is able to capture the excitonic features in solids, otherwise missing in simpler approximation for the kernel. For further details about the implementation in the code, see SAD-2009.

LRC kernels include the long-range component, which is missing at the level of the adiabatic local-density approximation (ALDA). The first model for a LRC kernel was proposed in REI-2002:

(1)
\begin{align} f_{\rm xc}(\mathbf{q},\omega) = - \frac{\alpha}{\mathbf{q}^2}\;. \end{align}

This kernel is static, non-local, and includes the long-range Coulomb tail. In this model, $\alpha$ is a material dependent parameter. An improvement of this kernel is given in BOT-2005, where a dynamical LRC kernel was developed, which explicitly presents a frequency dependence:

(2)
\begin{align} f_{\rm xc}(\mathbf{q},\omega)=-\frac{1}{\mathbf{q}^2}\left(\alpha+\beta \omega^2\right)\;. \end{align}

The values of $\alpha$ and $\beta$ are also material dependent and have to be tuned in order to correctly reproduce the experimental data for the excitons. For further details about the model we refer to the original paper BOT-2005.

The RBO kernel (see SHA-2011 and RIG-2015) shares similar characteristics to the LRC kernel: It is long-ranged and static. An important difference is, however, that it does not depend on adjustable parameters, like the parameter $\alpha$, but is fully ab-initio. It is explicitly given by (RIG-2015)

(3)
\begin{align} f_{\rm xc}^\textrm{RBO}\!(\mathbf{q}) = \frac{1}{\varepsilon_\textrm{M}^\textrm{RPA}\!(\mathbf{q},\omega=0)\:\overline{\chi}^\textrm{RPA}\!(\mathbf{q},\omega=0)}\;, \end{align}

where $\varepsilon_\textrm{M}^\textrm{RPA}$ and $\overline{\chi}^\textrm{RPA}$ are, respectively, the macroscopic dielectric function and the modified density-response function (ORR-2002) in the random-phase approximation (RPA).

### 3. TDDFT calculations using the BSE-derived xc-kernel

We are now ready to set up the TDDFT calculation of LiF, using a BSE-derived xc kernel. First of all, we move back to the parent directory. There, we create a new directory called BSE-kernel and we move into it:

$cd ..$ mkdir BSE-kernel
$cd BSE-kernel  We copy into it the files STATE.OUT and EFERMI.OUT, as well as the input file from the GS folder. $ cp ../GS/{input.xml,EFERMI.OUT,STATE.OUT} ./


Then, in the file input.xml the ground-state calculation can be skipped. To do so, set the attribute do = "skip" inside the groundstate element. Then, paste the following xs block:

...
<xs
xstype="TDDFT"
ngridk="4 4 4"
vkloff="0.097 0.273 0.493"
ngridq="4 4 4"
nempty="3"
gqmax="3.0"
scissor="0.2095"
tevout="true">

<energywindow
intv="0.0 1.0"
points="1200"/>

<screening
screentype="full"
nempty="100"/>

<tddft
fxctype="MB1"
aresdf="false"
aresfxc="false"/>

<BSE
bsetype="singlet"
nstlbse="1 5 1 4"
/>

<qpointset>
<qpoint>0.0 0.0 0.0</qpoint>
</qpointset>

</xs>
...


This block is very similar to the one presented in Excited States from BSE, to which we refer for an exhaustive description of the input attributes. In the following, we discuss only the relevant parameters for the TDDFT calculation with a BSE-derived kernel.

1. In the xs element we set the attribute xstype = "TDDFT", as we are performing a TDDFT calculation;
2. We decrease the number of empty states with respect to the example in Excited States from BSE by choosing (nempty = "3") to speed up the calculation;
3. Both tddft and bse elements appear in the input file;
4. Inside the tddft element, we specify the parameters related to TDDFT calculations. We refer to Excited States from TDDFT for additional details about this element. In this case we choose a many-body, BSE-derived kernel fxctype = "MB1" and we set to "false" the attributes aresdf and aresfxc to speed up the calculation (for further details see the Input Reference);
5. The bse element must be specified to generate the kernel.

To run the calculation, simply type:

$time excitingser &  Once the run is completed, we can analyze the results. As in any TDDFT calculation, a number of output files are created (see also Excited States from TDDFT for additional details). Here we are interested in the files named EPSILON_NAR_FXCMB1_OCYY_QMT001.OUT (YY = 11, 22, and 33) and specifically in the imaginary part of the macroscopic dielectric function, giving the optical absorption spectrum. In order to visualize, e.g., the "11" component of the spectrum (notice that, as LiF has a cubic crystal structure, the three diagonal components are the same!), type: $ PLOT-spectra.py EPSILON_NAR_FXCMB1_OC11_QMT001.OUT.xml


This will generate a file named PLOT.png as output, that you can visualize with any png visualizer (e.g., with "okular", "display" (ImageMagik), "gwenview", etc.). The resulting plot is the following:

The main features of the optical spectrum are clearly visible in the graph. The intense excitonic peak at about 13 eV dominates the low energy part of the spectrum and another strong peak is found above 20 eV. This result is in agreement with the spectrum obtained by solving the BSE equation (see Excited States from BSE). For further comparison with the literature, we refer to MDR-2003.

### 4. TDDFT calculations using the LRC and RBO Kernels

As usual, we recommend to move to the parent directory and create a new folder where running the TDDFT calculation with the LRC and RBO kernels:

$cd ..$ mkdir LRC-RBO-kernel
$cd LRC-RBO-kernel  We copy into it the files STATE.OUT and EFERMI.OUT, as well as the input file from the GS folder. $ cp ../GS/{input.xml,EFERMI.OUT,STATE.OUT} ./


Then, in the file input.xml the ground-state calculation can be skipped by setting the attribute do ="skip" inside the groundstate block. Furthermore, in the file input.xml the following xs block should be inserted:

...
<xs
xstype="TDDFT"
ngridk="4 4 4"
vkloff="0.097 0.273 0.493"
nempty="5"
gqmax="3.0"
scissor="0.2095"
tevout="true">

<energywindow
intv="0.0 1.0"
points="1200"/>

<tddft
fxctype="LRCstatic"
alphalrc="10.0"/>

<qpointset>
<qpoint>0.0 0.0 0.0</qpoint>
</qpointset>

</xs>
...


With respect to the previous example on BSE-derived kernel first of all we notice that the BSE element has disappeared: This is actually a purely TDDFT calculation. Inside the tddft element we have changed the attribute fxctype to LRCstatic: In this way we are setting the exchange-correlation kernel to be the LRC static one (REI-2002). The crucial parameter in this calculation is alphalrc, which determines the value of $\alpha$ in the expression of the kernel $f_{xc}(\mathbf{q},\omega) = - \alpha/\mathbf{q}^2$. This is a dimensionless number, which we choose here to be equal to 10.0. Run the calculation (note that it will take much shorter than the calculation with the BSE-derived kernel!!).

To plot the resulting optical absorption spectrum run

$PLOT-spectra.py EPSILON_FXCLRCstatic_OC11_QMT001.OUT.xml  The result should look like this: In the plot above, the strong excitonic peak at about 12.5 eV characterizing the spectrum of LiF is correctly reproduced by the TDDFT calculation with the static LRC kernel. Compared to the result obtained with the BSE-derived kernel the main differences appear in the higher energy region of the spectrum, above 20 eV. However, the purpose to correctly reproduce the intense bound exciton of LiF is fulfilled. Next, we now compare this result with that obtained with the dynamical LRC kernel. In this case, the xs block in the input file is the following: ... <xs xstype="TDDFT" ngridk="4 4 4" vkloff="0.097 0.273 0.493" nempty="5" gqmax="3.0" broad="0.007" scissor="0.2095" tevout="true"> <energywindow intv="0.0 1.0" points="1200"/> <tddft fxctype="LRCdyn" alphalrcdyn="2.0" betalrcdyn="35.6"/> <qpointset> <qpoint>0.0 0.0 0.0</qpoint> </qpointset> </xs> ...  Only the tddft element is modified compared to the static case. The attribute fxctype is set to "LRCdyn" to choose the dynamical$f_{xc}$quoted in BOT-2005. In this case, two parameters must be tuned, namely alphalrcdyn and betalrcdyn, which correspond respectively to$\alpha$(dimensionless) and$\beta(Ha-2) in the model kernel of Eq. (2). The parameters appearing in Eqs. (1) and (2) are chosen such that the calculated BSE energy peak coincides with the experimental one. For more details we refer to the original paper BOT-2005. As a rule of a thumb, the parameters alphalrcdyn and betalrcdyn can be set in order to mimic the behavior of alphalrc in the static LRC kernel: (4) \begin{align} \alpha_{\rm dyn} + \beta\; \omega_{\rm peak}^2 = \alpha_{\rm static}\;, \end{align} where\omega_{\rm peak}$indicates the position of the first excitonic peak for LiF at 12.9 eV. Choosing a value of 2.0 for alphalrcdyn, as suggested in BOT-2005 one obtains 35.6 for betalrcdyn. To plot the spectrum run $ PLOT-spectra.py EPSILON_FXCLRCdyn_OC11_QMT001.OUT.xml


The resulting optical spectrum is:

The first intense excitonic peak is again well reproduced by the LRC kernel.

Finally, we calculate the spectrum with the RBO kernel. In this case, the xs block in the input file is the following:

...
<xs
xstype="TDDFT"
ngridk="4 4 4"
vkloff="0.097 0.273 0.493"
nempty="5"
gqmax="3.0"
scissor="0.2095"
tevout="true">

<energywindow
intv="0.0 1.0"
points="1200"/>

<tddft
fxctype="RBO"
/>

<qpointset>
<qpoint>0.0 0.0 0.0</qpoint>
</qpointset>

</xs>
...


As you can notice, now the tddft element does not contain any empirical parameter, in accordance with the defining Eq. (3) above. The obtained spectrum shows a bound excitonic peak inside the band-gap, whose binding energy is in good agreement with experiment. In order to compare the RBO result with the LRC static result and the BSE-derived xc kernel, we can generate a plot containing all the spectra together by executing:

\$ PLOT-spectra.py EPSILON_FXCLRCstatic_OC11_QMT001.OUT.xml EPSILON_FXCRBO_OC11_QMT001.OUT.xml ../BSE-kernel/EPSILON_NAR_FXCMB1_OC11_QMT001.OUT.xml


The generated PLOT.png file looks like:

##### Exercises
• If you have already done the tutorial Excited states from TDDFT, calculate the optical absorption spectrum for LiF using RPA and ALDA kernels. What happens to the excitonic peak?
• Decrease the parameter alphalrc in the calculation with the static LRC kernel and check what happens to the excitonic peak. Compare your results with the onset of the spectrum obtained from the RPA calculation.
• Tune the parameters alphalrcdyn and betalrcdyn, following the rule of a thumb suggested above. What happens to the spectrum?

### Literature

• ORR-2002: G. Onida, L. Reining, and A. Rubio, Rev. Mod. Phys. 74, 601 (2002).
• MDR-2003: A. Marini, R. Del Sole, and A. Rubio, Phys. Rev. Lett. 91, 256402 (2003).
• SAD-2009: S. Sagmeister and C. Ambrosch-Draxl, Phys. Chem. Chem. Phys. 11, 4451 (2009).
• REI-2002: L. Reining, V. Olevano, A. Rubio, and G. Onida Phys. Rev. Lett. 88, 066404 (2002).
• BOT-2005: S. Botti, A. Fourreau, F. Nguyen, Y.-O. Renault, F. Sottile, and L. Reining, Phys. Rev. B 72, 125203 (2005).
• SHA-2011: S. Sharma, J. K. Dewhurst, A. Sanna, and E. K. U. Gross, Phys. Rev. Lett. 107, 186401 (2011).
• RIG-2015: S. Rigamonti, S. Botti, V. Veniard, C. Draxl, L. Reining, and F. Sottile, Phys. Rev. Lett. 114, 146402 (2015).