Excited States from TDDFT

for lithium version of Exciting

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).

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

Before starting, be sure that relevant shell variables are already defined and that the excitingscripts directory has already been downloaded, as specified in Tutorial scripts and environment variables. Here is a list of the scripts which are relevant for this tutorial with a short description.

• 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 silver_tddft and we move inside it. $ mkdir silver_tddft
$cd silver_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>

<!--
Below is the groundstate part of the input
-->
<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.

Be sure to set the correct path for the exciting root directory (indicated in this example by $EXCITINGROOT). This can be done using the following line command. $ SETUP-excitingroot.sh

##### ii) Running a groundstate calculation

Start the calculation by invoking the exciting executable (serial version) in background.


##### Output files

While the calculation is running you can check the info file for the excited states (INFOXS.OUT) and the terminal output file (outputXS.txt). 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.

### 3. Visualizing the output

From the XML data file for the loss function, we can generate the files PLOT.ps and PLOT.png typing the following command inside the directory silver_tddft.

$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 online entry which is the name of an xml output file. You can add as many files 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

In the figure below, thick lines show the result for a converged 30$\times$30$\times$30 k-mesh. Also shown is the result for the coarser mesh from the input above (thin lines). As you can see, the results for the coarser grid resemble the main features although it is not as smooth.

In the figure above, the labels LFE and no-LFE indicate calculation performed including and excluding local-field effects. 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. Phys. Rev. Lett. 88, 037601 (2002)].

### 4. Converging the results

Before starting a series of calculations to investigate the momentum-dependence of different optical quantities it is necessary to find a proper set of parameters leading to reliable results. The choice of the parameters listed below is crucial for the accuracy of the calculation.

• ngridk: 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.
• swidth: In particular for metals, the convergence test with respect to the smearing width swidth must go hand in hand with that of ngridk. The swidth parameter is necessary to smooth the convergence of the self-consistent DFT loop.
• rgkmax: It determines the size of the basis set and, therefore, influences the quality of the eigenvalues. Be careful when increasing it, since the computational time rises rapidly.
• nempty: It determines the number of empty states considered in the calculation and depends on the energy range of interest.
• gqmax: Local-field effects G-cutoff. If set to zero, local-field effects are neglected.

### 5. 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 approximation ALDA 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.

$excitingser >& outputXS.txt &  Plots can be generated from the LOSS*.xml files using PLOT-loss-function.py, as discussed for the RPA case, with the following command $ PLOT-loss-function.py LOSS_FXCRPA_OC11_QMT001.OUT.xml LOSS_FXCALDA_OC11_QMT001.OUT.xml

Here both results, from the RPA and the ALDA, are compared for a converged 30$\times$30$\times$30 k-mesh.

### Literature

• Tutorial talk by Stefan Sagmeister (PDF)
• 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)
• More details on the implementation of the TDDFT formalism within the LAPW method: S. Sagmeister, PhD thesis, University of Graz, August 2009 (PDF)