Studying Higher-Harmonic Generation using RT-TDDFT

by Ronaldo Rodrigues Pela for exciting neon

Purpose: In this tutorial, you will learn how to employ real-time time-dependent density-functional theory (RT-TDDFT) to study higher-harmonic generation. We consider two-dimensional MoS₂ to illustrate.

---

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-multitask.py: Python script for making plots from output files of RT-TDDFT.

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 RT-TDDFT calculation

Please refer to Ref.[2] and the tutorial Real-time TDDFT for more details about RT-TDDFT in exciting.

---

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, MoS2-hhg [1], and to enter it.

$ mkdir MoS2-hhg
$ cd MoS2-hhg

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>MoS2</title>
 
   <structure speciespath="$EXCITINGROOT/species/">
      <crystal scale="6.0274">
         <basevect>  0.5  0.86602540378   0.0 </basevect>
         <basevect> -0.5  0.86602540378   0.0 </basevect>
         <basevect>  0.0  0.00000000000  10.0 </basevect>
      </crystal>
      <species speciesfile="Mo.xml" rmt="2.3">
         <atom coord="0.666667  0.666667  0.0"/>
      </species>
      <species speciesfile="S.xml"  rmt="2.05">
         <atom coord="0.333333  0.333333 -0.04938"/>
         <atom coord="0.333333  0.333333  0.04938"/>
      </species>
   </structure>
 
   <groundstate
      do="fromscratch"
      xctype="LDA_PW"
      ngridk="12 12 1"
      rgkmax="5.0d0"
      nempty="5"
      epsengy="1.d-8">
   </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:

$ SETUP-excitingroot.sh

ii) Running the ground-state calculation

You can start the calculation by invoking the script EXECUTE-single.sh.

$ EXECUTE-single.sh GS

After the calculation is finished, move to the subdirectory GS.

$ cd GS

You can check the output files, especially the main output file INFO.OUT, for general information. In the case of a successfully finished calculation, the last lines of the INFO.OUT will contain the message

...
================================================================================
| EXCITING NEON stopped                                                        =
================================================================================

Check if this is indeed the case and if the EFERMI.OUT and STATE.OUT files are present. They 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 ground-state density and potential, it is necessary to perform 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. Performing the RT-TDDFT calculations

i) First amplitude

Move to the previous directory, create a new one named ampl-1, and enter it:

$ cd ../
$ mkdir ampl-1
$ cd ampl-1

You can prepare the input file (input.xml) of a RT-TDDFT calculation starting from the existing ground-state input file. Copy the following files to the current directory. You can perform all of this with the following commands:

$ cp ../GS/{input.xml,EFERMI.OUT,STATE.OUT} .

Edit the file input.xml and avoid the explicit ground-state calculation with the attribute do = "skip" inside the element groundstate.

...
   <groundstate 
      do="skip"
      ... >
      ...
   </groundstate>
...

Also, to perform the "RT-TDDFT" calculation, add the element xs (as shown below) the input file inside the input element.

...
   <xs
      xstype="RT-TDDFT"
      ngridk="6 6 1"
      vkloff="0.01 0.02 0.004"
      nempty="5"
      nosym="true"
      reducek="false">
 
      <realTimeTDDFT
         propagator="AETRS"
         timeStep="0.5d0"
         endTime="200.d0">
         <laser fieldType="total">
            <sinSq amplitude="20.d0" omega="0.25d0" phase="0.d0"
               t0="0.1d0" pulseLength="160.d0" direction="x" />
         </laser>
      </realTimeTDDFT>
   </xs>
...

The element laser describes the external field applied to the system. Here we consider a cosine modulated by a sin squared, as given in the element sinSq. The expression for the corresponding vector potential is

(1)

\begin{align} \mathbf{A}(t) = \mathbf{A}_0\, f(t) \cos ( \omega t + \phi ) \end{align}

where $f(t)$, the envelope function, is given as

(2)

\begin{align} 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{align}

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. In our example, we consider a vector field along the x-direction (direction="x") with an amplitude of 20 a.u. (amplitude="20.d0"). The cosine function has an angular frequency of $\omega$ = 0.25 a.u. (omega="0.25d0") and phase $\phi$ = 0.0 (phase="0.d0"). The envelope function has parameters: initial time $t_0$= 0.1 a.u. and width of $w$=160 a.u. (as given by: t0="0.1d0", pulseLength="160.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 inside the subdirectory amplitude-1 using

$ export OMP_STACKSIZE=2G
$ time exciting_smp

Be aware that, depending on the computer you are using, this calculation may take several minutes.

iii) Analyzing the polarization

The time-evolution of the polarization vector $P(t)$ is printed out in the file PVEC.OUT. With this file and AVEC.OUT, we can inspect how the polarization of the system reacts to the applied electric field. To do so, we first obtain the electric field from AVEC.OUT:

$ PLOT-multitask.py --preprocess --get_efield -f AVEC.OUT --x -o EFIELD.OUT

This command generates the output file EFIELD.OUT with the x component of the electric field, obtained from $A_x(t)$ in AVEC.OUT [3].

To better compare the $P(t)$ and $E(t)$, it is advisable to scale the electric field by a factor of 1/10. This can be done with:

$ PLOT-multitask.py --preprocess -f EFIELD.OUT --scale 1.0 0.1 -o EFIELD-SCALED.OUT

The previous command stores the scaled electric field in the file EFIELD-SCALED.OUT. Then, to plot $P(t)$ and the scaled electric field, $E(t)$/10:

$ PLOT-multitask.py -f PVEC.OUT EFIELD-SCALED.OUT -c '$P(t)$' '$E(t)/10$' --ylabel 'Field [a.u.]'

The resulting figure should look like

From the figure, although we see that there is a phase delay of almost 180° between $E(t)$ and $E(t)$, it is still not clear whether higher-harmonic components are present in $P(t)$. To verify it, we can employ the Fourier-Transform can help, which can be obtained through:

$ PLOT-multitask.py --preprocess --fourier -f PVEC.OUT -o P.OUT --wcut 0.025

In the previous command, we selected a broadening of $\eta$=0.025 a.u. This includes an exponential $\exp(-\eta t)$ in the definition of the Fourier Transform, in order to smooth it.

Before we plot the Fourier transform of $P(t)$, it is good to scale the angular frequencies $\omega$ with respect to 0.25 a.u. (the frequency of the applied laser), and also $P(t)$ with respect to the maximum of the electric field, $E_m$=0.03633 a.u. (see Note [4]).

$ PLOT-multitask.py --preprocess -f P.OUT --scale 4.0 27.52546105147261 --handle_complex abs -o P-SCALED.OUT

To plot the absolute value of the Fourier Transform in a semilog graph:

$ PLOT-multitask.py -f P-SCALED.OUT --xlabel '$\omega/\omega_0$' '$P(\omega)/E_m$' --xlim 0 7 --ylim 1e-5 1 --semilog

You should obtain a figure similar to the one bellow.

In this figure, $\omega$ is the usual variable employed for the frequency domain; and $\omega_0$ stands for the angular frequency of the electric field (here $\omega_0$=0.25 a.u.). Now, you can see a non-negligible component at $\omega/\omega_0$ = 3.

iv) Other amplitudes

The amplitude considered in the previous calculation is rather high. It corresponds to an intensity of 4.6×10¹³ W/cm² [5]. As an exercise, you can run calculations for other amplitudes. You see in the next Figure that when the amplitude is reduced by a factor of 10 or 100, the third harmonic component almost practically disappears from the spectrum.

You can obtain this figure with the command

$ PLOT-multitask.py -f ampl-1/P-SCALED.OUT ampl-2/P-SCALED.OUT ampl-3/P-SCALED.OUT --xlabel '$\omega/\omega_0$' --ylabel '$P(\omega)/E_m$' --xlim 0 7 --ylim 1e-5 1 --semilog

Have in mind that you should follow the same procedure for the new amplitudes, obtaining for each of them P-SCALED. We assume that the calculations have been performed in the subfolders ampl-1, ampl-2, and ampl-3. In every of the subfolders, run the commands

$ PLOT-multitask.py --preprocess --fourier -f PVEC.OUT -o P.OUT --wcut 0.025
$ PLOT-multitask.py --preprocess -f P.OUT --scale 4.0 27.52546105147261 --handle_complex abs -o P-SCALED.OUT

in folders ampl-2 and ampl-3 change the argument 27.52546105147261 in the second command above to 275.2546105147261 and 2752.546105147261, respectively, to have overlapping curves in the final plot. Then, you can use the previous command within the parent directory, MoS2-hhg.

---

Bibliography

1. Here hhg stands for higher-harmonic generation.

2. Ronaldo Rodrigues Pela, Claudia Draxl. "All-electron full-potential implementation of real-time TDDFT in exciting" (link).

3. More information about the script PLOT-multitask.py can be found here.

4. Note that dividing by 0.25 and 0.03633 is equivalent to multiplying by 4.0 and 27.52546105147261, respectively, the given scaling factors within the command.

5. The intensity of an electromagnetic wave is

(3)

\begin{align} I = \frac{1}{2}\varepsilon_0\, c\, E_m^2 \end{align}

where $\varepsilon_0$ is the dielectric constant, $c$ is the speed of light, and $E_m$ is the peak-value of the electric field. Given $E_m$ in atomic units, $I$ is found in W/cm² with $(\varepsilon_0\, c)/2$ = 3.50941×10¹⁶. Please note that, due to reflections on the interface, $E_m$ inside the considered material is not the same as that produced by a laser. This also means that $I$ calculated using $E_m$ inside the material is not expected to match the intensity of the light generated by the laser, rather the intensity which penetrates the material.