Is there a relation between transition density and density differences?

Question

When I have an excited state, is there a relation between the transition density associated with this electronic transition and the density difference between the excited state density and the ground state density? Can the transition density be equal to the difference density?

Answer 1

By definition, the transition density is not the same as the difference density. In the following derivation, the transition density is $\mathbf{T}^{\mathrm{CIS}}$ and the difference density is $\mathbf{P}^{\Delta}$. Although it may seem like they can be formed from the same quantities, the true difference density contains orbital relaxation terms and can be used to build a relaxed density, which is the correct excited state density. Neglecting those terms (the $P_{ia}^{\Delta}$) will give an unrelaxed density. I will also try to show that the unrelaxed difference density is still not equivalent to the transition density.

Derivation of individual CIS densities

Here is the derivation of the transition density for doubly-occupied spin-restricted MOs, with quoted blocks and equations directly taken from the original CIS paper.

The CIS approximation is identical to the TDA approximation in TD-DFT. Full TD-HF or TD-DFT follows the RPA equations, which include deexcitation terms in addition to excitation terms. The equations would be different for RPA, but conclusion would be identical.

The standard indexing convention is followed: $\mu, \nu, \lambda, \sigma$ for AOs, $i,j,k,l$ for occupied MOs, $a,b,c,d$ for virtual MOs, and $p,q,r,s$ for any MO. Define the Hartree-Fock (or SCF-type, valid for DFT as well) density as a sum over occupied MOs:

$$ P_{\mu\nu}^{\mathrm{HF}} = \sum_{i=1}^{n} c_{\mu i} c_{\nu i} $$

The two-particle CI-singles density matrix, $\mathbf{\Gamma}^{\mathrm{CIS}}$, can be written in terms of the HF ground-state density matrix and the ground-to-excited-state transition density matrix, $\mathbf{T}^{\mathrm{CIS}}$:

$$ \Gamma_{\mu\nu\lambda\sigma}^{\mathrm{CIS}} = \frac{1}{2} \left[ P_{\mu\nu}^{\mathrm{HF}} P_{\lambda\sigma}^{\mathrm{HF}} + 2 T_{\mu\nu}^{\mathrm{CIS}} T_{\lambda\sigma}^{\mathrm{CIS}} - P_{\mu\sigma}^{\mathrm{HF}} P_{\lambda\nu}^{\mathrm{HF}} - 2 T_{\mu\sigma}^{\mathrm{CIS}} T_{\lambda\nu}^{\mathrm{CIS}} \right] $$

where $\mathbf{P}^{\mathrm{HF}}$ is given (above) while $\mathbf{T}^{\mathrm{CIS}}$

$$ T_{\mu\nu}^{\mathrm{CIS}} = \sum_{ia} a_{ia} c_{\mu i} c_{\nu a} $$

I can't find the definition of $\{a\}$ in the paper, but it is the set of coefficients from a converged CIS calculation that describe the single-particle excitation from every occupied orbital to every unoccupied orbital. As an aside, they are the elements of the $\mathbf{X}$ vector in the RPA eigenvalue problem:

$$ \begin{pmatrix} \mathbf{A} & \mathbf{B} \\ \mathbf{B} & \mathbf{A} \end{pmatrix} \begin{pmatrix} \mathbf{X} \\ \mathbf{Y} \end{pmatrix} = \omega \begin{pmatrix} \mathbf{1} & \mathbf{0} \\ \mathbf{0} & \mathbf{-1} \end{pmatrix} \begin{pmatrix} \mathbf{X} \\ \mathbf{Y} \end{pmatrix} $$

$\mathbf{A}$ and $\mathbf{B}$ are parts of the orbital Hessian. If you diagonalize the whole thing, you get the set $\{\omega\}$ of all excitation energies back as eigenvalues. As mentioned above, setting $\mathbf{B} = \mathbf{0}$ gives the TDA or CIS approximation, so there is no longer a deexcitation vector $\mathbf{Y}$.

Importantly, we always start from the set of ground-state MO coefficients $\{c\}$ and describe excitations on top of them.

The CI-singles excited-state density matrix, $\mathbf{P}^{\mathrm{CIS}}$, is also constructed as a sum of HF and excited-state terms:

$$ P_{\mu\nu}^{\mathrm{CIS}} = P_{\mu\nu}^{\mathrm{HF}} + P_{\mu\nu}^{\Delta} $$

Here we have introduced $\mathbf{P}^{\Delta}$, the CI-singles $\Delta$ density matrix. This can also be called a "difference density matrix", since it represents the changes in electronic distribution upon excitation. It is not, however, the same as the transition density matrix $\mathbf{T}^{\mathrm{CIS}}$ defined above. As we shall demonstrate, it is the use of the true CI-singles density matrix required by (the CIS gradient expression that isn't important here) and not the simple one-particle density matrix (1PDM) which allows the realistic computation of charge distributions, orbital populations, and electronic moments of the excited state. To see this distinction, first consider the $\Delta$ density matrix which would be added to the HF density matrix in order to generate the 1PDM for an excited state. In the MO basis, it is a symmetric matrix with both occupied-occupied (OO) and virtual-virtual contributions

$$ \begin{align} P_{ij}^{\Delta} &= -\sum_{ab} a_{ia} a_{jb} \\ P_{ab}^{\Delta} &= +\sum_{ij} a_{ia} a_{jb} \end{align} $$

with the occupied-virtual (OV) elements all zero. The true CI-singles density matrix required in (the CIS gradient expression) will have exactly the same OO and VV contributions, but the OV terms are not all zero. The appearance of these off-diagonal block elements in the excited-state density matrix can be interpreted as orbital relaxation following the initial gross charge rearrangement due to excitation. That is to say, the CI coefficients will by themselves describe some of the gross features of charge redistribution in the excited state, but the response of the wave function to an external perturbation will account for further refinement in electronic properties. These OV terms can be found by solving a single set of CPHF equations:

$$ L_{ai} = \sum_{bj} \left[ (ij||ab) - (ib||ja) \right] P_{bj}^{\Delta} + (\epsilon_a - \epsilon_i) P_{ai}^{\Delta} $$

where the $\mathbf{L}$ vector is the CI-singles Lagrangian:

$$ \begin{align} L_{ai} &= C1_{ai} - C2_{ai} + \sum_{kl} P_{kl}^{\Delta} (al||ik) + \sum_{bc} P_{bc}^{\Delta} (ab||ic) \\ C1_{ci} &= -2 \sum_{jab} a_{ia} a_{jb} (cb||ja) \\ C2_{bk} &= -2 \sum_{ija} a_{ia} a_{jb} (ik||ja) \end{align} $$

The total CI-singles $\Delta$ density matrix required for $\mathbf{P}^{\mathrm{CIS}}$ can be generated by backtransforming the entire MO basis $\Delta$ density matrix defined by $P_{ij}^{\Delta}$, $P_{ab}^{\Delta}$, and $P_{ia}^{\Delta}$:

$$ P_{\mu\nu}^{\Delta} = \sum_{pq} P_{pq}^{\Delta} c_{\mu p} c_{\nu q} $$

Show that the transition density and unrelaxed difference density are not equivalent

We want to compare $T_{\mu\mu}^\mathrm{CIS}$ and $P_{pq}^{\Delta}$, but they are not in the same basis, so one must be transformed. Because I don't want to mess with the block-diagonal structure of $\mathbf{P}_{\mathrm{MO}}^{\Delta}$, I will transform $\mathbf{T}$ to the MO basis:

$$ \begin{align} T_{pq}^{\mathrm{CIS}} &= \sum_{\mu\nu} c_{\mu p} T_{\mu\nu}^{\mathrm{CIS}} c_{\nu q} \\ &= \sum_{\mu\nu} \sum_{ia} c_{\mu p} c_{\mu i} a_{ia} c_{\nu a} c_{\nu q}. \end{align} $$

Now consider the unrelaxed difference density, meaning the ov/vo blocks are zero due to neglecting orbital relaxation effects:

$$ P_{pq}^{\Delta} = \left( \begin{array}{c|c} -\sum_{ab} a_{ia} a_{jb} & \mathbf{0} \\ \hline \mathbf{0} & +\sum_{ij} a_{ia} a_{jb} \end{array} \right)_{pq}. $$

Here is a slightly hand-waving proof. I am not sure if the MO coefficients in the double sum can be simplified, but it doesn't matter; assume they are unity. $T_{pq}^{\mathrm{CIS}}$ contains terms with $a$, while $P_{pq}^{\Delta}$ contains terms with $|a|^2$. Unless there is an extra transformation that can be done, I don't see how they can be equivalent.

If you want a little test, generate some fake matrices and see for yourself. $\texttt{x} \equiv X_{ia} \equiv a_{ia}$.

#!/usr/bin/env python

from __future__ import print_function

import numpy as np


np.random.seed(42)

nocc = 2
nvirt = 2
norb = nocc + nvirt
nov = nocc * nvirt
nbasis = norb

x = np.random.rand(nocc, nvirt)
c = np.random.rand(nbasis, norb)
cocc = c[:, :nocc]
cvirt = c[:, nocc:]

p_delta_oo = -np.einsum('ia,jb->ij', x, x)
p_delta_vv = np.einsum('ia,jb->ab', x, x)

p_delta_ov = np.zeros(shape=(nocc, nvirt))
p_delta_vo = p_delta_ov.T

p_delta = np.asarray(np.bmat([[p_delta_oo, p_delta_ov],
                              [p_delta_vo, p_delta_vv]]))

# same as np.dot(cocc, np.dot(x, cvirt.T))
t_ao = np.einsum('ia,mi,na->mn', x, cocc, cvirt)
# same as np.dot(c.T, np.dot(t_ao, c))
t_mo = np.einsum('mp,mi,ia,na,nq->pq', c, cocc, x, cvirt, c)

print(t_mo)
print(p_delta)

which gives

[[ 1.82856579  1.89892758  0.59008502  3.09931192]
 [ 1.48134568  1.53535055  0.49193325  2.47904162]
 [ 0.54242862  0.56250364  0.17874155  0.9109361 ]
 [ 1.79242654  1.85791848  0.59456093  3.00118592]]
[[-1.75629929 -1.76345302  0.          0.        ]
 [-1.76345302 -1.77063588  0.          0.        ]
 [ 0.          0.          1.22441763  1.71443377]
 [ 0.          0.          1.71443377  2.40055604]]

Final comment

My original understanding of the term "difference density" comes from a coworker who was taking the real-space representation (cube files) of two state densities, say the ground state $0$ and an excited state $n$, and subtracting them.

$$ \Delta \mathbf{P}^{0\rightarrow n} = \mathbf{P}^{(n)} - \mathbf{P}^{(0)} $$

Is this equivalent to the definition of the difference density from above?

$$ \begin{align} \mathbf{P}^{(n)} - \mathbf{P}^{(0)} &= \mathbf{P}^{\mathrm{CIS}(n)} - \mathbf{P}^{\mathrm{HF}} \\ &= \left( \mathbf{P}^{\mathrm{HF}} + \mathbf{P}^{\Delta(n)} \right) - \mathbf{P}^{\mathrm{HF}} \\ &= \mathbf{P}^{\Delta(n)} \end{align} $$

Looks good!

References

• James B. Foresman, Martin Head-Gordon, John A. Pople, and Michael J. Frisch. Toward a systematic molecular orbital theory for excited states. The Journal of Physical Chemistry 1992 96 (1), 135-149, DOI: 10.1021/j100180a030

• Frank Neese. Prediction of molecular properties and molecular spectroscopy with density functional theory: From fundamental theory to exchange-coupling. Coordination Chemistry Reviews 2009 253 (5-6), 526-563, DOI: 10.1016/j.ccr.2008.05.014