1. Optical properties from Real-Time TD-DFT

In this tutorial we will use the real-time evolution of the electronic states to compute the optical properties of simple molecules. Different features of the TD-DFT implementation will be discussed, but the calculations here will assume that the atomic positions are fixed, and the time-evolution only affects the electronic wavefunctions.

[Additional information and context is provided in the video on TD-DFT.]

Note

This approach is different from the perturbation theory used for periodic crystals, which is covered in another tutorial.

The basic idea was proposed by Yabana and Bertsch [Phys. Rev. B 54, 4484 (1996)] and consists in following the evolution of the induced polarisation:

\[p(t)=\int d^3r\ r\cdot n(r,t)\]

when an external (time-dependent) electric field is applied to the system. The polarizability describes the response of the charge density to the field. The first-order contribution (linear polarizability) \(\alpha(\omega)\) is related to the dipole strength function, \(S(\omega)\), which is proportional to the photoabsorption cross section, allowing direct comparison with experiments.

We will take as example the water molecule. For the linear response calculation we will set the value of the field to 0.01 eV/Ang. The files in tddft-water contain the files required:

  • water.0.fdf

  • water.tddft.fdf

  • O.psf and H.psf

  • myfft.py

First we need to find the electronic solution in the presence of the external field. The input variables that we have to include in “water.0.fdf” are:

TDED.WF.Initialize      T

%block ExternalElectricField
0.000   0.010  0.000  V/Ang
%endblock ExternalElectricField

This will set the system under an electric field along the x-axis at t=0, find the ground-state using standard (time-independent) DFT, and finally save the wavefunctions into a file SystemLabel.TDWF.

Next, we switch the electric field off and propagate the electronic states using the Crank-Nicholson operator.

The basic flags for the input file (water.tddft.fdf) are:

MD.TypeOfRun         TDED
MD.FinalTimeStep     1
TDED.TimeStep       1.0E-03 fs
TDED.Nsteps          1000

In our example, the MD option is an “electronic dynamics only” (we have set MD.FinalTimeStep to 1, meaning that the atoms are kept fixed in their positions), which will use a time-step of 0.001fs, with a total time evolution of 3fs. We want to save the information of the dipole P(t), and the total energy E(t), using:

TDED.Write.Etot              T     # default value
TDED.Write.Dipole            T     # default value

After running the calculation, SIESTA will generate a SystemLabel.TDDIPOL with the 3 spatial components of the total dipole at each time step. The Fourier decomposition gives us the linear polarizability:

\[\alpha(\omega) \sim \int dt e^{i\omega t}p(t)\]

and the dipole strength function

\[S(\omega) \sim \omega \operatorname{Im}\alpha(\omega) $\]

The python script provided with the material (myfft.py) can be used to read the polarization vector from SystemLabel.TDDIPOL, FFT each component and print \(S(\omega)_{ij}\) (note that scaling factors are not taken into account in this simple version). As we have used an initial electric field pointing in the y-axis, our simulation will give us \(S(\omega)_{yj}\) with j=x,y,z. Frequencies are given in eV. You can plot the results with gnuplot or your favourite plotting software.

Our first simulation was relatively short to have enough resolution on the polarizability. If you wanted to continue with the simulation SIESTA had to be instructed to save the wavefunctions at the end of the simulation:

TDED.WF.Save     True

You can run the same calculation again adding this line. After the job is finished, copy the output files:

cp water.TDETOT water.TDETOT-lap1
cp water.TDDIPOL water.TDDIPOL-lap1

You can now resubmit the same job to have a total simulation run of 2 fs. Note however that the output files are overwritten, and the running time is also re-initialized. You will have to take this into account when analysis the results:

gnuplot> plot "water.TDETOT-lap1" w lp, "water.TDETOT" u ($1+1.0):2 w lp

The last command plots the 2 datasets from each simulation, shifting the time-frame for the second one to take into account that it corresponds to the second evolution.

Next, we check the effect of the time-step used during the simulation. Remember that you will have to run again the initialization state calculation in the presence of the external electric field (the SystemLabel.TDWF file was modified in the previous calculations). Once you have the initial state properly set, we can start the dynamic evolution. Take a longer 5 attosecond time step (5.0E-03 fs), and repeat the simulation with 1000 Nsteps. Plot the total energy as a function of time, and compare with the previous plot. You will see that the energy conservation is worse the larger the time-step. To improve the stability of the algorithm, you can use the hamiltonian extrapolation:

TDED.TimeStep               5.0E-03 fs
TDED.Extrapolate            True
TDED.Extrapolate.Substeps    3

You can try with different time-steps and extrapolation substeps to fine-tune the stability of the time-evolution.

Another important aspect to take into consideration is the quality of the basis set. You can play with the basis size, but also with the localization of the atomic orbitals. In general, it has been observed that confined basis (shorter cutoff radii) tend to shift the optical spectra to higher energies.

Finally remember that to get the optical susceptibility an average on the three spatial directions must be done. Place the different orientations of the electric field on x, y and z axis and compute:

\[\alpha(\omega) = Tr[\alpha_{ij}(\omega)]/3\]

where i and j correspond to x,y,z