by Stephan Sagmeister, Santiago Rigamonti & Ronaldo Rodrigues Pela for exciting neon
(Jupyter notebook by Megha Arya and Mara Voiculescu)
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 running any Jupyter tutorials, please refer to the 00_before_starting.md
document on how to correctly set up the environment. This only needs to be done once. After which, the venv can be (re)activated from exciting
's root directory:
source tools/excitingjupyter/venv/excitingvenv/bin/activate
Important note: All input parameters that will appear will be given in atomic units!
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(q,ω) is defined as
\begin{equation} \tag{1} L(\mathbf q, \mathit ω) = -Im \Big[\varepsilon_{M}(\mathbf q,\mathit ω)\Big]^{-1}, \end{equation}where $\varepsilon_{M}(\mathbf q,\mathit ω)$ is the macroscopic dielectric function which, in turn, is related to the inverse dielectric matrix by
\begin{equation} \tag{2} \Big[\varepsilon_{M}(\mathbf q,\mathit ω)\Big]^{-1} = \varepsilon_{00}^{-1}(\mathbf q,\mathit ω). \end{equation}The inverse dielectric matrix is given by $\varepsilon^{-1} = 1 + \nu \chi$, where $\nu$ 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
\begin{equation} \tag{3} \chi = \chi_{0} + \chi_{0} (\nu + f_{xc}) \chi , \end{equation}where $\chi_{0}$ is the Kohn-Sham independent-particle polarizability and $f_{xc}$ is the exchange-correlation kernel. The neglect of the $f_{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 *X*G,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).
The calculation of the different quantities involved determines the work-flow of a TDDFT linear-response calculation, as discussed in the following.
Eigenvalues $\varepsilon_{n \mathbf k}$ and wavefunctions $\psi_{n \mathbf k} (\mathbf r)$ of the Kohn-Sham equations are determined for a given k-mesh.
Momentum matrix elements, $\rho_{α,nm \mathbf k} = < \psi_{n \mathbf k} | \rho_{α} | \psi_{m \mathbf k} >$, are calculated.
q-dependent matrix elements of plane waves, $M_{nm \mathbf k}(\mathbf q,\mathbf G) = < \psi_{n \mathbf k} | exp[-i(\mathbf q + \mathbf G) \mathbf r] | \psi_{m \mathbf k+q} >$, are evaluated.
The Kohn-Sham response function $\chi_{0}$ is assembled including the matrix elements
In a last step the Dyson equation for the interacting response function $\chi$ is solved and the macroscopic dielectric function $\varepsilon_{M}$ is obtained therefrom. For q=0, all diagonal tensor components of the macroscopic dielectric function are calculated by default.
In the following flow chart, each level along a horizontal axis denotes quantities which are independent from each other.
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 run_Ag_tddft.
%%bash
mkdir -p run_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: Ground-State Calculation</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>
Do not forget to replace in the input.xml 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:
%%bash
cd run_Ag_tddft
python3 -m excitingscripts.setup.excitingroot
cd ..
ii) Running the Ground-State Calculation
You can start the calculation by invoking the script excitingscripts.execute.single
.
%%bash
python3 -m excitingscripts.execute.single -r run_Ag_tddft
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 NEON 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 the tutorial Simple convergence tests.
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
.
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:
xstype
is set to "TDDFT").fxctype
of the subelement tddft
is set to "RPA").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
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.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.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.broad
is a Lorentzian broadening for all spectra.tevout
sets the energy output to be given in electronvolt (eV).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.
%%bash
python3 -m excitingscripts.execute.single -r run_Ag_tddft
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 NEON 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. 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.
From the data file for the loss function, we can generate a PNG output by executing the following command:
%%bash
cd run_Ag_tddft
python3 -m excitingscripts.plot.files -f LOSS_FXCRPA_OC11_QMT001.OUT -lx 'Energy [eV]' -ly 'Loss function' -x 0 54 -y 0 1.4
cd ..
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. 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.
%%bash
cd run_Ag_tddft
python3 -m excitingscripts.plot.files -f LOSS_FXCRPA_OC11_QMT001.OUT LOSS_NLF_FXCRPA_OC11_QMT001.OUT -lx 'Energy [eV]' -ly 'Loss function' -x 0 54 -y 0 1.4 -ll 'RPA (LFE)' 'RPA (no-LFE)' -pn Loss_RPA
cd ..
The resulting plot is saved as Loss_RPA.png 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]).
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:
Attribute | Element | Description |
---|---|---|
ngridk |
xs |
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. The ngridk attribute inside the xs element is independent of the ngridk inside the groundstate element. By default, the TDDFT part of exciting triggers a one-step ground-state calculation from the saved electronic density and potentials (for other possibilities on how to trigger the ground-state calculation, see [1], but using its own parameters (e.g., ngridk , rgkmax , etc.) which do not conflict with those of the ground-state calculation. |
rgkmax |
xs |
It determines the size of the basis set and, therefore, influences the quality of the eigenvalues. When increasing it, the computational time rises rapidly. The rgkmax attribute inside the xs element is independent of the rgkmax inside the groundstate element. It is highly recommended to assign explicitly a value to rgkmax in any calculation. |
gqmax |
xs |
This is the value of the LFE cutoff Gmax. If it is set to zero, local-field effects are neglected. |
nempty |
xs |
It determines the number of empty states used in the calculation of the KS independent-particle response function $\chi_{0}$ (which is in turn used for the calculation of the interacting $\chi$ through the solution of the Dyson equation as explained above). It is to some extent related to the energy range covered by your spectra and also by the LFE cutoffs: In some case, a larger value of gqmax may require a larger value for nempty . Thus, you always have to converge with respect to nempty after selecting a new gqmax . |
For example, you could repeat the calculations done before, in a different folder, for a finer k-grid of 16×16×16 points, by changing the ngridk
attribute accordingly. In this case, the obtained loss spectra 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
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.
%%bash
python3 -m excitingscripts.execute.single -r run_Ag_tddft
Plots can be generated as discussed for the RPA case, with the following commands
%%bash
cd run_Ag_tddft
python3 -m excitingscripts.plot.files -f LOSS_FXCRPA_OC11_QMT001.OUT LOSS_FXCALDA_OC11_QMT001.OUT -lx 'Energy [eV]' -ly 'Loss function' -x 0 54 -y 0 1.4 -ll 'RPA (LFE)' 'ALDA (LFE)' -pn RPA_ALDA_Comparison
cd ..
Here both results, from RPA and ALDA, are compared for a coarse 8×8×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 the tutorial 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].
xs
element with the dogroundstate
attribute of xs
, this is indeed needed in spin polarized calculations. For an example refer to the tutorial Magneto-optical Kerr effect (MOKE).