(sec:mcrpa.detailed)= # 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.{cite}`Yeager1979,Joergensen1988,Helmich-Paris2019a` 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. ## 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{cite}`Christiansen1998` $$\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 $$\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}$$ (mcrpa:o1rsp) The first-order response equations [\[mcrpa:o1rsp\]](#mcrpa:o1rsp){reference-type="eqref" reference="mcrpa:o1rsp"} 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: $$ \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) $$ 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{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}$$ 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 {ref}`sec:mcrpa.ntos.detailed` the most important natural transition orbitals{cite}`Martin2003,Helmich-Paris2019a` and active natural orbitals of MC-RPA and SA-CASSCF are shown, respectively. (fig:MCRPA_ni-dmg-uv)= ```{figure} ../../images/MCRPA_ni-dmg-2-uv-vis.* Calculated UV/Vis spectra of Ni(dmg)$_2$. ``` (sec:casscf.instability)= ## 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. ```{literalinclude} ../../orca_working_input/mcrpa_instab_benzene.inp :language: orca ``` 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. ``` 18 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 -------- -------- -------- -------- -------- -------- 2 H s 11.6 0.0 0.0 0.0 0.0 12.8 3 H s 11.6 0.0 0.0 0.0 0.0 12.8 4 H s 11.6 0.0 0.0 0.0 0.0 12.8 5 H s 11.6 0.0 0.0 0.0 0.0 12.8 6 C pz 0.0 0.0 32.2 0.0 30.9 0.0 7 C pz 0.0 0.0 32.2 0.0 30.9 0.0 8 C pz 0.0 24.1 8.0 23.2 7.7 0.0 8 C px 6.7 0.0 0.0 0.0 0.0 3.0 9 C pz 0.0 24.1 8.0 23.2 7.7 0.0 9 C px 6.7 0.0 0.0 0.0 0.0 3.0 10 C pz 0.0 24.1 8.0 23.2 7.7 0.0 10 C px 6.7 0.0 0.0 0.0 0.0 3.0 11 C pz 0.0 24.1 8.0 23.2 7.7 0.0 11 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 ``` (sec:mcrpa.ntos.detailed)= ## Natural Transition Orbitals Natural transition orbitals{cite}`Martin2003,Helmich-Paris2019a` (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 {numref}`fig:MCRPA-nto-tdm`. (fig:MCRPA-nto-tdm)= ```{figure} ../../images/MCRPA_nto-tdm.* 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 77 -> 69 : n= 0.16680786 76 -> 70 : n= 0.06575768 75 -> 71 : n= 0.02841330 74 -> 72 : n= 0.01485889 73 -> 73 : n= 0.01138840 72 -> 74 : n= 0.01099324 71 -> 75 : n= 0.00899487 70 -> 76 : n= 0.00764206 69 -> 77 : n= 0.00546103 68 -> 78 : n= 0.00517046 67 -> 79 : n= 0.00511544 66 -> 80 : n= 0.00485839 65 -> 81 : n= 0.00405297 64 -> 82 : n= 0.00379270 63 -> 83 : n= 0.00325052 62 -> 84 : n= 0.00315477 61 -> 85 : n= 0.00297172 60 -> 86 : n= 0.00291081 59 -> 87 : n= 0.00268810 58 -> 88 : n= 0.00243609 57 -> 89 : n= 0.00240536 56 -> 90 : n= 0.00238574 55 -> 91 : n= 0.00183946 54 -> 92 : n= 0.00181492 53 -> 93 : n= 0.00165943 52 -> 94 : n= 0.00154428 51 -> 95 : n= 0.00146299 50 -> 96 : n= 0.00137434 49 -> 97 : n= 0.00136656 48 -> 98 : n= 0.00128274 47 -> 99 : n= 0.00121332 46 -> 100 : n= 0.00117752 45 -> 101 : n= 0.00107962 44 -> 102 : n= 0.00105733 43 -> 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. {ref}`sec:utilities.plot.detailed`) ``` 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 {numref}`fig:MCRPA-ni-dmg-2-ntos`. 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^*$. (fig:MCRPA-ni-dmg-2-ntos)= :::{subfigure} AB :subcaptions: below :class-grid: outline ![(a) SA natural orbitals](../../images/MCRPA_ni-dmg-2-nat-orb.*) ![(b) NTOs](../../images/MCRPA_ni-dmg-2-ntos.*) Calculated UV/Vis spectra of Ni(dmg)$_2$. ::: ## 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. {ref}`sec:casscf.general.detailed`). 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. ## 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 ```