by Ronaldo Rodrigues Pela for exciting neon
(Jupyter notebook by Mara Voiculescu and Martin Kuban)
Purpose: In this tutorial, you will learn how to perform an ilustrative real-time time-dependent density-functional theory (RT-TDDFT) calculation. We consider diamond as an example.
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! 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.
In TDDFT, we consider a system of electrons submitted to an external time-dependent (TD) perturbation, expressed usually in terms of an applied electric field $\mathbf{E}(t)$. This leads to a TD Kohn-Sham (KS) hamiltonian $\hat{H}(\mathbf{r},t)$:
\begin{equation} \tag{1} \hat{H}(\mathbf{r},t) = \frac{1}{2}\left(-\mathrm{i}\nabla + \frac{1}{c}\mathbf{A}(t)\right)^{\!\!2} + v^\phantom{I}_\textrm{KS}(\mathbf{r},t), \end{equation}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 TD-KS potential, a sum of the TD ionic, Hartree and exchange-correlation (XC) potentials. We assume in exciting
the adiabatic approximation for the TD-XC 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 ground-state is constructed from the ground-state density $n(\mathbf{r},0)$.
A KS wavefunction $\psi_{j\mathbf{k}}(\mathbf{r},t)$ labeled with index $j$ and associated to a wavevector k evolves as
\begin{equation} \tag{2} | \psi_{j\mathbf{k}}(t+\Delta t)\rangle = \hat{U}(t+\Delta t,t)\,| \psi_{j\mathbf{k}}(t)\rangle, \end{equation}where the propagator $\hat{U}(t+\Delta t,t)$ is given by:
\begin{equation} \tag{3} \hat{U}(t+\Delta t,t) = \hat{\mathcal{T}}\left[ \exp\left( -\mathrm{i}\int_t^{t+\Delta t} \mathrm{d}\tau\, \hat{H}(\tau) \right) \right], \end{equation}$\hat{\mathcal{T}}$ being the time-ordering operator. The following propagators are implemented:
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
\begin{equation} \tag{4} \mathbf{J}(t) = \frac{\mathrm{i}}{\Omega}\sum_{j\mathbf{k}}w_{\mathbf{k}}f_{j\mathbf{k}} \left\langle \psi_{j\mathbf{k}}(t) \big|\nabla \big|\psi_{j\mathbf{k}}(t)\right\rangle -\frac{N\mathbf{A}(t)}{c\,\Omega}, \end{equation}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 k-point, 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 k-point, the number of electrons that have been excited to an unoccupied KS state, labeled with $j$, is
\begin{equation} \tag{5} m^e_{j\mathbf{k}} (t)= \sum_{i} f_{i\mathbf{k}}\, | \langle \psi_{j\mathbf{k}}(0) | \psi_{i\mathbf{k}}(t)\rangle |^2, \end{equation}whereas the number of holes created in an occupied KS state $j′$ is
\begin{equation} \tag{6} m^h_{j'\mathbf{k}} (t)= f_{j'\mathbf{k}} - \sum_{i} f_{i\mathbf{k}}\, | \langle \psi_{j'\mathbf{k}}(0) | \psi_{i\mathbf{k}}(t)\rangle |^2. \end{equation}Thus, $N_\textrm{exc} (t)$ can be obtained by considering all the unoccupied states
\begin{equation} \tag{7} N_\textrm{exc}(t)= \sum_{j\mathbf{k}}^{j\, \textrm{unocc}} w_\mathbf{k}\, m^e_{j\mathbf{k}} (t) = \sum_{j'\mathbf{k}}^{j'\, \textrm{occ}} w_\mathbf{k}\, m^h_{j'\mathbf{k}} (t) . \end{equation}i) Preparation of the Input File
The first step is to create a directory for the system that we want to investigate, run_diamond_rt
%%bash
mkdir -p run_diamond_rt
As initial condition for the time evolution done in RT-TDDFT, we need the electron density and potential from a ground-state 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.d-8">
</groundstate>
</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_diamond_rt
python3 -m excitingscripts.setup.excitingroot
cd ..
%%bash
python3 -m excitingscripts.execute.single -r run_diamond_rt
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 RT-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.
i) Setting Elements and Attributes
You can prepare the input file (input.xml) of a RT-TDDFT calculation starting from the existing ground-state input file. The explicit ground-state calculation can be avoided with the attribute do="skip"
inside the element groundstate
.
...
<groundstate
do="skip"
... >
...
</groundstate>
...
In order to perform a "RT-TDDFT" calculation, the element xs
, must be added to the input file inside the input
element, as shown below.
...
<xs
xstype="RT-TDDFT"
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>
<pmat readFromFile="false" writeToFile="false" forceHermitian="false" />
</realTimeTDDFT>
</xs>
...
These settings in the input file correspond to a calculation with the following main features:
It is a RT-TDDFT calculation (the attribute xstype
is set to "RT-TDDFT"
).
The attribute ngridk
inside the xs
element determines the k-grid used for the KS wavefunctions $| \psi_{i\mathbf{k}}(t)\rangle$ evolved as given by Eq. (2). This k-grid 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 k-point reduction is to be used (Ref.[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 time-evolution in Eq. (2).
ngridk
, rgkmax
, vkloff
, nempty
, nosym
, and reducek
given in the xs
element are independent of their counterparts inside the groundstate
element.ngridk
, rgkmax
, vkloff
, nempty
, nosym
, and reducek
given inside xs
.ngridk
, …, reducek
inside the groundstate
element.The element realTimeTDDFT
determines the specific parameters for the RT-TDDFT calculation.
Here, we consider a propagator an expression which approximately enforces time-reversal symmetry (propagator="AETRS"
), with time-step 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 RT-TDDFT 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
\begin{equation} \tag{7} \mathbf{A}(t) = \mathbf{A}_0\, f(t) \cos ( \omega t + \phi ) \end{equation}
where $f(t)$ is the trapezoidal function:
\begin{equation} \tag{8} f(t) = \left\{ \begin{array}{ll} 0 & t \le t_0 \quad\mbox{ or }\quad t \ge t_0 + w + 2 t_r \\ 1 & t_0 + t_r \le t \le t_0 + t_r + w \\ (t-t_0)/t_r & t_0 < t < t_0 + t_r \\ (t_0 + w + 2 t_r - t )/t_r & t_0 + t_r + w < t < t_0 + w + 2 t_r \end{array} \right. \end{equation}
which is illustrated in the following figure.
amplitude="1.0d0"
) and is applied along the x-direction (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 "RT-TDDFT" calculation in Input Reference.
%%bash
python3 -m excitingscripts.execute.single -r run_diamond_rt
Be aware that, depending on the computer you are using, this calculation may take some minutes.
In our example, the one-step 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 self-consistent cycle triggered by the groundstate
element. The only exception is TIMING_RTTDDFT.OUT, which has the timings required to execute different subroutines invoked in the evolution of KS wavefunctions.
Other output files generated specifically by the RT-TDDFT calculations are:
To visualize the current density, type the following command inside the running directory:
%%bash
cd run_diamond_rt
python3 -m excitingscripts.plot.multitask --jind -f JIND.OUT --x --plot_name JIND
cd ..
This plots the x component of the current density [3]. The plot is saved in the running directory as JIND.png. Your figure should be similar to the following one:
By using the command
%%bash
cd run_diamond_rt
python3 -m excitingscripts.plot.multitask --nexc -f NEXC.OUT --plot_name NEXC
cd ..
you can plot the number of excited electrons per unit cell, $N_\textrm{exc}(t)$, as a function of time. The plot is saved in the running directory as NEXC.png. Bellow, we show how your plot is expected to be:
If you intend to obtain high-quality results, it is necessary to converge them with respect to a proper set of reliable parameters. In RT-TDDFT, typically they are:
ngridk
rgkmax
What regards the time-step, calculations are typically not too sensitive to it. Usually, there is a critical value $t_{cr}$ and when the employed time-step is larger than $t_{cr}$, a sort of divergence is observed. However, below tcr the results tend to be very similar and almost independent of the time-step (see discussion in Ref.[1]).
Below is, e.g., the current density for a finer k-grid of 16×16×16 points compared to 4×4×4 as used in this tutorial (see Ref.[5]):
Apart from a cosine function modulated by a trapezoid, as considered in this tutorial, it is also possible the following fields:
These other functions are explored in the tutorials Simulating pump-probe spectroscopy with RT-TDDFT and Studying higher-harmonic generation using RT-TDDFT.
Impulsive Function
In this case, the field is
\begin{equation} \tag{10} \mathbf{A}(t) = \left\{ \begin{array}{ll} 0 & t < t_0-w \\ -c\,\mathbf{E}_0 \frac{1}{16} \left( \frac{t+t_0}{w} +1 \right)^{\!\!3} \left[ 3\left( \frac{t-t_0}{w} \right)^{\!\!2} - 9\left( \frac{t-t_0}{w} \right) + 8 \right] \quad & t_0-w \le t \le t_0+w \\ -c\,\mathbf{E}_0 & t > t_0 + w \end{array} \right. \end{equation}The next figure depicts the vector potential as a function of time.
The electric field is
\begin{equation} \tag{11} \mathbf{E}(t) = -\frac{1}{c}\frac{d \mathbf{A} }{dt} = \left\{ \begin{array}{ll} \mathbf{E}_0 \frac{15}{16} \left( \frac{t-t_0}{w} +1 \right)^{\!\!2} \left( \frac{t-t_0}{w} -1 \right)^{\!\!2} & t_0 - w \le t_0 + w \\ 0 & t < t_0-w \mbox{ or } t > t_0+w \end{array} \right. \end{equation}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(t-t_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 t0. 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 time-step 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:
t0
;width
;amplitude
determines the amplitude $E_0$;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
\begin{equation} \tag{12} f(t) = \left\{ \begin{array}{ll} 0 , & t \le t_0 \quad\mbox{ or }\quad t \ge t_0 + w \\ \sin^2( \pi(t-t_0)/w)\quad & t_0 \le t \le t_0 + w \\ \end{array} \right. \end{equation}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:
t0
;omega
and phase
, respectively;pulseLength
;amplitude
determines the amplitude $A_0$, as given in Eq. (8);direction
defines the direction of the field.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.
The details on the implementation of real-time RT-TDDFT within the exciting code can be found in Ref.[1].
vkloff
attribute (see discussion in Ref.[1]). We also recommend to set nosym="true"
and reducek="false"
inside the element xs
.excitingscripts.plot.multitask
can be found here.where $\mathbf{A}_\textrm{ext}(t)$ and $\mathbf{A}_\textrm{ind}(t)$ are the external and induced vector potentials, respectively. On the one hand, the external is related to the displacement field $\mathbf{D}(t)$ as:
\begin{equation} \tag{14} \mathbf{A}_\textrm{ext}(t) = -c\int_0^t \mathbf{D}(t')\, dt' \end{equation}On the other hand, the induced vector potential, $\mathbf{A}_\textrm{ind}(t)$, is obtained from the current density as
\begin{equation} \tag{15} \frac{d^2\mathbf{A}_\textrm{ind}}{dt^2} = 4 \pi c\, \mathbf{J}(t). \end{equation}python3 -m excitingscripts.plot.multitask --jind -f JIND.OUT kpt16/JIND.OUT --x -c '4x4x4' '16x16x16'