**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

`executable (serial version) in the background.`

**exciting**`$ 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

`for warning or error messages and the`

**output.txt**`file if the calculation has finished properly.`

**INFO.OUT**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.

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

`element.`

**<input>**<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

`label. This stands for`

**_QMTxxx***Q momentum transfer*and corresponds to the label

`of the Q-points listed in the`

**xxx**`file.`

**QPOINTS.OUT**We will concentrate on the

`files and pick the XML output file:`

**LOSS_***`.`

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

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

`element of the input file:`

**<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

`inside the`

**fromkernel**`element.`

**<tddft>**<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:

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