---
authors: XG
---

# TDDFT (Casida)

## Time-Dependent Density Functional Theory, Casida's approach.

This tutorial aims at showing how to get the following physical properties, for finite systems:

  * Excitation energies
  * Associated oscillator strengths
  * Frequency-dependent polarizability and optical spectra

in the Casida approach, within Time-Dependent Density Functional Theory.

This tutorial should take about 30 minutes.

[TUTORIAL_README]

## Brief theoretical introduction

In order to perform Time-Dependent Density Functional Theory calculations (TDDFT)
of electronic excitations and oscillator strengths, in the Casida's approach,
you should first have some theoretical background.
TDDFT was first developed in the eighties, but the direct calculation of
electronic excitations was introduced much later, by Casida and co-workers.
A comprehensive description of the underlying formalism is given in

    "Time-Dependent Density Functional Response Theory of Molecular Systems:
     Theory, Computational Methods, and Functionals"
     M. E. Casida
     in Recent Development and Applications of Modern Density Functional Theory,
     edited by J.M. Seminario (Elsevier, Amsterdam, 1996), p. 391.
     <http://dx.doi.org/10.1016/S1380-7323(96)80093-8>

However this reference might be hard to get, that is why we have based the
tutorial instead on the following (also early) papers:
[[cite:Casida1998]], [[cite:Casida1998a]], and [[cite:Vasiliev1998]].

The first of these papers, [[cite:Casida1998]], will be used as main reference for our tutorial.

From these papers, you will learn that a TDDFT calculation of electronic
excitation energies start first from a usual ground-state calculation, with a
chosen XC functional. Such a calculation produces a spectrum of Kohn-Sham
electronic energies. It is widely known that differences between occupied and
unoccupied Kohn-Sham electronic energies resemble excitation energies (the
difference in energy between an excited state and the ground state), although
there is no real theoretical justification for this similarity.

These differences between Kohn-Sham electronic energies are the starting point
of Casida's approach: in the framework of TDDFT, their square give the main
contribution to the diagonal part of a matrix, whose eigenvalues will be the
square of the excitation energies. One has to add to the diagonal matrix made
from the squares of Kohn-Sham energy differences, a coupling matrix, whose
elements are four-wavefunction integrals of the Coulomb and exchange-correlation kernel.
The exchange-correlation kernel contribution will differ
in the spin-singlet and in the spin-triplet states, this being the only
difference between spin-singlet and spin-triplet states. See Eqs.(1.3) and
(1.4) of [[cite:Casida1998a]], and Eqs.(1-2) of [[cite:Vasiliev1998]].

The construction of the coupling matrix can be done on the basis of an
exchange-correlation kernel that is derived from the exchange-correlation
functional used for the ground-state, but this is not a requirement of the
theory, since such a correspondence only holds for the exact functional. In
practice, the approximation to the XC potential and the one to the XC kernel
are often different. See section III of [[cite:Casida1998]].

A big drawback of the currently known XC potentials and XC kernels is observed
when the system is infinite in at least one direction (e.g. polymers, slabs,
or solids). In this case, the addition of the coupling matrix is unable to
shift the edges of the Kohn-Sham band structure (each four-wavefunction
integral becomes too small). There is only a redistribution of the oscillator
strengths. In particular, the DFT band gap problem is NOT solved by TDDFT.
Also, the Casida's approach relies on the discreteness of the Kohn-Sham spectrum.

Thus, the TDDFT approach to electronic excitation energies in ABINIT is ONLY
valid for finite systems (atoms, molecules, clusters). Actually, only one
k-point can be used, and a "box center" must be defined, close to the center
of gravity of the system.

The Casida formalism also gives access to the oscillator strengths, needed to
obtain the frequency-dependent polarizability, and corresponding optical
spectrum. In the ABINIT implementation, the oscillators strengths are given as
a second-rank tensor, in cartesian coordinates, as well as the average over
all directions usually used for molecules and clusters. It is left to the user
to generate the polarisability spectrum, according to e.g. Eq.(1.2) of [[cite:Casida1998a]].

One can also combine the ground state total energy with the electronic
excitation energies to obtain Born-Oppenheimer potential energy curves for
excited states. This is illustrated for formaldehyde in [[cite:Casida1998a]].

Given its simplicity, and the relatively modest CPU cost of this type of
calculation, Casida's approach enjoys a wide popularity. There has been
hundreds of papers published on the basis of methodology. Still, its accuracy
might be below the expectations, as you will see. As often, test this method
to see if it suits your needs, and read the recent literature ...

## A first computation of electronic excitation energies and oscillator strengths, for N$_2$

We will now compute and analyse the excitation energies of the diatomic molecule N$_2$.
This is a rather simple system, with cylindrical symmetry,
allowing interesting understanding. Although we will suppose that you are
familiarized with quantum numbers for diatomic molecules, this should not play
an important role in the understanding of the way to use Abinit
implementation of Casida's formalism.

*Before beginning, you might consider to work in a different subdirectory as
for the other tutorials. Why not Work_tddft?*
Copy the file *ttddft_1.abi* in *Work_tddft*:

```sh
cd $ABI_TESTS/tutorial/Input
mkdir Work_tddft
cd Work_tddft
cp ../ttddft_1.abi .
```

So, issue now:

    abinit ttddft_1.abi > log &

The computation is quite fast: about 3 secs on a 2.8 GHz PC.
Let's examine the input file *ttddft_1.abi*.

{% dialog tests/tutorial/Input/ttddft_1.abi %}

There are two datasets: the first one corresponds to a typical ground-state
calculation, with only occupied bands. The density and wavefunctions are
written, for use in the second data set. The second dataset is the one where
the TDDFT calculation is done. Moreover, the non-self-consistent calculation
of the occupied eigenfunctions and corresponding eigenenergies is also
accomplished. This is obtained by setting [[iscf]] to -1.

Please, take now some time to read the information about this value of [[iscf]], and the few
input variables that acquire some meaning in this context (namely,
[[boxcenter]], [[td_mexcit]], and [[td_maxene]]). Actually, this is most of
the information that should be known to use the TDDFT in ABINIT!

You will note that we have 5 occupied bands (defined for dataset 1), and that
we add 7 unoccupied bands in the dataset 2, to obtain a total of 12 bands. The
box is not very large (6x5x5 Angstrom), the cutoff is quite reasonable, 25
Hartree), and as requested for the Casida's formalism, only one k point is
used. We are using the Perdew-Wang 92 LDA functional for both the self-
consistent and non-self-consistent calculations ([[ixc]] = -1012), as deduced
by abinit by looking at the pseudopotential file..

We can now examine the output file *ttddft_1.abo.*

{% dialog tests/tutorial/Refs/ttddft_1.abo %}

One can jump to the second dataset section, and skip a few non-interesting
information, in order to reach the following information:

     *** TDDFT : computation of excited states ***
     Splitting of  12 bands in   5 occupied bands, and   7 unoccupied bands,
     giving    35 excitations.

The matrix that is diagonalized, in the Casida's formalism, is thus a 35x35 matrix.
It will give 35 excitation energies.
Then, follows the list of excitation energies, obtained from the difference of
Kohn-Sham eigenvalues (occupied and unoccupied), for further reference. They
are ordered by increasing energy. In order to analyze the TDDFT as well as
experimental data in the next section, let us mention that the Kohn-Sham
eigenfunctions in this simulation have the following characteristics:

  * The first and fifth states are (non-degenerate) occupied $\sigma$ states (m=0), with even parity
  * The second state is a (non-degenerate) occupied $\sigma$ state (m=0), with odd parity
  * The third and fourth states are degenerate occupied $\pi$ states (m=+1,-1), with odd parity
  * The sixth and seventh states are degenerate unoccupied $\pi$ states (m=+1,-1), with even parity
  * The state 8 is a (non-degenerate) unoccupied $\sigma$ state (m=0), with even parity

Combining states 3,4 and 5 with 6, 7 and 8, give the first nine Kohn-Sham energy differences:

      Transition  (Ha)  and   (eV)   Tot. Ene. (Ha)  Aver     XX       YY       ZZ
       5->  6 3.06494E-01 8.34013E+00 -2.03684E+01 0.0000E+00 0.00E+00 0.00E+00 0.00E+00
       5->  7 3.06494E-01 8.34013E+00 -2.03684E+01 0.0000E+00 0.00E+00 0.00E+00 0.00E+00
       5->  8 3.50400E-01 9.53488E+00 -2.03245E+01 0.0000E+00 0.00E+00 0.00E+00 0.00E+00
       4->  6 3.60026E-01 9.79682E+00 -2.03149E+01 5.5534E-01 1.67E+00 0.00E+00 0.00E+00
       3->  6 3.60026E-01 9.79682E+00 -2.03149E+01 3.7821E-03 1.13E-02 0.00E+00 0.00E+00
       4->  7 3.60026E-01 9.79682E+00 -2.03149E+01 3.7822E-03 1.13E-02 0.00E+00 0.00E+00
       3->  7 3.60026E-01 9.79682E+00 -2.03149E+01 5.5534E-01 1.67E+00 0.00E+00 0.00E+00
       4->  8 4.03933E-01 1.09916E+01 -2.02710E+01 3.7976E-02 0.00E+00 1.14E-01 0.00E+00
       3->  8 4.03933E-01 1.09916E+01 -2.02710E+01 3.7976E-02 0.00E+00 0.00E+00 1.14E-01

Without the coupling matrix, these would be the excitation energies, for both
the spin-singlet and spin-triplet states. The coupling matrix modifies the
eigenenergies, by mixing different electronic excitations, and also lift some
degeneracies, e.g. the quadruplet formed by the combination of the degenerate
states 3-4 and 6-7 that gives the (four-fold degenerate) excitation energies 
with 3.60026E-01 Ha in the above table.

Indeed, concerning the spin-singlet, the following excitation energies are
obtained (see the next section of the output file):

      TDDFT singlet excitation energies (at most 20 of them are printed),
      and corresponding total energies.
      Excit#   (Ha)    and    (eV)    total energy (Ha)    major contributions
       1    3.45362E-01   9.39779E+00   -2.032952E+01    1.00(  5->  7)  0.00(  1->  7)
       2    3.45434E-01   9.39975E+00   -2.032945E+01    1.00(  5->  6)  0.00(  1->  6)
       3    3.60026E-01   9.79681E+00   -2.031486E+01    0.50(  3->  6)  0.50(  4->  7)
       4    3.68693E-01   1.00326E+01   -2.030619E+01    0.99(  5->  8)  0.00(  2-> 10)
       5    3.83765E-01   1.04428E+01   -2.029112E+01    0.50(  4->  7)  0.50(  3->  6)
       6    3.83798E-01   1.04437E+01   -2.029108E+01    0.50(  4->  6)  0.50(  3->  7)
       7    4.03285E-01   1.09740E+01   -2.027160E+01    0.99(  3->  8)  0.01(  4->  8)
       8    4.03304E-01   1.09745E+01   -2.027158E+01    0.99(  4->  8)  0.01(  3->  8)
       9    4.59051E-01   1.24914E+01   -2.021583E+01    0.91(  2->  8)  0.04(  3->  7)
       ...

The excitation energies are numbered according to increasing energies, in Ha
as well as in eV. The total energy is also given (adding excitation energy to
the ground-state energy), and finally, the two major contributions to each
of these excitations are mentioned (size of the contribution then identification).

It is seen that the first and second excitations are degenerate (numerical
inaccuracies accounts for the meV difference), and mainly comes from the first
and second Kohn-Sham energy differences (between occupied state 5 and
unoccupied states 6 and 7). This is also true for the fourth excitation, that
comes from the third Kohn-Sham energy difference (between occupied state 5 and
unoccupied state 8). The quadruplet of Kohn-Sham energy differences, that was
observed at 3.60026E-01 Ha, has been split into one doublet and two singlets,
with numbers 3 (the lowest singlet), 5-6 (the doublet) while the last singlet
is not present in the 20 lowest excitations.

The list of oscillator strength is then provided.

      Oscillator strengths :  (elements smaller than 1.e-6 are set to zero)
      Excit#   (Ha)   Average    XX        YY        ZZ         XY        XZ        YZ
       1 3.45362E-01 0.000E+00 0.000E+00 0.000E+00 0.000E+00  0.00E+00  0.00E+00  0.00E+00
       2 3.45434E-01 0.000E+00 0.000E+00 0.000E+00 0.000E+00  0.00E+00  0.00E+00  0.00E+00
       3 3.60026E-01 0.000E+00 0.000E+00 0.000E+00 0.000E+00  0.00E+00  0.00E+00  0.00E+00
       4 3.68693E-01 0.000E+00 0.000E+00 0.000E+00 0.000E+00  0.00E+00  0.00E+00  0.00E+00
       5 3.83765E-01 0.000E+00 0.000E+00 0.000E+00 0.000E+00  0.00E+00  0.00E+00  0.00E+00
       6 3.83798E-01 0.000E+00 0.000E+00 0.000E+00 0.000E+00  0.00E+00  0.00E+00  0.00E+00
       7 4.03285E-01 5.881E-02 0.000E+00 1.409E-03 1.750E-01  0.00E+00  0.00E+00 -1.57E-02
       8 4.03304E-01 5.685E-02 0.000E+00 1.692E-01 1.361E-03  0.00E+00  0.00E+00  1.52E-02
       9 4.59051E-01 8.613E-02 2.584E-01 0.000E+00 0.000E+00  0.00E+00  0.00E+00  0.00E+00
      10 4.61447E-01 0.000E+00 0.000E+00 0.000E+00 0.000E+00  0.00E+00  0.00E+00  0.00E+00
       ...

The first six transitions are forbidden, with zero oscillator strength. The
seventh and eighth transitions are allowed, with sizeable YY, YZ and ZZ components.

Next, one finds the excitation energies for the spin-triplet states:

      TDDFT triplet excitation energies (at most 20 of them are printed),
      and corresponding total energies.
      Excit#   (Ha)    and    (eV)    total energy (Ha)    major contributions
       1    2.84779E-01   7.74923E+00   -2.039010E+01    1.00(  5->  7)  0.00(  5->  6)
       2    2.84781E-01   7.74928E+00   -2.039010E+01    1.00(  5->  6)  0.00(  5->  7)
       3    2.97188E-01   8.08689E+00   -2.037770E+01    0.50(  3->  7)  0.50(  4->  6)
       4    3.30296E-01   8.98780E+00   -2.034459E+01    0.50(  4->  7)  0.50(  3->  6)
       5    3.30344E-01   8.98913E+00   -2.034454E+01    0.50(  4->  6)  0.50(  3->  7)
       6    3.43692E-01   9.35234E+00   -2.033119E+01    1.00(  5->  8)  0.00(  2-> 10)
       7    3.60026E-01   9.79681E+00   -2.031486E+01    0.50(  3->  6)  0.50(  4->  7)
       8    3.85278E-01   1.04839E+01   -2.028961E+01    0.91(  2->  7)  0.09(  3->  8)
       9    3.85310E-01   1.04848E+01   -2.028957E+01    0.91(  2->  6)  0.08(  4->  8)
       ...

Spin-triplet energies are markedly lower than the corresponding spin-singlet
energies. Also, the highest singlet derived from the Kohn-Sham quadruplet is
now the excitation number 7. The oscillator strengths also follow. At this
stage, we are in position to compare with experimental data, and try to
improve the quality of our calculation.

To summarize our results, we obtain the following five lowest-lying spin-
singlet excitation energies, with corresponding quantum numbers (that we
derive from the knowledge of the Kohn-Sham states quantum numbers):

    9.40 eV   m=+1,-1  even parity (Pi_g state)
    9.80 eV   m=0      odd parity  (Sigma_u state)
    10.03 eV  m=0      even parity (Sigma_g state)
    10.44 eV  m=+2,-2  odd parity  (Delta_u state)
    10.97 eV  m=+1,-1  odd parity  (Pi_u state)

and the following five lowest-lying spin-triplet excitations energies, with
corresponding quantum numbers:

    7.75 eV   m=+1,-1  even parity (Pi_g state)
    8.09 eV   m=0      odd parity  (Sigma_u state)
    8.99 eV   m=+2,-2  odd parity  (Delta_u state)
    9.35 eV   m=0      even parity (Sigma_g state)
    9.80 eV   m=0      odd parity  (Sigma_u state)

The quantum number related to the effect of a mirror plane, needed for $\Sigma$
states, could not be attributed on the sole basis of the knowledge of Kohn-
Sham orbitals quantum numbers.

The lowest-lying experimental spin-singlet excitation energies, see table III
of [[cite:Casida1998]], are as follows:

    9.31 eV   m=+1,-1  even parity (Pi_g state)
    9.92 eV   m=0      odd parity  (Sigma_u- state)
    10.27 eV  m=+2,-2  odd parity  (Delta_u state)

and the lowest-lying experimental spin-triplet excitations energies are:

    7.75 eV   m=0      odd parity  (Sigma_u+ state)
    8.04 eV   m=+1,-1  even parity (Pi_g state)
    8.88 eV   m=+2,-2  odd parity  (Delta_u state)
    9.67 eV   m=0      odd parity  (Sigma_u- state)

In several cases, the agreement is quite satisfactory, on the order of 0.1-0.2
eV. However, there are also noticeable discrepancies. Indeed, we have to understand, in our simulation:

  * The appearance of the spin-singlet $^1\Sigma_g$ state at 10.03 eV (Spin-singlet state 3)
  * The inversion between the spin-triplet $^3\Pi_g$ and $^3\Sigma_u$ states (Spin-triplet states 1 and 2)
  * The appearance of the spin-triplet $^3\Sigma_g$ state at 9.35 eV (Spin-triplet state 4)

Still, the agreement between these TDDFT values and the experimental values is
much better than anything that can be done on the sole basis of Kohn-Sham
energy differences, that are (for spin-singlet and -triplet):

    8.34 eV   m=+1,-1   even parity (Pi_g state)
    9.53 eV   m=0       odd parity  (Sigma_u state)
    9.80 eV   m=0(twice),+2,-2 odd parity  (Sigma_u and Delta_u states)
    10.99 eV  m=+1,-1   odd parity  (Pi_u state)

## Convergence studies

There are several parameters subject to convergence studies in this context:
the energy cut-off, the box size, and the number of unoccupied bands.

We will start with the number of unoccupied states. The only input parameter
to be changed in the input file is the value of nband2. The following results
are obtained, for nband2 = 12, 30, 60, 100 and 150 (Energies given in eV):

    Singlet 1 :  9.40   9.37   9.33   9.30   9.28
    Singlet 2 :  9.80   9.80   9.80   9.79   9.80
    Singlet 3 : 10.03   9.91   9.85   9.83   9.83
    Singlet 4 : 10.44  10.43  10.43  10.42  10.42
    Singlet 5 : 10.97  10.97  10.97  10.97  10.97
    Triplet 1 :  7.75   7.75   7.73   7.73   7.72
    Triplet 2 :  8.06   8.02   7.98   7.96   7.95
    Triplet 3 :  8.99   8.97   8.96   8.96   8.95
    Triplet 4 :  9.35   9.34   9.34   9.34   9.34
    Triplet 5 :  9.80   9.80   9.80   9.80   9.80

You might try to obtain one of these...

The computation with nband2 = 100 takes a bit more than 1 minute,
and gives a result likely converged within 0.01 eV.
Let us have a look at these data. Unfortunately, none of the above-
mentioned discrepancies with experimental data is resolved, although the
difference between the first and second spin-triplet states decreases
significantly. Although we see that at least 60 bands are needed to obtain
results converged within 0.05 eV, we will continue to rely on 12 bands to try
to understand the most important discrepancies, while keeping the CPU time to a low level.

We next try to increase the cut-off energy. Again, this is fairly easy. One
can e.g. set up a double dataset loop. The following results are obtained, for
ecut = 25, 30, 35 and 45 Ha:

    Singlet 1 :  9.40   9.37   9.36   9.36 
    Singlet 2 :  9.80   9.78   9.77   9.77 
    Singlet 3 : 10.03  10.03  10.04  10.04 
    Singlet 4 : 10.46  10.37  10.32  10.30
    Singlet 5 : 10.97  10.98  10.98  10.98
    Triplet 1 :  7.75   7.72   7.72   7.72 
    Triplet 2 :  8.09   8.06   8.06   8.06
    Triplet 3 :  8.99   8.97   8.96   8.96
    Triplet 4 :  9.35   9.35   9.35   9.35
    Triplet 5 :  9.80   9.78   9.77   9.77  

You might try to obtain one of these...

The computation with ecut=30 takes a couple of seconds
and gives a result likely converged within 0.01 eV. 
The modifications with respect to the results with ecut=25 Ha are quite small.

We finally examine the effect of the cell size. Again, this is fairly easy.
e keep 12 bands, but stick to ecut=30 Ha.
One can e.g. set up a double dataset loop. The following results are obtained,
for acell = (6 5 5), (8 7 7), (10 9 9), (12 11 11), (14 13 13), (16 15 15) and (20 19 19):

    Singlet 1 :  9.37   9.25   9.24   9.24   9.25   9.25   9.25
    Singlet 2 :  9.78   9.73   9.72   9.72   9.72   9.72   9.72
    Singlet 3 : 10.03  10.04  10.18  10.26  10.29  10.32  10.34
    Singlet 4 : 10.42  10.35  10.34  10.34  10.34  10.34  10.35
    Singlet 5 : 10.98  11.34  11.40  11.12  10.95  10.84  10.70
    Triplet 1 :  7.72   7.60   7.60   7.60   7.60   7.60   7.60
    Triplet 2 :  8.06   8.07   8.07   8.08   8.08   8.08   8.08
    Triplet 3 :  8.97   8.94   8.94   8.94   8.94   8.94   8.94
    Triplet 4 :  9.35   9.73   9.72   9.72   9.72   9.72   9.72
    Triplet 5 :  9.78   9.75  10.00  10.12  10.18  10.22  10.27

Obviously, the cell size plays an important role in the spurious appearance of
the states, that was remarked when comparing against experimental data.
Indeed, the singlet 3 and triplet 4 states energy increases strongly with the
cell size, while all other states quickly stabilize (except the still higher singlet 5 state).
Actually, the singlet 3 and 4 switch between (16 15 15) and (20 19 19),
and the triplet 4 and 5 switch between (6 5 5) and (8 7 7)  ...
The presence of the singlet 3 state (Sigma_g) in the three lowest ones,
an the triplet 4 state (Sigma_g) were two of our discrepancies.

There is one lesson to be learned from that convergence study: the
convergence of different states can be quite different. Usually, converging
the lower excited states do not require too much effort, while it is quite
difficult, especially concerning the supercell size, to converge higher states.

At this stage, we will simply stop this convergence study, and give the
results of an ABINIT calculation using ecut 30 Hartree, acell 10 9 9, and 100
bands, focusing only on those states already converged, then compare the results with other
LDA/TDLDA results (from [[cite:Casida1998]]) and experimental results:

                       present  Casida experimental
    Singlet Pi_g      :  9.20    9.05     9.31
    Singlet Sigma_u-  :  9.72    9.63     9.92
    Singlet Delta_u   : 10.33   10.22    10.27
    Triplet Sigma_u+  :  7.99    7.85     7.75
    Triplet Pi_g      :  7.59    7.54     8.04
    Triplet Delta_u   :  8.92    8.82     8.88
    Triplet Sigma_u-  :  9.72    9.63     9.67

Our calculation is based on pseudopotentials, while Casida's calculation is an
all-electron one. This fact might account for the discrepancy
between both calculations (maximal 0.15 eV - it is of course the user's responsibility to test
the influence of different pseudopotentials on his/her calculations). The
agreement with experimental data is on the order of 0.2 eV, with the exception
of the $^3\Pi_g$ state (0.45 eV). In particular, we note that LDA/TDLDA is
not able to get the correct ordering of the lowest two triplet states ... One
of our problems was intrinsic to the LDA/TDLDA approximation ...

## The choice of the exchange-correlation potential and kernel

As emphasized in [[cite:Casida1998]], choosing a different functional for the self-consistent part
(XC potential) and the generation of the coupling matrix (XC kernel) 
can give a better description of the higher-lying states. 
Indeed, a potential with a -1/r tail (unlike the LDA or GGA) like the van Leeuwen-Baerends one [[cite:VanLeeuwen1994]],
can reproduce fairly well the ionisation energy, giving a much
better description of the Rydberg states. Still, the LDA kernel works pretty well.

In order to activate this procedure, set the value of [[ixc]] in dataset 1 to the
SCF functional, and the value of ixc in dataset 2 to the XC functional to be
used for the kernel. Use pseudopotentials that agree with the SCF functional in dataset 1.
As of writing (end of 2020), the ABINIT implementation has a strong restriction on the XC kernels that can be used. They
must be of LDA-type. See the list of allowed [[ixc]] values in the description of [[iscf]]=-1.