Excited States from TDDFT

by Stephan Sagmeister & Santiago Rigamonti for exciting boron

Purpose: In this tutorial, you will learn how to perform a basic time-dependent density-functional theory (TDDFT) calculation. As an example, the loss function of bulk silver is studied for the optical case (long-wavelength limit q=0).

---

0. Define relevant environment variables

Read the following paragraphs before starting with the rest of this tutorial!

Before starting, be sure that relevant environment variables are already defined as specified in How to set environment variables for tutorials scripts. Here is a list of the scripts which are relevant for this tutorial with a short description.

• EXECUTE-single.sh: (Bash) shell script for running a single exciting calculation.
• PLOT-loss-function.py: Python script for generating plots of the loss function.

From now on the symbol $ will indicate the shell prompt.

Requirements: Bash shell. Python numpy, lxml, matplotlib.pyplot, and sys libraries.

Important note: All input parameters that will appear will be given in atomic units!

---

1. Groundstate calculation

i) Preparation of the input file

The first step is to create a directory for the system that we want to investigate. In this tutorial, we consider as an example the calculation of the loss function for silver in the fcc cubic structure. Thus, we will create a directory Ag-tddft and we move inside it.

$ mkdir Ag-tddft
$ cd Ag-tddft

As a starting point for the TDDFT calculation, we need converged electron density and potential. To this end, we perform a groundstate calculation. Therefore, we create (or copy from a previous calculation) the file input.xml that could look like the following.

<input>
 
   <title>Loss function of Ag</title>
 
   <structure speciespath="$EXCITINGROOT/species/">
      <crystal scale="7.72">
         <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="Ag.xml">
         <atom coord="0.0  0.0  0.0" />
      </species>
   </structure>
 
   <groundstate 
      ngridk="10  10  10" />
 
</input>

Please, remember that the input file for an exciting calculation must be always called input.xml. For further details see Input Reference and How to start an exciting calculation.

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

ii) Running the groundstate calculation

You can start the calculation by invoking the script EXECUTE-single.sh.

$ EXECUTE-single.sh test-1

After the execution, all results are written in the subdirectory test-1. Move inside it.

$ cd test-1

You can check the bunch of files created during the run, especially the main output file INFO.OUT, for convergence information. If the calculation of the ground state has been finished successfully, in the last lines of the INFO.OUT file you should find the message

...
================================================================================
| EXCITING BORON stopped                                                       =
================================================================================

Check if the calculation finishes gracefully and if the EFERMI.OUT and STATE.OUT files are present. These files contain the Fermi level and the converged electron density and potential, respectively, and are the starting point of the TDDFT calculation.

Please note: To obtain reliable results it is necessary to perform careful convergence tests with respect to the k-point mesh (parameter ngridk) and the size of the basis set (parameter rgkmax). For details see Simple convergence tests.

---

2. Theoretical background for TDDFT calculation

In this application of TDDFT, we want to calculate the loss function for the optical case (i.e., for q=0). The loss function is the quantity which is directly measurable using electron-energy loss spectroscopy (EELS). The loss function $L(\mathbf{q},\omega)$ is defined as

(1)

\begin{align} L(\mathbf{q},\omega)=-\text{Im}\,\left[\varepsilon_{\!M}\!(\mathbf{q},\omega)\right]^{-1}\;, \end{align}

where $\varepsilon_{\!M}\!(\mathbf{q},\omega)$ is the macroscopic dielectric function which, in turn, is related to the inverse dielectric matrix by

(2)

\begin{align} \left[\varepsilon_{\!M}\!(\mathbf{q},\omega)\right]^{-1}=\varepsilon_{00}^{-1}(\mathbf{q},\omega)\,. \end{align}

The inverse dielectric matrix is given by $\varepsilon^{-1}=1+v\,\chi$, where $v$ is the bare Coulomb kernel and $\chi$ is the reducible polarizability. The latter is obtained in exciting using the linear response to TDDFT through the solution of the Dyson equation

(3)

\begin{align} \chi=\chi_0+\chi_0\,(v+f_{\rm xc})\,\chi\;, \end{align}

where $\chi_0$ is the Kohn-Sham independent-particle polarizability and $f_{\rm xc}$ is the exchange-correlation kernel. The neglect of the $f_{\rm xc}$ term corresponds to the random-phase approximation (RPA). In actual calculations, the different response functions are represented by their Fourier expansion. The size of the matrices of Fourier components which represent the response functions is limited through the cutoff value G_max in such a way that only matrix elements X_G,G'(q) with both |q+G|, |q+G'| < G_max are taken into account (X = $\chi_0$, $\chi$, etc.). If G_max > 0 the calculation is said to include (fully or partially, depending on the size of G_max) local field effects (LFE).

Work-flow of a TDDFT linear-response calculation

[Show work-flow]

The calculation of the different quantities involved determines the work-flow of a TDDFT linear-response calculation, as discussed in the following.

1. Eigenvalues ɛ_nk and wavefunctions ψ_nk(r) of the Kohn-Sham equations are determined for a given k-mesh.
2. Momentum matrix elements, p_α,nmk = <ψ_nk|p_α|ψ_mk>, are calculated.
3. q-dependent matrix elements of plane waves, M_nmk(q,G) = <ψ_nk|exp[-i(q+G)r]|ψ_mk+q>, are evaluated.
4. The Kohn-Sham response function $\chi_0$ is assembled including the matrix elements
   • only p_α,nmk for q=0 without LFE,
   • both p_α,nmk and M_nmk(0,G) for q=0 with LFE,
   • only M_nmk(q,G) for finite momentum transfer (q different from 0).
5. In a last step the Dyson equation for the interacting response function $\chi$ is solved and the macroscopic dielectric function $\varepsilon_M$ is obtained therefrom. For q=0, all diagonal tensor components of the macroscopic dielectric function are calculated by default.

In this flow chart, each level along a horizontal axis denotes quantities which are independent from each other.

---

3. How to perform a TDDFT calculation

Setting elements and attributes

If you have already performed the (groundstate) calculation shown in Section 1, you can prepare the input file (input.xml) of a TDDFT calculation starting from the existing groundstate input file ([1]) according to the following: First, the explicit groundstate calculation can now be skipped by adding the do = "skip" attribute to the element groundstate.

...
   <groundstate 
      do="skip"
      ... >
      ...
   </groundstate>
...

In order to perform a calculation of the excited states, a new block, corresponding to the element xs, must be added to the input file inside the input element.

...
   <xs 
      xstype ="TDDFT"
      ngridk="8 8 8" 
      vkloff="0.097 0.273 0.493"
      nempty="30"
      gqmax="2.0"
      broad="0.004"
      tevout="true">
 
      <energywindow 
         intv="0.0 2.0" 
         points="500" />
 
      <tddft 
         fxctype="RPA"/>
 
      <qpointset>
          <qpoint> 0.0 0.0 0.0 </qpoint>
      </qpointset>
 
   </xs>
...

These settings in the input file correspond to a calculation with the following main features:

• It is a TDDFT calculation (the attribute xstype is set to "TDDFT").
• It is using the random-phase approximation (RPA) for the exchange-correlation kernel (the attribute fxctype of the subelement tddft is set to "RPA").
• Local-field effects are included (the value of the cutoff G_max, defined by the attribute gqmax of the element xs, is larger than zero). In order to neglect LFE in the calculation, gqmax should be set to zero. Notice that, for the sake of comparison, results without the inclusion of LFE are always stored and can be found in the output files flagged with the _NLF_ label (see TDDFT output files).

Now, make yourself familiar with the input parameters appearing in above (for details, see Input Reference) and check if and how the chosen values differ from the default values.

General comments on input parameters for TDDFT

• The k-point mesh (specified by the ngridk attribute) is crucial for a good resolution of the spectra, though choosing a dense mesh can result into a time-consuming calculation. Therefore, the spectra always have to be converged with respect to the choice of this mesh.
• In general, it is a good practice to shift the k-mesh off symmetry by modifying the vkloff attribute. Following this procedure, all points of the k mesh will be shifted by a small displacement. Such a displacement should break all symmetry relations among the k-points of the mesh. In this way, all the k-points in the mesh will be crystallographically inequivalent and there will be no redundant contribution to the spectrum.
• The q-vector (specified by the subelement qpoint of qpointset) defines whether an optical spectrum is calculated (q=0), or a finite momentum transfer is considered (q different from zero). Note that the q-vector is given in crystal coordinates.
• The value of the attribute broad is a Lorentzian broadening for all spectra.
• The attribute tevout sets the energy output to be given in electronvolt (eV).
• The element energywindow and its attributes allow to set the number and range of the energy values for which the spectra are calculated.

Running exciting

Having modified the input file according to the above prescription, you should now be able to launch the TDDFT calculation inside the subdirectory test-1.

$ time excitingser >& outputXS.txt &

You can check the progress of the TDDFT run by inspection of the file INFOXS.OUT. You will see detailed information about the several tasks executed. After every task (identified with a specific task number that is used for internal purposes of the code) is finished, the following message appears in the file

================================================================================
| EXCITING BORON stopped for task    XXX                                   =
================================================================================

This does not mean that the complete TDDFT calculation has finished, but only the specific task XXX. The last task of a TDDFT calculation consists on the solution of the Dyson equation to obtain the dielectric function, the calculation of different spectra and the printing of them to corresponding files. It is labeled with the number 350.

Output files

While the calculation is running you can check the file INFOXS.OUT and the terminal output file, outputXS.txt. The latter one should be empty at any stage of the calculation for a successful run. Once the calculation is finished, a bunch of output files (see TDDFT output files) will be present. Most of them containing a _QMTxxx label. This stands for q momentum transfer and corresponds to the label xxx of the q-points listed in the QPOINTS.OUT file. For our q=0 case we will only have files containing a _QMT001 label.
We will concentrate on the LOSS_* files and pick the XML output file LOSS_FXCRPA_OC11_QMT001.OUT.xml. This file contains the data for the loss function and the dynamical structure factor. Get familiar with the details on the output files described in TDDFT output files.

---

4. Visualizing the output

From the XML data file for the loss function, we can generate a PostScript (PLOT.ps) and a PNG (PLOT.png) output by typing the following command inside the subdirectory test-1:

$ PLOT-loss-function.py LOSS_FXCRPA_OC11_QMT001.OUT.xml

As you can see, we are executing the PLOT-loss-function.py script with one argument, which is the name of one XML output file. You can add as many arguments (XML file names) in the command line as you like, and the script will add the corresponding curves to the same plot. For example, to make a comparison between the loss function with and without local-field effects included, we can run the following command to generate a plot containing both loss functions.

$ PLOT-loss-function.py LOSS_FXCRPA_OC11_QMT001.OUT.xml LOSS_NLF_FXCRPA_OC11_QMT001.OUT.xml

You should obtain a result similar to what is shown below:

In the figure above, the labels LFE and no-LFE indicate calculations performed including and excluding local-field effects, respectively. One immediate conclusion from the result above is about the importance of the local-field effects in the high-energy range, especially in the range above 35 eV. They drastically reduce the oscillator strength. At higher energies, more localized states contribute to the determination of the dielectric response, so the loss function is more sensitive to the inclusion of LFE in this part of the spectrum (see, e.g., Ref.[2]).

---

5. Converging the results

If you intend to obtain high-quality spectra, it is necessary to find a proper set of parameters leading to reliable results. The choice of the values to be assigned to the parameters listed below is crucial for the accuracy of the calculation:

Attribute Element Description ngridk xs The quality of the spectra depends on the size of the k-mesh. A denser k-mesh results in a better resolution of the spectrum. Thus, always check for convergence with respect to the number of k points. The ngridk attribute inside the xs element is independent of the ngridk inside the groundstate element. By default, the TDDFT part of exciting triggers a one-step groundstate calculation from the saved electronic density and potentials (for other possibilities on how to trigger the groundstate calculation, see [1]), but using its own parameters (e.g., ngridk, rgkmax, etc.) which do not conflict with those of the groundstate calculation. rgkmax xs It determines the size of the basis set and, therefore, influences the quality of the eigenvalues. When increasing it, the computational time rises rapidly. The rgkmax attribute inside the xs element is independent of the rgkmax inside the groundstate element. It is highly recommended to assign explicitly a value to rgkmax in any calculation. gqmax xs This is the value of the LFE cutoff Gmax. If it is set to zero, local-field effects are neglected. nempty xs It determines the number of empty states used in the calculation of the KS independent-particle response function $\chi_0$ (which is in turn used for the calculation of the interacting $\chi$ through the solution of the Dyson equation as explained above). It is to some extent related to the energy range covered by your spectra and also by the LFE cutoffs: In some case, a larger value of gqmax may require a larger value for nempty. Thus, you always have to converge with respect to nempty after selecting a new gqmax.

For example, the loss spectra calculated for a finer k-grid of 16$\times$16$\times$16 points are shown below:

As you can see, while the position and relative intensity of the main peaks the do not change significantly using a finer grid, the whole spectrum is smoother and better resolved when obtained with a larger number of k-points.

Exercise

• Check the convergence of the loss spectra with respect to the other computational parameters listed above. Be aware that more accurate calculations are also more computationally demanding.

---

6. Try a different kernel

Alternatively, the loss function can be calculated using the adiabatic LDA (ALDA) exchange-correlation (xc) kernel, being the simplest non-trivial xc kernel based on an LDA potential for the time-dependent case. The ALDA approximation consists in calculating the xc potential by taking the static xc potential evaluated at the time-dependent density. Such a calculation can be performed from scratch or on top of the latter one by changing the fxctype attribute in the element tddft.

... 
   <tddft
      ...
      fxctype="ALDA"
      ...>
...

Moreover, the Kohn-Sham response function and the matrix elements do not have to be recalculated upon a change in the xc kernel. This can be avoided setting the do attribute to fromkernel inside the tddft element.

...
   <tddft
      do="fromkernel" 
      ...>
...

After doing these small changes to the input.xml file, lets run exciting again.

$ time excitingser >& outputXS.txt &

Plots can be generated from the LOSS_xxx.xml files using PLOT-loss-function.py, as discussed for the RPA case, with the following commands

$ PLOT-loss-function.py LOSS_FXCRPA_OC11_QMT001.OUT.xml LOSS_FXCALDA_OC11_QMT001.OUT.xml

Here both results, from RPA and ALDA, are compared for a coarse 8$\times$8$\times$8 k-mesh.

As one can see, in this case using a different xc kernel affects the loss spectra only in the high-energy region to a minor extent.

For more detailed information on the different kernels available in TDDFT, we refer the reader to the tutorial TDDFT Calculations with Different Kernels.

More details on the implementation of the TDDFT formalism within the LAPW method can be found in Refs.[3] and [4]. A detailed study of the loss function of Silver using exciting can be found in Ref.[5].

---

Bibliography

1. The groundstate calculation can be also triggered from scratch from the xs element with the dogroundstate attribute of xs, this is indeed needed in spin polarized calculations. For an example refer to Magneto-optical Kerr effect (MOKE).

2. "Local Field Effects in the Electron Energy Loss Spectra of Rutile TiO₂", Phys. Rev. Lett. 88, 037601 (2002).

3. S. Sagmeister, PhD thesis, University of Graz, August 2009 (PDF).

4. Tutorial talk by Stefan Sagmeister (PDF) at the HoW exciting! 2010 workshop in Lausanne.

5. Further information on the momentum-dependent loss function of Ag: A. Alkauskas, S. Schneider, S. Sagmeister, C. Ambrosch-Draxl, and C. Hèbert, Theoretical analysis of the momentum-dependent loss function of bulk Ag, Ultramicroscopy 110, 1081 (2010) (PDF).