Tutorial:Optical Spectra from Casida

From OctopusWiki

In this tutorial we will again calculate the absorption spectrum of methane, but this time using Casida's equations. Once again our first step will be the calculation of the ground state (see the previous tutorial).

Unoccupied States

The Casida equation is a (pseudo-)eigenvalue equation written in the basis of particle-hole states. This means that we need both the occupied states -- computed in the ground-state calculation -- as well as the unoccupied states, that we will now obtain. The input file we will use is

CalculationMode = unocc
Units = eV_angstrom

radius = 4
spacing = 0.175

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

NumberUnoccStates = 10
EigenSolverMaxIter = 500

By running Octopus, you will obtain the first 10 unoccupied states. Note that we increased the value of EigenSolverMaxIter to 500, as the default 50 are insufficient to converge the unoccupied states to the required precision. You can take a look at the eigenvalues of the unoccupied states in the file static/eigenvalues:

All unoccupied states converged.
Criterium =      0.100000E-05

Eigenvalues [eV]
 #st  Spin   Eigenvalue     Occupation       Error
   1   --   -16.159264       2.000000      (0.0E+00)
   2   --    -9.060724       2.000000      (0.0E+00)
   3   --    -9.060724       2.000000      (0.0E+00)
   4   --    -9.060724       2.000000      (0.0E+00)
   5   --    -0.109982       0.000000      (9.0E-07)
   6   --     1.638999       0.000000      (9.8E-07)
   7   --     1.638999       0.000000      (9.0E-07)
   8   --     1.638999       0.000000      (9.2E-07)
   9   --     2.327133       0.000000      (9.7E-07)
  10   --     2.327133       0.000000      (9.0E-07)
  11   --     2.327133       0.000000      (9.0E-07)
  12   --     3.162588       0.000000      (9.8E-07)
  13   --     3.162588       0.000000      (9.9E-07)
  14   --     5.836262       0.000000      (9.8E-07)

Note that only the first unoccupied state is below the continuum, which you will probably immediately recognize as a deficiency of the LDA.

Casida calculation

We can now modify the CalculationMode to casida and rerun Octopus. Note that by default Octopus will use all occupied and unoccupied states that it has available. Sometimes, it is useful not to use all states. For example, if you have a molecule with 200 atoms and 400 occupied states ranging from -50 to -2 eV, and you are interested in looking at excitations in the visible, you can try to use only the states that are within 10 eV from the Fermi energy. You can select the states to use with LinearResponseKohnShamStates.

A new directory will appear named linear, where you can find the file linear/casida:

              E             <x>             <y>             <z>             <f>
 9.20088336E+00 -2.08444419E-01 -2.12670629E-01  2.69122278E-01  1.29685197E-01
 9.20088349E+00 -2.97034838E-01  2.69407196E-01 -1.71678474E-02  1.29685194E-01
 9.20088357E+00  1.71539838E-01  2.08076059E-01  2.97293050E-01  1.29685238E-01
 1.07120722E+01  8.90212407E-08 -1.21594999E-07  1.07681966E-07  3.21506958E-14
 1.07120722E+01  8.14629741E-08 -3.14971969E-08 -1.02332063E-07  1.69632157E-14
 1.07120722E+01 -3.30195091E-10 -3.91013251E-08  4.84439645E-08  3.63238701E-15
...

The <x>, <y> and <z> have obvious meaning as transition dipole moments:

where Φ₀ is the ground state and Φ_I is the given excitation. The oscillator strength, in atomic units, is given by:

although typically one considers the average over the three directions. The optical absorption spectrum can be given as the "strength function", which is

Note that the excitations are degenerate with degeneracy 3. This could already be expected from the T_d symmetry of methane. Further information concerning the excitations can be found in the directory linear/excitations. For example, the first excitation at 9.2 eV corresponds to the file linear/excitations/00001:

# Energy [eV] =    9.20088E+00
# <x> [A] =   -2.08444E-01
# <y> [A] =   -2.12671E-01
# <z> [A] =    2.69122E-01
 1 5 1 -3.07533510648477177E-8
 2 5 1 -1.94830984025791965E-2
 3 5 1 -0.19126894910470688
 4 5 1 0.97565019205719927
 1 6 1 -5.53379579328759467E-3
...

These files contain basically the eigenvector of the Casida equation. The first two columns are respectively the index of the occupied and the index of the unoccupied state, and the fourth is the coefficient of that state in the Casida eigenvector. This eigenvector is normalized to one, so in this case one can say that 95% (0.975²) of the excitation is from state 4 to state 5 (HOMO->LUMO) with a 4% (0.191²) contribution of the transition from state 3 (degenerate with the HOMO) to state 5.

oct-broad

Absorption spectrum of CH₄ calculated with time-propagation and with the Casida equation.

To visualize the spectrum, we need to broaden these delta functions with the utility oct-broad (run in your working directory, not in linear). It convolves the delta functions with Lorentzian functions. The operation of oct-broad is controlled by the variables LinBroadening (the width of this Lorentzian), LinEnergyStep, LinMinEnergy, LinMaxEnergy (not yet documented, but with obvious meaning). If you run oct-broad you obtain the file linear/spectrum.casida. It contains all columns of casida broadened. If you are interested in the total absorption spectrum, then you should plot the first and fifth columns. You should obtain a picture like the one on the right.

Comparing the spectrum obtained with the time-propagation (previous tutorial) and the Casida approach, we can see that

• The peaks of the time-propagation are broader. This can be solved by either increasing the total propagation time, or by increasing the broadening in the Casida approach.
• The first two peaks are nearly the same. Probably also the third is OK, but the low resolution of the time-propagation does not allow to distinguish the two close peaks that compose it.
• For high energies the spectra differ a lot. The reason is that we only used 10 empty states in the Casida approach. In order to describe better this region of the spectrum we would need more.

You probably noticed that the Casida calculation took much less time than the time-propagation. This is clearly true for small or medium-sized systems. However, the implementation of Casida in Octopus has a much worse scaling with the size of the system than the time-propagation, so for larger systems the situation may be different. Note also that in Casida one needs a fair amount of unoccupied states which are fairly difficult to obtain.