Excited States from TDDFT

for helium 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 momentum-dependent loss function of bulk silver is studied.



0. Preparations…

Be sure to have at hand:

  • The exciting root directory, EXCITINGROOT, that is, simply the path to the exciting installation directory: (e.g. /home/lisa/programs/exciting …). Note that the root directory can also be a relative path (e.g. ../../).

The EXCITINGROOT variable has to be replaced in the subsequent input files by its actual value (the path itself).
The directories for examples, program binaries, species files and visualization templates are then located at

EXCITINGROOT/examples
EXCITINGROOT/bin
EXCITINGROOT/species
EXCITINGROOT/xml/visualizationtemplates

Note that the shell prompt is denoted by a Dollar sign

$

First the variable $EXCITINGROOT is defined

$ echo 'EXCITINGROOT=_exciting_root_directory_' >> ~/.bashrc; . ~/.bashrc

where _exciting_root_directory_ contains the actual path of the exciting installation directory.
For the execution of the binaries it is required to add the bin-directory to your path.
$ echo 'PATH="$PATH:$EXCITINGROOT/bin"' >> ~/.bashrc; . ~/.bashrc

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

1. Ground state calculation

As a starting point for the TDDFT calculation we need a converged density and potential. To this end we perform a groundstate calculation first considering the following input:

<?xml version="1.0" encoding="UTF-8"?>
<?xml-stylesheet href="inputtohtml.xsl" type="text/xsl"?>
 
<input xsi:noNamespaceSchemaLocation="EXCITINGROOT/xml/excitinginput.xsd"
    xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsltpath="EXCITINGROOT/xml/"
    scratchpath="/tmp">
 
    <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>

To perform the actual calculation, copy and paste the text from above into a file called

input.xml

within a directory of your choice and start the calculation by invoking the exciting executable (serial version) in the 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 ground state has finished, check the output.txt for warning or error messages and the INFO.OUT file if the calculation has finished properly.

Although the ground state run produced several output files, only two are needed as starting point for the TDDFT calculation. These files contain the Fermi level as well as the density and potential (the “state”):

EFERMI.OUT
STATE.OUT

2. TDDFT calculation

Symmetries

Before we actually start with the calculation, some initial knowledge about the structure of the dielectric function tensor - and hence the loss function or the dynamical structure factor - can be obtained from symmetry considerations. After the ground state calculation finishes, it leaves one file called

SYMT2.OUT

containing the structure of a general rank-2 tensor (without taking into account magnetism and spin) after a symmetry average with respect to the crystal symmetries of the system under study. It contains an upper bound for the number of independent components and its matrix structure is depicted. (The matrices listed below in the remaining part of this file shall not be of any importance for us now.)
For instance, this output could look like the following (for a lower-dimensional system)
(structure of a rank 2 tensor, e.g. the dielectric tensor, from symmetry considerations
 with respect to Cartesian coordinates)

 Upper limit for number of independent components :      5

 ( e_11   0    e_13 )
 (  0    e_22   0   )
 ( e_31   0    e_33 )

For systems with cubic symmetry, however, we expect only one independent component.

Work-flow

For the actual TDDFT application we want to calculate the loss function, related to a finite momentum transfer process.
Before going into detail, the work-flow of a TDDFT linear-response calculation is discussed.
The eigenvalues and wavefunctions are determined for a given k-mesh. The momentum matrix elements are then calculated. In a next step the ${\bf 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 ${\bf Q}=0$ all diagonal tensor components are calculated per default.

TDDFTsimple.png

In this flow chart, each level along a horizontal axis denotes quantities which are independent from one another.
Fortunately, this work-flow is already built into the program and the user can control the work-flow from within one input file.

The groundstate calculation can now be skipped by adding the do attribute setting it to skip:

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

Below a choice of parameters for the excited states, the <xs> element is given. This block has to be inserted inside the <input> element.

<xs xstype="TDDFT"
     ngridk="4 4 4" vkloff="0.05 0.15 0.25"
     nempty="60"
     lmaxapwwf="5" lmaxemat="5"
     gqmax="2.0"
     broad="0.004"
     tevout="true">
        <energywindow intv="0.0 2.5" points="500" />
        <tddft fxctype="RPA"/>
        <qpointset>
            <qpoint> 0.0 0.0 0.1 </qpoint>
        </qpointset>
</xs>

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

A few additional remarks

  • The k-point mesh 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 k-point mesh.
  • The symmetry relations between the k-points are not exploited in the code, but for the case of optics (Q=0) and the absence of local field effects. In any case it is a common practice to shift the k-mesh off symmetry. This means to shift all points of the mesh by a small displacement $k_\Delta$. Such a displacement should break all symmetry relations among the k-points of the mesh. This way all the k-points resulting from a regular mesh, which has been shifted, are crystallographically inequivalent and there are no redundant contributions to the spectrum.
  • The Q-vector is defined in the <qpointset> element in the input file and 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 lmaxapwwf parameter is reduced since the full lmaxapw value (the default) is rarely needed.
  • The lmaxemat parameter controls the Rayleigh expansion of the $e^{-i{\bf (q+G)r}}$ factor and is slightly increase for more stability in the high frequency region for small but finite Q-vectors.
  • The xc kernel is set to "RPA". It can be changed to "ALDA" at a later stage.

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 as well as the resume file which indicates if the program is currently running, and, finally, the terminal output file:

INFOXS.OUT
resume
outputXS.txt

Once the calculation is finished no resume file should be present in the directory. However, a bunch of 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.
We will concentrate on the LOSS_* files and pick the XML output file: LOSS_FXCRPA_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 here.

Visualizing the output

From the XML data file for the loss function we can generate xmgrace files using the xmlloss2agr.xsl template from the EXCITINGROOT/xml/visualizationtempaltes directory:

$ xsltproc --stringparam filename lossfunction $EXCITINGROOT/xml/visualizationtemplates/xmllossfunction2agr.xsl LOSS_FXCRPA_QMT001.OUT.xml

This template will create two xmgrace files, namely

lossfunction_Loss.agr
lossfunction_DynSfac.agr

Right now you have plotted the spectrum including local field effects. You can also apply this template to the loss function file without local field effects.

In order to view these files, invoke the xmgrace program (see xmgrace-quickstart for help on xmgrace)

$ xmgrace lossfunction_Loss.agr

Right now you have plotted the spectrum including local field effects. To plot the corresponding spectrum without

Here the result for a converged 30x30x30 k-mesh is displayed. The result obtained from the coarser mesh from the input above should, however, resemble the main features although it will be not as smooth.

Ag-RPA-LFE-NLF.png

In the figure above, LFE stands for local field effects and NLF for no local field (effects).
One immediate conclusion from the result is the importance of the local field effects, especially in the range above 35 eV. They drastically reduce the oscillator strength.

Converging the results

The spectra obtained from TDDFT for a given xc kernel depend on many parameters - on some of them although in a crucial way.

Description Parameter Info
k-mesh input/xs/@ngridk very crucial and should always be converged
basis set size input/xs/@rgkmax defines the quality of eigenvalues and wavefunctions
energy cutoff (number of empty states) input/xs/@nempty related to the energy range of the spectrum
local field effects G-cutoff input/xs/@gqmax can give rise to large effects depending on the system
broadening input/xs/@broad Lorentzian broadening, be careful it is in Hartrees

Changes of the input parameters give rise to changes in the computation time. A rough estimate for the scaling of the Kohn-Sham response function, which is the most time-consuming part of the calculation in general, is the following

(1)
\begin{align} T_{\rm TDDFT} \sim N_{\bf k}N_{\bf G}^2 N_{\omega}N_v N_c \nonumber \end{align}

where $N$ stands for the number of and the subscripts denote the k-points, G-vectors, frequencies, valence states and conduction states.

Excercise

Once you successfully run the previous TDDFT calculation you can play around with the input parameters to achieve a better quality for the loss function.

  • Try to gently increase parameters from the table above, except the basis set size which should be sufficiently large.
  • Be careful when increasing the local field effects cutoff, since the calculation might then quickly be very time-consuming.
  • It is therefore suggested to get a feeling for the k-mesh first.

A different kernel

Alternatively, the loss function can be calculated using the ALDA xc kernel, being the simplest non-trivial xc kernel based on an LDA potential for the time-dependent case, which is the static xc potential evaluated at the time-dependent density.
Such a calculation can be performed from scratch or on top the latter one by changing the fxctype attribute in <tddft> element of the input file:

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

The xmgrace file is again generated from the LOSS*.xml files via the XSL template discussed for the RPA case.
Here both results, from the RPA and the ALDA, are compared for a converged 30x30x30 k-mesh:

Ag-RPA-ALDA.png

Excercise

Once you came that far it should be no problem to do the same type of calculation for a different systems. With a tiny set of changes to the input we can perform a calculation for bulk gold (Au) if we know:

  • the lattice constant (can be found e.g. from http://www.webelements.com), use the scale attribute and note that inputs are in bohr!
  • that the spacegroup is the same as for Ag
  • the species file of Au

Try to do the calculation for Au and do some convergence tests on your own. Compare the results for the spectra to those of Ag.

Literature

Tutorial talk 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