by Stephan Sagmeister, Santiago Rigamonti, & Ronaldo Rodrigues Pela for exciting oxygen
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).
Table of Contents
|
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. 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)where $\varepsilon^\phantom{I}_{\hspace{-0.4mm}\textrm{M}}\hspace{-0.4mm}(\mathbf{q},\omega)$ is the macroscopic dielectric function which, in turn, is related to the inverse dielectric matrix by
(2)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)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 Gmax in such a way that only matrix elements XG,G'(q) with both |q+G|, |q+G'| < Gmax are taken into account (X = $\chi_0$, $\chi$, etc.). If Gmax > 0 the calculation is said to include (fully or partially, depending on the size of Gmax) local field effects (LFE).
2. Ground-state 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 ground-state 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 xctype="GGA_PBE_SOL" 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 ground-state 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 OXYGEN 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.
3. How to perform a TDDFT calculation
Setting elements and attributes
If you have already performed the (ground-state) calculation shown in Section 2, you can prepare the input file (input.xml) of a TDDFT calculation starting from the existing ground-state input file ([1]) according to the following: First, the explicit ground-state 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 Gmax, 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 exciting_smp >& 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 OXYGEN 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 by 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.
Now, we will concentrate on the LOSS_ files and pick the output file LOSS_FXCRPA_OC11_QMT001.OUT. This file contains the data for the loss function and the dynamical structure factor. Details about the output files are described in TDDFT output files.
4. Visualizing the output
From the 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
As you can see, we are executing the PLOT-loss-function.py script with one argument, which is the name of one output file. You can add as many arguments (file names) in the command line as you like, and the script will add the corresponding curves to the same plot. If you wish to plot more than a file, do not forget to add in the command line the option —kernel. 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 --kernel LOSS_FXCRPA_OC11_QMT001.OUT LOSS_NLF_FXCRPA_OC11_QMT001.OUT
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:
|
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 exciting_smp >& outputXS.txt &
Plots can be generated from the LOSS_xxx.OUT files using PLOT-loss-function.py, as discussed for the RPA case, with the following commands
$ PLOT-loss-function.py --kernel LOSS_FXCRPA_OC11_QMT001.OUT LOSS_FXCALDA_OC11_QMT001.OUT
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 Many-body kernels for TDDFT calculations.
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].