To illustrate the TDDFT capabilities of DIRAC we will consider excited states of the carbon monoxide molecule. The experimental spectrum can be found here . Note that this source gives only adiabatic excitation energies \(T_e\), that is, the energy difference between the minima of the potential curve of the ground and excited state, whereas the TDDFT calculations give vertical excitation energies \(T_v\). The experimental adiabatic excitation energies has been converted to vertical ones by [Nielsen_JCP1980]. We will begin this tutorial in spinfree mode to faciliate the connection to the non-relativistic domain. In the following we employ the molecular input file CO.mol

```
INTGRL
Carbon monoxide
Dunning aug-cc-pVDZ
C 2 A
6. 1
C 0.0000000000000 0.0000000000000 0.0000000000000
LARGE BASIS aug-cc-pVDZ
8. 1
O 0.0000000000000 0.0000000000000 1.128323
LARGE BASIS aug-cc-pVDZ
FINISH
```

Here we do not provide any symmetry information, meaning that we ask DIRAC to detect it. DIRAC will find that the full group is \(C_{\infty v}\). With spin-orbit coupling DIRAC will then activate linear supersymmetry, but in the spin-orbit free case it will simply use the highest Abelian single point group, that is \(C_{2v}\). Irreps of \(C_{\infty v}\) and \(C_{2v}\) correlate as follows:

\(C_{\infty v}\) | \(C_{2v}\) |
---|---|

\(\Sigma^{+}\) | \(A1\) |

\(\Sigma^{-}\) | \(A2\) |

\(\Pi\) | \(B_1, B_2\) |

\(\Delta\) | \(A_1, A_2\) |

Before actually doing TDDFT let us first have a look at the electronic structure of the ground state of CO. We therefore carry out Kohn-Sham DFT using the PBE functional and the eXact 2-Component (X2C) Hamiltonian followed by Mulliken population analysis using the input file PBE.inp

```
**DIRAC
.TITLE
CO
.WAVE FUNCTIONS
.ANALYZE
**HAMILTONIAN
.SPINFREE
.X2C
.DFTAUTO ! uses xcfun
PBE
**INTEGRALS
*TWOINT
.SCREEN
1.0D-12
*READIN
.UNCONT
**WAVE FUNCTIONS
.SCF
*SCF
.CLOSED SHELL
14
**ANALYZE
.MULPOP
*MULPOP
.VECPOP
1..10
*END OF
```

Based on orbital energies and the Mulliken population analysis we can set up the following MO diagram:

(Orbital energies for the atoms were obtained from atomic calculations using fractional occupation)

According to the above MO diagram the valence electron configuration of carbon monoxide is
\(1\sigma^{2}2\sigma^{2}1\pi^{4}3\sigma^{2}\). From projection analysis we find that the electron configurations
of the carbon and oxygen atoms in the molecule are \(1s^{2.0}2s^{2.0}2p^{2.0}\) and \(1s^{2.0}2s^{1.7}2p^{4.3}\) , respectively,
with a slight negative charge of -0.3 on oxygen. We note that the HOMO is a \(\sigma\) orbital dominated (~90 %) by carbon (2s,2p),
whereas the doubly degenerate LUMO is a \(\pi\) orbital with about 75% C *2p* and 25% O *2p*. The LUMO+1 is \(\sigma\) orbital
with a combination of carbon and oxygen *3s* dominated by the former atom.

The excited states will be calculated in a “bottoms-up” fashion for each irrep by expansion of the solution vector in trial vector [Bast2009].
For each irrep the user specifies the number of desired excitations. It may be a good idea of including some extra excitations in case of
root flipping during the iterative solution of the TDDFT equations. We are working within the adiabatic approximation of TDDFT
and consider single excitations from a closed-shell ground state. The final states will be therefore be singlet or triplets in the spin-orbit
free case. We need to tell DIRAC how these states are distributed on the four irreps of \(C_{2v}\).
At this point we should emphasize that DIRAC will employ the *total* symmetry of the final states,
that is the combined spin and spatial symmetry. A singlet spin function is totally symmetric (\(A_1\)), whereas the
triplet spin functions transform as rotations. We know the triplet functions as:

\[T_{-1} = \alpha_1\alpha_2;\quad T_0=\frac{1}{\sqrt{2}}\left(\alpha_1\beta_2-\beta_1\alpha_2\right);\quad T_{+1}=\beta_1\beta_2\]

but for our purposes it will be more convenient to form the combinations:

\[T_x=\frac{1}{\sqrt{2}}\left(T_{-1}-T_{+1}\right);\quad T_y=\frac{i}{\sqrt{2}}\left(T_{-1}+T_{+1}\right);\quad T_z=T_0\]

which transform as rotations \(R_x\left(B_{2}\right)\) , \(R_y\left(B_{1}\right)\) and \(R_z\left(A_{2}\right)\) , respectively. For a given point group you find the symmetry of the rotations in most standard character tables (for instance here ).

As an example of how to determine total symmetry of the final states let us consider single excitations HOMO into the LUMO (\(\pi\)) and LUMO+1 (\(\sigma\)) orbitals. This leads to \(2 \times\ 6 = 12\) determinants which translate into the following final states:

Configuration | States |
---|---|

\(3\sigma^{-1}2\pi^{1}\) | \(^{1,3}\Pi\) |

\(3\sigma^{-1}4\sigma^{1}\) | \(^{1,3}\Sigma\) |

We can now set up the following correlation of states:

State | Spin | Spatial | Spin \(\otimes\) Spatial |
---|---|---|---|

\(^1\Sigma\) | \(A_1\) | \(A_1\) | \(A_1\) |

\(^1\Pi_{x,y}\) | \(A_1\) | \(B_{1}, B_{2}\) | \(B_{1}, B_{2}\) |

\(^3\Sigma\) | \(B_{2}, B_{1}, A_{2}\) | \(A_1\) | \(B_{2}, B_{1}, A_{2}\) |

\(^3\Pi_{x,y}\) | \(B_{2}, B_{1}, A_{2}\) | \(B_{1}, B_{2}\) | \(\left(A_2, A_{1}, B_2\right), \left(A_1, A_2, B_2\right)\) |

(The group multiplication table is found in the DIRAC output.)

Counting total symmetries we find that the 12 microstates are evenly distributed amongst the four irreps of \(C_{2v}\):

Irrep | Valence-excited state | |
---|---|---|

1 | \(A_1\) | \(^1\Sigma, ^{3(y)}\Pi_{x}, ^{3(x)}\Pi_{y}\) |

2 | \(B_1\) | \(^{1}\Pi_{x}, ^{3(y)}\Sigma, ^{3(z)}\Pi_{y}\) |

3 | \(B_2\) | \(^{1}\Pi_{y}, ^{3(x)}\Sigma, ^{3(z)}\Pi_{x}\) |

4 | \(A_2\) | \(^{3(z)}\Sigma, ^{3(x)}\Pi_{x}, ^{3(y)}\Pi_{y}\) |

Note that the order of the irreps follows the output of DIRAC.

We now set up the following menu file for our calculation:

```
**DIRAC
.TITLE
CO
.WAVE FUNCTIONS
.ANALYZE
.PROPERTIES
**HAMILTONIAN
.SPINFREE
.X2C
.DFT
CAMB3LYP
**INTEGRALS
*READIN
.UNCONT
**WAVE FUNCTIONS
.SCF
*SCF
.CLOSED SHELL
14
**ANALYZE
.MULPOP
*MULPOP
.VECPOP
1..10
**PROPERTIES
*EXCITATION
#A1 :
.EXCITA
1 10
#B1 :
.EXCITA
2 10
#B2 :
#.EXCITA
#3 10
#A2 :
.EXCITA
4 10
.ANALYZE
.OPERATOR
XDIPLEN
.OPERATOR
YDIPLEN
.OPERATOR
ZDIPLEN
*END OF
```

Note that we skip excitations in \(B_{2}\) since they are related by symmetry to the excitations of \(B_{1}\). On the other hand we ask for 10 excitations in the other symmetries. We base our TDDFT calculation on the long-range corrected CAMB3LYP functional [Yanai_CPL2004] which provides a better description of Rydberg and charge-transfer excitations than standard continuum functionals. We furthermore ask for transition moments to be calculated with respect to the component of the dipole moment operator. These will be non-zero only for excitations in irreps \(B_{1}, B_{2}\) and \(A_{1}\). Finally we ask for analysis of what orbitals contribute to the various excitations. For this the Mulliken population analysis may come in handy as reference.

After running the calculation, let us now look at the output. Let us take the final output of excitations of \(A_1\) symmetry as an example:

```
================================================================================
*** Excitations of boson symmetry 1 : A1
--------------------------------------------------------------------------------
Excitation no. 1 excitation energy 0.21796849 a.u. 2.89E-06 (converged)
Nonrel. sym.: TxB2 5.9312 eV
*** Transition moments <0|A|n> for this excitation:
A - ZDIPLEN A1 T+ 0.00000000D+00 a.u.
--------------------------------------------------------------------------------
Excitation no. 2 excitation energy 0.21796849 a.u. 6.32E-06 (converged)
Nonrel. sym.: TyB1 5.9312 eV
*** Transition moments <0|A|n> for this excitation:
A - ZDIPLEN A1 T+ 0.00000000D+00 a.u.
--------------------------------------------------------------------------------
Excitation no. 3 excitation energy 0.31960825 a.u. 5.31E-06 (converged)
Nonrel. sym.: TzA2 8.6970 eV
*** Transition moments <0|A|n> for this excitation:
A - ZDIPLEN A1 T+ 0.00000000D+00 a.u.
--------------------------------------------------------------------------------
Excitation no. 4 excitation energy 0.35719034 a.u. 9.84E-06 (converged)
Nonrel. sym.: TzA2 9.7196 eV
*** Transition moments <0|A|n> for this excitation:
A - ZDIPLEN A1 T+ 0.00000000D+00 a.u.
--------------------------------------------------------------------------------
Excitation no. 5 excitation energy 0.37137212 a.u. 2.71E-06 (converged)
Nonrel. sym.: S0A1 10.1055 eV
*** Transition moments <0|A|n> for this excitation:
A - ZDIPLEN A1 T+ 0.75578153D-10 a.u.
--------------------------------------------------------------------------------
Excitation no. 6 excitation energy 0.40056879 a.u. 3.42E-06 (converged)
Nonrel. sym.: S0A1 10.9000 eV
*** Transition moments <0|A|n> for this excitation:
A - ZDIPLEN A1 T+ -0.100967888969 a.u.
-> dipole oscillator strength (length) : 0.00272240
-> dipole oscillator strength (velocity) : 0.00000000
Total dipole radiation rate (length) 1.40352E+07 s-1
- corresponding radiation life time 7.12496E-08 s
--------------------------------------------------------------------------------
Excitation no. 7 excitation energy 0.41485654 a.u. 4.09E-06 (converged)
Nonrel. sym.: S0A1 11.2888 eV
*** Transition moments <0|A|n> for this excitation:
A - ZDIPLEN A1 T+ 0.850545002869 a.u.
-> dipole oscillator strength (length) : 0.20007889
-> dipole oscillator strength (velocity) : 0.00000000
Total dipole radiation rate (length) 1.10639E+09 s-1
- corresponding radiation life time 9.03841E-10 s
--------------------------------------------------------------------------------
Excitation no. 8 excitation energy 0.41547529 a.u. 2.39E-06 (converged)
Nonrel. sym.: TxB2 11.3057 eV
*** Transition moments <0|A|n> for this excitation:
A - ZDIPLEN A1 T+ 0.00000000D+00 a.u.
--------------------------------------------------------------------------------
Excitation no. 9 excitation energy 0.41547529 a.u. 3.46E-06 (converged)
Nonrel. sym.: TyB1 11.3057 eV
*** Transition moments <0|A|n> for this excitation:
A - ZDIPLEN A1 T+ 0.00000000D+00 a.u.
--------------------------------------------------------------------------------
Excitation no. 10 excitation energy 0.44215019 a.u. 3.79E-06 (converged)
Nonrel. sym.: TxB2 12.0315 eV
*** Transition moments <0|A|n> for this excitation:
A - ZDIPLEN A1 T+ 0.00000000D+00 a.u.
================================================================================
```

There is sufficient symmetry in the calculation (symmetry distinct rotations) to allow DIRAC to pinpoint the symmetry of the valence-excited states (see the entry Nonrel. sym.). The triplet states have strictly zero transition moment (and thus intensity), whereas the sixth and seventh states are singlet states with non-zero transition moments. We can also look at a snippet from the excitation analysis:

```
** E solution vectors : PP EXCITATIONA1 Irrep: A1 Trev: 1 Length: 441
Freq.: 0.217968 Norm: 0.71160560E+00 Residual norm: 0.29E-05
Dominant inactive orbitals:
7(E1 ) 99.56%
Dominant virtual orbitals:
8(E1 ) 87.92%
13(E1 ) 10.80%
Excitation amplitudes larger than threshold 0.34E-01
7(i:E1 ) ---> 8(v:E1 ) 0.66634909E+00
7(i:E1 ) ---> 13(v:E1 ) 0.23316934E+00
7(i:E1 ) ---> 15(v:E1 ) 0.55294899E-01
7(i:E1 ) ---> 24(v:E1 ) 0.41863354E-01
4(i:E1 ) ---> 8(v:E1 ) 0.34733952E-01
```

We see that for each excitation we get information about the dominant occupied and virtual orbitals. The character of each orbital can be found from Mulliken population analysis or from projection analysis.

In order to properly identify the excited states we collect data from the four irreps and list the states in order of increasing excitation energy. This gives the following table:

Level (eV) | Symmetry | Dominant orbitals | |
---|---|---|---|

1 | 5.9312 | A1(TxB2), A1(TyB1), B1(TzB2), A2(TyB2), A2(TxB1) | \(7\rightarrow (8,9)\) |

2 | 7.9482 | B1(TyA1), A2(TzA1) | \((5,6)\rightarrow (8,9)\) |

3 | 8.4859 | B1(S0B1) | \(7 \rightarrow (8,9)\) |

4 | 8.6970 | A1(TzA2), B1(TxA2), B1(TyA1), A2(TzA1) | \((5,6)\rightarrow (8,9)\) |

5 | 9.7196 | A1(TzA2), B1(TxA2), A2(S0A2) | \((5,6)\rightarrow (8,9)\) |

6 | 10.1055 | A1(S0A1), A2(S0A2) | \((5,6)\rightarrow (8,9)\) |

7 | 10.2650 | B1(TyA1), A2(TzA1) | \(7\rightarrow 10,11\) |

8 | 10.9000 | A1(S0A1) | \(7\rightarrow 10\) |

9 | 11.0348 | B1(TyA1), A2(TzA1) | \(7\rightarrow 11\) |

10 | 11.2888 | A1(S0A1) | \(7\rightarrow 11\) |

11 | 11.3057 | A1 TxB2, TyB1, TzB2 | \(7\rightarrow (12,13)\) |

12 | 11.4598 | B1 S0B1 | \(7\rightarrow (12,13)\) |

13 | 12.0315 | A1 TxB2 | \(4\rightarrow (8,9)\) |

from which we identify 13 energy levels. Keeping in mind that for every excitation of total symmetry \(B_1\) there is a \(B_2\) partner we now consider each energy level:

- The first level has six-fold degeneracy and is a triplet. The space symmetries are \(B_1\) and \(B_2\) and from the correspondence table between \(C_{2v}\) and \(C_{\infty v}\) we identify it as \(^{3}\Pi\). The dominant excitation is \(3\sigma \rightarrow 2\pi\).
- The second level is a triplet of three-fold degeneracy identified as \(^{3}\Sigma^{+}\). The dominant excitation is \(1\pi \rightarrow 2\pi\).
- The third level is a singlet of two-fold degeneray identified as \(^{1}\Pi\). The dominant excitation is \(3\sigma \rightarrow 2\pi\).
- The first level has six-fold degeneracy and is a triplet. The space symmetries are \(A_1\) and \(A_2\) and from the correspondence table between \(C_{2v}\) and \(C_{\infty v}\) we identify it as \(^{3}\Delta\). The dominant excitation is \(1\pi \rightarrow 2\pi\).
- The fifth level consists of a singlet and a triplet, so two separate electronic states, namely \(^{1}\Sigma^{-}\) and \(^{3}\Sigma^{-}\). The dominant excitation is \(1\pi \rightarrow 2\pi\).
- The sixth level is a doubly-degenerate singlet identified as \(^{1}\Delta\). The dominant excitation is \(1\pi \rightarrow 2\pi\).
- \(^{3}\Sigma^{+}\). The dominant excitation is \(3\sigma \rightarrow 4\sigma\). The \(4\sigma\) orbital is dominated by
*3s*from both atoms, but with some*3s*character as well. - \(^{1}\Sigma^{+}\). The dominant excitation is \(3\sigma \rightarrow 4\sigma\).
- \(^{3}\Sigma^{+}\). The dominant excitation is \(3\sigma \rightarrow 5\sigma\). The \(5\sigma\) orbital is dominated by
*3p*from both atoms, but with some*3p*character as well. - \(^{1}\Sigma^{+}\). The dominant excitation is \(3\sigma \rightarrow 5\sigma\).
- \(^{3}\Pi\). The dominant excitation is \(3\sigma \rightarrow 3\pi\), the latter dominated by carbon
*3p*. - \(^{1}\Pi\).The dominant excitation is \(3\sigma \rightarrow 3\pi\).

The assignment and energy of states may be compared to [Tozer_JCP1998]. Below we give a summary of results obtained with a variety of functionals available in DIRAC.

\(T_e\) | \(T_v\) | CAMB3LYP | LDA | PBE | PBE0 | BLYP | B3LYP | |
---|---|---|---|---|---|---|---|---|

\(a\ ^{3}\Pi\) | 6.04 | 6.32 | 5.93 | 5.98 | 5.77 | 5.76 | 5.85 | 5.88 |

\(a\ ^{3}\Sigma^{+}\) | 6.92 | 8.51 | 7.95 | 8.43 | 8.12 | 7.85 | 8.10 | 7.93 |

\(A\ ^{1}\Pi\) | 8.07 | 8.51 | 8.49 | 8.19 | 8.26 | 8.45 | 8.25 | 8.41 |

\(d\ ^{3}\Delta\) | 7.58 | 9.36 | 8.70 | 9.21 | 8.78 | 8.64 | 8.72 | 8.66 |

\(e\ ^{3}\Sigma^{-}\) | 7.96 | 9.88 | 9.72 | 9.89 | 9.87 | 9.79 | 9.79 | 9.73 |

\(I\ ^{1}\Sigma^{-}\) | 8.07 | 9.88 | 9.72 | 9.89 | 9.87 | 9.79 | 9.79 | 9.73 |

\(D\ ^{1}\Delta\) | 8.17 | 10.23 | 10.11 | 10.36 | 10.21 | 10.21 | 10.03 | 10.05 |

\(b\ ^{3}\Sigma^{+}\) | 10.39 | 10.4 | 10.27 | 9.59 | 9.35 | 10.12 | 9.29 | 9.95 |

\(B\ ^{1}\Sigma^{+}\) | 10.78 | 10.78 | 10.90 | 9.95 | 9.78 | 10.72 | 9.65 | 10.43 |

\(j\ ^{3}\Sigma^{+}\) | 11.28 | 11.3 | 10.27 | 10.37 | 10.09 | 10.84 | 10.09 | 10.71 |

\(C\ ^{1}\Sigma^{+}\) | 11.40 | 11.4 | 11.29 | 10.68 | 10.56 | 11.25 | 10.49 | 11.05 |

\(E\ ^{1}\Pi\) | 11.52 | 11.53 | 11.46 | 10.70 | 10.58 | 11.36 | 10.49 | 11.13 |

\(c\ ^{3}\Pi\) | 11.55 | 11.55 | 11.31 | 10.31 | 10.31 | 11.12 | 10.30 | 10.97 |

MAE (all states) | 0.29 | 0.49 | 0.62 | 0.29 | 0.68 | 0.39 | ||

MAE (valence) | 0.30 | 0.15 | 0.26 | 0.31 | 0.31 | 0.33 | ||

MAE (Rydberg) | 0.29 | 0.89 | 1.05 | 0.26 | 1.11 | 0.46 |

For all functionals we provide the mean absolute error (MAE), and it can indeed be seen that the performance of CAMB3LYP is quite good. PBE0 is also display nice performance, whereas its GGA equivalent comes out significantly worse. However, if we do separate error analysis of the valence states (the first seven states) and the Rydberg states (the last six states), we see that PBE performs well for the valence states, but not for the Rydberg states. This is connected to the wrong asymptotic behaviour of LDA/GGA functionals.

The performance of standard LDA/GGA functionals can be improved by invoking asymptotic corrections. One option is the statistical average of orbital potentials (SAOP). Below we show the result of SAOP-corrected functionals.

SAOP-corrected | \(T_e\) | \(T_v\) | SAOP | LDA | PBE | PBE0 | BLYP | B3LYP |
---|---|---|---|---|---|---|---|---|

\(a\ ^{3}\Pi\) | 6.04 | 6.32 | 6.33 | 6.11 | 5.76 | 5.76 | 5.84 | 5.87 |

\(a\ ^{3}\Sigma^{+}\) | 6.92 | 8.51 | 8.65 | 8.47 | 7.88 | 7.88 | 8.19 | 8.00 |

\(A\ ^{1}\Pi\) | 8.07 | 8.51 | 8.60 | 8.47 | 8.57 | 8.57 | 8.61 | 8.64 |

\(d\ ^{3}\Delta\) | 7.58 | 9.36 | 9.39 | 9.30 | 8.64 | 8.64 | 8.79 | 8.71 |

\(e\ ^{3}\Sigma^{-}\) | 7.96 | 9.88 | 10.04 | 10.00 | 9.77 | 9.77 | 9.84 | 9.76 |

\(I\ ^{1}\Sigma^{-}\) | 8.07 | 9.88 | 10.04 | 10.00 | 9.77 | 9.77 | 9.84 | 9.76 |

\(D\ ^{1}\Delta\) | 8.17 | 10.23 | 10.50 | 10.49 | 10.27 | 10.27 | 10.39 | 10.31 |

\(b\ ^{3}\Sigma^{+}\) | 10.39 | 10.4 | 10.48 | 10.87 | 10.49 | 10.49 | 10.21 | 10.52 |

\(B\ ^{1}\Sigma^{+}\) | 10.78 | 10.78 | 10.97 | 11.25 | 11.21 | 11.21 | 10.83 | 11.22 |

\(j\ ^{3}\Sigma^{+}\) | 11.28 | 11.3 | 11.42 | 11.70 | 11.34 | 11.34 | 11.20 | 11.43 |

\(C\ ^{1}\Sigma^{+}\) | 11.40 | 11.4 | 11.73 | 11.97 | 11.80 | 11.80 | ||

\(E\ ^{1}\Pi\) | 11.52 | 11.53 | 11.89 | 11.97 | 11.97 | |||

\(c\ ^{3}\Pi\) | 11.55 | 11.55 | 11.76 | 11.73 | 11.68 | 11.68 | 11.42 | 11.70 |

MAE (all states) | 0.17 | 0.24 | 0.29 | 0.29 | 0.20 | 0.26 | ||

MAE (valence) | 0.12 | 0.12 | 0.32 | 0.32 | 0.24 | 0.29 | ||

MAE (Rydberg) | 0.21 | 0.42 | 0.25 | 0.25 | 0.12 | 0.21 |