Among the most useful class of properties to be able to predict for a molecule or crystal are its optical properties. These can be used to find the frequencies of optical radiation that will be absorbed or emitted, based on its electronic structure. This can be used to find, for example, what colour a dye molecule will have in a solvent.
There are several ways these properties could be calculated starting from density functional theory. The method we’ll look at in this lab is time-dependent density functional theory (TDDFT). A TDDFT code could be used to calculate the evolution of the electrons under the effect of the oscillating electric field associated with the presence of a photon. This would be a fairly intensive calculation, and would need to be repeated at many different energies to build up the full optical spectrum.
The set of codes that comes with the quantum espresso package, turboTDDFT, actually use time-dependent density functional perturbation theory (TDDFPT) to calculate optical spectra of molecules. This calculates the response of the system to the oscillating electric field associated with an incoming photon using perturbation theory (polarizability), in a similar way to how DFPT finds the response of the system to perturbations of the atomic positions. It then uses this to find the susceptibility as a function of photon energy. The code also uses an approach known as the Liouville-Lanczos method which allows it to calculate the full optical spectrum for a wide energy range, at a very moderate cost.
Typically if you want to understand what energy of photon will be absorbed,
you might expect you would need to know the difference between the energy of
the occupied and empty states: when a photon is absorbed it excites an
electron from an occupied state to an empty state higher in energy by an
amount equal to the photon energy. The turboTDDFT codes use a clever approach
to avoid the need for including a large number of additional empty states in
the calculation. This involves instead using the projection operator on to the
occupied states. Details of this are given in a paper outlining the methods
used in the turboTDDFT code. This is distributed with the quantum espresso
package. You can find it on the mt-student server in
/opt/share/quantum-espresso/doc-6.3/turboTDDFT-CPC.pdf
. There are also
slides on the quantum espresso website outlining the approach and making a
comparison to the more basic TDDFT method at
https://www.quantum-espresso.org/resources/tutorials/shanghai-2013/TDDFT_Talk_Gebauer.handouts.pdf
Despite the various advantages offered by the approaches used in this code, this will still be the most intensive calculation you’ll have done so far. However it is much faster than would be the case if we needed to manually look at many photon energies and including many empty states in our DFT calculation.
You should also keep in mind that this approach is for fixed nuclear positions and does not include the effect of interatomic vibrations which may also couple to an incoming photon, either directly if the energy range is similar, or indirectly if, for example, a photon creates an excited electron and vibration simultaneously. Generally by looking in a particular energy range, the dominant effects will be from one mechanism or another. For example, the homo-lumo gap in methane is around 10 eV in LDA-DFT, while the highest energy atomic vibrations are around 0.4 eV. This means that if we’re looking at the effect of photons in the 10-50 eV range (vacuum UV to extreme UV), the effect of the atomic vibrations will be minor. But we don’t immediately have the spectrum across all energy ranges.
As with some of the other more advanced calculations we looked at, we’ll be doing this calculation in several stages:
- A self consistent calculation of the molecule is performed. We can do this in the same way as previously. No special inputs are needed.
- We’ll use the
turbo_lanczos.x
code which does the TDDFPT calculation using the aforementioned Liouville-Lanczos method. This is the first, and more computationally intensive step in calculating the polarizability. It calculates the components of a tridiagonal matrix that will be used in a subsequent code, for each optical polarization direction. - We’ll use the
turbo_spectrum.x
code to take the tridiagonal matrix generated in the previous step, generate the polarizability and from it find both the susceptibility tensor and the oscillator strength as a function of photon energy.- The susceptibility tensor relates the electric field to the induced electric polarization
- the oscillator strength gives the probability of absorption or emission of photons.
Optical Spectrum of Methane
The directory CH4
contains a set of example inputs to perform
this calculation for a methane molecule.
- Start by examining the
pw.x
input file01_CH4_scf.in
. You’ll see that we have set this up as usual. We have specified both theprefix = CH4
andoutdir = './out'
as this calculation generates many intermediate files, so it will be easier to keep them all together. We have also set the positions such that the C atom is at the centre of the cell. This will make some subsequent visualization a little easier. - Run
pw.x
with this input file and examine the output to ensure it worked as expected. - Next take a look at the
turbo_lanczos.x
input file02_CH4_tl.in
. This is fairly short. As usual, we have accepted the default values for most of the input parameters. You can see the full list of inputs in theINPUT_LANCZOS.txt
file in the quantum espresso documentation directory. We have specified two sections in the input file:LR_INPUT
where we set theprefix
andoutdir
to match the scf calculation, andLR_CONTROL
where we setitermax = 400
which tells it how many elements of the Lanczos chain to calculate (elements of the tridiagonal matrix). While in principle more elements will give us a more accurate result, the subsequent code can extrapolate these to a very high number. Usually on the order of 500 to 1000 is sufficient to calculate explicitly, though we’ll look at how this converges later. It is connected to how well the spectral features can be resolved in your calculation.ipol = 4
which tells it to calculate the response in all 3 directions: x, y and z.
- Run
turbo_lanczos.x
with this input file. It will likely take somewhere around five minutes on the mt-student server. If you follow the output as it is generated, you will see it calculating each of the 400 requested iterations in turn for each of the three polarization directions. - Finally take a look at the
turbo_spectrum.x
input file03_CH4_ts.in
. We need to set a few more things here. As before, you can see a full description of all the input parameters that can be used with this code inINPUT_SPECTRUM.txt
in the quantum espresso documentation directory. It has a single section:LR_INPUT
. Again we set the prefix and outdir to match the previous calculations. Then we have:itermax0
which needs to match the number of iterations calculated explicitly byturbo_lanczos.x
.itermax
which says how many iterations in total to calculate. Those not explicitly calculated are generated using an extrapolation scheme.extrapolation
sets the approach used in the extrapolation.espil
sets a broadening to use (in Ry). Very small wiggles will appear in the output spectrum if this is too small, and if it is too large it may obscure important features.units
sets the output units. We use1
here we means the output energies are in eV (rather than Ry by default).start
,end
andincrement
specify the range of energies to output and the spacing between them (in the units specified by theunits
input).ipol
specifies the polarization direction. We set this to4
again for all three cartesian directions.
- Run
turbo_spectrum.x
with this input file. This is quite a quick calculation, and should finish in 10 to 20 seconds. The output file mainly contains the stages of the calculation. - The main output is instead in the file
CH4.plot_chi.dat
which will be generated. This lists each energy value and its corresponding susceptibility tensor and oscillator strength.
Take a look at the file CH4.plot_chi.dat
now. You’ll see the file is not
really amenable for plotting.
- First we can extract the oscillator strength with the following awk command:
awk '/S\(E\)=/{print $2, $3}' CH4.plot_chi.dat > ostrength.dat
This will print the second and third character on every line containing the
text S(E)=
(we need to use a \
to escape the parentheses). This extracts
the oscillator strength as a function of energy (in eV since we set units =
1
in the turbo_spectrum.x
input file).
- Plot the file containing the oscillator strength as a function of energy in gnuplot so you can see where the peaks in absorption are.
It’s also interesting to look at the susceptibility. As the system is
isotropic, the polarizability, and hence susceptibility in each of the three
Cartesian directions is the same. This means looking at the _11
component
alone is sufficient. This is a complex number, and we should look at both
the real and imaginary parts.
- Parse the energy, the real and the imaginary part of the susceptibility
(from the 11 component of chi in the
plot_chi.dat
file) and plot both the real and imaginary parts together in gnuplot.- If you parse this to a file with three columns: x, y1, y2, you can plot
them together in gnuplot as
plot "file" with lines, "file" using 1:3 with lines
. - You can plot several plots together by separating them with commas, and
you can use the keyword
using
to specify the columns to plot (it is1:2
by default).
- If you parse this to a file with three columns: x, y1, y2, you can plot
them together in gnuplot as
Charge Density Response
It can also be interesting to examine how the charge density actually varies in response to a photon at a given energy (and hence frequency) and polarization. In particular it’s useful to look at the response to photons at the energy where the highest peak in the oscillator strength is.
This can be done using an additional flag and input section in a
turbo_lanczos.x
calculation, following a previous turbo_lanczos.x
calculation for all polarization directions.
Take a look at the input file 04_CH4_tl.in
. In this
file we have set the LR_CONTROL
section inputs:
itermax = 400
which matches what we had previously.charge_response = 1
which turns on the calculation of the charge response.ipol = 1
which set the polarization of the calculated charge response to be along the x-direction.
Then we have added an addition section LR_POST
, which is needed whenever
charge_response
is set to 1. Here we have set the following:
omeg = 1.08
which set’s the energy of the photon for which we want to calculate the response. This is in Ry, and corresponds to 14.7 eV which is roughly where the peak of the oscillator strength is for this level of (under)convergence.epsil = 0.02
this is an energy broadening to use when calculating the response (in Ry).w_T_npol = 3
which says that we calculated the response to all three polarization directions in our previous calculation.plot_type = 3
which selects the output file to be in the gaussian cube format. This can be opened with xcrysden.- Note: it seems there may a bug with this output for the version compiled on the mt-student server such that for particular choices of energy cut-off, box-size and output format, the output file is cut off before the end and can’t be plot. For the values chosen here, it should work hopefully work correctly.
Run turbo_lanczos.x
with this input file. It’ll take a couple of minutes to
finish. Once completed you’ll see either a pair of files named
CH4-absorbtive-pol1.cube
and CH4-dispersive-pol1.cube
have been created,
or a single file named CH4-summed-rho-pol1.cube
. The former are generated
where the code determines the energy omeg
(including the broadening epsil
)
you are calculating at corresponds to a resonance where there is a peak in the
spectrum, and the latter is output otherwise. In our case we should be within
the broadening of a resonance peak and so two files should be output.
Open xcrysden
, which we used to look at crystal and molecular structures in
the extra section at the end of lab 2, (don’t forget to load the module). Use
this to open one of the generated files by navigating through File -> Open
Structure -> Gaussian98 Cube File, and selecting the file. You can look at the
regions where the charge density is increased and decreased in response to the
oscillating electric field of the photon by going to Tools -> Data Grid, click
OK. Then select an appropriate Isovalue in the menu that appears, tick the
Render +/- isovalue
box, and click Submit
. To find a good isovalue you can
repeat this process with different values within the listed maximum and
minimum on the grid till you see something interesting; the menu will stay
visible after you click submit.
Convergence
The calculated spectrum is quite sensitive to both the box size, the
plane-wave energy cut off, and the number of elements of the Lanczos chain
explicitly calculated, as set by the itermax
variable in the
turbo_lanczos.x
input file. (And subsequently used as itermax0
in the
turbo_spectrum.x
input file). These are all somewhat somewhat underconverged
in the input file given to you earlier so that the calculation would complete
in a reasonable time. And are somewhat connected to each other: e.g. the
spectrum may be converged at a higher value of itermax
for a higher value of
ecutwfc
so that testing convergence can be time consuming.
Underconverging your spectra can lead to the appearance of spurious peaks and shoulders and to changes in the locations of different peaks. So you will always need to explicitly check the convergence of your spectra with respect to these parameters.
As it would take quite a few heavy calculations to for you to test this, the
oscillator strengths for methane for various box sizes, energy cut-offs and
values of itermax
have been pre-calculated in the CH4_convergence
directory. For the various files here, in the filename the box dimension in
Angstrom is indicated by the number following A
, the plane-wave energy
cut-off in Ry is indicated by the number following e
, and the value of
itermax is indicated by the number following i
.
Try plotting and comparing combinations of these in gnuplot to see how the
spectrum changes as these parameters are changed.
Note: if you are testing the convergence with respect to
itermax
, you can generate output fromturbo_spectrum.x
whereitermax0
is any value up to the value chosen foritermax
in theturbo_lanczos.x
calculation. So for example, if you have calculated a spectrum withitermax = 500
inturbo_lanczos.x
, you can quickly generate output for say 400 by settingitermax0 = 400
inturbo_spectrum.x
, rather than doing a full recalculation with a lower value ofitermax
inturbo_lanczos.x
which would be far more time consuming.
Summary
In this lab you have seen how to use the turboTDDFT codes for performing and analysing TDDFT from the quantum espresso package to
- calculate the optical spectrum of a molecule by
- Performing a self-consistent DFT calculation with
pw.x
. - using the
turbo_lanczos.x
code to perform a TDDFPT calculation as the first step in calculating the polarizability. - using the
turbo_spectrum.x
code to generate the polarizability and find both the susceptibility tensor and oscillator strength as a function of photon energy.
- Performing a self-consistent DFT calculation with
- You have also used
turbo_spectrum.x
to calculate the charge density response of a molecule and visualised it withxcrysden
.
Extra - Li2
If you want to practice this type of calculation a bit more, the directory
Li2
contains a pseudopotential for lithium. You can use this
to do the following set of calculations:
- Find the optimal bond length for an Li2 molecule by relaxing the atomic positions as shown in previous labs.
- Calculate the oscillator strength as a function of energy.
- You will likely need to play around with the values of the energy range to
use in the
turbo_spectrum.x
calculation to find the best values to show the main peaks.