Introduction

In this section, we will describe the ultafast spectroscopy calculations with PERTURBO. We will compute separately the electron and hole dynamics (independent particle picture) and will model an optical pump pulse:

\[\frac{\partial f_{n\mathbf{k}}(t)}{\partial t} = \mathcal{I}^{e-\mathrm{ph}}[f_{n\mathbf{k}}(t)] + F^{\mathrm{pump}}_{n\mathbf{k}}(t),\]

where \(f_{n\mathbf{k}}(t)\) is the time-dependent electron or hole occupation of a Bloch state with band index \(n\) and crystal momentum \(\mathbf{k},\) \(\mathcal{I}^{e-\mathrm{ph}}[f_{n\mathbf{k}}(t)]\) is the electron-phonon (e-ph) scattering integral computed from ab initio, and \(F^{\mathrm{pump}}_{n\mathbf{k}}(t)\) describes the pump pulse, accounting for the finite pulse duration and the energy broadening.

Having both electron \(f^e_{n\mathbf{k}}(t)\) and hole \(f^h_{n\mathbf{k}}(t)\) occupations as a function of time on the momentum grid, we will compute the transient absorption as a function of probe energy \(E\) and time \(t\):

\[\Delta A\big( E,t \big) \propto \sum_{nm\mathbf{k}} \big[ f^e_{n\mathbf{k}}(t) + f^h_{m\mathbf{k}}(t) \big] \times \delta(E - (\varepsilon^e_{n\mathbf{k}}-\varepsilon^h_{m\mathbf{k}})).\]

In this tutorial, we will consider the GaAs direct gap system. We will use the gaas_epr.h5 file from the example04-gaas folder. The epr file can be downloaded using this this link.

Overview of the workflow

The computational workflow will consist of the following stages:

  • run the bands calculation (calc_mode = 'bands') to determine the appropriate energy windows for electrons and holes
  • setup the ultrafast dynamics calculations for electrons and holes (calc_mode = 'setup'), generate the grid HDF5 file
  • run the electron and hole dynamics for 1 time step to compute the e-ph coupling elements and generate one time snapshot for carriers
  • generate the pump pulse HDF5 file with Perturbopy
  • run the ultafast dynamics simulation for 500 fs
  • postprocess simulation with Perturbopy: compute and plot transient absorption

Bands calculation

Here, we perform a typical bands calculation in GaAs. Details about this calculation mode of PERTURBO can be found here.

This calculation will be run in the bands folder:

mkdir bands

In the bands folder, we link the gaas_epr.h5 file; setup a gaas_band.kpt file:

6
  0.500 0.500 0.500  50 !L
  0.000 0.000 0.000  50 !G
  0.500 0.000 0.500  20 !X
  0.500 0.250 0.750  20 !W
  0.375 0.375 0.750  50 !K
  0.000 0.000 0.000  1  !G

and pert.in:

&perturbo
 prefix      = 'gaas'
 calc_mode   = 'bands'
 fklist   = 'gaas_band.kpt'
/

and run perturbo.x.

We then postprocess the bands calculation with a simple Perturbopy script (named plot_bands.py in the bands folder):

#!/usr/bin/env python3
import numpy as np
import perturbopy.postproc as ppy
import matplotlib.pyplot as plt
plt.rcParams.update(ppy.plot_tools.plotparams)

gaas_bands = ppy.Bands.from_yaml('gaas_bands.yml')

fig, ax  = plt.subplots(figsize=(10, 8))
labels_dict = {'L': [0.5, 0.5, 0.5], 'G': [0.0, 0.0, 0.0], 'X': [0.5, 0.0, 0.5], 'W': [0.5, 0.25, 0.75], 'K': [0.375, 0.375, 0.75], 'G': [0.0, 0.0, 0.0]}
gaas_bands.kpt.add_labels(labels_dict)
gaas_bands.plot_bands(ax, show_kpoint_labels=True)#, energy_window=[0,12])

emin_holes, emax_holes = 5.05, 5.3
emin_electrons, emax_electrons = 5.9, 6.7

xmin, xmax = ax.get_xlim()
ax.fill_between([xmin, xmax], [emin_holes, emin_holes],
                [emax_holes, emax_holes], color='blue', alpha=0.2)
ax.fill_between([xmin, xmax], [emin_electrons, emin_electrons],
                [emax_electrons, emax_electrons], color='red', alpha=0.2)
ax.set_xlim(xmin, xmax)

plt.show()

Here, we initialize the gaas_bands Perturbopy object from the YAML output file from PERTURBO, setup the k-point labels, and plot bands. The selected electron and hole energy windows are highlighted with red and blue shaded regions. More about the bands postprocessing can be found on the Perturbpy website.

gaas_bands

Energy windows for electron and hole bands depend on the pump energy of the optical excitations: the higher the energy, the larger the energy windows should be.

Setup the electron and hole grids

Here, we will set up the \(\mathbf{k}\)- and \(\mathbf{q}\)-grids for further electron and hole dynamics, using the calc_mode = 'setup', more about this calculation mode can be found here.

We will perform the separate electron and hole simulations in the cdyna-elec and cdyna-hole:

mkdir cdyna-elec
mkdir cdyna-hole

Link the epr file to both electron and hole directories:

ln -s <path to gaas_epr.h5>/gaas_epr.h5 cdyna-elec
ln -s <path to gaas_epr.h5>/gaas_epr.h5 cdyna-hole

First, for electrons (in the cdyna-elec folder), we create a pert_setup.in file (one can generate it using the interactive workflow).

&perturbo
! ***Mandatory parameters***
 calc_mode           = 'setup'
 prefix              = 'gaas'
 boltz_kdim(1)       = 120
 boltz_kdim(2)       = 120
 boltz_kdim(3)       = 120
 ftemper             = 'gaas.temper'

! ***Optional parameters***
 band_min           = 5
 band_max           = 5
 boltz_emin         = 5.9       ! eV
 boltz_emax         = 6.7       ! eV
 find_efermi        = .true.
 hole               = .false.
/

The top of the energy is chosen from the bands calculation mode. The bottom of the window is in the gap. We choose the first conduction band, with index 5. For the transient absorption calculations, the dense grid is particularly important. Here we use a quite dense grid of 120x120x120. However, for the converged calculations, even denser grids are required.

We then setup the gaas.temper file:

    1
  300.00          6.0456674805           0.1012649E+18

We then do similar calculations for holes (see the cdyna-hole folder for reference), where we select the corresponding energy window, select all the valence bands within the energy window (2-4), and set hole = .true. in pert_setup.in.

Dynamics initialization run

Now, we will run the electron and hole dynamics only for one step to generate the electron-phonon coupling elements (in the dyna-elec/tmp and dyna-hole/tmp folders), and to generate the gaas_cdyna.h5 files that will be later used by Perturbopy to generate the pump pulse excitation data.

For electrons, we create a cdyna-elec/pert_dyna_init.in input file:

&perturbo
  prefix      = 'gaas'

  calc_mode   = 'dynamics-run'

  yaml_fname = 'gaas_dyna_init.yml'

  ! grids
  boltz_kdim(1) = 120
  boltz_kdim(2) = 120
  boltz_kdim(3) = 120

  boltz_qdim(1) = 40
  boltz_qdim(2) = 40
  boltz_qdim(3) = 40

  ! band- and energy- windows
  boltz_emin         = 5.9         ! eV
  boltz_emax         = 6.7         ! eV
  band_min = 5
  band_max = 5

  ! carrier dynamics
  time_step = 1  ! femto-second
  boltz_nstep = 1
  output_nstep = 1
  solver = 'euler' ! default is Euler

  ! initial excitation (is irrelevant since it will be modified later)
  boltz_init_dist = 'gaussian'
  boltz_init_e0 = 6.4  ! eV
  boltz_init_smear = 40 ! meV
  boltz_init_ampl = 0

  boltz_de = 1

  ftemper  = 'gaas.temper'

  tmp_dir = "./tmp"
  load_scatter_eph = .false.

  !phonon energy cutoff
  phfreq_cutoff = 1 ! meV

  !smear parameter for delta function
  delta_smear = 8 ! meV
/

The input file is typical for calc_mode = 'dynamics-run' (see here). However, we pay special attention to the following parameters:

  • time_step: the time step is set to 1 fs and will be consistently used for all the dynamics calculations for electrons and holes
  • boltz_nstep: here, we run the dynamics only for one step
  • boltz_emin, boltz_emax, band_min, band_max: we use the same bands parameters as in the calc_mode = 'setup' for electrons
  • boltz_init_ampl: currently, we are not interested in the carrier initial state, therefore, we set the zero amplitude for excitation with boltz_init_ampl = 0
  • load_scatter_eph: we will compute the e-ph coupling elements. Therefore, we set this parameter to .false.
  • boltz_qdim: for this tutorial, we use an unconverged grid for phonons, 40x40x40. The \(\mathbf{q}\)-grid is commensurate with the \(\mathbf{k}\)-grid
  • delta_smear: the Gaussian smearing parameters should be cross-converged with the momenta grid convergence

For holes, see a similar file, pert_dyna_init.in in the cdyna-hole folder.

We then run perturbo.x for both electrons and holes (indicating -i pert_dyna_init.in). This is a relatively expensive calculation. Therefore, using HPC resources is recommended. The number of MPI tasks used for this calculation must be kept the same for the following dynamics runs. The calculation for holes is more expensive because of the three bands near the valence band maximum. Therefore, we are particularly cautious about the energy window for holes.

Generation of the pump pulse

Once we computed one step of the real-time dynamics, we have all the necessary components to generate the additional occupations due to a pump pulse, \(F^{\mathrm{pump}}_{n\mathbf{k}}(t)\). We use Perturbopy utils.spectra_generate_pulse.setup_pump_pulse function to generate the pulse, (see the file generate_pump_pulse.py in the example10-spectroscopy folder):

#!/usr/bin/env python3
import perturbopy.postproc as ppy

prefix = 'gaas'

# Define electron paths and DynaRun object
cdyna_elec_path = f'cdyna-elec/{prefix}_cdyna.h5'
elec_tet_path = f'cdyna-elec/{prefix}_tet.h5'
elec_yaml_path = f'cdyna-elec/{prefix}_dyna_init.yml'
elec_dyna_run = ppy.DynaRun.from_hdf5_yaml(cdyna_elec_path, elec_tet_path, elec_yaml_path)

# Define hole paths and DynaRun object
cdyna_hole_path = f'cdyna-hole/{prefix}_cdyna.h5'
hole_tet_path = f'cdyna-hole/{prefix}_tet.h5'
hole_yaml_path = f'cdyna-hole/{prefix}_dyna_init.yml'
hole_dyna_run = ppy.DynaRun.from_hdf5_yaml(cdyna_hole_path, hole_tet_path, hole_yaml_path)

# Print DynaRun info
print(elec_dyna_run)
print(hole_dyna_run)

# Generate the pump pulse HDF5 files in the
# elec_pump_pulse_path and hole_pump_pulse_path folders
pump_energy = 1.5
elec_pump_pulse_path = f'cdyna-elec/pump_pulse_elec_Epump_{pump_energy:05.2f}.h5'
hole_pump_pulse_path = f'cdyna-hole/pump_pulse_hole_Epump_{pump_energy:05.2f}.h5'

ppy.utils.spectra_generate_pulse.setup_pump_pulse(
    elec_pump_pulse_path,
    hole_pump_pulse_path,
    elec_dyna_run,
    hole_dyna_run,
    pump_energy=pump_energy,
    pump_time_step=1.0,
    pump_fwhm=20.0,
    pump_energy_broadening=0.040,
    pump_time_window=50.0,
    finite_width=True,
    animate=True)

elec_dyna_run.close_hdf5_files()
hole_dyna_run.close_hdf5_files()

In this script, we specify the paths to the prefix_cdyna.h5, prefix_tet.h5, prefix_dyna_init.yml files for both electrons and holes. The pulse files will be written to the cdyna-elec/pump_pulse_elec_Epump_01.50.h5 and cdyna-hole/pump_pulse_hole_Epump_01.50.h5 HDF5 files.

Here are the parameters for the setup_pump_pulse() Perturbopy function:

  • elec_pump_pulse_path: Path to the HDF5 file for the pump pulse excitation for electrons that will be created.
  • hole_pump_pulse_path: Path to the HDF5 file for the pump pulse excitation for holes
  • elec_dyna_run: The DynaRun object for electrons from the HDF5 and YAML files
  • hole_dyna_run: The DynaRun object for holes from the HDF5 and YAML files
  • pump_energy: Energy of the pump pulse in eV
  • pump_time_step: Time step in fs for the pump pulse generation. This must match the time step used in the PERTURBO dynamics run
  • pump_fwhm: Full width at half maximum (FWHM) of the pump pulse in fs
  • pump_energy_broadening: Energy broadening of the pump pulse in eV
  • pump_time_window: Time window for the pump pulse in fs. Determines the duration of the pump pulse activity
  • finite_width: If True, the pulse is finite in time. If False, the pulse is a step function
  • animate: If True, creates animation of the pump pulse.

After running this script, it generates the pump pulse HDF5 files and generates the pump pulse animation: pump_pulse

Ultrafast dynamics simulation with the pump pulse

Once we have the pump pulse HDF5 files generated, we are ready to run a long ultrafast dynamics. We run the dynamics separately for electrons and holes with the calc_mode = 'dynamics-run' PERTURBO calculation. We create pert_pump.in input files in both cdyna-elec and cdyna-hole folders.

For electrons, the input file is as follows:

&perturbo
  prefix      = 'gaas'

  calc_mode   = 'dynamics-run'

  boltz_kdim(1) = 120
  boltz_kdim(2) = 120
  boltz_kdim(3) = 120

  boltz_qdim(1) = 40
  boltz_qdim(2) = 40
  boltz_qdim(3) = 40

  ! band- and energy- windows
  boltz_emin         = 5.9         ! eV
  boltz_emax         = 6.7         ! eV
  band_min = 5
  band_max = 5

  ! carrier dynamics
  time_step = 1  ! femto-second
  boltz_nstep = 500
  output_nstep = 10
  solver = 'euler' ! default is Euler

  ! initial excitation (is irrelevant since it will be modified later)
  boltz_init_dist = 'gaussian'
  boltz_init_e0 = 6.4  ! eV
  boltz_init_smear = 40 ! meV
  boltz_init_ampl = 0

  ! pump pulse
  pump_pulse = .true.
  pump_pulse_fname = 'pump_pulse_elec_Epump_01.50.h5'

  boltz_de = 1

  ftemper  = 'gaas.temper'

  tmp_dir = "./tmp"
  load_scatter_eph = .true.

  !phonon energy cutoff
  phfreq_cutoff = 1 ! meV

  !smear parameter for delta function
  delta_smear = 8 ! meV
/

Here we outline the important parameters of the input file:

  • time_step: we consistently use 1 fs time step for the dynamics for electrons and holes
  • boltz_nstep: we will run dynamics for 500 steps
  • output_nstep: we output every 10th step of the simulation to reduce the prefix_cdyna.h5 file sizes. It is recommended to keep all of the snapshots for better temporal resolution, if storage allows.
  • boltz_init_ampl: we keep the zero-amplitude Gaussian excitation, as the actual excitation will be read from the pump pulse
  • pump_pulse: we set this parameter to .true. for the pump pulse excitation
  • pump_pulse_fname: we specify the name for the HDF5 pump pulse file previously generated with Perturbopy
  • load_scatter_eph: we load the precomputed el-ph scattering elements

For holes, see the hole input file in the cdyna-hole folder. After running perturbo.x -i pert_pump.in (using the MPI parallelization if it was used for the dynamics initialization run with pert_init_dyna.in), we obtain the YAML (gaas_dyanmics-run.yml) and HDF5 (gaas_cdyna.h5) files that will be used for the postprocessing.

Computing and ploting the transient absorption spectra

Computing transient absorption

At this final step, we will use Perturbopy to compute and plot transient absorption maps.

We compute the transient absorption with the .spectra_trans_abs.compute_trans_abs() function from Perturbopy (see the file compute_trans_abs.py in the example10 folder):

#!/usr/bin/env python3
import perturbopy.postproc as ppy

prefix = 'gaas'

# Define electron paths and DynaRun object
cdyna_elec_path = f'cdyna-elec/{prefix}_cdyna.h5'
elec_tet_path = f'cdyna-elec/{prefix}_tet.h5'
elec_yaml_path = f'cdyna-elec/{prefix}_dynamics-run.yml'
elec_dyna_run = ppy.DynaRun.from_hdf5_yaml(cdyna_elec_path, elec_tet_path, elec_yaml_path)

# Define hole paths and DynaRun object
cdyna_hole_path = f'cdyna-hole/{prefix}_cdyna.h5'
hole_tet_path = f'cdyna-hole/{prefix}_tet.h5'
hole_yaml_path = f'cdyna-hole/{prefix}_dynamics-run.yml'
hole_dyna_run = ppy.DynaRun.from_hdf5_yaml(cdyna_hole_path, hole_tet_path, hole_yaml_path)

# Print DynaRun info
print(elec_dyna_run)
print(hole_dyna_run)

# Print the PumpPulse info
print(elec_dyna_run.pump_pulse)

ppy.utils.spectra_trans_abs.compute_trans_abs(elec_dyna_run, hole_dyna_run)

elec_dyna_run.close_hdf5_files()
hole_dyna_run.close_hdf5_files()

The function compute_trans_abs() requires as input only the electron and hole DynaRun objects, defined in the script as elec_dyna_run and hole_dyna_run. Once the run is complete, the script generates several numpy binary files:

  • trans_abs_dA_Epump_1.500.npy: total transient absorption, 2D array with energy grid index as the first axis and time grid index as the second axis
  • trans_abs_dA_elec_Epump_1.500.npy: electron contribution to the transient absorption spectrum
  • trans_abs_dA_hole_Epump_1.500.npy: hole contribution to the transient absorption spectrum
  • trans_abs_T_Epump_1.500.npy: time grid array
  • trans_abs_E_Epump_1.500.npy: energy grid array

The suffix 1.500 in the filenames is the pump energy in eV.

Plotting transient absorption

Having the above-mentioned .npy files, we are ready for the final step of the tutorial: plotting the transient absorption. It can be done with a few lines of Python code:

#!/usr/bin/env python3
import numpy as np
import perturbopy.postproc as ppy
import matplotlib.pyplot as plt

pump_energy = 1.5

# Load the numpy binary files
dA_elec = 'trans_abs_dA_elec_Epump_1.500.npy'
dA_hole = 'trans_abs_dA_hole_Epump_1.500.npy'
time_grid = 'trans_abs_T_Epump_1.500.npy'
energy_grid = 'trans_abs_E_Epump_1.500.npy'

dA_elec = np.load('trans_abs_dA_elec_Epump_1.500.npy')
dA_hole = np.load('trans_abs_dA_hole_Epump_1.500.npy')
time_grid = np.load('trans_abs_T_Epump_1.500.npy')
energy_grid = np.load('trans_abs_E_Epump_1.500.npy')

# Plot the total transient absorption
fig, ax = plt.subplots(1, 1, figsize=(8, 6))
plt.subplots_adjust(left=0.12, right=0.88, top=0.9, bottom=0.12)
ppy.utils.spectra_plots.plot_trans_abs_map(ax,  time_grid, energy_grid, dA_elec+dA_hole, vmin=-2.5, vmax=0.0)
ax.set_title(f'Transient absorption in GaAs, Epump={pump_energy:.1f} eV')
ax.set_ylim(ymax=2.3)
fig.savefig('dA_total.png')
plt.show()

The script produces the following plot: dA

One can run the plot_trans_abs_map() function with the dA_elec or dA_hole arrays to plot the electron and hole contributions to the spectrum. See the file plot_trans_abs.py in the example.