# Tutorial:Optical spectra from time-propagation

In this tutorial we will learn how to obtain the absorption spectrum of a molecule from the explicit solution of the time-dependent Kohn-Sham equations. We choose as a test case methane (CH4).

## Ground state

Before starting out time-dependent simulations, we need to obtain the ground state of the system. For this we use basically the same inp file as in a previous tutorial:

CalculationMode = gs
Units = eV_angstrom

Radius = 4
Spacing = 0.18

CH = 1.087
%Coordinates
"C" |           0 |          0 |           0
"H" |  CH/sqrt(3) | CH/sqrt(3) |  CH/sqrt(3)
"H" | -CH/sqrt(3) |-CH/sqrt(3) |  CH/sqrt(3)
"H" |  CH/sqrt(3) |-CH/sqrt(3) | -CH/sqrt(3)
"H" | -CH/sqrt(3) | CH/sqrt(3) | -CH/sqrt(3)
%


After running Octopus, we will have the Kohn-Sham wave-functions of the ground-state in the directory restart/gs. As we are going to propagate these wave-functions, they have to be well converged. It is important not only to converge the energy (that is relatively easy to converge), but the density. The default ConvRelDens = 1e-05 is usually enough, though.

## Time-dependent run

To calculate absorption, we excite the system with an infinitesimal electric-field pulse, and then propagate the time-dependent Kohn-Sham equations for a certain time T. The spectrum can then be evaluated from the time-dependent dipole moment.

This is what you need to add to the input file (don't forget to change the CalculationMode to "td"):

TDDeltaStrength = 0.01
TDPolarizationDirection = 1

TDPropagator = aetrs
TDTimeStep = 0.0023
TDMaxSteps = 4350  # ~ 10.0/TDTimeStep


First we define our perturbation. Its strength is given by TDDeltaStrength. This number should be small to keep the response linear, but should be sufficiently large to avoid numerical problems. The variable TDPolarizationDirection sets the polarization of our perturbation to be on the first axis (x).

We then set up the time-evolution. Octopus has a large number of possible propagators that you can use (see J. Chem. Phys 121 3425-3433 (2004) for an overview). Here we use the simple Approximately Enforced Time-Reversal Symmetry (aetrs). Then we need a time step (variable TDTimeStep). The choice of this number depends on the propagation method, and it is a quite crucial parameter. You should have learned how to set it up in the previous tutorial. Finally, we set the number of time steps with the variable TDMaxSteps. To have a maximum propagation time of 10 we will need around 4350 iterations.

Note that you will be calculating the singlet dipole spectrum. You can also obtain the triplet by using TDDeltaStrengthMode, and other multipole responses by using TDKickFunction. For details on triplet calculations see Tutorial:Triplet Excitations.

You can now start Octopus and go for a quick coffee (this should take around 5 minutes depending on your machine). Propagations are slow, but the good news is that they scale very well with the size of the system. This means that even if methane is very slow, a molecule with 200 atoms can still be calculated without big problems.

## Symmetries

The proper way to get the spectrum from fewer than three runs is to take advantage of the symmetries of the system. However, this is fairly complicated to understand for the purposes of this introductory tutorial. To learn more about symmetry, look at this page: Tutorial:Optical Spectra from TD:Symmetries. In this case, we can actually just do one run, and plot the absorption spectrum for x-polarization, and that will in fact be equivalent to the spectrum averaged over the three directions.

## Optical spectra

When the run is finished, you should find the file td.general/multipoles that contains all the information needed to compute the spectrum of methane. For a general system we need three time-dependent runs to get the multipoles files called td.general/multipoles.1, td.general/multipoles.2, and td.general/multipoles.3.

The beginning of the multipoles file will look similar to this:

################################################################################
# nspin         1
# lmax          1
# kick mode     0
# kick strength    0.005291772086
# pol(1)           1.000000000000    0.000000000000    0.000000000000
# pol(2)           0.000000000000    1.000000000000    0.000000000000
# pol(3)           0.000000000000    0.000000000000    1.000000000000
# direction     1
# Equiv. axes   0
# wprime           0.000000000000    0.000000000000    1.000000000000
# kick time        0.000000000000
# Iter             t          Electronic charge(1)       <x>(1)              <y>(1)              <z>(1)
#[Iter n.]     [hbar/eV]           Electrons              [A]                 [A]                 [A]
################################################################################
0  0.000000000000e+00  8.000000000000e+00 -5.105069841261e-11 -2.458728679159e-11 -6.191777548908e-11
1  2.300000000000e-03  7.999999999898e+00 -1.441204189916e-03 -2.454483489431e-11 -6.178536414960e-11
2  4.600000000000e-03  7.999999999794e+00 -2.849596623682e-03 -2.441311503813e-11 -6.139129111463e-11
3  6.900000000000e-03  7.999999999692e+00 -4.212595751069e-03 -2.420204397245e-11 -6.075589410104e-11


To obtain the spectrum, in your working directory run the utility oct-propagation_spectrum. Don't worry about the warnings generated, we know what we are doing!

File "multipoles" found. This will be the only file to be processed.
(If more than one file is to be used, the files should be called
"multipoles.1", "multipoles.2", etc.

** Warning:
**   The file "multipoles" tells me that the system has no usable symmetry.
**   However, I am only using this file; cannot calculate the full tensor.
**   A file "XXXX_vector" will be generated instead.

Octopus emitted      1 warning(s).


You will notice that the file cross_section_vector is created. If you have all the information (either a symmetric molecule and one multipole file, or a less symmetric molecule and multiple multipole files), oct-propagation_spectrum will generate the whole polarizability tensor, cross_section_tensor. If not enough multipole files are available, it can only generate the response to the particular polarization you chose in you time-dependent run, cross_section_vector. Let us look first at cross_section_vector:

# nspin         1
# kick mode    0
# kick strength    0.005291772086
# pol(1)           1.000000000000    0.000000000000    0.000000000000
# pol(2)           0.000000000000    1.000000000000    0.000000000000
# pol(3)           0.000000000000    0.000000000000    1.000000000000
# direction    1
# Equiv. axes  0
# wprime           0.000000000000    0.000000000000    1.000000000000
# kick time        0.000000000000
#%


The beginning of the file just repeats some information concerning the run.

# Number of time steps =     4350
# PropagationSpectrumDampMode   =    2
# PropagationSpectrumDampFactor =     4.0817
# PropagationSpectrumStartTime  =     0.0000
# PropagationSpectrumEndTime    =    10.0050
# PropagationSpectrumMaxEnergy  =    20.0000
# PropagationSpectrumEnergyStep =     0.0100


Now comes a summary of the variables used to calculate the spectra. Note that you can change all these settings just by adding these variables to the input file before running oct-propagation_spectrum. Of special importance are perhaps PropagationSpectrumMaxEnergy that sets the maximum energy that will be calculated, and PropagationSpectrumEnergyStep that determines how many points your spectrum will contain. To have smoother curves, you should reduce this last variable.

# Electronic sum rule       =         3.946570
# Static polarizability (from sum rule) =         2.145862 A



Now comes some information from the sum rules. The first is just the f-sum rule, which should yield the number of active electrons in the calculations. We have 8 valence electrons (4 from carbon and 4 from hydrogen), but the sum rule gives a number that is much smaller than 8! The reason is that we are just summing our spectrum up to 20 eV (see PropagationSpectrumMaxEnergy), but for the sum rule to be fulfilled, we should go to infinity. Of course infinity is a bit too large, but increasing 20 eV to a somewhat larger number will improve dramatically the f-sum rule. The second number is the static polarizability also calculated from the sum rule. Again, do not forget to converge this number with PropagationSpectrumMaxEnergy.

#       Energy        sigma(1, nspin=1)   sigma(2, nspin=1)   sigma(3, nspin=1)   StrengthFunction(1)
#        [eV]               [A^2]               [A^2]               [A^2]               [1/eV]
0.00000000E+00     -0.00000000E+00     -0.00000000E+00     -0.00000000E+00     -0.00000000E+00
0.10000000E-01     -0.13419076E-08     -0.33805821E-12     -0.85198023E-12     -0.12225726E-08
0.20000000E-01     -0.53784954E-08     -0.13506268E-11     -0.34038646E-11     -0.49001891E-08
...


Finally comes the spectrum. The first column is the energy (frequency), the next three columns are a row of the cross-section tensor, and the last one is the strength function for this run.

The dynamic polarizability is related to optical absorption cross-section via in atomic units, or more generally (where is the fine-structure constant) or . The cross-section is related to the strength function by .

If all three directions were done, we would also have a second file, cross_section_tensor:

# nspin         1
# kick mode    0
# kick strength    0.010000000000
# pol(1)           1.000000000000    0.000000000000    0.000000000000
# pol(2)           0.000000000000    1.000000000000    0.000000000000
# pol(3)           0.000000000000    0.000000000000    1.000000000000
# direction    1
# Equiv. axis  3
# wprime           0.000000000000    0.000000000000    1.000000000000
# kick time        0.000000000000
#       Energy         (1/3)*Tr[sigma]    Anisotropy[sigma]      sigma(1,1,1)        sigma(1,2,1)        sigma(1,3,1)         ...
#        [eV]               [A^2]               [A^2]               [A^2]               [A^2]               [A^2]             ...
0.00000000E+00      0.00000000E+00      0.00000000E+00      0.00000000E+00      0.00000000E+00      0.00000000E+00      ...
0.10000000E-01      0.51922187E-09      0.12093134E-10      0.51922187E-09      0.46057845E-12     -0.60377899E-11      ...
0.20000000E-01      0.20467565E-08      0.48314996E-10      0.20467565E-08      0.18401322E-11      0.20467565E-08      ...
...


The columns are now the energy, the average absorption coefficient (Tr ), the anisotropy, and then all the 9 components of the tensor. The third number is the spin component, just 1 here since it is unpolarized. The anisotropy is defined as The reason for this definition is that it is identically equal to: where are the eigenvalues of . An "isotropic" tensor is characterized by having three equal eigenvalues, which leads to zero anisotropy. The more different that the eigenvalues are, the larger the anisotropy is.

If you now plot the absorption spectrum (column 5 vs 1 in cross_section_vector), you should obtain the plot shown on the right. Of course, you should now try to converge this spectrum with respect to the calculation parameters. In particular:

• Increasing the total propagation time will reduce the width of the peaks. In fact, the width of the peaks in this methods is absolutely artificial, and is inversely proportional to the total propagation time. Do not forget, however, that the area under the peaks has a physical meaning: it is the oscillator strength of the transition.
• As you are running in a box with zero boundary conditions, you will see the box states in the spectrum. This can be improved by increasing the radius of the box. However, as one is usually interested in bound-bound transitions in the optical or near-ultraviolet regimes, a fairly small box (like the one used in this tutorial) is often enough.

Some questions to think about:

• What is the equation for the peak width in terms of the propagation time? How does the observed width compare to your expectation?
• What is the highest energy excitation obtainable with the calculation parameters here? Which is the key one controlling the maximum energy?
• Why does the spectrum go below zero? Is this physical? What calculation parameters might change this?

Back to [[Tutorial Series:{{{series}}}|{{{series}}}]]