by Ronaldo Rodrigues Pela for exciting oxygen
Purpose: In this tutorial, you will learn how to perform an ilustrative realtime timedependent densityfunctional theory (RTTDDFT) calculation. We consider diamond as an example.
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.
 EXECUTEsingle.sh: (Bash) shell script for running a single exciting calculation.
 PLOTmultitask.py: Python script for making plots from output files of RTTDDFT.
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! For this tutorial, it can be useful to remind the conversion between a.u. and the SI for times: 1 a.u. = 0.02418884254 fs.
1. Theoretical background for RTTDDFT calculation
In TDDFT, we consider a system of electrons submitted to an external timedependent (TD) perturbation, expressed usually in terms of an applied electric field $\mathbf{E}(t)$. This leads to a TD KohnSham (KS) hamiltonian $\hat{H}(\mathbf{r},t)$:
(1)where $\mathbf{A}(t)$ is the vector potential $\mathbf{A}(t)=c\int_0^t \mathbf{E}(t')\,dt'$, $c$ is the speed of light, and $v^\phantom{I}_\textrm{KS}(\mathbf{r},t)$ is the TDKS potential, a sum of the TD ionic, Hartree and exchangecorrelation (XC) potentials. We assume in exciting the adiabatic approximation for the TDXC potential: $v^\phantom{I}_\textrm{XC}(\mathbf{r},t)$ is constructed using $n(\mathbf{r},t)$, the TD electronic density, the same way that $v^\phantom{I}_\textrm{XC}(\mathbf{r},0)$ in the groundstate is constructed from the groundstate density $n(\mathbf{r},0)$.
A KS wavefunction $\psi_{j\mathbf{k}}(\mathbf{r},t)$ labeled with index $j$ and associated to a wavevector $\mathbf{k}$ evolves as
(2)where the propagator $\hat{U}(t+\Delta t,t)$ is given by:
(3)$\hat{\mathcal{T}}$ being the timeordering operator. The following propagators are implemented:
 simple exponential,
 exponential at midpoint rule,
 approximate enforced timereversal symmetry,
 CommutatorFree Magnus expansion of 4th order,
 exponential using a basis of the hamiltonianeigenvectors, and
 RungeKutta of 4th order.
Please refer to Ref.[1] for the expressions of $\hat{U}(t+\Delta t,t)$.
As consequence of the exciting electric field, electrons are moved away from their equilibrium position. This gives origin to a macroscopic current density $\mathbf{J}(t)$, which for local and semilocal KS functionals, can be obtained as
(4)where $N$ stands for the number of valence electrons in the unit cell with volume $\Omega$, $w_{\mathbf{k}}$ is the weight of the considered kpoint, and $f_{j\mathbf{k}}$ is the occupation number of the corresponding KS state.
Another quantity of interest is the number of excited electrons (per unit cell) $N_\textrm{exc} (t)$ at a given time. It is possible to describe them by projecting $ \psi_{i\mathbf{k}}(t)\rangle$ onto the ground state wavefunctions. For a given kpoint, the number of electrons that have been excited to an unoccupied KS state, labeled with $j$, is
(5)whereas the number of holes created in an occupied KS state $j'$ is
(6)Thus, $N_\textrm{exc} (t)$ can be obtained by considering all the unoccupied states
(7)2. Groundstate calculation
i) Preparation of the input file
The first step is to create a directory for the system that we want to investigate, diamondrt, and to enter it.
$ mkdir diamondrt
$ cd diamondrt
As initial condition for the time evolution done in RTTDDFT, we need the electron density and potential from a groundstate calculation. Therefore, we create the file input.xml that could look like the following:
<input> <title>Diamond</title> <structure speciespath="$EXCITINGROOT/species/"> <crystal scale="6.7407"> <basevect>0.0 0.5 0.5</basevect> <basevect>0.5 0.0 0.5</basevect> <basevect>0.5 0.5 0.0</basevect> </crystal> <species speciesfile="C.xml" rmt="1.4"> <atom coord="0.00 0.00 0.00"/> <atom coord="0.25 0.25 0.25"/> </species> </structure> <groundstate do="fromscratch" xctype="LDA_PW" ngridk="8 8 8" rgkmax="5.0d0" nempty="5" epsengy="1.d8"> </groundstate> </input>
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:
$ SETUPexcitingroot.sh
ii) Running the groundstate calculation
You can start the calculation by invoking the script EXECUTEsingle.sh.
$ EXECUTEsingle.sh test1
After the calculation is finished, move to the subdirectory test1, into where all the output files are written.
$ cd test1
You can check the 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 RTTDDFT calculation.
Please note: To obtain reliable groundstate density and potential, it is necessary to perform convergence tests with respect to the kpoint mesh (parameter ngridk) and the size of the basis set (parameter rgkmax). For details see Simple convergence tests.
3. Performing a RTTDDFT calculation
i) Setting elements and attributes
You can prepare the input file (input.xml) of a RTTDDFT calculation starting from the existing groundstate input file. The explicit groundstate calculation can be avoided with the attribute do = "skip" inside the element groundstate.
... <groundstate do="skip" ... > ... </groundstate> ...
In order to perform a "RTTDDFT" calculation, the element xs must be added to the input file inside the input element, as shown below.
... <xs xstype="RTTDDFT" ngridk="4 4 4" rgkmax="5.0d0" vkloff="0.01 0.02 0.004" nempty="5" nosym="true" reducek="false"> <realTimeTDDFT propagator="AETRS" timeStep="0.25d0" endTime="50.d0" printTimingGeneral="true" calculateNExcitedElectrons="true" printAfterIterations="1"> <laser> <trapCos amplitude="1.0d0" omega="1.0d0" phase="0.d0" t0="0.25d0" riseTime="5.d0" width="30.d0" direction="x" /> </laser> </realTimeTDDFT> </xs> ...
These settings in the input file correspond to a calculation with the following main features:
 It is a RTTDDFT calculation (the attribute xstype is set to "RTTDDFT").
 The attribute ngridk inside the xs element determines the kgrid used for the KS wavefunctions $ \psi_{i\mathbf{k}}(t)\rangle$ evolved as given by Eq. (2). This kgrid is also affected by vkloff inside xs, which, in turn, specifies an offset (given in lattice coordinates) for the grid generation. Additionally, nosym=true and reducek=false inside xs enforce that no symmetry and no kpoint reduction is to be used ([2]).
 rgkmax inside xs determines the size of the basis for the expansion used for $ \psi_{i\mathbf{k}}(t)\rangle$.
 The attribute nempty inside xs specifies the number of unoccupied states also considered in the timeevolution in Eq. (2).
 The attributes ngridk, rgkmax, vkloff, nempty, nosym, and reducek given in the xs element are independent of their counterparts inside the groundstate element.
 Before evolving KS wavefunctions, a onestep groundstate calculation is carried out using ngridk, rgkmax, vkloff, nempty, nosym, and reducek given inside xs.
 However, the saved electronic density and the KS potential from STATE.OUT, used for this onestep calculation, had been obtained beforehand using the parameters ngridk, …, reducek inside the groundstate element.
 The element realTimeTDDFT determines the specific parameters for the RTTDDFT calculation.
 Here, we consider a propagator an expression which approximately enforces timereversal symmetry (propagator="AETRS"), with timestep of 0.25 atomic units (a.u.), timeStep="0.25d0", and the time evolution is carried out up to the end time of 50 a.u. (endTime="50.d0").
 The meaning of calculateNExcitedElectrons="true" is to trigger the calculation of the number of the excited electrons. They are printed out in the file NEXC.OUT.
 printTimingGeneral="true" triggers the measurement of the timings during the execution of the subroutines called in the RTTDDFT calculation. The file TIMING_RTTDDFT.OUT stores such information. More detailed information is available adding printTimingDetailed="true".
 The element laser describes the external field applied to the system. Here we consider a cosine modulated by a trapezoid, as given in the element trapCos. The expression for the corresponding vector potential is
where $f(t)$ is the trapezoidal function:
(9)which is illustrated in the following figure.
In our example, the field has an amplitude of 1.0 a.u. (amplitude="1.0d0") and is applied along the xdirection (direction="x"). The cosine function has an angular frequency of $\omega = 1.0$ a.u. (omega="1.0d0") and phase $\phi=0.0$ (phase="0.d0"). The trapezoid function has parameters: initial time $t_0=0.25$ a.u., rise time of $t_r=5.0$ a.u., and width of $w=30$ a.u. (as given by: t0="0.25d0", riseTime="5.d0", width="30.d0").
You can check all relevant parameters for a "RTTDDFT" calculation in Input Reference.
ii) Running exciting
Having modified the input file as described above, you can launch the RTTDDFT calculation inside the subdirectory test1 using.
$ time exciting_smp
Be aware that, depending on the computer you are using, this calculation may take some minutes.
iii) Output files
In our example, the onestep calculation performed before the evolution of the KS wavefunctions generates output files ending with "_RTTDDFT.OUT. They contain the same information as their counterparts generated within a usual selfconsistent cycle triggered by the groundstate element. The only exception is "TIMING_RTTDDFT.OUT, which brings timings required to execute different subroutines invoked in the evolution of KS wavefunctions.
Other output files generated specifically by the RTTDDFT calculations are:
 JIND.OUT: with the three components of the current density. The time is written in the first column, and the next three columns contain the $x$, $y$ and $z$ components of the current density.
 PVEC.OUT, into where the polarization vector is written. The convention for the four columns is the same of JIND.OUT.
 AVEC.OUT, which contains the vector potential. The time is written in the first column, and the next six columns contain the $x$, $y$ and $z$ components of the induced and the total vector fields in the following order: $A_{\textrm{ind},x}$, $A_{x}$, $A_{\textrm{ind},y}$, $A_{y}$, $A_{\textrm{ind},z}$, $A_{z}$ (see Ref.[4]).
 NEXC.OUT, where you find the number of excited electrons per unit cell. The first column contains the time; the second, the number of electrons per unit cell on the ground state; the third, the number of excited electrons per unit cell; the fourth, the sum of the previous two columns.
4. Visualizing the outputs
To visualize the current density, type the following command inside the subdirectory test1:
$ PLOTmultitask.py jind f JIND.OUT x
This plots the x component of the current density [3]. Your figure should be similar to the following one:
By using the command
$ PLOTmultitask.py nexc f NEXC.OUT
you can plot the number of excited electrons per unit cell, $N_\textrm{exc}(t)$, as a function of time. Bellow, we show how your plot is expected to be:
5. Converging the results
If you intend to obtain highquality results, it is necessary to converge them with respect to a proper set of reliable parameters. In RTTDDFT, typically they are:
What regards the timestep, calculations are typically not too sensitive to it. Usually, there is a critical value $t_{cr}$ and when the employed timestep is larger than $t_{cr}$, a sort of divergence is observed. However, below $t_{cr}$ the results tend to be very similar and almost independent of the timestep (see discussion in Ref.[1]).
Below is, e.g., the current density for a finer kgrid of 16$\times$16$\times$16 points compared to 4$\times$4$\times$4 as used in this tutorial (see [5]):
6. Other expressions for the laser field
Apart from a cosine function modulated by a trapezoid, as considered in this tutorial, it is also possible the following fields:
 an impulse (deltakick)
 a cosine with a sin squared as envelope (which mimics a laser with a Gaussian envelope, employed in many experiments).
These other functions are explored in the tutorials Simulating pumpprobe spectroscopy with RTTDDFT and Studying higherharmonic generation using RTTDDFT.
Impulsive function
In this case, the field is
(10)The next figure depicts the vector potential as a function of time.
The electric field is
(11)This is schematically shown in the figure bellow.
With the figure, it is ease to verify that this function converges to $\mathbf{E}_0 \,\delta(tt_0)$ in the limit $w \to 0$. It is also possible to specify $w=0$. In this case, the vector potential is equal to 0, for times not larger than $t_0$. And for the smallest time larger than $t_0$, the vector potential assumes a value of $c\,\mathbf{E}_0$ (this smallest time depends clearly on the timestep employed in the calculation).
The corresponding entry for this field in the input.xml file is
... <laser> <kick t0="1.d0" width="0.1d0" amplitude="0.01d0" direction="z"/> </laser> ...
The meaning of the parameters is:
 $t_0$ is specified by the attribute t0;
 the width of the impulse function, $w$, is specified by the attribute width;
 the attribute amplitude determines the amplitude $E_0$;
 the attribute direction defines the direction of the field.
Cosine modulated by a sin squared
In this case, the field has the same expression as in Eq. (8), but with the envelope given as
(12)The following figure illustrates the considered field
The parameter $w$ determines the width of the pulse, whereas $t_0$, the instant at which the pulse is applied.
The entry for this field in input.xml is
... <laser> <sinSq t0="0.5d0" omega="1.d0" phase="0.d0" amplitude="1.0d0" pulseLength="5.d0" direction="y"/> </laser> ...
The meaning of the parameters is:
 $t_0$ is specified by the attribute t0;
 the angular frequency $\omega$ and phase, given in Eq. (8) are determined by the attributes omega and phase, respectively;
 the width of the pulse, $w$, is specified by the attribute pulseLength;
 the attribute amplitude determines the amplitude $A_0$, as given in Eq. (8);
 the attribute direction defines the direction of the field.
Combining fields
It is possible to combine fields by adding more elements to laser element. Actually, all the entries in this element are added up to build the resulting vector potential. To use the example considered in this tutorial, this entry:
... <laser> <trapCos amplitude="1.0d0" omega="1.0d0" phase="0.d0" t0="0.25d0" riseTime="5.d0" width="30.d0" direction="x" /> <trapCos amplitude="1.0d0" omega="1.0d0" phase="0.d0" t0="0.25d0" riseTime="5.d0" width="30.d0" direction="y" /> </laser> ...
will produce a resulting field $\mathbf{A}(t) = A_\textrm{tutorial}(t)\,\hat{x}+A_\textrm{tutorial}(t)\,\hat{y}$, where we are using here $A_\textrm{tutorial}$ to denote the field employed in this tutorial. Here we are using the same field and different directions. But it is also possible to combine different kinds of fields.
7. Reference
The details on the implementation of realtime RTTDDFT within the exciting code can be found in Ref.[1].
$ PLOTmultitask.py jind f test1/JIND.OUT kpt16/JIND.OUT x c '4x4x4' '16x16x16'