by Ronaldo Rodrigues Pela for exciting fluorine
Purpose: In this tutorial, you will learn how to employ realtime timedependent densityfunctional theory (RTTDDFT) to study higherharmonic generation. We consider twodimensional MoS_{2} to illustrate.
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
Please refer to Ref.[2] and the tutorial Realtime TDDFT for more details about RTTDDFT in exciting.
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, MoS2hhg [1], and to enter it.
$ mkdir MoS2hhg
$ cd MoS2hhg
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>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.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 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 FLUORINE 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 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 the RTTDDFT calculations
i) First amplitude
Move to the previous directory, create a new one named ampl1, and enter it:
$ cd ../
$ mkdir ampl1
$ cd ampl1
You can prepare the input file (input.xml) of a RTTDDFT 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:
$ cp ../GS/{input.xml,EFERMI.OUT,STATE.OUT} .
Edit the file input.xml and avoid the explicit groundstate calculation with the attribute do = "skip" inside the element groundstate.
... <groundstate do="skip" ... > ... </groundstate> ...
Also, to perform the "RTTDDFT" calculation, add the element xs (as shown below) the input file inside the input element.
... <xs xstype="RTTDDFT" 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)where $f(t)$, the envelope function, is given as
(2)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 xdirection (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 "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 amplitude1 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 timeevolution 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:
$ PLOTmultitask.py ylabel "" 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:
$ PLOTmultitask.py ylabel "" preprocess f EFIELD.OUT scale 1.0 0.1 o EFIELDSCALED.OUT
The previous command stores the scaled electric field in the file EFIELDSCALED.OUT. Then, to plot $P(t)$ and the scaled electric field, $E(t)$/10:
$ PLOTmultitask.py f PVEC.OUT EFIELDSCALED.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 higherharmonic components are present in $P(t)$. To verify it, we can employ the FourierTransform can help, which can be obtained through:
$ PLOTmultitask.py ylabel "" 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]).
$ PLOTmultitask.py ylabel "" preprocess f P.OUT scale 4.0 27.52546105147261 handle_complex abs o PSCALED.OUT
To plot the absolute value of the Fourier Transform in a semilog graph:
$ PLOTmultitask.py f PSCALED.OUT xlabel '$\omega/\omega_0$' ylabel '$P(\omega)/E_m$' xlim 0 7 ylim 1e5 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 nonnegligible 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^{13} W/cm^{2} [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
$ PLOTmultitask.py f ampl1/PSCALED.OUT ampl2/PSCALED.OUT ampl3/PSCALED.OUT xlabel '$\omega/\omega_0$' ylabel '$P(\omega)/E_m$' xlim 0 7 ylim 1e5 1 semilog
Have in mind that you should follow the same procedure for the new amplitudes, obtaining for each of them PSCALED. We assume that the calculations have been performed in the subfolders ampl1, ampl2, and ampl3. In every of the subfolders, run the commands
$ PLOTmultitask.py ylabel "" preprocess fourier f PVEC.OUT o P.OUT wcut 0.025
$ PLOTmultitask.py ylabel "" preprocess f P.OUT scale 4.0 27.52546105147261 handle_complex abs o PSCALED.OUT
in folders ampl2 and ampl3 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, MoS2hhg.