Studying Higher-Harmonic Generation using RT-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 employ real-time time-dependent density-functional theory (RT-TDDFT) to study higher-harmonic generation. We consider two-dimensional MoS2 to illustrate it.
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.
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_MoS2_hhg and a subdirectory for the ground-state calculation, GS.
%%bash
mkdir -p run_MoS2_hhg && cd run_MoS2_hhg
mkdir -p GS
cd ..
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>
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_MoS2_hhg/GS
python3 -m excitingscripts.setup.excitingroot
cd ../..
%%bash
cd run_MoS2_hhg
python3 -m excitingscripts.execute.single -r GS
cd ..
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 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 the RT-TDDFT Calculations¶
i) First Amplitude
Create a new subdirectory ampl_1 and enter it. You can prepare the input file (input.xml) of a RT-TDDFT calculation starting from the existing groundstate input file. Copy the following files to the current directory. You can perform all of this with the following commands:
%%bash
cd run_MoS2_hhg
mkdir -p ampl_1 && cd ampl_1
cp ../GS/{input.xml,EFERMI.OUT,STATE.OUT} .
cd ../..
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="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>
<pmat readFromFile="false" writeToFile="false" forceHermitian="false" />
</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
where $f(t)$, the envelope function, is given as
\begin{equation} \tag{2} 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.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.
%%bash
cd run_MoS2_hhg
python3 -m excitingscripts.execute.single -r ampl_1
cd ..
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:
%%bash
cd run_MoS2_hhg/ampl_1
python3 -m excitingscripts.plot.multitask --preprocess --get_efield -f AVEC.OUT --x -o EFIELD.OUT
cd ../..
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:
%%bash
cd run_MoS2_hhg/ampl_1
python3 -m excitingscripts.plot.multitask --preprocess -f EFIELD.OUT --scale 1.0 0.1 -o EFIELD-SCALED.OUT
cd ../..
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:
%%bash
cd run_MoS2_hhg/ampl_1
python3 -m excitingscripts.plot.multitask -f PVEC.OUT EFIELD-SCALED.OUT -c '$P(t)$' '$E(t)/10$' --ylabel 'Field [a.u.]' --plot_name 'PLOT_01'
cd ../..
The resulting figure should look like
From the figure, although we see that there is a phase delay of almost 180° between $P(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:
%%bash
cd run_MoS2_hhg/ampl_1
python3 -m excitingscripts.plot.multitask --preprocess --fourier -f PVEC.OUT -o P.OUT --wcut 0.025
cd ../..
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]).
%%bash
cd run_MoS2_hhg/ampl_1
python3 -m excitingscripts.plot.multitask --preprocess -f P.OUT --scale 4.0 27.52546105147261 --handle_complex abs -o P-SCALED.OUT
cd ../..
To plot the absolute value of the Fourier Transform in a semilog graph:
%%bash
cd run_MoS2_hhg/ampl_1
python3 -m excitingscripts.plot.multitask -f P-SCALED.OUT --xlabel '$\omega/\omega_0$' --ylabel '$P(\omega)/E_m$' --xlim 0 7 --ylim 1e-5 1 --semilog --plot_name 'PLOT_02'
cd ../..
You should obtain a figure similar to the one below.
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.
The amplitude considered in the previous calculation is rather high. It corresponds to an intensity of 4.6×1013 W/cm2 [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
python3 -m excitingscripts.plot.multitask -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
python3 -m excitingscripts.plot.multitask --preprocess --fourier -f PVEC.OUT -o P.OUT --wcut 0.025
python3 -m excitingscripts.plot.multitask --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, run_MoS2_hhg.
Bibliography¶
- Here hhg stands for higher-harmonic generation.
- Ronaldo Rodrigues Pela, Claudia Draxl. "All-electron full-potential implementation of real-time TDDFT in exciting" (link).
- More information about the script
excitingscripts.plot.multitaskcan be found here. - 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.
- The intensity of an electromagnetic wave is
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/cm2 with $(\varepsilon_0\, c)/2$ = 3.50941×1016. 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.