# Tutorial:Optical Spectra from Casida

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, via a non-self-consistent calculation using the density computed in `gs`. The input file we will use is

`CalculationMode`

= unocc`UnitsOutput`

= eV_angstrom`Radius`

= 3.5*angstrom`Spacing`

= 0.18*angstrom CH = 1.097*angstrom %`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) %`ExtraStates`

= 10

By running Octopus, you will obtain the first 10 unoccupied states (you need a ground-state calculation first). The solution of the unoccupied states is controlled by the variables in section `SCF::Eigensolver`

. You can take a look at the eigenvalues of the unoccupied states in the file ** static/eigenvalues**:

All states converged. Criterion = 0.100000E-05 Eigenvalues [eV] #st Spin Eigenvalue Occupation Error 1 -- -16.865648 2.000000 (5.1E-07) 2 -- -9.433544 2.000000 (4.1E-07) 3 -- -9.433544 2.000000 (4.4E-07) 4 -- -9.433544 2.000000 (3.8E-07) 5 -- 0.687525 0.000000 (8.7E-07) 6 -- 2.568509 0.000000 (8.1E-07) 7 -- 2.568509 0.000000 (9.3E-07) 8 -- 2.568509 0.000000 (8.0E-07) 9 -- 3.186757 0.000000 (8.1E-07) 10 -- 3.186757 0.000000 (7.9E-07) 11 -- 3.186757 0.000000 (9.3E-07) 12 -- 4.065409 0.000000 (3.9E-07) 13 -- 4.065409 0.000000 (5.9E-07) 14 -- 7.707535 0.000000 (5.8E-07)

## Casida calculation

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 could select the states to use with `CasidaKohnShamStates`

or `CasidaKSEnergyWindow`

.

A new directory will appear named ** casida**, where you can find the file

**:**

`casida/casida`E [eV] <x> [A] <y> [A] <z> [A] <f> 1 1.03574479E+01 -4.58729232E-02 -4.70248327E-02 -3.92558757E-01 1.43553318E-01 2 1.03574479E+01 -2.65971628E-01 2.96070779E-01 -4.38604830E-03 1.43553330E-01 3 1.03574480E+01 -2.92528292E-01 -2.61818268E-01 6.55471171E-02 1.43553314E-01 4 1.20129521E+01 3.30909578E-08 4.43526414E-08 -9.50853878E-08 1.27207569E-14 5 1.20129521E+01 1.67462244E-08 7.03225255E-08 -3.81171337E-08 7.01926516E-15 6 1.20129522E+01 4.29169979E-08 -1.77764987E-07 -9.41377619E-08 4.44619857E-14 7 1.21052958E+01 -1.81367255E-07 -1.96694143E-07 1.77922847E-07 1.09339205E-13 8 1.21052959E+01 -6.47696451E-07 -7.08655209E-08 1.59024918E-07 4.76399730E-13 ...

The <x>, <y> and <z> are the transition dipole moments:

where is the ground state and is the given excitation. The oscillator strength is given by:

as 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 symmetry of methane.

*Note that within the degenerate subspaces, there is some arbitrariness (possibly dependent on the details of your compilation and machine) in the linear combinations of transitions. Therefore, you should not be concerned if you do not have the same results for the components of the transition dipole moments (above) and analysis of the excitations (below). However, the energies and resulting spectra should agree.*

Further information concerning the excitations can be found in the directory ** casida/casida_excitations**. For example, the first excitation at 9.96 eV is analyzed in the file

**:**

`casida/casida_excitations/00001`# Energy [eV] = 9.96919E+00 # <x> [A] = -1.84120E-01 # <y> [A] = -1.90185E-01 # <z> [A] = -2.91670E-01 1 5 1 -3.1386000894250019E-008 2 5 1 -7.6787866663895837E-002 3 5 1 0.82480869066203977 4 5 1 0.55798340346785968 1 6 1 1.9670247480165549E-003 ...

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, the third is the spin index (always 1 when spin-unpolarized), 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 99.6% (-0.998^{2}) of the excitation is from state 3 to state 5 (one of 3 HOMO orbitals->LUMO) with small contribution from some other transitions.

## Absorption spectrum

To visualize the spectrum, we need to broaden these delta functions with the utility ** oct-casida_spectrum** (run in your working directory, not in

**). It convolves the delta functions with Lorentzian functions. The operation of**

`casida`**is controlled by the variables**

`oct-casida_spectrum``CasidaSpectrumBroadening`

(the width of this Lorentzian), `CasidaSpectrumEnergyStep`

, `CasidaSpectrumMinEnergy`

, and `CasidaSpectrumMaxEnergy`

. If you run **you obtain the file**

`oct-casida_spectrum`**. It contains all columns of**

`casida/spectrum.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.**

`casida`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.

If you are interested, you may also compare the Casida results against the Kohn-Sham eigenvalue differences in ** casida/spectrum.eps_diff** and the Petersilka approximation to Casida in

**.**

`casida/spectrum.petersilka`