Simulating Pump-Probe Spectroscopy with 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 simulate a simple pump-probe experiment. We consider two-dimensional MoS2 to illustrate.
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¶
In a pump-probe experiment, there are basically two laser pulses that shine a sample with different goals. The first pulse (pump) excites the sample during a specific time, generating a non-equilibrium state. The second pulse, usually much weaker than the first, probes specific properties, such as optical constants. With this technique, it is possible to analyze the impact of excitations on the optical properties, and how the systems relax back to the equilibrium.
In exciting, it is possible to simulate this kind of experiment, employing RT-TDDFT. In one possible scenario, we first consider a pump field described by a cosine modulated by a sin squared. This field mimics a gaussian envelope, which is quite common in experiments. The expression for the vector potential $\mathbf{A}_\textrm{pump}(t)$ 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
After some time $t$ when the pump is not acting any longer, we can consider a delta-kick as the probe
\begin{equation} \tag{3} \mathbf{A}_\textrm{probe}(t) = -c\,\mathbf{E}'_0\, u( t-t_\textrm{probe} ), \end{equation}where $u( t-t_\textrm{probe} )$ is a step-function applied at time $t=t_\textrm{probe}$. The electric field of the probe is
\begin{equation} \tag{4} \mathbf{E}_\textrm{probe}(t)=-\frac{1}{c}\frac{d}{dt}\mathbf{A}_\textrm{probe}(t) = \mathbf{E}'_0\, \delta( t-t_\textrm{probe} )\:. \end{equation}In this tutorial, we will analyze how the imaginary part of the dielectric function (related to the absorption spectrum) of the sample excited by the pump differs from its ordinary counterpart (evaluated starting with the sample on the ground state). With this purpose, we consider the general relationship between the Fourier-transforms of the current density $J$ and the electric field $E$:
\begin{equation} \tag{5} \sigma_{\alpha\beta}(\omega) = \frac{J_\alpha(\omega)}{E_\beta(\omega)}, \quad \quad\varepsilon_{\alpha\beta}(\omega) = \delta_ {\alpha\beta}+ \frac{4\pi \mathrm{i}\,\sigma_{\alpha\beta}(\omega)}{\omega}, \end{equation}in which $\sigma$ and $\varepsilon$ are the optical conductivity and the dielectric tensors, respectively; the indexes $\alpha$, $\beta$ denote the cartesian directions x, y, or z [1].
To simulate the pump-probe experiment, we can evaluate the dielectric function as indicated in Eq. (5), but taking the current density as the difference between the values after the pump and the probe ($J_\textrm{pump-probe}$) and the pump pulse $J_\textrm{pump}$, as obtained from two separate calculations. And as comparison, we consider a third calculation just with the probe pulse, from which we can extract the dielectric function from the system starting on the groundstate.
Please refer to Ref [2] for one different example of simulation of pump-probe spectroscopy with 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, run_MoS2_pump_probe and a subdirectory for the ground-state calculation, GS.
%%bash
mkdir -p run_MoS2_pump_probe && cd run_MoS2_pump_probe
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_pump_probe/GS
python3 -m excitingscripts.setup.excitingroot
cd ../..
%%bash
cd run_MoS2_pump_probe
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 Pump-Probe Calculations¶
i) First Case: Pump-Probe
Create a new subdirectory pump_probe. 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_pump_probe
mkdir -p pump_probe && cd pump_probe
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="500.d0"
vectorPotentialSolver="improvedeuler">
<laser fieldType="total">
<sinSq amplitude="20.d0" omega="0.5d0" phase="0.d0"
t0="0.1d0" pulseLength="80.d0" direction="x" />
<kick amplitude="0.1d0" t0="100.d0" width="0.d0" direction="x" />
</laser>
<pmat readFromFile="false" writeToFile="false" forceHermitian="false" />
</realTimeTDDFT>
</xs>
...
Please check all relevant parameters for a "RT-TDDFT" calculation in Input Reference.
Having modified the input file as described above, you can launch the RT-TDDFT calculation inside the subdirectory pump_probe using
%%bash
cd run_MoS2_pump_probe
python3 -m excitingscripts.execute.single -r pump_probe
cd ..
Be aware that, depending on the computer you are using, this calculation may take several minutes up to some hours.
ii) Second Case: Pump
Now, we can study the effect of the pump alone on MoS2. Create a new subdirectory named pump, and enter it. Copy the files input.xml, EFERMI.OUT and STATE.OUT from the ground-state calculation to the current directory:
%%bash
cd run_MoS2_pump_probe
mkdir -p pump && cd pump
cp ../GS/{input.xml,EFERMI.OUT,STATE.OUT} .
cd ../..
As done before, edit the file input.xml to avoid the explicit ground-state calculation by setting do="skip" inside the element groundstate.
Now you need to add the element xs 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="500.d0"
vectorPotentialSolver="improvedeuler">
<laser fieldType="total">
<sinSq amplitude="20.d0" omega="0.5d0" phase="0.d0"
t0="0.1d0" pulseLength="80.d0" direction="x" />
</laser>
<pmat readFromFile="false" writeToFile="false" forceHermitian="false" />
</realTimeTDDFT>
</xs>
...
Now you are ready to execute the RT-TDDFT calculation inside the subdirectory pump using
%%bash
cd run_MoS2_pump_probe
python3 -m excitingscripts.execute.single -r pump
cd ..
Have in mind that depending on your computational resources, this calculation may take some minutes or even hours.
iii) Third Case: Probe
Finally, we can carry out a calculation to obtain the dielectric function of MoS2 starting from the groundstate. Similar as you did before, you need to go to the previous directory, create a new one named probe, and enter it. You can again copy the files input.xml, EFERMI.OUT and STATE.OUT from the ground-state calculation to the current directory:
%%bash
cd run_MoS2_pump_probe
mkdir -p probe && cd probe
cp ../GS/{input.xml,EFERMI.OUT,STATE.OUT} .
cd ../..
As done before, edit the file input.xml to skip the ground-state calculation and add the element xs 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="500.d0"
vectorPotentialSolver="improvedeuler">
<laser fieldType="total">
<kick amplitude="0.1d0" t0="100.d0" width="0.d0" direction="x" />
</laser>
<pmat readFromFile="false" writeToFile="false" forceHermitian="false" />
</realTimeTDDFT>
</xs>
...
Then, perform the RT-TDDFT calculation (inside the subdirectory probe) using
%%bash
cd run_MoS2_pump_probe
python3 -m excitingscripts.execute.single -r probe
cd ..
As warned before, this calculation may take more some minutes or hours, depending on the computational power that you have.
iv) Analyzing the Dielectric Function
First, it is convenient to obtain the difference $J_\textrm{pump-probe}-J_\textrm{pump}$. We first move to the parent directory and then preprocess the files JIND.OUT in the folders pump_probe and pump with the command:
%%bash
cd run_MoS2_pump_probe
python3 -m excitingscripts.plot.multitask --preprocess --sub -f pump_probe/JIND.OUT pump/JIND.OUT -o DIFF.OUT
cd ..
With the previous command, the difference $J_\textrm{pump-probe}-J_\textrm{pump}$ is stored in the file DIFF.OUT. It assumes that we take the x components, otherwise, we would need to add --y or --z.
Now, you can plot the current densities $J_\textrm{pump}$, $J_\textrm{probe}$, $J_\textrm{pump-probe}$ and the difference $J_\textrm{pump-probe}-J_\textrm{pump}$ with [3]:
%%bash
cd run_MoS2_pump_probe
python3 -m excitingscripts.plot.multitask --jind -f pump/JIND.OUT probe/JIND.OUT pump_probe/JIND.OUT DIFF.OUT --legend 'lower right' --plot_name 'JIND'
cd ..
Your figure should be similar to the one given bellow:
Now, we can compare the imaginary part dielectric function in the pump-probe case, and compare it with the probe case. First, we determine the dielectric function for the probe case:
%%bash
cd run_MoS2_pump_probe/probe
python3 -m excitingscripts.plot.multitask --preprocess --get_eps --xx -f AVEC.OUT -f JIND.OUT --wcut 0.05 -o EPSILON.OUT
cd ../..
Now, we move back to parent directory, and obtain the dielectric function for the pump-probe:
%%bash
cd run_MoS2_pump_probe
python3 -m excitingscripts.plot.multitask --preprocess --get_eps --xx -f probe/AVEC.OUT -f DIFF.OUT --wcut 0.05 -o EPSILON-DIFF.OUT
cd ..
And, finally, we can plot the imaginary part of both dielectric functions in the same figure using:
%%bash
cd run_MoS2_pump_probe
python3 -m excitingscripts.plot.multitask --imag_eps -f probe/EPSILON.OUT EPSILON-DIFF.OUT --xlim 5 30 --ylim 0 1 --plot_name 'EPSILON'
cd ..
Your figure should be similar to the one bellow:
Be aware that the calculations shown here employ parameters which are not converged, and should be considered only for introductory purposes. To improve the quality of the calculation, special attention needs to be paid to the following parameters:
ngridkrgkmaxtimeStepandendTime
Bibliography¶
- Note that to obtain $\sigma_{\beta\alpha}$ it is convenient to take electric field as a delta-kick $E_0\,\delta(t)$ in a specific direction $\alpha$ such that $E_\alpha(\omega) = E_0$.
- 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