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

0. Define relevant shell variables and download scripts

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.

$ excitingser >& output.txt &

Meanwhile, you can check the bunch of files created during the run, especially the INFO.OUT file for convergence information. After the groundstate calculation has finished check the last lines of the INFO.OUT file. If the self-consistent cycle was successfully terminated, you should find the message

EXCITING Lithium stopped

Check the output.txt file for warning or error messages which should be dumped there. Note: In the INFO.OUT file you can verify that the convergence critetria have been met by looking at the output lines just before the block regarding the last iteration.

Although the groundstate run produce several output files, only two files are needed as a starting point for the TDDFT calculation. These are the files EFERMI.OUT and STATE.OUT containing the Fermi level and the converged electron density and potential, respectively.

Please notice: 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. 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 $L(\mathbf{q},\omega)$ is defined as the imaginary part of the inverse dielectric tensor.

(1)
\begin{align} L(\mathbf{q},\omega)=-\Im\left(\frac{1}{\varepsilon(\mathbf{q},\omega)}\right) \end{align}

Before going into detail, we discuss in the next the work-flow of a TDDFT linear-response calculation.

Work-flow of a TDDFT linear-response calculation

The eigenvalues and wavefunctions of the Kohn-Sham equations are determined for a given k-mesh. Momentum matrix elements are then calculated. In the next step, the Q-dependent matrix elements are calculated. Next, the Kohn-Sham response function is assembled with the help of the latter matrix elements. In a last step the Dyson equation for the interacting response function is solved and the macroscopic dielectric function is obtained therefrom. For Q=0, all diagonal tensor components are calculated by default.

TDDFTsimple.png

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

Setting elements and attributes

The groundstate calculation can now be skipped by adding the do="skip" attribute to the element groundstate.

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

In order to perform the calculation for 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.05 0.15 0.25"
     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>
...

The input file corresponds to a calculation including local-field effects. To neglect local-field effects in the calculation, the gqmax attribute should be set to zero
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.
  • The symmetry relations between the k-points are not exploited in the code, except for the optical case (Q=0) in absence of local field effects. 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 0). Note that the Q-vector is given in crystal coordinates.
  • The exchange correlation kernel (specified by the fxctype attribute) is set to RPA (random-phase approximation). It will be changed to ALDA (adiabatic local-density approximation) at a later stage.
Running exciting

Having modified the input file according to the above prescription, you should now be able to launch the TDDFT calculation.

$ excitingser >& outputXS.txt &
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.
PLOT-loss-comp.png

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.
PLOT-RPA-ALDA.png

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)
Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License