10.2 Calculation of Excitation Energies

With the ricc2 program excitation energies can presently be calculated with the RI variants of the methods CCS/CIS, CIS(D), CIS(D\(_\infty\)), ADC(2) and CC2. The CC2 excitation energies are obtained by standard coupled-cluster linear response theory as eigenvalues of the Jacobian, defined as derivative of the vector function with respect to the cluster amplitudes. \[\begin{equation} {\bf A}^{\text{CC2}}_{\mu\nu} = \frac{d \Omega_{\mu}}{d t_{\nu}} = \left( \begin{array}{cc} \langle \mu_{1} | [[\hat{H}+[\hat{H},T_2],\tau_{\nu_1}]|\text{HF} \rangle & \langle \mu_{1} | [\hat{H},\tau_{\nu_2}] |\text{HF} \rangle \\ \langle \mu_{2} | [\hat{H},\tau_{\nu_1}] |\text{HF} \rangle & \langle \mu_{2} | [F,\tau_{\nu_2}] |\text{HF} \rangle \end{array} \right) \end{equation}\](10.8) Since the CC2 Jacobian is a non-symmetric matrix, left and right eigenvectors are different and the right (left) eigenvectors \(E^i_{\nu}\) (\(\bar{E}^i_{\mu}\)) are not orthogonal among themselves, but form a biorthonormal basis (if properly normalized): \[\begin{equation} \bar{E}^{i} E^{j} = \bar{E}^i_{\mu_1} E^j_{\nu_1} + \bar{E}^i_{\mu_2} E^j_{\nu_2} = \delta_{ij} ~~. \end{equation}\](10.9) To obtain excitation energies only the right or the left eigenvalue problem needs to be solved, but for the calculation of transition strengths and first-order properties both, left and right, eigenvectors are needed (see below). A second complication that arises from the non-symmetric eigenvalue problem is that in the case of close degeneracies within the same irreducible representation (symmetry) it can happen that instead of two close lying real roots a degenerate complex conjugated pair of excitation energies and eigenvectors is obtained. CC2 (and also other standard coupled-cluster response methods) are thus not suited for the description of conical intersections etc. For the general theory behind coupled cluster response calculations see e.g. ref. [268, 269] or other reviews.

The ricc2 program exploits that the doubles/doubles block of the CC2 Jacobian is diagonal and the (linear) eigenvalue problem in the singles and doubles space can be reformulated as a (non-linear) eigenvalue problem in single-substitution space only: \[\begin{align*} {\bf A}^{eff}_{\mu_1\nu_1}(t,\omega) & = {\bf A}^{\text{CC2}}_{\mu_1\nu_1}(t) - {\bf A}^{\text{CC2}}_{\mu_1\gamma_2}(t) ({\bf A}_{\gamma_2\gamma_2}-\omega) {\bf A}^{\text{CC2}}_{\gamma2\nu_1}(t) \\ {\bf A}^{eff}_{\mu_1\nu_1}(t^{\text{CC2}},\omega^{\text{CC2}}) E_{\nu_1} & = \omega^{\text{CC2}} E_{\nu_1} \end{align*}\] This allows to avoid the storage of the double-substitution part of the eigen- or excitation vectors \(E_{\nu_2}\), \(\bar{E}_{\nu_2}\). The algorithms are described in refs. [12, 13], about the RI error see ref. [267].

The solution of the CC2 eigenvalue problem can be started from the solutions of the CCS eigenvalue problem (see below) or the trial vectors or solutions of a previous CC2 excitation energy calculation. The operation count per transformed trial vector for one iteration for the CC2 eigenvalue problem is about 1.3–1.7 times the operation count for one iteration for the cluster equations in the ground-state calculation—depending on the number of vectors transformed simultaneously. The disk space requirements are about \(O(V+N)N_x\) double precision words per vector in addition to the disk space required for the ground state calculation.

CCS excitation energies are obtained by the same approach, but here double-sub­sti­tutions are excluded from the expansion of the excitation or eigenvectors and the ground-state amplitudes are zero. Therefore the CCS Jacobian, \[\begin{equation} {\bf A}^{\text{CCS}}_{\mu\nu} = \frac{d \Omega_{\mu}}{d t_{\nu}} = \langle\mu_{1} | [H,\tau_{\nu_1}]|\text{HF}\rangle ~~, \end{equation}\](10.10) is a symmetric matrix and left and right eigenvectors are identical and form an orthonormal basis. The configuration interaction singles (CIS) excitation energies are identical to the CCS excitation energies. The operation count for a RI-CIS calculation is \({\cal O}(ON^2N_x)\) per iteration and transformed trial vector.

The second-order perturbative correction CIS(D) to the CIS excitation energies is calculated from the expression \[\begin{equation} \omega^{\text{CIS(D)}} = \omega^{\text{CIS}}+\omega^{\text{(D)}} = {\bf E}^{\text{CIS}} {\bf A}^{eff}(t^{\text{MP1}},\omega^{\text{CIS}}) {\bf E}^{\text{CIS}} \end{equation}\](10.11) (Note that \(t^{\text{MP1}}\) are the first-order double-substitution amplitudes from which also the MP2 ground-state energy is calculated; the first-order single-substitution amplitudes vanish for a Hartree–Fock reference due to the Brillouin theorem.) The operation count for a RI-CIS(D) calculation is similar to that of a single iteration for the CC2 eigenvalue problem. Also disk space requirements are similar.

Running excitation energy calculations:

The calculation of excitation energies is initiated by the data group $excitations in which at least the symmetries (irreducible representations) and the number of the excited states must be given (for other options see Section 25.2.23). With the following input the ricc2 program will calculate the lowest two roots (states) for the symmetries \(A_1\) and \(B_1\) of singlet multiplicity

¹

at the CIS, CIS(D) and CC2 level with default convergence thresholds. Ground-state calculations will be carried out for MP2 (needed for the CIS(D) model and used as start guess for CC2) and CC2.

$ricc2
  cis
  cis(d)
  cc2
$excitations
  irrep=a1 nexc=2
  irrep=b1 nexc=2 

The single-substitution parts of the right eigenvectors are stored in files named CCRE0-s--m-xxx, where s is the number of the symmetry class (irreducible representation), m is the multiplicity, and xxx the number of the excitation within the symmetry class. For the left eigenvectors the single-substitution parts are stored in files named CCLE0-s--m-xxx. These files can be kept for later restarts.

Trouble shooting:

For the iterative second-order methods CIS(D\(_\infty\)), ADC(2), and CC2 the solution of the nonlinear partitioned eigenvalue problem proceeds usually in three steps:

1. solution of the CCS/CIS eigenvalue problem to generate reasonable start vectors; the eigenvectors are converged in this step only to a remaining residual norm \(<\) preopt

2. pre-optimization of the eigenvectors by a robust modified Davidson algorithm (see ref. [12]) using the LINEAR CC RESPONSE SOLVER until the norm of all residuals are below preopt, combined with a DIIS extrapolation for roots assumed to be converged below the threshold thrdiis.

3. solution of the nonlinear eigenvalue problem with a DIIS algorithm using the DIIS CC RESPONSE SOLVER until the norm of the residuals are below the required threshold conv

This procedure is usually fairly stable and efficient with the default values for the thresholds. But for difficult cases it can be necessary to select tighter thresholds. In case of convergence problems the first thing do is to verify that the ground state is not a multireference case by checking the D1 diagnostic. If this is not the case the following situations can cause problems in the calculation of excitation energies:

• almost degenerate roots in the same symmetry class

• complex roots (break down of the CC approximation close to conical intersections)

• large contributions from double excitations

The first two reasons can be identified by running the program with a print level \(\le 3\). It will then print in each iteration the actual estimates for the eigenvalues. If some of these are very close or if complex roots appear, you should make sure that the DIIS procedure is not switched on before the residuals of the eigenvectors are small compared to the differences in the eigenvalues. For this, thrdiis (controlling the DIIS extrapolation in the linear solver) should be set about one order of magnitude smaller than the smallest difference between two eigenvalues and preopt (controlling the switch to the DIIS solver) again about one order of magnitude smaller then thrdiis.

Tighter thresholds or difficult situations can make it necessary to increase the limit for the number of iterations maxiter.

In rare cases complex roots might persist even with tight convergence thresholds. This can happen for CC2 and CIS(D\(_\infty\)) close to conical intersections between two states of the same symmetry, where CC response can fail due to its non-symmetric Jacobian. In this case one can try to use instead the ADC(2) model.

Diagnostics for double excitations:

As pointed out in ref. [14], the \(\%T_1\) diagnostic (or \(\%T_2 = 100 - \%T_1\)) which is evaluated directly from the squared norm of the single and double excitation part of the eigenvectors \(\%T_1 = 100 \cdot T_1 /(T_1+T_2)\) with \(T_i = \sum_{\mu_i} E_{\mu_i}^2\) where the excitation amplitudes are for spin-free calculations in a corresponding spin-adapted basis (which is not necessarily normalized) has the disadvantage that the results depend on the parameterization of the (spin-adapted) excitation operators. This prevents in particular a simple comparison of the results for singlet and triplet excited states if the calculations are carried out in a spin-free basis. With the biorthogonal representation for singlet spin-coupled double excitations [268] results for \(\%T_1\) also differ largely between the left and right eigenvectors and are not invariant with respect to unitary transformations of the occupied or the virtual orbitals.

The ricc2 module therefore uses since release 6.5 an alternative double excitation diagnostic, which is defined by \(\%{\cal T}_1 = 100 * {\cal T}_1 /({\cal T}_1+{\cal T}_2)\) with \({\cal T}_1 = \sum_{ai} E_{ai}^2\) and \({\cal T}_2 = \sum_{i>j}\sum_{a>b} E_{aibj}^2\) with \(E_{ai}\) and \(E_{aibj}\) in the spin-orbital basis. They are printed in the summaries for excitation energies under the headings %t1 and %t2. For spin-adapted excitation amplitudes \({\cal T}_1\) and \({\cal T}_2\) have to be computed from respective linear combinations for the amplitudes which reproduce the values in the spin-orbital basis. For ADC(2), which has a symmetric secular matrix with identical left and right normalized eigenvectors \({\cal T}_1\) and \({\cal T}_2\) are identical with the contributions from the singles and doubles parts for the eigenvectors to the trace of the occupied or virtual block of the (orbital unrelaxed) difference density between the ground and the excited state, i.e. the criterium proposed in ref. [14]. Compared to the suggestion from ref. [14] \({\cal T}_1\) and \({\cal T}_2\) have the additional advantage of that they are for all methods guaranteed to be positive and can be evaluated with the same insignificantly low costs as \(T_1\) and \(T_2\). They are invariant with respect to unitary transformations of the occupied or the virtual orbitals and give by construction identical results in spin-orbital and spin-free calculations. For CC2 and CIS(D\(_\infty\)) the diagnostics \({\cal T}_1\) and \({\cal T}_2\) agree for left and right eigenvectors usually within a few \(0.01\%\), for CIS(D) and ADC(2) they are exactly identical. For singlet excitations in spin-free calculations, \(\%{\cal T}_2\) is typically by a factors of 1.5–2 larger than \(\%T_2\). The second-order methods CC2, ADC(2), CIS(D\(_\infty\)) and CIS(D) can usually be trusted for \(\%{\cal T}_2 \le 15\%\).

For compatibility, the program can be switched to use of the old \(\%{\cal T}_1\) and \(\%{\cal T}_2\) diagnostics (printed with the headers ||T1|| and ||T2||) by setting the flag oldnorm in the data group $excitations. Note that the choice of the norm effects the individual results left and right one- and two-photon transition moments, while transition strengths and all other observable properties independent of the individual normalization of the right and left eigenvectors.

The \(\%{\cal T}_2\) and \(\%T_2\) diagnostics can not be monitored in the output of the (quasi-) linear solver. But it is possible to do in advance a CIS(D) calculation. The CIS(D) results for the \(\%{\cal T}_2\) and \(\%T_2\) correlate usually well with the results for this diagnostic from the iterativ second-order models, as long as there is clear correspondence between the singles parts of the eigenvectors. Else the DIIS solver will print the doubles diagnostics in each iteration if the print level is set \(> 3\). States with large double excitation contributions converge notoriously slow (a consequence of the partitioned formulation used in the ricc2 program). However, the results obtained with second-order methods for doubly excited states will anyway be poor. It is strongly recommended to use in such situations a higher-level method.

Visualization of excitations:

An easy way to visualize single excitations is to plot the natural transition orbitals that can be obtained from a singular value decomposition of the excitation amplitudes. See Sec. 22.1 for further details.

Another, but computational more involved possibility is to plot the difference density between the ground and the respective excited state. This requires, however, a first-order property or gradient calculation for the excited state to obtain the difference density. For further details see Sec. 10.3.3.

10.2.1 Core-Valence Separation (CVS) Approximation for Core Spectra

Core excited states are high in energy and, at the CC2 or ADC(2) level, embedded in a dense spectrum of doubly excited states for which the eigenvectors can not easily be converged in the doubles-direct implementation that is used in the ricc2 program. Two alternatives implemented in ricc2 for core spectra are damped response (see Secs. 10.5.1 and 10.6.1) and the core-valence separation (CVS) approximation[270, 5].

The CVS approximation decouples the core excited from the valence states by neglecting in the second-quantized Hamiltonian the coupling terms that change the number of electrons in the core orbitals. For the excitation energies for valence states and transition moments between valence states, the CVS approximation is equivalent to the frozen-core approximation. The excitation energies and amplitudes for singly core excited states are then obtained as eigenpairs of the block of the Jacobian matrix \(\mathbf{A}\) or, for ADC(2), the secular matrix with one core hole index (\(1ch\)): \[\begin{align} \left( \begin{array}{cc} \mathbf{A}_{val,val} & \mathbf{0} \\ \mathbf{0} & \mathbf{A}_{1ch,1ch} \end{array} \right) \end{align}\](10.12) The off-diagonal blocks \(\mathbf{A}_{val,1ch}\) and \(\mathbf{A}_{1ch,val}\) vanish in the CVS approximation. The calculation of core excited states needs as prerequisite an additional data group $core_excitations that specifies the subset of MOs from which core excitations are allowed. These need to be a subset of the core orbitals specified in the data group $freeze, e.g. for indole (C\(_8\)NH\(_7\)) in \(C_S\) symmetry at the carbon K-edge:

$freeze
  a' 1-9
$core_excitations
  a' 2-9

In addition the number of core holes has to be specified in the irrep option in data group $excitations. Possible values are 0 for valence states, which is the default, and 1 for states with one core hole:

$excitations
  irrep= a' multiplicity=1 nexc=4 npre=8 nstart=16 ncore=1
  irrep= a' multiplicity=1 nexc=2 npre=4 nstart=6  ncore=0
  spectrum states=all operators=diplen

The CVS approximation is available in the ricc2 program for excitation energies and transition strengths for between core-excited states and the ground-state or other core- or valence-excited states with CIS, ADC(2), and CC2. First-order properties, gradients, or non-linear spectra (two-photon, phosphorescence, MCD, etc.) for core excited states are not yet available.