:orphan:
====================================
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 :cite:`Bast2009`.
Excitation energies :math:`\omega_m` are found by solving the generalized eigenvalue problem
.. math::
\left(E_0^{[2]}-\hbar\lambda S^{[2]}\right)X_m=0;\qquad (1)
where the matrices :math:`E_0^{[2]}` and :math:`S^{[2]}` are the electronic Hessian and generalized metric, respectively. From their structure, it can be shown that solution vectors
come in pairs
.. math::
\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].
Transition moments associated with operator :math:`\hat{h}_{A}` for transitions :math:`n\leftarrow0` are generally calculated as
.. math::
\langle n|\hat{h}_{A}|0\rangle=\boldsymbol{X}_{+;n}^{\dagger}\boldsymbol{E}_{A}^{[1]};\qquad (2)
where appears solution vector :math:`\boldsymbol{X}_{+;n}^{\dagger}` and
property gradient :math:`\boldsymbol{E}_{A}^{[1]}`. The property gradient
is given by
.. math::
\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}.
An important generalization above is that, in addition to Hermitian
operators :math:`\hat{h}_{B}` (:math:`\Theta_{hB}=+1`), imposed by the tenets
of quantum mechanics, we also allow anti-Hermitian ones (:math:`\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 :math:`\boldsymbol{X}_{+}`
and :math:`\boldsymbol{X}_{-}`, we may form Hermitian and anti-Hermitian
combinations
.. math::
\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}
The inverse relations therefore provide a separation of solution vectors
into Hermitian and anti-Hermitian contributions
.. math::
\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
.. math::
\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
.. math::
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 :cite:`Saue2002a`
.. math::
\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
.. math::
\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)
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 :math:`1s_{1/2}\rightarrow 3p_{3/2}\,(\left|m_{j}\right|=1/2)`
of the magnesium atom. We limit attention to boson irrep :math:`B_{1u}`, corresponding to the symmetry of the :math:`z` - coordinate (:math:`\Gamma_z`).
We work within the electric dipole approximation and from a character
table of :math:`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`
.. literalinclude:: ed.inp
and the xyz-file `Mg.xyz`
.. literalinclude:: Mg.xyz
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 (:math:`1s_{1/2}`) to the sixth *ungerade* one (:math:`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
.. math::
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 :math:`g_{ai}` may be written as
.. math::
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 :math:`g_{a\overline{i}}`. For present purposes it suffices to think in terms of 2-component
spinors. Since :math:`s_{1/2}` and :math:`p_{3/2}` have opposite parity, they
will in the conventions of DIRAC have different symmetry structures
.. math::
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].
In addition we have to remember the relation between Kramers partners
.. math::
\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].
In terms of symmetry the orbital distributions are
.. math::
\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 :math:`\Gamma_{z}`. Since
the integrand of the matrix element has
to be totally symmetric, we see that only the imaginary part of :math:`g_{ai}`
will be non-zero; the partner element :math:`g_{a\overline{i}}` is strictly
zero. The generally complex gradient elements are combined into a
quaternion one
.. math::
^{Q}g_{AI}=g_{ai}+g_{a\overline{i}}\check{j}
and, indeed, when we look in the output we find
.. literalinclude:: propgrad.txt
Here the notation *EpolLnlXnxYnyZnz* refers to a Cartesian electric multipole operator :math:`\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
.. literalinclude:: E+_solvec.txt
Searching in the output, we find the corresponding calculated transition moment
.. literalinclude:: tmom.txt
which is easily seen to be :math:`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
.. math::
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
.. math::
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
.. literalinclude:: rspread.txt
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
.. literalinclude:: length.inp
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
.. literalinclude:: radint.txt
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
.. image:: ztm.png
To get some idea of the extent of the z-electric dipole moment transition density associated with the :math:`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 :math:`^{1/2}=\left(++\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. | :math:`` | :math:`` | :math:`` | :math:`^{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 :math:`a_0` which is slightly outside the size of the :math:`1s_{1/2}` orbital.