Real-Time TDDFT

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.

Table of Contents

0. Before Starting

1. Theoretical Background for RT-TDDFT Calculation

2. Ground-State Calculation

3. Performing a RT-TDDFT Calculation

4. Visualizing the Output

5. Converging the Results

6. Other Expressions for the Laser Field

Reference

0. Before Starting

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.

1. Theoretical Background for RT-TDDFT Calculation

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:

  • simple exponential,
  • exponential at midpoint rule,
  • approximate enforced time-reversal symmetry,
  • Commutator-Free Magnus expansion of 4th order,
  • exponential using a basis of the hamiltonian-eigenvectors, and
  • Runge-Kutta 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

\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}

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, run_diamond_rt

In [1]:
%%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:

In [3]:
%%bash
cd run_diamond_rt
python3 -m excitingscripts.setup.excitingroot
cd ..

ii) Running the Ground-State Calculation

You can start the calculation by invoking the script excitingscripts.execute.single.

In [4]:
%%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.

3. Performing a RT-TDDFT Calculation

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

    • 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 one-step ground-state 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 one-step calculation, had been obtained beforehand using the parameters 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.

  • In our example, the field has an amplitude of 1.0 a.u. (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.

ii) Running exciting

Having modified the input file as described above, you can launch the RT-TDDFT calculation using the following command.

In [ ]:
%%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:

  • 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 running directory:

In [12]:
%%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

In [14]:
%%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:

5. Converging the results

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
  • local-orbitals

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

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 (delta-kick)
  • 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 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:

  • $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

\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:

  • $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.

Reference

The details on the implementation of real-time RT-TDDFT within the exciting code can be found in Ref.[1].

Bibliography

  1. Ronaldo Rodrigues Pela, Claudia Draxl. "All-electron full-potential implementation of real-time TDDFT in exciting" (link).
  2. In general, it is a good practice break the k-mesh off symmetry by modifying the vkloff attribute (see discussion in Ref.[1]). We also recommend to set nosym="true" and reducek="false" inside the element xs.
  3. More information about the script excitingscripts.plot.multitask can be found here.
  4. The (total) vector potential $\mathbf{A}(t)$, appearing in the Hamiltonian Eq. (1), is a sum
\begin{equation} \tag{13} \mathbf{A}(t) = \mathbf{A}_\textrm{ext}(t)+\mathbf{A}_\textrm{ind}(t), \end{equation}

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}
  1. Assuming that you performed a new calculation for a k-grid of 16×16×16 in a subfolder kpt16 inside the folder run_diamond_rt, and that your current directory is run_diamond_rt, then the following command can be used to plot two current densities
python3 -m excitingscripts.plot.multitask --jind -f JIND.OUT kpt16/JIND.OUT --x -c '4x4x4' '16x16x16'