7.34. Excited States via STEOM-CCSD¶

The EOM-CCSD approach for excitation energies becomes prohibitively costly for large systems because of its O(\(N^{6}\)) scaling. Therefore, one needs a more compact form of the wave-function ansatz. A second similarity transformation can compress the final matrix diagonalization step to the CIS space only. The resulting STEOM-CCSD method of Marcel Nooijen and co-workers [638] is an efficient way for accurate calculations of excitation energies.

7.34.1. General Description¶

In the standard EOM-CC method, the transformed Hamiltonian is diagonalized over a singles and doubles space to obtain ionized, attached, or excited states of the reference state. In STEOM-CC, one performs a second similarity transformation

\[\begin{aligned} \hat{G}=\{e^{\hat{S} }\}^{-1}\hat{\bar{H} }\{e^{\hat{S} }\}\end{aligned}\]

The transformation operator \(\hat{S}\), including singles and doubles, is defined as

\[\begin{aligned} \hat{S}=\hat{S}^{IP}+\hat{S}^{EA},\end{aligned}\]

\[\begin{aligned} \hat{S}^{IP}=S_{i^{\prime} }^{m}\hat{E}^{m}_{i^{\prime} }+\dfrac{1}{2}S^{mb}_{ij}\hat{E}^{mb}_{ij},\end{aligned}\]

\[\begin{aligned} \hat{S}^{EA}=S^{a^{\prime} }_{e}\hat{E}_{e}^{a^{\prime} }+\dfrac{1}{2}S^{ab}_{ej}\hat{E}^{ab}_{ej}.\end{aligned}\]

In the above equations, \(m\) and \(e\) denote active indices of the hole and particle type respectively, while a prime denotes a restriction to orbitals that are not active. The amplitudes of the operator \(\hat{S}\) are defined in such a way that matrix elements of the transformed Hamiltonian, in second quantized notation, become equal to zero.

\[\begin{aligned} g_{i^{\prime} }^{m}=g_{ij}^{mb}=g^{a^{\prime} }_{e}=g^{ab}_{ej}=0\end{aligned}\]

In addition, the zeros which pre-existed in \(\bar{H}\), after solving the CCSD equations, remain preserved. The above equations are linear in \(\hat{S}\) and are equivalent to the Fock space multireference coupled cluster equations for the one valence problem. However, to ensure numerical stability, the equations are re-casted as matrix diagonalization problem and solved as IP-EOM-CCSD and EA-EOM-CCSD problems. The \(\hat{S}^{IP}\) and \(\hat{S}^{EA}\) are then extracted from the converged previous calculations, respectively, by invoking intermediate normalization on the suitably chosen eigenvectors corresponding to active holes and active particles. The total process can be described as following

Solution of the ground state coupled cluster equations
Construct the first similarity transformed Hamiltonian as \(\hat{\bar{H} }=e^{-\hat{T} }\hat{H}e^{\hat{T} }\)
Solution of the IP-EOM and EA-EOM equations
Extraction of the \(\hat{S}\) amplitudes
Construct the second similarity transformed Hamiltonian as \(\hat{\bar{G} }=e^{-\hat{S} }\hat{\bar{H} }e^{\hat{S} }\)
Diagonalization of \(\hat{\bar{G} }\) in CIS space

The advantage of the above method is that, instead of one iterative O(\(N^{6}\)) scaling diagonalization step, it requires two iterative O(\(N^{5}\)) scaling steps, one non-iterative O(\(N^{5}\)) scaling step and one iterative O(\(N^{4}\)) scaling matrix diagonalization step. The presence of so-called ‘implicit triples excitation’ term ensures the charge transfer separability of the excited states, which is absent in EOM-CCSD. In addition, since the final diagonalization step is performed in a CIS space, the spin adaption is trivial and excited states of triplet multiplicity can be obtained without going through the complications of a spin orbital based implementation.

The STEOMCC approach has also recently been extended for applications to open-shell systems within the UHF formalism [403]. In this case, the expressions for the operators \(\hat{S}^{IP}\) and \(\hat{S}^{EA}\) take the form,

\[\begin{split}\begin{aligned} \hat{S}^{IP}&=\frac{1}{2}\sum_{i,e,a,b}s_{ie}^{ab}\big\{\hat{a}^{\dagger}\,\hat{b}^{\dagger}\,\hat{e}\,\hat{i}\big\}+\sum_{\bar{i},e,\bar{a},b}s_{\bar{i}e}^{\bar{a}b}\big\{\hat{\bar{a} }^{\dagger}\,\hat{b}^{\dagger}\,\hat{e}\,\hat{\bar{i} }\big\}\nonumber\\ &+\frac{1}{2}\sum_{\bar{i},\bar{e},\bar{a},\bar{b} }s_{\bar{i}\bar{e} }^{\bar{a}\bar{b} }\big\{\hat{\bar{a} }^{\dagger}\,\hat{\bar{b} }^{\dagger}\,\hat{\bar{e} }\,\hat{\bar{i} }\big\}+\sum_{i,\bar{e},a,\bar{b} }s_{i\bar{e} }^{a\bar{b} }\big\{\hat{a}^{\dagger}\,\hat{\bar{b} }^{\dagger}\,\hat{\bar{e} }\,\hat{i}\big\},\\ \widehat{S}_{-}&=\frac{1}{2}\sum_{i,j,a,m}s_{ij}^{am}\big\{\hat{a}^{\dagger}\,\hat{m}^{\dagger}\,\hat{j}\,\hat{i}\big\}+\sum_{\bar{i},j,\bar{a},m}s_{\bar{i}j}^{\bar{a}m}\big\{\hat{\bar{a} }^{\dagger}\,\hat{m}^{\dagger}\,\hat{j}\,\hat{\bar{i} }\big\}\nonumber\\ &+\frac{1}{2}\sum_{\bar{i},\bar{j},\bar{a},\bar{m} }s_{\bar{i}\bar{j} }^{\bar{a}\bar{m} }\big\{\hat{\bar{a} }^{\dagger}\,\hat{\bar{m} }^{\dagger}\,\hat{\bar{j} }\,\hat{\bar{i} }\big\}+\sum_{i,\bar{j},a,\bar{m} }s_{i\bar{j} }^{a\bar{m} }\big\{\hat{a}^{\dagger}\,\hat{\bar{m} }^{\dagger}\,\hat{\bar{j} }\,\hat{i}\big\}.\end{aligned}\end{split}\]

where we use overbars to distinguish the \(\beta\) orbitals from the \(\alpha\) orbitals. The amplitudes \(\big\{s_{ie}^{ab}, s_{\bar{i}e}^{\bar{a}b}\big\}\) are determined by solving the UHF EA-EOM-CCSD equations for the attachment of an \(\alpha\) electron, while the \(\big\{s_{\bar{i}\bar{e} }^{\bar{a}\bar{b} }, s_{i\bar{e} }^{a\bar{b} }\big\}\) amplitudes are extracted from a UHF EA-EOM-CCSD calculation for the attachment of a \(\beta\) electron. Similarly, the sets of amplitudes \(\big\{s_{ij}^{am}, s_{\bar{i}j}^{\bar{a}m}\}\) and \(\big\{s_{\bar{i}\bar{j} }^{\bar{a}\bar{b} }, s_{i\bar{j} }^{a\bar{m} }\big\}\) are determined by solving the decoupled UHF IP-EOM-CCSD problems for the ionization of an \(\alpha\) electron and the ionization of a \(\beta\) electron, respectively. Hence, an UHF STEOMCC calculation involves two separate IP calculations (O(\(N^{5}\)) scaling) and two separate EA calculations (O(\(N^{5}\)) scaling steps).

All the speed up options, including CCSD(2) (only available in RHF implementation) and COSX, which are available for EOM-CCSD are also available for STEOMCC. The most important steps in a STEOMCC calculation are the IP-EOM and EA-EOM calculations. These steps are performed using the EOM-CCSD module and the relevant keywords are the same as that described in Excited States via EOM-CCSD. The keywords which are exclusive to the RHF STEOM module are:

%mdci
#RHF STEOM parameters - defaults displayed
DoCISNat true        # automatic selection of active space       
NActIP  3            # number of states defined as active in the IP calculation
NActEA  2            # number of states defined as active in the EA calculation
DoTriplet false      # target state of triplet multiplicity
DoDbFilter true      # filters out states with doubles excitation  character 
DoNewActSch true    # new active space selection scheme for STEOM-CCSD                             
DoSOLV               # perturbative correction for solvation effects (experimental) 
#Default values for automatic active space selection scheme
OThresh 0.001        # Cut off occupation of CIS natural orbitals in IP calculation 
VThresh 0.001        # Cut off occupation of CIS natural orbitals in EA calculation 
IPSThrs 80            #  The percentage singles threshold for the IP calculation
EASThrs 80           #  The percentage singles threshold for the EA calculation
end

The keywords pertaining to the UHF STEOM module are:

%mdci
#UHF STEOM parameters - defaults displayed
DoCISNat true        # automatic selection of active space       
NActIP_a  3          # number of states defined as active in the IP calculation 
                     # for the removal of an α electron
NActIP_b  3          # number of states defined as active in the IP calculation 
                     # for the removal of a β electron
NActEA_a  2          # number of states defined as active in the EA calculation 
                     # for the attachment of an α electron
NActEA_b 2           # number of states defined as active in the EA calculation 
                     # for the attachment of a β electron
DoDbFilter true      # filters out states with doubles excitation  character 
UseQROs  false       # use QROs or not
DoNewActSch true     # new active space selection scheme for STEOM-CCSD                             
#Default values for automatic active space selection scheme
OThresh 0.001        # Cut off occupation of CIS natural orbitals in tIP calculations
VThresh 0.001        # Cut off occupation of CIS natural orbitals in EA calculations 
IPSThrs 80           #  The percentage singles threshold for the IP calculations
EASThrs 80           #  The percentage singles threshold for the EA calculations
end

7.34.2. Selection of Active space¶

The results of a STEOM-CC calculation depend upon the number of roots selected as active in the EOMIP and EOMEA calculations. In ORCA, they are chosen automatically, by using state-averaged CIS natural transition orbitals (NTO). By default, the number of roots included in this initial CIS computation is equal to the number of roots requested in STEOM (NRoots). However, this can be modified setting NRootsCISNAT to higher values. The orbitals up to a predefined occupation are then chosen to be active in the EOMIP and EOMEA calculations, and this is controlled by the keywords OThresh and VThresh respectively. Now, there are two possible ways to chose active space. One is to use the criteria of percentage occupation of NTO’s as described in ref [236]. However, a newer and more robust approach is to use the criteria of absolute occupation, which is default in the current implementation. One can switch on the old percentage occupation based active space selection by setting DoNewActSch to false (not recommended).

One can also select the active spaces manually by turning the DoCISNat to false and setting the NActIP and NActEA (RHF STEOM calculation) or the NActIP_a, NActIP_b, NActEA_a and NActEA_b (UHF STEOM calculation) to desired values. However, this is not recommended for general uses. The following snippet shows the output of the active orbital selection procedure on a closed-shell molecule:

------------------------------------------
STATE AVERAGED NATURAL ORBITALS FOR ACTIVE SPACE SELECTION
------------------------------------------

Solving eigenvalue problem for the occupied space  ... Occupied block occupation :  
 0.000478 
 0.002266 
 0.169928 
 0.171663 
 0.310125 
 0.345541 
Orbital taken as active for IP roots:  
 0.345541 
 0.310125 
 0.171663 
 0.169928 
done
Solving eigenvalue problem for the virtual space   ... Virtual block occupation :  
 0.640886 
 0.332262 
 0.017272 
 0.005326 
 0.001752 
 0.000667 
 0.000574 
 0.000540 
 0.000160 
 0.000150 
 0.000139 
 0.000086 
 0.000082 
 0.000037 
 0.000023 
 0.000016 
 0.000013 
 0.000008 
 0.000003 
 0.000002 
 0.000001 
 0.000000 
 0.000000 
 0.000000 
 0.000000 
 0.000000 
 0.000000 
 0.000000 
 0.000000 
-0.000000 
Orbital taken as active for EA roots :  
 0.640886 
 0.332262 
 0.017272 
done
 No of roots active in IP calculation:    4  
 No of roots active in EA calculation:    3

7.34.3. Active space selection using TD-DFT densities¶

Instead of using a CIS calculation for selected the Active Space roots, a TD-DFT based one can also be considered. Be aware that using DFT Kohn-Sham orbitals for computing the CCSD GS energy can lead to some instabilities and give incorrect results.

The main interest of this approach is to start the STEOM-CCSD calculation with TD-DFT electronic densities which are in general better than the CIS one, especially for some specific compounds (metallic complexes for example). The computed TD-DFT densities are also often more stable than the CIS one. It will however slow down the calculation.

The input has to be written like this:

!RHF/UHF BHANDHLYP STEOM-CCSD TZVP 

%mdci
nroots 10
tddftguess true
end

%tddft
nroots 10
end

*xyz 0 1

Any DFT functional can be used but we recommend one with a decent amount of HF exchange. On top of this, the keyword TDDFTGuess has to be set to true in mdci block and the tddft has to be added together with the NRoots keyword. In both input blocks (%mdci and %tddft) the same number of roots has to be given. Starting from ORCA 6, using TD-DFT guess with an UHF reference is also possible.

7.34.4. The reliability of the calculated excitation energy¶

The excitation energy for any states calculated in STEOM-CC are only reliable when the dominant excitation for that states are confined within the active space. This can be verified from the percentage active character of the calculated states, an a posteriori diagnostic which is defined as

\[\begin{aligned} \% active character= \frac{\sum\limits_{m,e} C(m,e)*C(m,e) }{\sum\limits_{i,a} C(i,a)*C(i,a) }*100\end{aligned}\]

for closed-shell systems and takes the form,

\[\begin{aligned} \% active character=\frac{\sum\limits_{m,e}C(m,e)*C(m,e)+\sum\limits_{\bar{m},\bar{e} }C(\bar{m},\bar{ e})*C(\bar{m},\bar{ e}) }{\sum\limits_{i,a}C(i,a)*C(i,a)+\sum\limits_{\bar{i},\bar{a} }C(\bar{i},\bar{a})*C(\bar{i},\bar{a}) }*100.\end{aligned}\]

within the UHF formalism. The roots which have \(\% active character\) higher than 98.0 are considered to be converged with respect to the active space.

------------------
STEOM-CCSD RESULTS
------------------

IROOT=  1:  0.145412 au     3.957 eV   31914.3 cm**-1
  Amplitude    Excitation
  -0.169361     4 ->   8
  -0.984822     7 ->   8

Percentage Active Character     99.86

  Amplitude    Excitation in Canonical Basis
  -0.166580     4 ->   8
  -0.975432     7 ->   8
  -0.124356     7 ->  13

IROOT=  2:  0.309409 au     8.419 eV   67907.5 cm**-1
  Amplitude    Excitation
   0.994141     7 ->   9

Percentage Active Character     99.78

  Amplitude    Excitation in Canonical Basis
  -0.990029     7 ->   9

IROOT=  3:  0.336993 au     9.170 eV   73961.4 cm**-1
  Amplitude    Excitation
  -0.994078     5 ->   8

Percentage Active Character     99.10

  Amplitude    Excitation in Canonical Basis
  -0.984116     5 ->   8
  -0.136769     5 ->  13

IROOT=  4:  0.357473 au     9.727 eV   78456.2 cm**-1
  Amplitude    Excitation
   0.181761     4 ->  10
   0.728209     6 ->   8
   0.611668     7 ->  10
  -0.191540     7 ->  12

Percentage Active Character     94.10

Warning:: the state may have not converged with respect to active space 
-------------------- Handle with Care -------------------- 

  Amplitude    Excitation in Canonical Basis
  -0.184144     4 ->  10
  -0.725183     6 ->   8
  -0.633718     7 ->  10

IROOT=  5:  0.386654 au    10.521 eV   84860.8 cm**-1
  Amplitude    Excitation
   0.980406     4 ->   8
  -0.178551     7 ->   8

Percentage Active Character     99.79

  Amplitude    Excitation in Canonical Basis
   0.971678     4 ->   8
   0.122877     4 ->  13
  -0.179242     7 ->   8

IROOT=  6:  0.444881 au    12.106 eV   97640.1 cm**-1
  Amplitude    Excitation
  -0.995150     6 ->   9

Percentage Active Character     99.69

  Amplitude    Excitation in Canonical Basis
  -0.989966     6 ->   9

If the \(\% active character\) for any calculated state is less than 98, that state may have not converged with respect to active space and the excitation energy for that particular state is less reliable. The user should request a larger number of roots under those conditions.

7.34.5. Removal of IP and EA states with double excitation character¶

To obtain accurate results with STEOM-CCSD, only the \(\hat{S}\) amplitudes corresponding to the states dominated by single excitations should be included in the second similarity transformation. This is ensured in ORCA in two ways. First, the root following (FollowCIS) is activated by default so that it converges to the states dominated by singly excited guess vectors. This avoids the calculation of so called ‘satellite states’, which are of double excitation character with respect to the ground state. Secondly, among the converged IP and EA roots, the states which have %singles character below a certain predefined threshold (i.e. controlled by the keywords IPThresh and EAThresh) are automatically excluded from the second similarity transformation.

EOM-CCSD RESULTS
----------------

IROOT=  1:  0.105316 au     2.866 eV   23114.2 cm**-1
  Amplitude    Excitation
   0.697547     x ->   8
IROOT=  2:  0.217925 au     5.930 eV   47829.1 cm**-1
  Amplitude    Excitation
  -0.701454     x ->   9
IROOT=  3:  0.304098 au     8.275 eV   66741.8 cm**-1
  Amplitude    Excitation
  -0.700458     x ->  10
IROOT=  4:  0.350387 au     9.535 eV   76901.1 cm**-1
  Amplitude    Excitation
   0.702705     x ->  11
IROOT=  5:  0.651462 au    17.727 eV  142979.4 cm**-1
  Amplitude    Excitation
   0.637352     x ->  12
   0.121747     x ->   8     4 ->  10
   0.177039     x ->   8     5 ->   9
   0.109987     x ->   9     5 ->   8
  -0.206789     x ->   8     7 ->  10
  -0.109870     x ->  10     7 ->   8

EA STATE=  0:  percentage singles     95.282 
EA STATE=  1:  percentage singles     96.981 
EA STATE=  2:  percentage singles     96.540 
EA STATE=  3:  percentage singles     97.844 
EA STATE=  4:  percentage singles     68.884 

Warning: high double excitation character, excluding from the STEOM transformation
 Final no active EA roots:    4

Note that the use of CIS natural transition orbitals can lead to convergence issues for the IP and EA states which are dominated by double excitation character. This can be remedied by setting DoDbFilter to true.

7.34.6. Transition and difference densities¶

At the end of a STEOM computation, it is possible to store the final eigenvectors in a file “job.cis”, in analogy with what is done for CIS and TD-DFT computations. This file can be obtained by setting DoStoreSTEOM true in the input. This file can then be processed by orca_plot to obtain the difference and / or the transition densities.

An Natural Transition Orbitals analysis can also be performed within the STEOM-CCSD scheme, as described in Natural Transition Orbitals. It can be performed by setting the keyword DoSTEOMNatTransOrb to true.

7.34.7. Properties¶

The dipolar and transition moments (as well as the oscillator strength) can be computed within the STEOM module using different kinds of approximations. Please cite our paper on these corrected STEOM transition densities [297]! Starting from ORCA 5, new defaults (DoSimpleDens false) are used that are much better than the previous CIS-like approximation, and the full option is of CC3-like quality.

%mdci

DosimpleDens false      # Default, using the STEOM-CCSD density + some doubles effect.
AddL2term true

DosimpleDens false      # using the STEOM-CCSD density + some doubles effect.
AddL2term true          # + neglected GS double
UpdateL1 true

DosimpleDens false      # using the STEOM-CCSD density + some doubles effect.
AddL2term true          # + neglected GS double + doubles from EOM-CCSD
UpdateL1 true           # (expensive, but of CC3 quality - see reference)
AddDDTerm true

end

By default, the STEOM-CCSD densities with AddL2term true should be used for all calculation as discussed in ref. [297].

7.34.8. Solvation (Experimental)¶

In STEOM-CCSD, the excitation energies and densities can be corrected using the CPCM solvation scheme in ORCA.

To use it, the keyword DoSolv has to be set to true in the %mdci block and the simple keyword CPCM (or SMD) + name of the solvent has to be given.

!CPCM(ethanol) STEOM-CCSD TightSCF def2-TZVP def2-TZVP/C def2/J 

%mdci
Nroots 5
DoSolv true
end

7.34.9. Spin-Orbit Coupling (Experimental)¶

You can compute the spin-orbit coupling between singlet and triplets states in STEOM-CCSD using the keyword STEOMSOC true. Please note that all SOC matrix elements and properties are currently computed from the right vector only!

7.34.10. Core excitation¶

The STEOM-CCSD (and bt-PNO-STEOM-CCSD) method can also be used to compute the K-edge core-excitation energy of molecules. See Core-Excitation for more details.

7.34.11. Transient absorption¶

Transient absorption spectra can be computed using the keyword DoTrans true. The IRoot keyword will select the targeted excited state.