7.32. Excited States via MC-RPA¶

MC-RPA excitation energies and transition moments are computed from the poles and residues of the linear response function of CASSCF a wave function.[380, 417, 902] By following similar lines, it is in principle possible to compute any kind of static and dynamic molecular property that is based on analytic derivatives of the CASSCF energy, which may be available in future releases of ORCA.

7.32.1. General Description¶

The starting point of response theory for variational wave functions like CASSCF is the time-dependent (TD) Schrödinger equation in its phase-isolated form[174]

\[\begin{aligned} \Big(\hat{H} - i\frac{\partial}{\partial t} - Q \Big)| \tilde{0} \rangle &= 0\end{aligned}\]

with the TD quasi energy

\[\begin{aligned} Q(t) &= \langle \tilde{0} | \Big(\hat{H} - i \frac{\partial}{\partial t} \Big)| \tilde{0} \rangle \text{.}\end{aligned}\]

The Hamiltonian \(\hat{H} = \hat{H}_0 + \hat{V}^t\) consists of the unperturbed time-independent Hamiltonian \(\hat{H}_0\) and a TD perturbation

\[\begin{aligned} \hat{V}^t &= \sum_k e^{-i \, \omega_k\, t} \sum_x \varepsilon_x( \omega_k )\, \hat{X}\end{aligned}\]

which is described as a sum of periodic perturbations, i.e. a Fourier series. TD molecular properties are obtained by applying the TD variational principal

\[\begin{aligned} \delta \{ Q(t) \}_T &= 0\end{aligned}\]

up to a certain order in the time-averaged quasi energy

\[\begin{aligned} \{ Q(t) \}_T &= \frac{1}{T} \int_{-T/2}^{T/2} Q(t) \,dt\end{aligned}\]

while \(\{ Q(t) \}_T\) is expanded by the perturbation strengths \(\varepsilon_X\) at vanishing frequencies \(\omega_k = 0\). Applying the TD variational principle for the second-order quasi energy leads to

(7.239)¶\[\begin{aligned} & 0 = \delta \{ Q(t) \}_T^{(2) } \Longrightarrow \frac{ \partial Q^{XY}(-\omega_Y, \omega_Y) }{ \partial \boldsymbol{\lambda}_X(-\omega_Y) } = \Big(\mathbf{E}^{(2) } - \omega_Y \, \mathbf{S}^{(2) } \Big)\boldsymbol{\lambda}^Y - \mathbf{V}^{Y} = 0 \text{.}\end{aligned}\]

The first-order response equations (7.239) become singular if the perturbation frequency \(\omega_Y\) approaches any eigenvalue

\[\begin{aligned} \mathbf{E}^{(2) } \, \boldsymbol{\lambda} = \mathbf{S}^{(2) } \, \boldsymbol{\lambda}\, \omega \text{,}\end{aligned}\]

of second-derivative matrices \(\mathbf{E}^{(2) }\) and \(\mathbf{S}^{(2) }\). The eigenvalues \(\omega\) correspond to the electronic excitation energies. The second-derivative matrices \(\mathbf{E}^{(2) }\) and \(\mathbf{S}^{(2) }\) have a paired structure as both kind of operators that express orbital excitation and de-excitations are involved:

\[\begin{split} \left[ \left( \begin{array}{cc} \mathbf{A} & \mathbf{B} \\ \mathbf{B}^* & {\mathbf{A}}^* \end{array} \right) - \omega \left( \begin{array}{cc} \boldsymbol{\Sigma} & \boldsymbol{\Delta} \\ -\boldsymbol{\Delta}^* & -\boldsymbol{\Sigma}^* \end{array} \right)\right] \left( \begin{array}{c} \mathbf{X} \\ \mathbf{Y}^* \end{array} \right) = \left( \begin{array}{c} \mathbf{0} \\ \mathbf{0} \end{array} \right) \end{split}\]

The eigenvalue equations above are valid for all variational wave functions methods, e.g. DFT, HF, CASSCF etc. The only difference is the operator manifold and the unperturbed Hamiltonian \(\hat{H}_0\) that is used. For the CASSCF linear response and eigen value equations, the super matrices \(\mathbf{A}\), \(\mathbf{B}\), \(\boldsymbol{\Sigma}\), and \(\boldsymbol{\Delta}\) have the following structure:

\[\begin{split}\begin{aligned} & \begin{array}{cc} \mathbf{A} = \left(\begin{array}{cc} \bra{0} [ q_i, [\hat{H}_0, q_j^{\dagger} ] ] \ket{0} & \bra{0} [ [ q_i, \hat{H}_0 ], R_j^{\dagger} ] \ket{0} \\ \bra{0} [ R_i , [ \hat{H}_0, q_j^{\dagger} ] ] \ket{0} & \bra{0} [ R_i, [ \hat{H}_0, R_j^{\dagger} ] ] \ket{0} \end{array} \right)% & \boldsymbol\Sigma = \left(\begin{array}{cc} \bra{0} [ q_i, q_j^{\dagger} ] \ket{0} & \bra{0} [ q_i, R_j^{\dagger} ] \ket{0} \\ \bra{0} [ R_i, q_j^{\dagger} ] \ket{0} & \bra{0} [ R_i, R_j^{\dagger} ] \ket{0} \end{array} \right) \\[2.0em] \mathbf{B} = \left(\begin{array}{cc} \bra{0} [ q_i, [\hat{H}_0, q_j ] ] \ket{0} & \bra{0} [ [ q_i, \hat{H}_0 ], R_j ] \ket{0} \\ \bra{0} [ R_i , [ \hat{H}_0, q_j ] ] \ket{0} & \bra{0} [ R_i, [ \hat{H}_0, R_j ] ] \ket{0} \end{array} \right)% & \boldsymbol\Delta = \left(\begin{array}{cc} \bra{0} [ q_i, q_j ] \ket{0} & \bra{0} [ q_i, R_j ] \ket{0} \\ \bra{0} [ R_i, q_j ] \ket{0} & \bra{0} [ R_i, R_j ] \ket{0} \end{array} \right) \end{array} \end{aligned}\end{split}\]

The TD CASSCF wave function is expressed in terms of orbital excitation \(q_i^{\dagger}\) and de-excitation operators \(q_i\),

\[\begin{aligned} & \begin{array}{lcr} q_i^{\dagger} = E_{pq} = a^{\dagger}_{p\alpha} a_{q\alpha} + a^{\dagger}_{p\beta} a_{q\beta} \text{,} & q_i = E_{qp} \text{,} & \text{ with } \quad p > q \end{array}\end{aligned}\]

as well as so called state transfer operators

\[\begin{aligned} & \begin{array}{lcr} R_i^{\dagger} = \ket{i}\bra{0} \text{,} & R_i = \ket{0}\bra{i} \text{,} & \text{ with } i \ne 0 \end{array}\end{aligned}\]

that account for relaxation of orbitals and CI coefficients when perturbed by an electromagnetic field, respectively.

The eigenvalue (and response) equations are solved iteratively by a customized version of the Davidson algorithm that simultaneously determines the N lowest lying roots. The most time-consuming step is the transformation of the trial vectors that contain an orbital and CI coefficient part with the electronic Hessian matrix \(\mathbf{E}^{(2) }\). The working equations are very similar to those of the CASSCF electronic gradient that is computed when minimizing the CASSCF ground state energy.

As a show case example the UV/Vis spectrum of a Nickel dimethylglyoximato complex (Ni(dmg) \(_2\)) was simulated with both SA-CASSCF and MC-RPA. A CAS (12/9) with the 3d electrons on Ni and 4 \(\pi\) orbitals and electrons from the ligands was selected; the def2-SVP basis set was used. For SA-CASSCF we have averaged over 21 states while for MC-RPA the 20 lowest roots were determined. Though both UV/Vis spectra have two intense peaks, their excitation energies and oscillator strengths differ quite substantially. This can be attributed to the lack of state-specific orbital relaxation that is only available in MC-RPA. In subsection Natural Transition Orbitals the most important natural transition orbitals[380, 562] and active natural orbitals of MC-RPA and SA-CASSCF are shown, respectively.

../../_images/MCRPA_ni-dmg-2-uv-vis.svg — Fig. 7.35 Calculated UV/Vis spectra of Ni(dmg)\(_2\).¶

7.32.2. Detecting CASSCF Instabilities¶

Selecting the right orbitals for the active space is not always an easy task. A wrong selection may lead to convergence to excited states or saddle points when minimizing the CASSCF energy. Such an instability in the wave function can be detected by computing the lowest excitation energy, i.e. the lowest root of the electronic Hessian with MC-RPA.

Instabilities may occur even for the simplest cases if the starting orbitals for CASSCF energy calculation were inappropriate. Let us look at a benzene CAS(6/6) calculation where we started from the model potential initial guess for the MOs.

! TZVP Def2-TZVP/C VeryTightSCF

%pal
 nprocs = 20
end

%casscf
 nel           6
 norb          6
 mult          1
 nroots        1
 TrafoStep     RI
 gtol  1e-8
 etol  1e-12
end

%mcrpa
 nroots        1
 TolR        1e-4
 MaxRed        80
end

*xyz 0 1
H           0.000000       2.484212       0.000000
H           0.000000      -2.484212       0.000000
H           2.151390       1.242106       0.000000
H          -2.151390      -1.242106       0.000000
H          -2.151390       1.242106       0.000000
H           2.151390      -1.242106       0.000000
C           0.000000       1.396792       0.000000
C           0.000000      -1.396792       0.000000
C           1.209657       0.698396       0.000000
C          -1.209657      -0.698396       0.000000
C          -1.209657       0.698396       0.000000
C           1.209657      -0.698396       0.000000
*

The energy converges smoothly

E(CAS)=  -230.560657053 Eh DE=     0.000000000
E(CAS)=  -230.767928524 Eh DE=    -0.207271471
E(CAS)=  -230.810472828 Eh DE=    -0.042544304
E(CAS)=  -230.811818980 Eh DE=    -0.001346152
E(CAS)=  -230.812866285 Eh DE=    -0.001047305
E(CAS)=  -230.812930081 Eh DE=    -0.000063796
E(CAS)=  -230.812944302 Eh DE=    -0.000014221
E(CAS)=  -230.812944635 Eh DE=    -0.000000333
E(CAS)=  -230.812943979 Eh DE=     0.000000656
E(CAS)=  -230.812944834 Eh DE=    -0.000000856
E(CAS)=  -230.812944944 Eh DE=    -0.000000110
E(CAS)=  -230.812944943 Eh DE=     0.000000001
E(CAS)=  -230.812944952 Eh DE=    -0.000000009
E(CAS)=  -230.812944953 Eh DE=    -0.000000000

as the gradient norm does

||g|| =      2.538097040 Max(G)=    -0.486818498 Rot=144,0
||g|| =      0.850498225 Max(G)=     0.219916319 Rot=23,11
||g|| =      0.192712320 Max(G)=    -0.055714161 Rot=23,11
||g|| =      0.144524323 Max(G)=     0.035741221 Rot=23,11
||g|| =      0.034101354 Max(G)=    -0.011346113 Rot=23,17
||g|| =      0.016336825 Max(G)=     0.005232633 Rot=122,17
||g|| =      0.003090776 Max(G)=     0.000935457 Rot=122,17
||g|| =      0.002539517 Max(G)=    -0.000597246 Rot=21,16
||g|| =      0.004134603 Max(G)=    -0.000980627 Rot=21,16
||g|| =      0.001410672 Max(G)=     0.000353494 Rot=21,16
||g|| =      0.000506790 Max(G)=    -0.000172518 Rot=23,17
||g|| =      0.000408528 Max(G)=    -0.000191970 Rot=122,17
||g|| =      0.000087961 Max(G)=    -0.000042271 Rot=23,17
||g|| =      0.000107447 Max(G)=    -0.000041858 Rot=23,17
||g|| =      0.000139470 Max(G)=    -0.000058693 Rot=23,11
||g|| =      0.000094867 Max(G)=     0.000020639 Rot=23,11
||g|| =      0.000033056 Max(G)=    -0.000011307 Rot=23,11
||g|| =      0.000009908 Max(G)=    -0.000005192 Rot=23,11
||g|| =      0.000013678 Max(G)=    -0.000007007 Rot=23,11
||g|| =      0.000011838 Max(G)=     0.000004991 Rot=23,17
||g|| =      0.000006431 Max(G)=    -0.000002076 Rot=122,17
||g|| =      0.000008817 Max(G)=    -0.000002828 Rot=122,17
||g|| =      0.000012070 Max(G)=     0.000004189 Rot=23,17
||g|| =      0.000001891 Max(G)=     0.000001284 Rot=23,17
||g|| =      0.000003505 Max(G)=     0.000001187 Rot=23,17
||g|| =      0.000002121 Max(G)=    -0.000000447 Rot=122,17
||g|| =      0.000002233 Max(G)=    -0.000000508 Rot=122,17
||g|| =      0.000000933 Max(G)=     0.000000494 Rot=23,17
||g|| =      0.000000711 Max(G)=     0.000000369 Rot=23,17
||g|| =      0.000000430 Max(G)=     0.000000230 Rot=23,17
||g|| =      0.000000200 Max(G)=     0.000000047 Rot=23,11
||g|| =      0.000000103 Max(G)=     0.000000030 Rot=23,11
||g|| =      0.000000025 Max(G)=     0.000000005 Rot=20,16

Though we have reached convergence for a CASSCF ground state energy calculation, the MC-RPA calculation however detects an instability

Davidson Eigenvalue solver (Iteration 10)

 State      Eigenvalue     RMSD error        Converged
   0          0.2405996888      1.4767724243e-01 F

 WARNING
 1 null space vectors in reduced space Hessians!
 This indicates an instability in your reference wave function!

 Davidson Eigenvalue solver (Iteration 11)

 State      Eigenvalue     RMSD error        Converged
   0          0.0000000000      0.0000000000e+00 T

by finding positive-indefiniteness by a Cholesky decomposition of the reduced space Hessians.

Instabilities in the CASSCF wavefunction can usually be avoided by carefully monitoring the active space orbitals in the

----------------------------
LOEWDIN REDUCED ACTIVE MOs  
----------------------------

section of the CASSCF output.

      19        20        21        22        23   
                  -0.62274  -0.32864  -0.32863   0.15983   0.16011   0.77971
                   1.99910   1.93810   1.93808   0.06203   0.06195   0.00075
                  --------  --------  --------  --------  --------  --------
H  s              11.6       0.0       0.0       0.0       0.0      12.8
H  s              11.6       0.0       0.0       0.0       0.0      12.8
H  s              11.6       0.0       0.0       0.0       0.0      12.8
H  s              11.6       0.0       0.0       0.0       0.0      12.8
C  pz              0.0       0.0      32.2       0.0      30.9       0.0
C  pz              0.0       0.0      32.2       0.0      30.9       0.0
C  pz              0.0      24.1       8.0      23.2       7.7       0.0
C  px              6.7       0.0       0.0       0.0       0.0       3.0
C  pz              0.0      24.1       8.0      23.2       7.7       0.0
C  px              6.7       0.0       0.0       0.0       0.0       3.0
C  pz              0.0      24.1       8.0      23.2       7.7       0.0
C  px              6.7       0.0       0.0       0.0       0.0       3.0
C  pz              0.0      24.1       8.0      23.2       7.7       0.0
C  px              6.7       0.0       0.0       0.0       0.0       3.0

In this particular example, MOs 18 and 23 are not part of the \(\pi\) system and have to be rotated with orbitals 16 and 31, respectively. After rotating all \(\pi\) orbitals into the active space, the CASSCF converges to a lower energy.

--------------
CASSCF RESULTS
--------------

Final CASSCF energy       : -230.844448647 Eh   -6281.5968 eV

The electronic CASSCF Hessian is now positive definite and the lowest MC-RPA excitation energy becomes

STATE  1:  E=   0.171023 au      4.654 eV    37535.3 cm**-1

7.32.3. Natural Transition Orbitals¶

Natural transition orbitals[380, 562] (NTO) are obtained from a singular value decomposition of the MC-RPA ground-to-excited state (f) transition density matrices \(\rho^{0\to f}_{pq}\). As for TD-DFT and ROCIS one obtains two sets of orbitals for each state that describe the donation (occupied and active) and acceptance (active and virtual) of an electron in the electronic transition. The orbital structure of \(\rho^{0\to f}_{pq}\) for CASSCF wave functions is illustrated in Fig. 7.36.

../../_images/MCRPA_nto-tdm.svg — Fig. 7.36 Structure of MC-RPA transition density matrix \(\rho^{0\to f}_{pq}\)¶

The compute NTOS only the following flag in the input has to switched ON:

%mcrpa
 nroots 20
 DoNTO true
 NTOThresh 5e-5
end

This will compute all NTOs with a singular value larger then the NTOThresh threshold for ALL roots.

------------------------------------------
NATURAL TRANSITION ORBITALS FOR STATE   12
------------------------------------------

Natural Transition Orbitals were saved in ni-dmg-2-svp-cas-12-9-mcrpa.s12.nto-donor
Threshold for printing occupation numbers 1.0000e-03
Natural Transition Orbitals were saved in ni-dmg-2-svp-cas-12-9-mcrpa.s12.nto-acceptor
Threshold for printing occupation numbers 1.0000e-03

STATE 13:  E=   0.214726 au      5.843 eV    47126.9 cm**-1

->  69  : n=  0.16680786
->  70  : n=  0.06575768
->  71  : n=  0.02841330
->  72  : n=  0.01485889
->  73  : n=  0.01138840
->  74  : n=  0.01099324
->  75  : n=  0.00899487
->  76  : n=  0.00764206
->  77  : n=  0.00546103
->  78  : n=  0.00517046
->  79  : n=  0.00511544
->  80  : n=  0.00485839
->  81  : n=  0.00405297
->  82  : n=  0.00379270
->  83  : n=  0.00325052
->  84  : n=  0.00315477
->  85  : n=  0.00297172
->  86  : n=  0.00291081
->  87  : n=  0.00268810
->  88  : n=  0.00243609
->  89  : n=  0.00240536
->  90  : n=  0.00238574
->  91  : n=  0.00183946
->  92  : n=  0.00181492
->  93  : n=  0.00165943
->  94  : n=  0.00154428
->  95  : n=  0.00146299
->  96  : n=  0.00137434
->  97  : n=  0.00136656
->  98  : n=  0.00128274
->  99  : n=  0.00121332
-> 100  : n=  0.00117752
-> 101  : n=  0.00107962
-> 102  : n=  0.00105733
-> 103  : n=  0.00104762

For the above example, the most important (controlled by NTOThresh) donating and accepting NTOs of state 13 are written to the gbw-type files

ni-dmg-2-svp-cas-12-9-mcrpa.s12.nto-donor
ni-dmg-2-svp-cas-12-9-mcrpa.s12.nto-acceptor

and can be plotted with the orca_plot program (see Sec. orca_plot)

orca_plot ni-dmg-2-svp-cas-12-9-mcrpa.s12.nto-donor -i

Please be aware of the different indices for states in the in and output!

To compute less or more NTOs the threshold NTOThresh can be adapted accordingly.

Let us come back to the UV/Vis spectra of Ni(dmg)\(_2\). For the two most intense peaks the natural orbitals and NTOs of MC-RPA and SA-CASSCF, respectively, are shown in Fig. 7.37. While the most intense peak in each spectrum (b and A) correspond to the same \(\pi \to \pi^*\) excitation, transition a and B are complete different, i.e. \(d \to \pi^*\) and \(\pi \to \pi^*\).

(a) SA natural orbitals — Fig. 7.37 Calculated UV/Vis spectra of Ni(dmg)\(_2\).¶

(b) NTOs — Fig. 7.37 Calculated UV/Vis spectra of Ni(dmg)\(_2\).¶

7.32.4. Computational Aspects¶

The code is intended to be used for medium-sized and larger open-shell molecules. It has the same scaling as ORCA’s first-order CASSCF energy implementation though a larger pre-factor as the computational cost grow “in principle” linearly with the number of roots.

The implementation is AO-driven meaning that the computational bottleneck is the Fock matrix construction for the several state-specific pseudo AO densities. Note that there are up to 6 pseudo AO densities for each state. The computational costs can be reduced significantly if the RIJCOSX approximation is employed, which is highly recommended.

The second most expensive part of the MC-RPA computation are the two-electron integrals with 3 active indices \(g_{ptuv}\). As we aim for running calculations on larger systems, there is only an implementation of the integral transformation that uses the resolution-of-the-identity (RI) approximation.

The restrictions on the auxiliary basis sets are the same as for the CASSCF code (Sec. General Description). That is

If the Fock matrices are constructed in Direct or Conventional mode, the /C bases are used for the RI approximation of the \(g_{ptuv}\) integrals.
If the RIJCOSX approximation for the Fock matrices is employed, the /JK bases are used for both the Fock matrices and the \(g_{ptuv}\) integrals.

Note that MC-RPA implementation can be run in parallel with MPI which allows for computing UV/Vis and ECD spectra large open-shell molecules in a limited amount of time.

Before starting running MC-RPA, it is recommended to converge the state specific CASSCF energy calculation until you hit the point of stagnating convergence. Note that property calculations in general assume vanishing electronic gradients otherwise numerical issues in the eigenvalue / response equations may occur.

7.32.5. Keyword List¶

%mcrpa

       NRoots    0         # The number of desired roots

       TolR      1e-5      # Convergence threshold for residual norm
       TolLDP    1e-5      # linear dependency threshold for generalized eigenvalue problem
       PrintWF   CFG       # print CI part in (CFG (default), CSF, or DET basis)
       TolPrint  1e-2      # Threshold for printing elements of excitation vector

       MaxRed    200       # maximum size of reduces space for ALL VECTORS IN TOTAL
       MaxIter   100       # maximum number of (Davidson) iterations
                           
       TDA       false     # Switch off full TD-CASSCF (Tamm-Dancoff approximation)
       
       DoNTO     false     # Generate Natural Transition Orbitals
       NTOThresh 1e-3      # Threshold for printing occupation numbers

       DoCD      true      #Request circular dichroism calculation
       DoDipoleLength true #Request the use of electric moments in a length formulation
       DoDipoleVelocity true #Request the use of electric moments in a velocity formulation
       DoHigherMoments true  #Request the calculation of electric quadrupole and magnetic dipole moments contributions
       DoFullSemiclassical true #Request the calculation of complete semiclassical multipolar moments
       DecomposeFoscLength true #Request the decomposition of the oscillator strengths in a multipolar expansion under a length formulation
       DecomposeFoscVelocity true #Request the decomposition of the oscillator strengths in a multipolar expansion under a velocity formulation