A TDDFT study of carbon monoxide¶
Introduction¶
To illustrate the TDDFT capabilities of DIRAC we will consider excited states of the carbon monoxide molecule.
The experimental spectrum can be found at this page .
Note that this source gives only adiabatic excitation energies download
):
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
Ground state electronic structure¶
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 (see *DFT) using the PBE functional and the eXact 2-Component (.X2C) Hamiltonian
followed by Mulliken population analysis using the input file SCF.inp (download
):
**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
Excited states¶
Preparing the input files¶
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
is totally symmetric (
but for our purposes it will be more convenient to form the combinations:
which transform as rotations
As an example of how to determine total symmetry of the final states let us consider single excitations
HOMO into the LUMO (
Configuration |
States |
---|---|
We can now set up the following correlation of states:
State |
Spin |
Spatial |
Spin |
---|---|---|---|
(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
Irrep |
Valence-excited state |
|
---|---|---|
1 |
||
2 |
||
3 |
||
4 |
Note that the order of the irreps follows the output of DIRAC.
We now set up the following menu file TDDFT.inp (download
) 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
.INTENS
0
*END OF
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 oscillator strengths within the electric dipole approximation. 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.
Looking at the output¶
Energy levels¶
After running the calculation, let us now look at the output. There is output for each symmetry, but we go directly to the final output at the end where we first find a list of levels
Level Rel eigenvalue Abs eigenvalue Total Energy Degeneracy
0 0.0000000000 0.000000000000 -113.360751212612 ( 1 * )
1 0.2179684486 0.217968448567 -113.142782764044 ( 6 * )
2 0.2920913701 0.292091370135 -113.068659842477 ( 3 * )
3 0.3118503049 0.311850304930 -113.048900907681 ( 2 * )
4 0.3196082449 0.319608244945 -113.041142967667 ( 3 * )
5 0.3196082723 0.319608272266 -113.041142940346 ( 3 * )
6 0.3571903374 0.357190337380 -113.003560875232 ( 4 * )
7 0.3713720977 0.371372097692 -112.989379114920 ( 1 * )
8 0.3713721204 0.371372120381 -112.989379092230 ( 1 * )
9 0.3772336605 0.377233660477 -112.983517552135 ( 3 * )
10 0.4005688235 0.400568823512 -112.960182389100 ( 1 * )
11 0.4055211550 0.405521155029 -112.955230057582 ( 3 * )
12 0.4148565628 0.414856562810 -112.945894649802 ( 1 * )
13 0.4154753129 0.415475312869 -112.945275899743 ( 6 * )
14 0.4211386643 0.421138664269 -112.939612548342 ( 2 * )
15 0.4421501163 0.442150116274 -112.918601096338 ( 1 * )
The degeneracies indicated in parenthesis for each level generally arise from both spatial and spin degrees of freedom.
DIRAC assumes degeneracy when energies are within
To learn more about the nature of the excitated states we look at the following list of energy levels, now in other units
(eV and cm
Symmetry Classification in the Abelian subgroup
Level eigenvalue (eV) Eigenvalue (cm-1) 1A1 |1B1 |1B2 |1A2 |3A1 |3B1 |3B2 |3A2 |
0 0.000000000 0.000000 1| 0| 0| 0| 0| 0| 0| 0|
1 5.931223636 47838.544869 0| 0| 0| 0| 0| 3| 3| 0|
2 7.948211082 64106.645746 0| 0| 0| 0| 3| 0| 0| 0|
3 8.485879088 68443.230673 0| 1| 1| 0| 0| 0| 0| 0|
4 8.696983390 70145.901697 0| 0| 0| 0| 0| 0| 0| 3|
5 8.696984133 70145.907693 0| 0| 0| 0| 3| 0| 0| 0|
6 9.719644222 78394.217575 0| 0| 0| 1| 0| 0| 0| 3|
7 10.105549579 81506.754190 0| 0| 0| 1| 0| 0| 0| 0|
8 10.105550196 81506.759170 1| 0| 0| 0| 0| 0| 0| 0|
9 10.265050828 82793.218521 0| 0| 0| 0| 3| 0| 0| 0|
10 10.900032961 87914.694822 1| 0| 0| 0| 0| 0| 0| 0|
11 11.034792767 89001.605956 0| 0| 0| 0| 3| 0| 0| 0|
12 11.288822154 91050.491136 1| 0| 0| 0| 0| 0| 0| 0|
13 11.305659200 91186.291077 0| 0| 0| 0| 0| 3| 3| 0|
14 11.459766843 92429.253037 0| 1| 1| 0| 0| 0| 0| 0|
15 12.031517578 97040.733717 0| 0| 0| 0| 0| 0| 1| 0|
The first level is simply the ground state. The reference state for calculation excitation energies is presently limited
to closed-shell and so this is
Below we give a summary of results obtained with a variety of functionals available in DIRAC. The assignment and energy of states may be compared to [Tozer_JCP1998].
CAMB3LYP |
LDA |
PBE |
PBE0 |
BLYP |
B3LYP |
|||
---|---|---|---|---|---|---|---|---|
6.04 |
6.32 |
5.93 |
5.98 |
5.77 |
5.76 |
5.85 |
5.88 |
|
6.92 |
8.51 |
7.95 |
8.43 |
8.12 |
7.85 |
8.10 |
7.93 |
|
8.07 |
8.51 |
8.49 |
8.19 |
8.26 |
8.45 |
8.25 |
8.41 |
|
7.58 |
9.36 |
8.70 |
9.21 |
8.78 |
8.64 |
8.72 |
8.66 |
|
7.96 |
9.88 |
9.72 |
9.89 |
9.87 |
9.79 |
9.79 |
9.73 |
|
8.07 |
9.88 |
9.72 |
9.89 |
9.87 |
9.79 |
9.79 |
9.73 |
|
8.17 |
10.23 |
10.11 |
10.36 |
10.21 |
10.21 |
10.03 |
10.05 |
|
10.39 |
10.4 |
10.27 |
9.59 |
9.35 |
10.12 |
9.29 |
9.95 |
|
10.78 |
10.78 |
10.90 |
9.95 |
9.78 |
10.72 |
9.65 |
10.43 |
|
11.28 |
11.3 |
10.27 |
10.37 |
10.09 |
10.84 |
10.09 |
10.71 |
|
11.40 |
11.4 |
11.29 |
10.68 |
10.56 |
11.25 |
10.49 |
11.05 |
|
11.52 |
11.53 |
11.46 |
10.70 |
10.58 |
11.36 |
10.49 |
11.13 |
|
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 |
SAOP |
LDA |
PBE |
PBE0 |
BLYP |
B3LYP |
||
---|---|---|---|---|---|---|---|---|
6.04 |
6.32 |
6.33 |
6.11 |
5.76 |
5.76 |
5.84 |
5.87 |
|
6.92 |
8.51 |
8.65 |
8.47 |
7.88 |
7.88 |
8.19 |
8.00 |
|
8.07 |
8.51 |
8.60 |
8.47 |
8.57 |
8.57 |
8.61 |
8.64 |
|
7.58 |
9.36 |
9.39 |
9.30 |
8.64 |
8.64 |
8.79 |
8.71 |
|
7.96 |
9.88 |
10.04 |
10.00 |
9.77 |
9.77 |
9.84 |
9.76 |
|
8.07 |
9.88 |
10.04 |
10.00 |
9.77 |
9.77 |
9.84 |
9.76 |
|
8.17 |
10.23 |
10.50 |
10.49 |
10.27 |
10.27 |
10.39 |
10.31 |
|
10.39 |
10.4 |
10.48 |
10.87 |
10.49 |
10.49 |
10.21 |
10.52 |
|
10.78 |
10.78 |
10.97 |
11.25 |
11.21 |
11.21 |
10.83 |
11.22 |
|
11.28 |
11.3 |
11.42 |
11.70 |
11.34 |
11.34 |
11.20 |
11.43 |
|
11.40 |
11.4 |
11.73 |
11.97 |
11.80 |
11.80 |
|||
11.52 |
11.53 |
11.89 |
11.97 |
11.97 |
||||
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 |
Orbital contributions¶
We can also get information about dominant (single) excitations leading to each excited state. We then go to the section starting with:
Analysis of response solution vectors
-------------------------------------
For the first excitation of
** 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%
12(E1 ) 10.80%
Excitation amplitudes larger than threshold 0.34E-01
7(i:E1 ) ---> 8(v:E1 ) 0.66634917E+00
7(i:E1 ) ---> 12(v:E1 ) 0.23316913E+00
7(i:E1 ) ---> 15(v:E1 ) 0.55294831E-01
7(i:E1 ) ---> 24(v:E1 ) 0.41863343E-01
4(i:E1 ) ---> 8(v:E1 ) 0.34734025E-01
This is an excitation to the first
* Electronic eigenvalue no. 7: -0.4549728114484 (Occupation : f = 1.0000) sym= A1
==========================================================================================
* Gross populations greater than 0.00010
Gross Total | L A1 C s L A1 C pz L A1 C dxx L A1 C dyy L A1 C dzz L A1 O s L A1 O pz
--------------------------------------------------------------------------------------------------------------------------------
alpha 1.0000 | 0.5192 0.3567 0.0081 0.0081 -0.0122 0.0290 0.0872
beta 0.0000 | 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
showing that this is the HOMO, which is a
* Electronic eigenvalue no. 8: 0.116603897E-01 (Occupation : f = 0.0000) sym= B1
==========================================================================================
* Gross populations greater than 0.00010
Gross Total | L B1 C px L B1 C dxz L B1 O px L B1 O dxz
-----------------------------------------------------------------------------------
alpha 0.0000 | 0.0000 0.0000 0.0000 0.0000
beta 1.0000 | 0.9302 0.0039 0.0622 0.0037
which is clearly a
Oscillator strenghts¶
DIRAC also allows provides oscillator strenghts corresponding to the calculated (single) excitations. The general definition of an oscillator strenght is
where
.INTENS
0
which means that we ask for oscillator strengths to order zero in the wave vector. This corresponds to the electric dipole approximation. Starting from DIRAC21, oscillator strengths have been implemented to arbitrary order in the wave vector, in both the length and velocity representation, see [List_JCP2020]. In addition DIRAC also allows the use of the full semi-classical light-matter interaction. Presently, we limit ourselves to the electric-dipole approximation. The isotropic dipole length-dipole length oscillator strenghts are given as
and are found to be
* Isotropic oscillator strengths (generalized length gauge) above threshold :1.00E-08
- [f( 0) corresponds to the electric dipole approximation]
================================================================
Level Frequency (eV) Symmetry f(total) f( 0)
----------------------------------------------------------------
10 8.48588 B1 1B1 8.693603E-02 8.693603E-02
11 8.48588 B2 1B2 8.693603E-02 8.693603E-02
27 10.90003 A1 1A1 2.722369E-03 2.722369E-03
31 11.28882 A1 1A1 2.000788E-01 2.000788E-01
38 11.45977 B2 1B2 5.287744E-02 5.287744E-02
39 11.45977 B1 1B1 5.287744E-02 5.287744E-02
----------------------------------------------------------------
Sum of oscillator strengths (general length) : 4.82428E-01
Corresponding rates of spontaneous emission are given by
This is a spin-orbit free calculation and so we only see oscillator strength to
singlet states (clearly dominated by
Energy eigenvalues in atomic units
Level Rel eigenvalue Abs eigenvalue Total Energy Degeneracy
0 0.0000000000 0.000000000000 -113.360753486552 ( 1 * )
1 0.2177547139 0.217754713945 -113.142998772607 ( 1 * )
2 0.2177548426 0.217754842630 -113.142998643923 ( 1 * )
3 0.2179681827 0.217968182656 -113.142785303897 ( 2 * )
4 0.2181824349 0.218182434936 -113.142571051617 ( 2 * )
5 0.2920874557 0.292087455726 -113.068666030827 ( 2 * )
6 0.2920912039 0.292091203926 -113.068662282626 ( 1 * )
7 0.3118487168 0.311848716812 -113.048904769740 ( 2 * )
8 0.3194494091 0.319449409064 -113.041304077488 ( 2 * )
9 0.3196035534 0.319603553362 -113.041149933191 ( 1 * )
10 0.3196035760 0.319603576039 -113.041149910513 ( 1 * )
11 0.3197697403 0.319769740262 -113.040983746290 ( 2 * )
12 0.3571899414 0.357189941378 -113.003563545175 ( 1 * )
13 0.3571902624 0.357190262395 -113.003563224157 ( 1 * )
14 0.3571939414 0.357193941428 -113.003559545124 ( 2 * )
15 0.3713761759 0.371376175880 -112.989377310672 ( 1 * )
16 0.3713761940 0.371376193965 -112.989377292587 ( 1 * )
17 0.3772335791 0.377233579067 -112.983519907486 ( 1 * )
18 0.3772335955 0.377233595510 -112.983519891042 ( 2 * )
19 0.4005686571 0.400568657147 -112.960184829405 ( 1 * )
20 0.4055210630 0.405521062959 -112.955232423594 ( 1 * )
21 0.4055210674 0.405521067358 -112.955232419194 ( 2 * )
22 0.4148563889 0.414856388866 -112.945897097687 ( 1 * )
23 0.4154637983 0.415463798333 -112.945289688219 ( 1 * )
24 0.4154638705 0.415463870460 -112.945289616092 ( 1 * )
25 0.4154752272 0.415475227205 -112.945278259348 ( 2 * )
26 0.4154867237 0.415486723685 -112.945266762867 ( 2 * )
27 0.4211386126 0.421138612595 -112.939614873957 ( 2 * )
28 0.4420085576 0.442008557560 -112.918744928993 ( 1 * )
Total average: -113.0232469843
Relative real eigenvalues in other units;
Symmetry Classification in the Abelian subgroup
Level eigenvalue (eV) Eigenvalue (cm-1) 0+ | 0- | 1 | 2 | 3 |
0 0.000000000 0.000000 1| 0| 0| 0| 0|
1 5.925407621 47791.635542 0| 1| 0| 0| 0|
2 5.925411122 47791.663785 1| 0| 0| 0| 0|
3 5.931216400 47838.486508 0| 0| 2| 0| 0|
4 5.937046502 47885.509448 0| 0| 0| 2| 0|
5 7.948104566 64105.786632 0| 0| 2| 0| 0|
6 7.948206559 64106.609267 0| 1| 0| 0| 0|
7 8.485835873 68442.882122 0| 0| 2| 0| 0|
8 8.692661245 70111.041251 0| 0| 0| 0| 2|
9 8.696855725 70144.872014 0| 0| 0| 1| 0|
10 8.696856343 70144.876991 0| 0| 0| 1| 0|
11 8.701377901 70181.345822 0| 0| 2| 0| 0|
12 9.719633446 78394.130663 1| 0| 0| 0| 0|
13 9.719642181 78394.201118 0| 1| 0| 0| 0|
14 9.719742293 78395.008572 0| 0| 2| 0| 0|
15 10.105660552 81507.649249 0| 0| 0| 1| 0|
16 10.105661044 81507.653218 0| 0| 0| 1| 0|
17 10.265048612 82793.200653 0| 1| 0| 0| 0|
18 10.265049060 82793.204262 0| 0| 2| 0| 0|
19 10.900028434 87914.658309 1| 0| 0| 0| 0|
20 11.034790262 89001.585749 0| 1| 0| 0| 0|
21 11.034790381 89001.586714 0| 0| 2| 0| 0|
22 11.288817420 91050.452960 1| 0| 0| 0| 0|
23 11.305345874 91183.763928 0| 1| 0| 0| 0|
24 11.305347837 91183.779759 1| 0| 0| 0| 0|
25 11.305656869 91186.272276 0| 0| 2| 0| 0|
26 11.305969704 91188.795462 0| 0| 0| 2| 0|
27 11.459765436 92429.241696 0| 0| 2| 0| 0|
28 12.027665569 97009.665171 1| 0| 0| 0| 0|
and that oscillator strength goes to more states
* Isotropic oscillator strengths (generalized length gauge) above threshold :1.00E-08
- [f( 0) corresponds to the electric dipole approximation]
================================================================
Level Frequency (eV) Symmetry f(total) f( 0)
----------------------------------------------------------------
3 5.93122 B2 1 2.425375E-07 2.425375E-07
4 5.93122 B1 1 2.425375E-07 2.425375E-07
7 7.94810 B2 1 1.611899E-06 1.611899E-06
8 7.94810 B1 1 1.611899E-06 1.611899E-06
10 8.48584 B1 1 8.690290E-02 8.690290E-02
11 8.48584 B2 1 8.690290E-02 8.690290E-02
16 8.70138 B2 1 3.112435E-05 3.112435E-05
17 8.70138 B1 1 3.112435E-05 3.112435E-05
18 9.71963 A1 0+ 8.239305E-07 8.239305E-07
20 9.71974 B2 1 1.270263E-06 1.270263E-06
21 9.71974 B1 1 1.270263E-06 1.270263E-06
25 10.26505 B2 1 1.962960E-08 1.962960E-08
26 10.26505 B1 1 1.962960E-08 1.962960E-08
27 10.90003 A1 0+ 2.722523E-03 2.722523E-03
29 11.03479 B2 1 1.125638E-07 1.125638E-07
30 11.03479 B1 1 1.125638E-07 1.125638E-07
31 11.28882 A1 0+ 2.000395E-01 2.000395E-01
33 11.30535 A1 0+ 4.016589E-05 4.016589E-05
34 11.30566 B2 1 3.321212E-07 3.321212E-07
35 11.30566 B1 1 3.321212E-07 3.321212E-07
38 11.45977 B2 1 5.287687E-02 5.287687E-02
39 11.45977 B1 1 5.287687E-02 5.287687E-02
40 12.02767 A1 0+ 4.211753E-06 4.211753E-06
----------------------------------------------------------------
Sum of oscillator strengths (general length) : 4.82436E-01
reflecting that triplet states “steal” intensity from singlet ones.