(sec:mreom.typical)= # MR-EOM-CC: Multireference Equation of Motion Coupled-Cluster The Multireference Equation of Motion Coupled-Cluster (MR-EOM-CC) methodology {cite}`datta2011pICMRCC,datta2012mreom,demel2013mreom,nooijen2014mreom,huntington2015mreom_orb_select,huntington2015mreom_bench` has been implemented in ORCA. The strength of the MR-EOM-CC methodology lies in its ability to calculate many excited states from a single state-averaged CASSCF solution, for which only a single set of amplitudes needs to be solved and the final transformed Hamiltonian is diagonalized over a small manifold of excited states only through an uncontracted MRCI problem. Hence, a given MR-EOM calculation involves three steps, performed by three separate modules in ORCA : 1. a state-averaged CASSCF calculation (CASSCF module), 2. the solution of amplitude equations and the calculation of the elements of the similarity transformed Hamiltonians (MDCI module), 3. and the uncontracted MRCI diagonalization of the final similarity transformed Hamiltonian (MRCI module). The current implementation allows for MR-EOM-T|T$^{\dagger}$-h-v, MR-EOM-T|T$^{\dagger}$|SXD-h-v and MR-EOM-T|T$^{\dagger}$|SXD|U-h-v calculations. A more detailed description of these methods and the available input parameters will be given in Sec. {ref}`sec:mreom.detailed`. We also note that the theoretical details underlying these methods can be found in Ref. {cite}`huntington2015mreom_bench`. In Sec. {ref}`sec:mreom.detailed`, we will discuss a strategy for the selection of the state-averaged CAS and other steps for setting up an MR-EOM calculation in detail. Furthermore, we will discuss how spin-orbit coupling effects can be included in MR-EOM calculations, a projection scheme to aid with convergence difficulties in the iteration of the $T$ amplitude equations, an orbital selection scheme to reduce the size of the inactive core and virtual subspaces in the calculation of excitation energies and a strategy for obtaining nearly size-consistent results in MR-EOM. The purpose of this section is simply to provide a simple example which illustrates the most basic usage of the MR-EOM implementation in ORCA. (mreom.calc.typical)= ## A Simple MR-EOM Calculation Let us consider an MR-EOM-T|T$^{\dagger}$|SXD|U-h-v calculation on formaldehyde. An MR-EOM-T|T$^{\dagger}$|SXD|U-h-v calculation is specified via the `MR-EOM` keyword along with the specification of a state-averaged CASSCF calculation (i.e. CASSCF(nel, norb) calculation with the number of roots of each multiplicity to be included in the state-averaging for the reference state) and the number of desired roots in each multiplicity block for the final MRCI diagonalization. We note that the CASSCF module is described in sections {ref}`sec:energygradients.casscf.typical` and {ref}`sec:casscf.detailed` and that a description of the MRCI module is given in sections {ref}`sec:mrci.typical` and {ref}`sec:mrci.detailed`. Here, we have a state-averaged CAS(6,4) calculation, comprised of 3 singlets and 3 triplets and we request 6 singlet roots and 6 triplet roots in our final MRCI diagonalization (i.e. the roots to be computed in the MR-EOM-T|T$^{\dagger}$|SXD|U-h-v calculation): ```{literalinclude} ../../orca_working_input/C05S10_195.inp :language: orca ``` One can alternatively perform an MR-EOM-T|T$^{\dagger}$-h-v or MR-EOM-T|T$^{\dagger}$|SXD-h-v calculation by replacing the `MR-EOM` keyword, in the first line of the input above, by `MR-EOM-T|Td` or `MR-EOM-T|Td|SXD`, respectively. Namely, replacing the first line of the input above with ```orca !MR-EOM-T|Td def2-TZVP VeryTightSCF ``` runs the MR-EOM-T|T$^{\dagger}$-h-v calculation, while ```orca !MR-EOM-T|Td|SXD def2-TZVP VeryTightSCF ``` runs the MR-EOM-T|T$^{\dagger}$|SXD-h-v calculation. The final MRCI diagonalization manifold includes 2h1p, 1h1p, 2h, 1h and 1p excitations in MR-EOM-T|T$^{\dagger}$-h-v calculations, 2h, 1p and 1h excitations in MR-EOM-T|T$^{\dagger}$|SXD-h-v calculations and 1h and 1p excitations in MR-EOM-T|T$^{\dagger}$|SXD|U-h-v calculations. Note that in the `%mdci` block, we have set the convergence tolerance (`STol`) for the residual equations for the amplitudes to $10^{-7}$, as this default value is overwritten with the usage of the `TightSCF`, `VeryTightSCF`, etc. keywords. It is always important to inspect the values of the largest $T$, $S$ (here, we use $S$ to denote the entire set of $S$, $X$ and $D$ amplitudes) and $U$ amplitudes. If there are amplitudes that are large (absolute values $> 0.15$), the calculated results should be regarded with suspicion. For the above calculation, we obtain: ``` -------------------- LARGEST T AMPLITUDES -------------------- 8-> 13 8-> 13 0.060331 4-> 17 4-> 17 0.029905 8-> 9 8-> 9 0.028160 8-> 16 8-> 16 0.027266 6-> 20 6-> 20 0.025885 8-> 21 8-> 21 0.025308 4-> 16 4-> 16 0.024803 8-> 12 8-> 12 0.023915 5-> 18 5-> 18 0.023553 8-> 23 8-> 23 0.023384 3-> 16 3-> 16 0.023182 7-> 19 7-> 19 0.023043 8-> 13 4-> 11 0.022010 3-> 19 3-> 19 0.021987 8-> 16 8-> 9 0.021230 8-> 9 8-> 16 0.021230 ``` for the $T$ amplitudes, ``` -------------------- LARGEST S AMPLITUDES -------------------- 4-> 8 8-> 11 0.074048 3-> 8 8-> 9 0.064886 4-> 5 5-> 11 0.045479 3-> 8 8-> 16 0.042657 4-> 7 7-> 11 0.042598 4-> 5 5-> 17 0.042076 4-> 5 8-> 11 0.039958 4-> 8 8-> 17 0.037532 3-> 5 8-> 9 0.035907 4-> 7 7-> 17 0.035767 2-> 6 6-> 19 0.034148 3-> 5 5-> 10 0.033339 2-> 6 6-> 10 0.032691 4-> 6 6-> 11 0.032181 8-> 8 3-> 16 0.031775 2-> 7 7-> 22 0.031238 ``` for the $S$ amplitudes, and ``` -------------------- LARGEST U AMPLITUDES -------------------- 3-> 8 3-> 8 0.026128 3-> 8 3-> 5 0.007683 2-> 8 2-> 8 0.006182 3-> 8 2-> 5 0.006154 2-> 8 3-> 5 0.004954 3-> 5 3-> 5 0.004677 4-> 8 4-> 8 0.003989 3-> 8 2-> 8 0.002040 2-> 8 3-> 8 0.002040 2-> 8 2-> 5 0.001818 4-> 8 4-> 5 0.001173 2-> 5 2-> 5 0.001107 4-> 5 4-> 5 0.000714 3-> 7 3-> 7 0.000607 3-> 6 3-> 6 0.000521 2-> 5 3-> 5 0.000365 ``` for the $U$ amplitudes. Hence, one can see that there are no unusually large amplitudes for this calculation. **We note that there can be convergence issues with the $\boldsymbol{T}$ amplitude iterations and that in such cases, the flag:** ```orca DoSingularPT true ``` **should be added to the `%mdci` block.** The convergence issues are caused by the presence of nearly singular $T_2$ amplitudes and setting the `DoSingularPT` flag to `true` activates a procedure which projects out the offending amplitudes (in each iteration) and replaces them by suitable perturbative amplitudes. For more information, see the examples in section {ref}`sec:mreom.singularpt.detailed`. After the computation of the amplitudes and the elements of the similarity transformed Hamiltonians, within the MDCI module, the calculation enters the MRCI module. For a complete, step by step description of the output of an MRCI calculation, we refer the reader to the example described in section {ref}`sec:mrci.example.typical`. Let us first focus on the results for the singlet states (`CI-BLOCK 1`). Following the convergence of the Davidson diagonalization (default) or DIIS procedure, the following results of the MRCI calculation for the singlet states are printed: ``` ---------- CI-RESULTS ---------- The threshold for printing is 0.30 percent The weights of configurations will be printed. The weights are summed over all CSFs that belong to a given configuration before printing STATE 0: Energy= -114.321368425 Eh RefWeight= 0.9781 0.00 eV 0.0 cm**-1 0.0137 : h---h---[0222] 0.0756 : h---h---[1221] 0.8879 : h---h---[2220] STATE 1: Energy= -114.176866027 Eh RefWeight= 0.9765 3.93 eV 31714.6 cm**-1 0.0039 : h---h---[1122] 0.9726 : h---h---[2121] 0.0071 : h---h 4[1222] 0.0085 : h---h 4[2221] STATE 2: Energy= -113.988050555 Eh RefWeight= 0.9774 9.07 eV 73154.8 cm**-1 0.0044 : h---h---[1212] 0.9730 : h---h---[2211] 0.0063 : h---h 3[1222] 0.0041 : h---h 3[2221] STATE 3: Energy= -113.963862283 Eh RefWeight= 0.8810 9.73 eV 78463.5 cm**-1 0.7459 : h---h---[1221] 0.0807 : h---h---[2022] 0.0533 : h---h---[2220] 0.0228 : h---h 4[2122] 0.0034 : h---h---[1220]p13 0.0072 : h---h---[1220]p18 0.0236 : h---h---[2120]p11 0.0148 : h---h---[2120]p14 0.0069 : h---h---[2120]p17 0.0056 : h---h---[2120]p20 0.0098 : h---h---[2210]p19 STATE 4: Energy= -113.931144468 Eh RefWeight= 0.0003 10.62 eV 85644.3 cm**-1 0.0045 : h---h---[0122]p9 0.0089 : h---h---[1121]p9 0.9333 : h---h---[2120]p9 0.0243 : h---h---[2120]p10 0.0080 : h---h---[2120]p12 0.0113 : h---h---[2120]p16 STATE 5: Energy= -113.929056780 Eh RefWeight= 0.6857 10.68 eV 86102.5 cm**-1 0.0061 : h---h---[0222] 0.0918 : h---h---[1221] 0.5784 : h---h---[2022] 0.0048 : h---h---[2202] 0.0047 : h---h---[2220] 0.2905 : h---h 4[2122] 0.0045 : h---h---[2021]p13 ``` For each state, the total energy is given in $E_\text{h}$; the weight of the reference configurations (`RefWeight`) in the given state is provided, and the energy differences from the lowest lying state are given in eV and cm$^{-1}$. Also, in each case, the weights and a description of the configurations which contribute most strongly to the given state are also provided. See section {ref}`sec:mrci.example.typical` for a discussion of the notation that is used for the description of the various configurations. To avoid confusion, we note that in the literature concerning the MR-EOM methodology {cite}`datta2012mreom,demel2013mreom,nooijen2014mreom,huntington2015mreom_orb_select,huntington2015mreom_bench,liu2015mreom_bench,liu2015mreom_bench_SOC`, the term "%active" is used to denote the reference weight multiplied by 100%. In general, `RefWeight ` should be $> 0.9$, such that the states are dominated by reference space configurations. This criterion is satisfied for the first three states and the reference weight of the fourth state is sufficiently close to $0.9$. **However, the reference weights of the two higher lying states (especially state 4) are too small and these states should be discarded as the resulting energies will be inaccurate (i.e. states with significant contributions from configurations outside the reference space cannot be treated accurately)** . In the case of the triplet states (`CI-BLOCK 2`), we obtain the following results: ``` ---------- CI-RESULTS ---------- The threshold for printing is 0.30 percent The weights of configurations will be printed. The weights are summed over all CSFs that belong to a given configuration before printing STATE 0: Energy= -114.190840989 Eh RefWeight= 0.9693 0.00 eV 0.0 cm**-1 0.9691 : h---h---[2121] 0.0079 : h---h 4[1222] 0.0115 : h---h 4[2221] STATE 1: Energy= -114.106733017 Eh RefWeight= 0.9941 2.29 eV 18459.6 cm**-1 0.9941 : h---h---[1221] STATE 2: Energy= -114.015150051 Eh RefWeight= 0.9787 4.78 eV 38559.7 cm**-1 0.9786 : h---h---[2211] 0.0050 : h---h 3[1222] STATE 3: Energy= -113.939299674 Eh RefWeight= 0.0006 6.84 eV 55206.9 cm**-1 0.0044 : h---h---[0122]p9 0.0084 : h---h---[1121]p9 0.9419 : h---h---[2120]p9 0.0131 : h---h---[2120]p10 0.0043 : h---h---[2120]p12 0.0173 : h---h---[2120]p16 STATE 4: Energy= -113.925571130 Eh RefWeight= 0.4017 7.22 eV 58220.0 cm**-1 0.3863 : h---h---[1122] 0.0154 : h---h---[2121] 0.1722 : h---h 4[1222] 0.4098 : h---h 4[2221] 0.0045 : h---h---[2120]p13 STATE 5: Energy= -113.910479339 Eh RefWeight= 0.0009 7.63 eV 61532.3 cm**-1 0.0088 : h---h---[0122]p10 0.0030 : h---h---[1121]p10 0.0120 : h---h---[2120]p9 0.9408 : h---h---[2120]p10 0.0106 : h---h---[2120]p16 0.0112 : h---h---[2120]p19 ``` Here, we see that the first three states have reference weights which are $> 0.9$, while the reference weights of the final three states are well below that threshold. Hence, the latter three states should be discarded from any meaningful analysis. Following the printing of the CI results for the final CI block, the states are ordered according to increasing energy and the vertical transition energies are printed: ``` ------------------- TRANSITION ENERGIES ------------------- The lowest energy is -114.321368425 Eh State Mult Irrep Root Block mEh eV 1/cm 0 1 -1 0 0 0.000 0.000 0.0 1 3 -1 0 1 130.527 3.552 28647.5 2 1 -1 1 0 144.502 3.932 31714.6 3 3 -1 1 1 214.635 5.841 47107.0 4 3 -1 2 1 306.218 8.333 67207.2 5 1 -1 2 0 333.318 9.070 73154.8 6 1 -1 3 0 357.506 9.728 78463.5 7 3 -1 3 1 382.069 10.397 83854.4 8 1 -1 4 0 390.224 10.619 85644.3 9 1 -1 5 0 392.312 10.675 86102.5 10 3 -1 4 1 395.797 10.770 86867.5 11 3 -1 5 1 410.889 11.181 90179.7 ``` Furthermore, following the generation of the (approximate) densities, the absorption and CD spectra are printed: ``` ========================================== MR-EOM Non Relativistic Properties ========================================== ---------------------------------------------------------------------------------------------------- ABSORPTION SPECTRUM VIA TRANSITION ELECTRIC DIPOLE MOMENTS ---------------------------------------------------------------------------------------------------- Transition Energy Energy Wavelength fosc(D2) D2 DX DY DZ (eV) (cm-1) (nm) (au**2) (au) (au) (au) ---------------------------------------------------------------------------------------------------- 0-1A -> 1-1A 3.932110 31714.6 315.3 0.000000000 0.00000 -0.00000 -0.00000 0.00000 0-1A -> 2-1A 9.070040 73154.8 136.7 0.002137450 0.00962 0.09808 -0.00000 -0.00000 0-1A -> 3-1A 9.728237 78463.5 127.4 0.157495738 0.66081 -0.00000 0.00000 -0.81290 0-1A -> 4-1A 10.618534 85644.3 116.8 0.025353906 0.09746 -0.00000 -0.31218 -0.00000 0-1A -> 5-1A 10.675343 86102.5 116.1 0.024673667 0.09434 -0.00000 -0.00000 0.30715 ------------------------------------------------------------------------------------------ CD SPECTRUM VIA TRANSITION ELECTRIC DIPOLE MOMENTS ------------------------------------------------------------------------------------------ Transition Energy Energy Wavelength R MX MY MZ (eV) (cm-1) (nm) (1e40*cgs) (au) (au) (au) ------------------------------------------------------------------------------------------ 0-1A -> 1-1A 3.932110 31714.6 315.3 -0.00000 0.00000 0.00000 0.56273 0-1A -> 2-1A 9.070040 73154.8 136.7 -0.00000 0.00000 -0.74486 0.00000 0-1A -> 3-1A 9.728237 78463.5 127.4 -0.00000 -0.00000 -0.00000 -0.00000 0-1A -> 4-1A 10.618534 85644.3 116.8 0.00000 0.35898 0.00000 -0.00000 0-1A -> 5-1A 10.675343 86102.5 116.1 0.00000 0.00000 -0.00000 -0.00000 ``` :::{Warning} - It is important to note that the transition moments and oscillator strengths (and state dipole moments) have been blindly computed by the MRCI module and currently, no effort has been made to include the effects of the various similarity transformations in the evaluation of these quantities. Hence these quantities are only approximate and should only be used as a qualitative aid to determine which states are dipole allowed or forbidden. Furthermore, since the calculated densities are approximate, so are the results of the population analysis that are printed before the absorption and CD spectra. - While both the CASSCF and MRCI modules can make use of spatial point-group symmetry to some extent, the MR-EOM implementation is currently limited to calculations in $C_1$ symmetry. ::: (mreom.capability.typical)= ## Capabilities The MR-EOM methodology can be used to calculate a desired number of states for both closed- and open-shell systems from a single state-averaged CASSCF solution. Currently, the approach is limited to serial calculations and to smaller systems in smaller active spaces. One should be aware that in the most cost-effective MR-EOM-T|T$^{\dagger}$|SXD|U-h-v approach (i.e. the smallest diagonalization manifold), an MRCI diagonalization is performed over all 1h and 1p excited configurations out of the CAS, which will inevitably limit the size of the initial CAS which can be used. We have also implemented an orbital selection scheme which can be used to reduce the size of the inactive core and virtual subspaces in the calculation of excitation energies, and this can be employed to extend the applicability of the approach to larger systems. The current implementation can also be used in conjunction with the spin-orbit coupling submodule ({ref}`sec:mrci.general.detailed`) of the MRCI module to calculate spin-orbit coupling effects in MR-EOM calculations to first order. These and other features of the current implementation will be discussed in {ref}`sec:mreom.detailed`. (mreom.perturbative.typical)= ## Perturbative MR-EOM-PT The MR-EOM family of methods now also features an almost fully perturbative approach called MR-EOMPT {cite}`Lechner2021`. This method shares the features of the MR-EOMCC parent method while using non-iterative perturbative estimates for the $\hat T$ and $\hat S,\hat X,\hat D$ amplitudes. This slightly reduces the accuracy compared to iterative MR-EOMCC while reducing runtime. Furthermore, convergence issues due to nearly singular $\hat T$ and $\hat S,\hat X,\hat D$ amplitudes cannot occur anymore. This method can be invoked by adding the keyword `DoMREOM_MRPT True` to the `%mdci` block.