Excitation energies at the SCF-level

In this tutorial we will look at the theory behind the calculation of excitation energies at the SCF-level, that is, either by time-dependent Hartree-Fock (TD-HF) or TD-DFT. The implementations in DIRAC are presented in [Bast2009].

Excitation energies \(\omega_m\) are found by solving the generalized eigenvalue problem

\[\left(E_0^{[2]}-\hbar\lambda S^{[2]}\right)X_m=0;\qquad (1)\]

where the matrices \(E_0^{[2]}\) and \(S^{[2]}\) are the electronic Hessian and generalized metric, respectively. From their structure, it can be shown that solution vectors come in pairs

\[\begin{split}\lambda_{+;n}=+\left|\omega_{n}\right|,\,\boldsymbol{X}_{+;n}=\left[\begin{array}{c} \boldsymbol{Z}_{n}\\ \boldsymbol{Y}_{n}^{\ast}\end{array}\right],\quad \lambda_{-;n}=-\left|\omega_{n}\right|,\,\boldsymbol{X}_{-;n}=\left[\begin{array}{c} \boldsymbol{Y}_{n}\\\boldsymbol{Z}_{n}^{\ast}\end{array}\right].\end{split}\]

Transition moments associated with operator \(\hat{h}_{A}\) for transitions \(n\leftarrow0\) are generally calculated as

\[\langle n|\hat{h}_{A}|0\rangle=\boldsymbol{X}_{+;n}^{\dagger}\boldsymbol{E}_{A}^{[1]};\qquad (2)\]

where appears solution vector \(\boldsymbol{X}_{+;n}^{\dagger}\) and property gradient \(\boldsymbol{E}_{A}^{[1]}\). The property gradient is given by

\[\begin{split}\boldsymbol{E}_{B}^{[1]}=\left[\begin{array}{c}\boldsymbol{g}_{B}\\\Theta_{hB}\boldsymbol{g}_{B}^{\ast}\end{array}\right];\quad g_{B;ai}=-h_{B;ai}.\end{split}\]

An important generalization above is that, in addition to Hermitian operators \(\hat{h}_{B}\) (\(\Theta_{hB}=+1\)), imposed by the tenets of quantum mechanics, we also allow anti-Hermitian ones (\(\Theta_{hB}=-1\)).

A key to computational efficiency is to consider a decomposition of solution vectors in terms of components of well-defined Hermiticity and time reversal symmetry. Using a pair of solution vectors \(\boldsymbol{X}_{+}\) and \(\boldsymbol{X}_{-}\), we may form Hermitian and anti-Hermitian combinations

\[\begin{split}\begin{array}{lcl} \boldsymbol{X}_{h} & = &\frac{1}{2}\left(\boldsymbol{X}_{+}+\boldsymbol{X}_{-}\right)=\left[\begin{array}{c} \boldsymbol{Z}+\boldsymbol{Y}\\ \boldsymbol{Y}^{\ast}+\boldsymbol{Z}^{\ast} \end{array}\right]=\left[\begin{array}{c} \boldsymbol{h}\\ \boldsymbol{h}^{\ast} \end{array}\right]\\ \boldsymbol{X}_{a} & = & \frac{1}{2}\left(\boldsymbol{X}_{+}-\boldsymbol{X}_{-}\right)=\left[\begin{array}{c} \boldsymbol{Z}-\boldsymbol{Y}\\ \boldsymbol{Y}^{\ast}-\boldsymbol{Z}^{\ast} \end{array}\right]=\left[\begin{array}{c} \boldsymbol{a}\\ -\boldsymbol{a}^{\ast} \end{array}\right]. \end{array}\end{split}\]

The inverse relations therefore provide a separation of solution vectors into Hermitian and anti-Hermitian contributions

\[\boldsymbol{X}_{+}=\boldsymbol{X}_{h}+\boldsymbol{X}_{a};\quad\boldsymbol{X}_{-}=\boldsymbol{X}_{h}-\boldsymbol{X}_{a}.\]

Further decomposition of each contribution into time-symmetric and time-antisymmetric parts give vectors that are well-defined with respect to both Hermiticity and time reversal symmetry

\[\boldsymbol{U}^{\dagger}\left(\Theta_{h},\Theta_{t}\right)=\left[\begin{array}{cccccccc} \boldsymbol{c}^{\dagger} & \boldsymbol{d}^{\dagger} & \Theta_{t}\boldsymbol{c}^{T} & \Theta_{t}\boldsymbol{d}^{T} & \Theta_{h}\boldsymbol{c}^{T} & \Theta_{h}\boldsymbol{d}^{T} & \Theta_{h}\Theta_{t}\boldsymbol{c}^{\dagger} & \Theta_{h}\Theta_{t}\boldsymbol{d}^{\dagger}\end{array}\right]\]

where

\[c_{ai} = x_{ai}\quad\mbox{and}\quad d_{ai} = x_{a\overline{i}}\]

and the index overbar refers to a Kramers’ partner in a Kramers-restricted orbital set. The scalar product of such vectors is given by [Saue2002a]

\[\boldsymbol{U}_{1}^{\dagger}\left(\Theta_{h1},\Theta_{t1}\right)\boldsymbol{U}_{2}\left(\Theta_{h2},\Theta_{t2}\right)= \left(1+\Theta_{h1}\Theta_{h2}\Theta_{t1}\Theta_{t2}\right)\left[z+\Theta_{h1}\Theta_{h2}z^{\ast}\right];\quad z=\left(\boldsymbol{c}_{1}^{\dagger}\boldsymbol{c}_{2}+\boldsymbol{d}_{1}^{\dagger}\boldsymbol{d}_{2}\right),\]

and one may therefore distinguish three cases

\[\begin{split}\boldsymbol{U}_{1}^{\dagger}\left(\Theta_{h1},\Theta_{t1}\right)\boldsymbol{U}_{2}\left(\Theta_{h2},\Theta_{t2}\right)=\begin{cases} 0 & ;\quad\Theta_{h1}\Theta_{h2}=-\Theta_{t1}\Theta_{t2}\\ 4Re\left[z\right] & ;\quad\Theta_{h1}\Theta_{h2}=\Theta_{t1}\Theta_{t2}=+1\\ 4iIm\left[z\right] & ;\quad\Theta_{h1}\Theta_{h2}=\Theta_{t1}\Theta_{t2}=-1 \end{cases}.\quad (3)\end{split}\]

One may show that Hermiticity is conserved when multiplying such a trial vector by the electronic Hessian, whereas it is reversed by the generalized metric. However, both the electronic Hessian and the generalized metric conserves time reversal symmetry. The implication is that the time-symmetric and time-antisymmetric components of a solution vector does not mix upon solving the generalized eigenvalue problem, and one can dispense with one of them. From a physical point of view, this can be understood from the observation that excited states can be reached through both time-symmetric and time-antisymmetric operators. In order to employ the quaternion symmetry scheme, we choose to work with the time-symmetric vectors. It follows from the above that their scalar product are either zero or real. In practice, a property gradient is therefore always contracted with the component of the solution vector having the same hermiticity, so that all transition moments are real.

Practical example

We consider the excitation \(1s_{1/2}\rightarrow 3p_{3/2}\,(\left|m_{j}\right|=1/2)\) of the magnesium atom. We limit attention to boson irrep \(B_{1u}\), corresponding to the symmetry of the \(z\) - coordinate (\(\Gamma_z\)). We work within the electric dipole approximation and from a character table of \(D_{2h}\) it is clear that we will only get non-zero transition moments from the z-component of the electric dipole operator.

We start from the menu file ed.inp

**DIRAC                                                                                             
.WAVE FUNCTION                                                                                      
.ANALYZE                                                                                            
.PROPERTIES                                                                                         
**INTEGRALS                                                                                         
*READIN                                                                                             
.UNCONT                                                                                             
**WAVE FUNCTION                                                                                     
.SCF                                                                                                
*SCF                                                                                                
.CLOSED SHELL                                                                                       
6 6                                                                                                 
**PROPERTIES                                                                                        
*EXCITATION ENERGIES                                                                                
.PRINT                                                                                              
10                                                                                                  
.OCCUP                                                                                              
1                                                                                                   
                                                                                                    
.VIRTUAL                                                                                            
                                                                                                    
6                                                                                                   
.ANALYZE                                                                                            
! Excitations s -> p                                                                                
.EXCITATION ! B1u                                                                                   
5 1                                                                                                 
**ANALYZE                                                                                           
.MULPOP                                                                                             
*MULPOP                                                                                             
.VECPOP                                                                                             
1..oo                                                                                               
1..oo                                                                                               
**MOLECULE                                                                                          
*BASIS                                                                                              
.DEFAULT                                                                                            
dyall.ae2z                                                                                          
**END OF                                                                                            

and the xyz-file Mg.xyz

1                                                                                                   
                                                                                                    
Mg  0.0D0 0.0D0 0.0D0                                                                              

and run our calculation using:

pam --inp=ed --mol=Mg.xyz --get "DFCOEF PAMXVC"

We are keeping some file for later use, see below. It can be seen in the menu file ed.inp that excitations are only allowed from the first gerade orbital (\(1s_{1/2}\)) to the sixth ungerade one (\(3p_{3/2,1/2}\)). The nature of these orbitals can be seen from the Mulliken population analysis in the output.

We have increased the print level significantly in order to be able to look a bit under the hood. Let us first look at the property gradient

\[g_{ai}\rightarrow-\langle3p_{3/2;1/2}|\hat{\mu}_{z}|1s_{1/2;{\color{red}+1/2}}\rangle;\quad g_{a\overline{i}}\rightarrow-\langle3p_{3/2;1/2}|\hat{\mu}_{z}|1s_{1/2;{\color{red}-1/2}}\rangle;\quad\hat{\mu}_{z}=-ez.\]

The matrix element \(g_{ai}\) may be written as

\[g_{ai}=-\int\Omega_{ai}\left(\boldsymbol{r}\right)\hat{\mu}_{z}\left(\boldsymbol{r}\right)d^{3}\boldsymbol{r}; \quad\Omega_{ai}\left(\boldsymbol{r}\right)=\phi_{a}^{\dagger}\left(\boldsymbol{r}\right)\phi_{i}\left(\boldsymbol{r}\right)\]

and similar for \(g_{a\overline{i}}\). For present purposes it suffices to think in terms of 2-component spinors. Since \(s_{1/2}\) and \(p_{3/2}\) have opposite parity, they will in the conventions of DIRAC have different symmetry structures

\[\begin{split}s_{1/2}\rightarrow\left[\begin{array}{c} \left(\Gamma_{0},\Gamma_{R_{z}}\right)\\ \left(\Gamma_{R_{y}},\Gamma_{R_{x}}\right) \end{array}\right];\quad p_{3/2}\rightarrow\left[\begin{array}{c} \left(\Gamma_{xyz},\Gamma_{z}\right)\\ \left(\Gamma_{y},\Gamma_{x}\right) \end{array}\right].\end{split}\]

In addition we have to remember the relation between Kramers partners

\[\begin{split}\phi=\left[\begin{array}{c} \phi^{\alpha}\\ \phi^{\beta} \end{array}\right];\quad\overline{\phi}=\left[\begin{array}{c} -\phi^{\beta\ast}\\ \phi^{\alpha\ast} \end{array}\right].\end{split}\]

In terms of symmetry the orbital distributions are

\[\Omega_{\left(p_{3/2};s_{1/2}\right)}=\left(\Gamma_{xyz},\Gamma_{z}\right);\quad\Omega_{\left(p_{3/2};\overline{s_{1/2}}\right)}=\left(\Gamma_{y},\Gamma_{x}\right).\]

In the present case the operator symmetry is \(\Gamma_{z}\). Since the integrand of the matrix element has to be totally symmetric, we see that only the imaginary part of \(g_{ai}\) will be non-zero; the partner element \(g_{a\overline{i}}\) is strictly zero. The generally complex gradient elements are combined into a quaternion one

\[^{Q}g_{AI}=g_{ai}+g_{a\overline{i}}\check{j}\]

and, indeed, when we look in the output we find

                                Output from GPGET
                                -----------------

* Gradient (e-e) of property EpolL01X00Y00Z01

*** A imag part ***; NROW,NCOL,LRQ,LCQ:   1   1   1   1


               Column   1
       1      -0.00682451
    ==== End of matrix output ====

Here the notation EpolLnlXnxYnyZnz refers to a Cartesian electric multipole operator \(\hat{Q}^{[nl]}=-ex^{nx}y^{ny}z^{nz};\quad nx+ny+nz=nl\). Since this is a time symmetric, Hermitian operator, it will be contracted with the Hermitian part of the solution vector. In this particular case we have only a single orbital rotation, and between two orbitals of positive energy, so the relevant part of the solution vector is

 E+ part of solution vectors for  1 excitations:
* Vector nr.    1

*** A imag part ***; NROW,NCOL,LRQ,LCQ:   1   1   1   1


               Column   1
       1       0.35355084
    ==== End of matrix output ====

Searching in the output, we find the corresponding calculated transition moment

< 0 |EpolL01X00Y00Z01|  1 > = -4.8256194928E-03 a.u.,

which is easily seen to be \(2*0.35355084*-0.00682451\). (With respect to Eq.(3) a factor two seems to be missing; this is appears to be connected to the scaling of eigenvectors in subroutine XPPORD.)

We may proceed to visualize the transition moment density

\[T(\boldsymbol{r})=-\sum_{ai}x^{\ast}_{+;n;ai}\Omega_{ai}\left(\boldsymbol{r}\right)\hat{\mu}_{z}\left(\boldsymbol{r}\right)\]

as a radial distribution, that is

\[T(r) = \int_{0}^{2\pi}\int_{0}^{\pi}T(\mathbf{r})r^2\sin\theta d\theta d\phi\]

As preparation we first inspect the file PAMXVC of solution vectors

dirac 2074>rspread.x PAMXVC

 Solution vectors found on file : PAMXVC

** Solution vector : PP EXCITATIONB1u Irrep:   5 Trev:   1
  1 Type : E+ Freq.:    48.935572 Rnorm:   0.16E-13 Length:         1
  2 Type : E- Freq.:    48.935572 Rnorm:   0.16E-13 Length:         1

We see that the Hermitian e-e part of the solution vector E+ is the first vector on the file, which therefore leads us to set up the input file length.inp as

**DIRAC
.PROPERTIES
**INTEGRALS
*READIN
.UNCONT
**WAVE FUNCTION
.SCF
*SCF
.CLOSED SHELL
6 6
**GRID
.ULTRAFINE
**PROPERTIES
*EXPECTATION VALUES
.ORBANA
.OPERATOR
 'XXSECMOM'
.OPERATOR
 'YYSECMOM'
.OPERATOR
 'ZZSECMOM'
**VISUAL
.EDIPZ
 PAMXVC 1
.RADIAL
0.0D0 0.0D0 0.0D0
10.0
500
.3D_INT
**MOLECULE
*BASIS
.DEFAULT
dyall.ae2z
**END OF

We now generate the radial distribution using:

pam --inp=length --mol=Mg.xyz --put "DFCOEF PAMXVC" --get "plot.radial.scalar=length.dat"

First we may note at the end of the output

   The radial distribution integrates to:    4.8280356220858738E-003

and so we see that the radial distribution integrates to the value obtained in our previous calculation.

Next, with a suitable plotting program (I used xmgrace) we can visualize the radial distribution

To get some idea of the extent of the z-electric dipole moment transition density associated with the \(1s_{1/2}\rightarrow 3p_{3/2}\,(\left|m_{j}\right|=1/2)\) transition, I have included in the input the calculation of orbital sizes as \(<r^2>^{1/2}=\left(<x^2>+<y^2>+<z^2>\right)^{1/2}\). Extracting and using the data from the output I obtain the following table (orbital labels are taken from the Mulliken analysis of the first calculation):

Sym.	Orb.	\(<x^2>\)	\(<y^2>\)	\(<z^2>\)	\(<r^2>^{1/2}\)
E1g	1s 1/2; 1/2	0.008	0.008	0.008	0.151
E1g	2s 1/2; 1/2	0.190	0.190	0.190	0.754
E1g	3s 1/2; 1/2	4.106	4.106	4.106	3.510
E1u	2p 1/2; 1/2	0.198	0.198	0.198	0.772
E1u	2p 3/2; -3/2	0.240	0.240	0.120	0.774
E1u	3p 3/2; 1/2	0.160	0.160	0.279	0.774

We see that the radial distribution peaks around 0.25 \(a_0\) which is slightly outside the size of the \(1s_{1/2}\) orbital.