7.40. Multireference Equation of Motion Coupled-Cluster (MR-EOM-CC) Theory¶

In analogy with single reference EOM-CC (see sections Excited States with EOM-CCSD and Excited States via EOM-CCSD) and STEOM-CC (see sections Excited States with STEOM-CCSD and Excited States via STEOM-CCSD), Multireference Equation of Motion Coupled-Cluster (MR-EOM-CC) theory [193, 194, 203, 406, 407, 639] can be viewed a transform and diagonalize approach to molecular electronic structure theory. An MR-EOM calculation involves a single state-averaged CASSCF calculation, incorporating a few low-lying states and the solution of a single set of cluster amplitudes, which define a sequence of similarity-transformed Hamiltonians. The ultimate goal of these many-body transformations is to effectively decouple the CAS configurations from important excited configurations (e.g., 2p2h, 2p1h, 1p1h, etc.) which comprise the first-order interacting space. Through the definition of suitable cluster operators, in each of the transformations, most of these excitations can be included in an internally contracted fashion. Hence, the resulting final transformed Hamiltonian can be diagonalized over a small subspace of the original first-order interacting space to gain access to many electronic states. As discussed in section MR-EOM-CC: Multireference Equation of Motion Coupled-Cluster, the MR-EOM implementation in ORCA therefore makes use of the CASSCF module (to obtain the state-averaged CASSCF reference), the MDCI module for the solution of the amplitude equations and the calculation of the elements of the various similarity transformed Hamiltonians and the MRCI module for the diagonalization of the final transformed Hamiltonian. Some desirable features of this methodology are:

Many states can be obtained through the diagonalization of a similarity transformed Hamiltonian over a compact diagonalization manifold (e.g. the final diagonalization space in MR-EOM-T|T\(^\dagger\)|SXD|U only includes the CAS configurations and 1h and 1p configurations).
Only a single state-averaged CASSCF calculation and the solution of a single set of amplitudes is required to define the final similarity transformed Hamiltonian and the results are typically quite insensitive to the precise definition of the CAS (only a few low-lying multiplets need to be included in the state-averaging)
The MR-EOM approach is rigorously invariant to rotations of the orbitals in the inactive, active and virtual subspaces, and it preserves both spin and spatial symmetry.

Table 7.22 The details of the various MR-EOM transformations that are considered in the ORCA implementation of MR-EOM. The equations for the operator components and the residual equations which determine the corresponding amplitudes also appear in the Table. Note that we use the usual (Einstein) convention that repeated indices are summed over.¶
Name	Transformation	Operators	Operator Components	Residual Equation
T	\(\widehat{\xbar{H}}=e^{-\widehat{T} }\widehat{H}e^{\widehat{T} }\)	\(\widehat{T}=\widehat{T}_1+\widehat{T}_2\)	\(\widehat{T}_1=t_i^a\widehat{E}_i^a\)	\(R^{a}_i=\sum_{\mathfrak{m} }w_{\mathfrak{m} }\left\langle\Phi_{\mathfrak{m} }\right\vert \widehat{E}^{i}_{a}\widehat{\xbar{H}}\left\vert\Phi_{\mathfrak{m} }\right\rangle\)
T	\(=\xbar{h}_0+\xbar{h}^p_q\big\{\widehat{E}^{p}_{q}\big\}+\xbar{h}^{pq}_{rs}\big\{\widehat{E}^{pq}_{rs}\big\}+\ldots\)		\(\widehat{T}_2=\frac{1}{2}t_{ij}^{ab}\widehat{E}_{ij}^{ab}\)	\(R_{ij}^{ab}=\xbar{h}^{ab}_{ij}\)
T\(^{\dagger}\)	\(\widehat{\xbar{H}}_2e^{-\widehat{T}^{\dagger} }\)	\(\widehat{T}^{\dagger}=\widehat{T}^{\dagger}_1+\widehat{T}^{\dagger}_2\)	\(\widehat{T}^{\dagger}_1=t_a^i\widehat{E}_a^i\)	None (i.e. set \(t^{i}_{a}\approx t^a_i\))
T\(^{\dagger}\)	\(=\widetilde{h}_0+\widetilde{h}\,^{p}_{q}\big\{\widehat{E}^{p}_{q}\big\}+\widetilde{h}\,^{pq}_{rs}\big\{\widehat{E}^{pq}_{rs}\big\}+\ldots\)		\(\widehat{T}^{\dagger}_2=\frac{1}{2}t_{ab}^{ij}\widehat{E}_{ab}^{ij}\)	None (i.e. set \(t^{ij}_{ab}\approx t^{ab}_{ij}\))
SXD	\(\widehat{\xbar{F}}=\big\{e^{\widehat{S}+\widehat{X}+\widehat{D} }\big\}^{-1}\widehat{\mathcal{H} }_2\big\{e^{\widehat{S}+\widehat{X}+\widehat{D} }\big\}\)	\(\widehat{S}=\widehat{S}_2\)	\(\widehat{S}_2=s_{i'j'}^{aw}\widehat{E}_{i'j'}^{aw}\)	\(R^{aw}_{i'j'}=\xbar{f}\,^{aw}_{i'j'}\)
	\(=\xbar{f}_0+\xbar{f}\,^{p}_{q}\big\{\widehat{E}^{p}_{q}\big\}+\xbar{f}\,^{pq}_{rs}\big\{\widehat{E}^{pq}_{rs}\big\}+\ldots\)	\(\widehat{X}=\widehat{X}_2\)	\(\widehat{X}_2=x_{i'x}^{aw}\widehat{E}_{i'x}^{aw}\)	\(R^{aw}_{i'x}=\xbar{f},^{aw}_{i'x}\)
		\(\widehat{D}=\widehat{D}_2\)	\(\widehat{D}_2=d_{xi'}^{aw}\widehat{E}_{xi'}^{aw}\)	\(R^{aw}_{xi'}=\xbar{f}\,^{aw}_{xi'}\)
U	\(\widehat{G}=e^{-\widehat{U} }\widehat{\xbar{F}}_2e^{\widehat{U} }\)	\(\widehat{U}=\widehat{U}_2\)	\(\widehat{U}_2=u_{i'j'}^{wx}\widehat{E}_{i'j'}^{wx}\)	\(R^{wx}_{i'j'}=g^{wx}_{i'j'}\)
	\(=g_0+g^{p}_{q}\big\{\widehat{E}^{p}_{q}\big\}+g^{pq}_{rs}\big\{\widehat{E}^{pq}_{rs}\big\}+\ldots\)

As the details concerning the MR-EOM methodology are rather involved, we refer the interested reader to Refs. [193, 194, 203, 406, 407, 639] for a more detailed discussion. Note that the details concerning the implementation of MR-EOM in ORCA can be found in Refs. [406] and [407]. In the following discussion, we note that general spatial orbitals \(p,q,r,s\), which comprise the molecular orbital basis, are partitioned into (doubly occupied) inactive core orbitals \(i',j',k',l'\), occupied orbitals \(i,j,k,l\) (i.e. the union of the inactive core and active orbital subspaces), active orbitals \(w,x,y,z\) and virtual orbitals \(a,b,c,d\). In general, the many-body similarity transformations assume the general form

\[\begin{aligned} \widehat{G}=\big\{e^{\widehat{Y} }\big\}^{-1}\widehat{\xbar{H}}_{2}\big\{e^{\widehat{Y} }\big\}=g_0+\sum_{p,q}g^{p}_{q}\big\{\widehat{E}^{p}_{q}\big\}+\sum_{p,q,r,s}g^{pq}_{rs}\big\{\widehat{E}^{pq}_{rs}\big\}+\ldots,\end{aligned}\]

in which \(\widehat{Y}\) is a cluster operator and \(\widehat{\xbar{H}}_2\) is the bare Hamiltonian or a similarity transformed Hamiltonian truncated up to two-body operators. The braces indicate Kutzelnigg-Mukherjee normal ordering [485, 596], which is used extensively in the definition of the MR-EOM formalism. The various transformations which need to be considered in the ORCA implementation of MR-EOM are summarized in Table Table 7.22. The table also includes the expressions for the operator components of the various internally contracted cluster operators and the residual equations that must be solved for the various amplitudes. Note that the residual equations are typically of the many-body type (i.e. obtained by setting the corresponding elements of the similarity transformed Hamiltonian to zero). The only exception is the residual equation which defines the \(t_i^a\) amplitudes, which is a projected equation of the form

\[R_i^a=\sum_{\mathfrak{m} }w_{\mathfrak{m} }\big\langle\Phi_{\mathfrak{m} }\vert \widehat{E}^{i}_a\widehat{H}\vert\Phi_{\mathfrak{m} }\big\rangle.\]

Here, \(\vert \Phi_{\mathfrak{m} }\big\rangle\) is the \(\mathfrak{m}^{\text{th} }\) state included in the state averaged CAS, with weight \(w_{\mathfrak{m} }\). The reason the equation for the singles is of the projected form is that it satisfies the Brillouin theorem (i.e. the first order singles vanish for all \(i\) and \(a\)), whereas the corresponding many-body equation (\(\bar{h}_i^a=0)\) does not.

Table 7.23 Details of the three MR-EOM approaches implemented in ORCA¶
Method	Input Keyword	Operators	Diagonalization Manifold
MR-EOM-T\|T\(^\dagger\)-h-v	`MR-EOM-T\|Td`	\(\widehat{T}\); \(\widehat{T}^\dagger\)	CAS, 2h1p, 1h1p, 2h, 1h, 1p
MR-EOM-T\|T\(^\dagger\)\|SXD-h-v	`MR-EOM-T\|Td\|SXD`	\(\widehat{T}\); \(\widehat{T}^\dagger\); \(\widehat{S}+\widehat{X}+\widehat{D}\)	CAS, 2h, 1h, 1p
MR-EOM-T\|T\(^\dagger\)\|SXD\|U-h-v	`MR-EOM`	\(\widehat{T}\); \(\widehat{T}^\dagger\); \(\widehat{S}+\widehat{X}+\widehat{D}\); \(\widehat{U}\)	CAS, 1h, 1p

Note that there are three different MR-EOM approaches which have been implemented in ORCA. Namely, 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. At this point it is useful to discuss the naming convention used for these approaches. We use a vertical line to separate each transformation involved in the sequence of transformations defining the given MR-EOM approach. For example, T|T\(^\dagger\)|SXD indicates that a T transformation, is followed by a T\(^\dagger\) transformation, which is then followed by an SXD transformation. The h-v indicates that the elements of the transformed Hamiltonian have been hermitized (h) and vertex symmetrized (v) before entering the MRCI diagonalization (see Ref. [407] for more information). Essentially, this means that the full eightfold symmetry of the two-electron integrals (and hermiticity of the one-body elements) have been enforced upon the elements of the transformed Hamiltonian. The details of the three MR-EOM approaches are summarized in Table Table 7.23. This table includes the keyword (in the first line of input) used to initiate the calculation in ORCA, the various operators involved, and the configurations included in the final diagonalization manifold. One can clearly see that the MR-EOM-T|T\(^\dagger\)|SXD|U-h-v approach is the most cost effective, as it only includes the 1h and 1p configurations, beyond the CAS, in the final diagonalization manifold.

The various %mdci keywords, which are important for controlling MR-EOM calculations are (i.e. default values are given here):

%mdci
  STol              1e-7        #Convergence Tolerance on Residual Equations
  MaxIter           100         #Maximum Number of Iterations
  DoSingularPT      false       #Activate the Singular PT/Projection Procedure
  SingularPTThresh  0.01        #Threshold for the Singular PT/Projection
                                #Procedure
  PrintOrbSelect    false       #Print the Eigenvalues of the Orbital Selection
                                #Densities (and R_core and R_virt values)
                                #and Terminate the Calculation
  CoreThresh        0.0         #Core Orbital Selection Threshold
  VirtualThresh     1.0         #Virtual Orbital Selection Threshold
end

As discussed below, the orbital selection scheme is activated by adding !OrbitalSelection to the simple input line. Keywords that are specific to the CASSCF and MRCI modules are discussed in sections The Complete Active Space Self-Consistent Field (CASSCF) Module and The Multireference Correlation Module, respectively. We note that in MR-EOM-T|T\(^\dagger\)-h-v and MR-EOM-T|T\(^\dagger\)|SXD-h-v calculations, it is possible to override the default excitation classes in the final MRCI diagonalization. This is done by specifying excitations none and then explicitly setting the excitation flags within a given multiplicity block. For example, if we wanted to have 1h, 1p, 1h1p, 2h and 2h1p excitations in the final diagonalization manifold, we would specify (i.e. here we have requested 6 singlets and have a CAS(6,4) reference):

%mrci
    newblock 1 *
        excitations none
        Flags[is]   1
        Flags[sa]   1
        Flags[ia]   1
        Flags[ijss] 1
        Flags[ijsa] 1
        nroots      6
        refs
            cas(6,4)
        end
    end
end

7.40.1. The Steps Required to Run an MR-EOM Calculation¶

To illustrate the various steps required in a typical MR-EOM calculation, we will consider the calculation of the excitation energies of the neutral Fe atom at the MR-EOM-T|T\(^{\dagger}\)|SXD|U-h-v level of theory.

7.40.1.1. State-Averaged CASSCF Calculation¶

Evidently, the first step is to determine a suitable state-averaged CASSCF reference for the subsequent MR-EOM calculation. In choosing the state-averaged CAS for an MR-EOM calculation, we typically include a few of the low-lying multiplets that have the same character as the (much larger number of) states that we wish to compute in the final MR-EOM calculation. For the neutral Fe atom, we typically have electronic states which have either 4s\(^2\)3d\(^6\) character or 4s\(^1\)3d\(^7\) character. From the NIST atomic spectra database [1, 605], we find that the lowest lying a\(^5\)D multiplet is of 4s\(^2\)3d\(^6\) character and the higher lying a\(^5\)F multiplet is of 4s\(^1\)3d\(^7\) character. Hence, we can set up a state-averaged CASSCF(8,6) calculation (i.e. 8 electrons in 6 orbitals (4s and 3d)) which includes the \(^5\)D and \(^5\)F states and choose the weights such that the average occupation of the 4s orbital is 1.5. As discussed in Ref. [529], this is done to avoid a preference toward either of the 4s configurations in the state-averaging. We will run the state-averaged CASSCF calculation, making use of the second order DKH (see The Douglas-Kroll-Hess Method) method for the inclusion of relativistic effects in a Def2-TZVPP basis (i.e. the DKH-Def2-TZVPP relativistically recontracted basis, listed in section Choice of Basis Set). The input file for the state-averaged CASSCF(8,6) calculation takes the form:

!CASSCF DKH-Def2-TZVPP VeryTightSCF DKH

%casscf
  nel 8
  norb 6
  mult 5
  nroots 12
  weights[0] = 0.7, 0.7, 0.7, 0.7, 0.7, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5
end

* xyz 0 5
  Fe 0.000000 0.000000 0.000000
end

Here, we have requested 12 quintet states (the lowest lying \(^5\)D and \(^5\)F multiplets) and we have chosen the weights to be 0.7 for the five \(^5\)D states and 0.5 for the seven \(^5\)F states, such that the overall occupation of the 4s orbital will be 1.5. Once the calculation has converged, it is important to inspect the results printed in the final macro-iteration of the CASSCF calculation (macro-iteration 8 in this case). In this case, we have:

MACRO-ITERATION   8:
   --- Inactive Energy E0 = -1249.82392764 Eh
   --- All densities will be recomputed
CI-ITERATION   0:
 -1271.258899450   0.000000000001 (    0.00)
 -1271.258899450   0.000000000001
 -1271.258899450   0.000000000001
 -1271.258899450   0.000000000001
 -1271.258899450   0.000000000001
 -1271.186288591   0.000000000001
 -1271.186288591   0.000000000001
 -1271.186288591   0.000000000002
 -1271.186288591   0.000000000001
 -1271.186288591   0.000000000001
 -1271.186288591   0.000000000001
 -1271.186288591   0.000000000001
   CI-PROBLEM SOLVED
   DENSITIES MADE
   E(CAS)= -1271.222594021 Eh DE=   -1.591616e-12
   --- Energy gap subspaces: Ext-Act = 0.276   Act-Int = 2.469
   N(occ)=  1.50000 1.30000 1.30000 1.30000 1.30000 1.30000
   ||g|| =     4.669702e-05 Max(G)=    1.104493e-05 Rot=68,1

Directly below CI-ITERATION 0, the final CAS-CI energies are printed, and one observes that they follow the correct degeneracy pattern (i.e. 5 states with energy -1271.258899450 and 7 states with energy -1271.186288591). Furthermore, the final state-averaged CASSCF energy (E(CAS)= -1271.222594021) and occupation numbers (N(occ)= 1.50000 1.30000 1.30000 1.30000 1.30000 1.30000) are also printed. As expected, the occupation number of the 4s orbital is indeed 1.5, while the 3d orbitals each have an occupation of 1.3.

7.40.1.2. Selection of the States to Include in the MR-EOM Calculation¶

Once a satisfactory CASSCF reference has been obtained, the next step is to determine the number of states to include in the MR-EOM calculation. From the NIST atomic spectra database, one finds that the higher lying states of 4s\(^2\)3d\(^6\) and 4s\(^1\)3d\(^7\) character are either singlets, triplets, or quintets. To figure out how many states should be included in each multiplicity block, one can perform an inexpensive CAS-CI calculation. This is done by reading in the orbitals from the previous CASSCF calculation (here they are stored in the file CAS.gbw) and requesting a single iteration (i.e. using the NoIter keyword) of a state-averaged CASSCF calculation:

!CASSCF DKH-Def2-TZVPP ExtremeSCF DKH NoIter

!MOREAD
%moinp "CAS.gbw"

%casscf
  nel 8
  norb 6
  mult 5,3,1
  nroots 15,90,55
end

* xyz 0 5
  Fe 0.000000 0.000000 0.000000
end

Here, after some experimentation, we have chosen 15 quintets, 90 triplets and 55 singlets. It is important that we calculate states up to sufficiently high energy (i.e. all the states that are of interest) and it is imperative that we have complete multiplets. Hence, several iterations of this procedure might be required to choose the proper number of states for each multiplet. The relevant section of the output file which should be analyzed is the SA-CASSCF TRANSITION ENERGIES. For the above calculation, we obtain (i.e. only the CAS-CI energies for the first 33 roots are shown here):

-----------------------------
SA-CASSCF TRANSITION ENERGIES
------------------------------

LOWEST ROOT (ROOT 0 ,MULT 5) =  -1271.258899450 Eh -34592.713 eV

STATE   ROOT MULT  DE/a.u.    DE/eV    DE/cm**-1
   1    5   0.000000     0.000      0.0
   2    5   0.000000     0.000      0.0
   3    5   0.000000     0.000      0.0
   4    5   0.000000     0.000      0.0
   5    5   0.072611     1.976  15936.2
   6    5   0.072611     1.976  15936.2
   7    5   0.072611     1.976  15936.2
   8    5   0.072611     1.976  15936.2
   9    5   0.072611     1.976  15936.2
  10    5   0.072611     1.976  15936.2
  11    5   0.072611     1.976  15936.2
   0    3   0.092859     2.527  20380.2
   1    3   0.092859     2.527  20380.2
   2    3   0.092859     2.527  20380.2
   3    3   0.092859     2.527  20380.2
   4    3   0.092859     2.527  20380.2
   5    3   0.092859     2.527  20380.2
   6    3   0.092859     2.527  20380.2
   7    3   0.092859     2.527  20380.2
   8    3   0.092859     2.527  20380.2
   9    3   0.092859     2.527  20380.2
  10    3   0.092859     2.527  20380.2
  11    3   0.101848     2.771  22353.1
  12    3   0.101848     2.771  22353.1
  13    3   0.101848     2.771  22353.1
  14    3   0.101848     2.771  22353.1
  15    3   0.101848     2.771  22353.1
  16    3   0.101848     2.771  22353.1
  17    3   0.101848     2.771  22353.1
  18    3   0.102559     2.791  22509.1
  19    3   0.102559     2.791  22509.1
  20    3   0.102559     2.791  22509.1

7.40.1.3. Running the MR-EOM Calculation¶

Now that we have chosen a suitable CASSCF reference and the states that we wish to calculate, we can finally proceed with the MR-EOM calculation. The following input file runs an MR-EOM-T|T\(^\dagger\)|SXD|U-h-v calculation for 15 quintet, 90 triplet and 55 singlet states of the neutral Fe atom (i.e. the CASSCF orbitals are read from CAS.gbw):

!MR-EOM DKH-Def2-TZVPP VeryTightSCF DKH

!MOREAD
%moinp "CAS.gbw"

%method frozencore fc_ewin end

%casscf
  nel 8
  norb 6
  mult 5
  nroots 12
  weights[0] = 0.7, 0.7, 0.7, 0.7, 0.7, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5
end

%mdci
  ewin -6, 10000
  STol 1e-7
end

%mrci
  ewin -6, 10000
  MaxIter 200
  newblock 5 *
    nroots 15
    refs cas(8,6) end
  end
  newblock 3 *
    nroots 90
    refs cas(8,6) end
  end
  newblock 1 *
    nroots 55
    refs cas(8,6) end
  end
end

* xyz 0 5
  Fe 0.000000 0.000000 0.000000
end

Note that since the default frozen core settings exclude the 3p orbitals from the correlation treatment, we have used an energy window (i.e. the line ewin -6, 10000 in both the %mdci and %mrci blocks) such that they are included in the current calculation. We note that a detailed discussion of the input and output of an MR-EOM calculation has already been given in section MR-EOM-CC: Multireference Equation of Motion Coupled-Cluster and thus, we do not repeat it here. It is important to reiterate that one should always inspect the values of the largest (T, S and U) amplitudes. Ideally, the largest amplitudes should be smaller than 0.1 and should not exceed 0.15. If some amplitudes are larger than 0.15, it might be necessary to revisit the definition of the CAS and the weights used. For the T amplitudes, an alternative solution is to use the projection/singular PT scheme discussed in section A Projection/Singular PT Scheme to Overcome Convergence Issues in the T Amplitude Iterations below.

As discussed in section MR-EOM-CC: Multireference Equation of Motion Coupled-Cluster, the excitation energies are printed under the heading TRANSITION ENERGIES. For the current calculation, we obtain the following results (only the results for 33 states are shown here):

-------------------
TRANSITION ENERGIES
-------------------

The lowest energy is  -1271.833871822 Eh

State Mult Irrep Root Block    mEh          eV      1/cm
  5    -1    0    0      0.000      0.000       0.0
  5    -1    1    0      0.000      0.000       0.0
  5    -1    2    0      0.000      0.000       0.0
  5    -1    3    0      0.000      0.000       0.0
  5    -1    4    0      0.000      0.000       0.0
  5    -1    5    0     33.901      0.922    7440.4
  5    -1    6    0     33.901      0.922    7440.4
  5    -1    7    0     33.901      0.922    7440.4
  5    -1    8    0     33.901      0.922    7440.4
  5    -1    9    0     33.901      0.922    7440.4
  5    -1   10    0     33.901      0.922    7440.4
  5    -1   11    0     33.901      0.922    7440.4
  3    -1    0    1     54.743      1.490   12014.8
  3    -1    1    1     54.743      1.490   12014.8
  3    -1    2    1     54.743      1.490   12014.8
  3    -1    3    1     54.743      1.490   12014.8
  3    -1    4    1     54.743      1.490   12014.8
  3    -1    5    1     54.743      1.490   12014.8
  3    -1    6    1     54.743      1.490   12014.8
  5    -1   12    0     78.790      2.144   17292.5
  5    -1   13    0     78.790      2.144   17292.5
  5    -1   14    0     78.790      2.144   17292.5
  3    -1    7    1     95.413      2.596   20940.8
  3    -1    8    1     95.413      2.596   20940.8
  3    -1    9    1     95.413      2.596   20940.8
  3    -1   10    1     95.413      2.596   20940.8
  3    -1   11    1     95.413      2.596   20940.8
  3    -1   12    1     95.413      2.596   20940.8
  3    -1   13    1     95.413      2.596   20940.8
  3    -1   14    1     95.413      2.596   20940.8
  3    -1   15    1     95.413      2.596   20940.8
  3    -1   16    1     95.413      2.596   20940.8
  3    -1   17    1     95.413      2.596   20940.8

It is also important to recall that one should always inspect the reference weights for each state, as only states which are dominated by reference space configurations can be treated accurately at the MR-EOM level of theory. Generally, the reference weights should be larger than (or close to) 0.9. In each multiplicity block, the individual state energies and reference weights can be found following convergence of the MRCI procedure, under the heading CI-RESULTS (see section MR-EOM-CC: Multireference Equation of Motion Coupled-Cluster for a more detailed discussion).

7.40.2. Approximate Inclusion of Spin-Orbit Coupling Effects in MR-EOM Calculations¶

The effects of spin-orbit coupling can approximately be included in MR-EOM calculations using the SOC submodule of the MRCI module, as outlined in section Properties Calculation Using the SOC Submodule. This can be viewed as a first order approximation to the inclusion of spin-orbit coupling effects in MR-EOM. In a more rigorous formulation, one would have to consider the various similarity transformations of the spin-orbit coupling operator. The details of the SOC submodule of the MRCI module have already been discussed in detail in Sec. Properties Calculation Using the SOC Submodule and its usage within the MR-EOM formalism is identical to that discussed therein. Let us consider the calculation of spin-orbit coupling effects in the excitation spectrum of the neutral Fe atom considered in the previous section. The input file for this calculation is:

!MR-EOM DKH-Def2-TZVPP ExtremeSCF DKH

%method frozencore fc_ewin end

%casscf
  nel 8
  norb 6
  mult 5
  nroots 12
  weights[0] = 0.7, 0.7, 0.7, 0.7, 0.7, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5
  etol 1e-12
  gtol 1e-12
end

%mdci
  ewin -6, 10000
  MaxIter 300
  STol 1e-12
end

%mrci
  ewin -6, 10000
  MaxIter 200
  newblock 5 *
    nroots 15
    refs cas(8,6) end
  end
  newblock 3 *
    nroots 90
    refs cas(8,6) end
  end
  newblock 1 *
    nroots 55
    refs cas(8,6) end
  end
  soc
    dosoc true  #include spin-orbit coupling effects
  end
end

* xyz 0 5
  Fe 0.000000 0.000000 0.000000
end

In contrast with the calculation performed in section The Steps Required to Run an MR-EOM Calculation, the convergence thresholds have been tightened in all aspects of the calculation (i.e. the use of the ExtremeSCF keyword, etol and gtol (CASSCF energy and orbital gradient convergence tolerance) are set to \(1 \times 10^{-12}\) and the convergence tolerance for the residuals in the MR-EOM amplitude iterations have been set to \(1 \times 10^{-12}\)). We note that with the use of the ExtremeSCF keyword, the convergence tolerance on the energy (Etol) and residual (Rtol) in the MRCI portion of the calculation are also set to \(1 \times 10^{-12}\). Although it is not absolutely necessary, we have used very strict convergence thresholds to preserve the degeneracies of the various multiplets as much as possible. The output of spin-orbit corrected MR-EOM spectrum appears under the heading SOC CORRECTED ABSORPTION SPECTRUM VIA TRANSITION DIPOLE MOMENTS:

--------------------------------------------------------------------------------------------------------
      SOC CORRECTED ABSORPTION SPECTRUM VIA TRANSITION ELECTRIC DIPOLE MOMENTS    
--------------------------------------------------------------------------------------------------------
      Transition         Energy     Energy  Wavelength fosc(D2)      D2       |DX|      |DY|      |DZ|  
                          (eV)      (cm-1)    (nm)  (*population)  (au**2)    (au)      (au)      (au)  
--------------------------------------------------------------------------------------------------------

It is possible to obtain more accurate results by performing an MR-EOM-T|T\(^\dagger\)|SXD-h-v calculation and including the 1h1p excitations. It is important to note that these calculations are significantly more expensive. As discussed above, to run an MR-EOM-T|T\(^\dagger\)|SXD-h-v calculation, the keyword MR-EOM-T|Td|SXD must appear in the first line of input and, in order to activate the 1h1p excitations in each multiplicity block of the MRCI calculation, the %mrci block takes the form:

%mrci
  ewin -6, 10000
  MaxIter 200
  newblock 5 *
    nroots 15
    excitations none
    Flags[is] 1
    Flags[sa] 1
    Flags[ia] 1
    Flags[ijss] 1
    refs cas(8,6) end
  end
  newblock 3 *
    nroots 90
    excitations none
    Flags[is] 1
    Flags[sa] 1
    Flags[ia] 1
    Flags[ijss] 1
    refs cas(8,6) end
  end
  newblock 1 *
    nroots 55
    excitations none
    Flags[is] 1
    Flags[sa] 1
    Flags[ia] 1
    Flags[ijss] 1
    refs cas(8,6) end
  end
  soc
    dosoc true
  end
end

We use excitations none to set the default excitation flags to false and then manually set the 1h (Flags[is]), 1p (Flags[sa]), 1h1p (Flags[ia]) and 2h (Flags[ijss]) excitation flags to true.

WARNINGS

Currently, MR-EOM-T|T\(^\dagger\)|SXD|U-h-v calculations can only be run with the default excitation classes in the final MRCI (i.e. 1h and 1p). Any other input options for the excitation flags will automatically be overwritten and set to the default values.
Only the inclusion of spin-orbit coupling effects has been tested for MR-EOM calculations. Other features that are available in the MRCI module (e.g. spin-spin coupling, magnetic property calculations, etc.) have not been tested for use within MR-EOM calculations.

7.40.3. A Projection/Singular PT Scheme to Overcome Convergence Issues in the T Amplitude Iterations¶

In certain cases, there may be nearly singular T\(_2\) amplitudes (often, but not always large in magnitude), which can cause convergence issues in the solution of the T amplitude equations. Hence, it is sometimes necessary to discard some of the amplitudes to remedy these convergence problems. The nearly singular T\(_2\) amplitudes are of the form \(t^{ab}_{wx}\), where (\(w,x\)) is a pair of active orbitals which corresponds to a small eigenvalue (pair occupation number \(n_{wx}\)) of the two-body reduced density matrix (RDM). When nearly singular amplitudes are present, it is possible to employ a singular PT/projection scheme (i.e. Scheme I described in Ref. [193]), using the two-body RDM as the metric matrix, to discard these nearly singular amplitudes and replace them with suitable perturbative estimates. As a first example, let us consider the following calculation on the cyclopentadiene molecule:

!MR-EOM def2-SVP VeryTightSCF

%casscf
  nel 4
  norb 4
  nroots 2
  mult 3
end

%mdci
  STol 1e-7
  MaxIter 60
end

%mrci
  newblock 1 *
    nroots 3
    refs cas(4,4) end
  end
  newblock 3 *
    nroots 3
    refs cas(4,4) end
  end
end

* xyz 0 3
  H       -0.879859        0.000000        1.874608
  H        0.879859        0.000000        1.874608
  H        0.000000        2.211693        0.612518
  H        0.000000       -2.211693        0.612518
  H        0.000000        1.349811       -1.886050
  H        0.000000       -1.349811       -1.886050
  C        0.000000        0.000000        1.215652
  C        0.000000       -1.177731        0.285415
  C        0.000000        1.177731        0.285415
  C        0.000000       -0.732372       -0.993420
  C        0.000000        0.732372       -0.993420
*

The T amplitude iterations do not converge after 60 iterations and show no signs of convergence (i.e. final largest residual of 0.000458135 and oscillatory behavior over a significant portion of the iterations). If we inspect the largest T amplitudes,

--------------------
LARGEST T AMPLITUDES
--------------------
  19-> 24  19-> 24       0.043103
  19-> 23  19-> 23       0.031162
  11-> 25  11-> 25       0.028458
  19-> 41  19-> 41       0.027950
  11-> 47  11-> 47       0.027026
  19-> 22  19-> 22       0.025163
  19-> 21  19-> 21       0.022167
  15-> 26  15-> 26       0.022084
  11-> 47  11-> 25       0.022033
  11-> 25  11-> 47       0.022033
  19-> 24  19-> 29       0.021769
  19-> 29  19-> 24       0.021769
  19-> 36  19-> 36       0.020986
  17-> 38  17-> 38       0.019743
  19-> 41  16-> 36       0.019107
  18-> 40  18-> 40       0.017949

one can see that there are no unusually large amplitudes. If we turn on the singular PT/projection scheme by adding the line DoSingularPT true to the %mdci block:

%mdci
  STol 1e-7
  MaxIter 60
  DoSingularPT true
end

and rerun the calculation, we find that the T amplitude iterations now successfully converge in 22 iterations. If we look at the largest T amplitudes:

--------------------
LARGEST T AMPLITUDES
--------------------
  11-> 25  11-> 25       0.028440
  11-> 47  11-> 47       0.027027
  15-> 26  15-> 26       0.022069
  11-> 47  11-> 25       0.022031
  11-> 25  11-> 47       0.022031
  19-> 41  19-> 41       0.020463
  17-> 38  17-> 38       0.018288
  11-> 43  11-> 43       0.017250
  11-> 39  11-> 39       0.016838
  15-> 27  15-> 27       0.016001
  13-> 26  13-> 26       0.015985
  16-> 36  16-> 36       0.015759
  19-> 41  16-> 36       0.015697
  18-> 40  18-> 40       0.015376
  17-> 31  17-> 31       0.015074
  18-> 40  17-> 38       0.014470

most of the amplitudes corresponding to the active pair \((w,x)=(19,19)\) no longer appear in the list (i.e. they are nearly singular amplitudes which have been projected out). The only one that does appear in the list, corresponds to a projected perturbative estimate (e.g. 19-> 41 19-> 41 0.020463).

By default, when the singular PT/projection scheme is active, the amplitudes \(t^{ab}_{wx}\) for which the pair occupation numbers satisfy \(n_{wx}<0.01\) (i.e. SingularPTThresh = 0.01), are replaced by perturbative amplitudes in the procedure. However, in some cases, it might be necessary to increase the SingularPTThresh threshold beyond the default value to achieve convergence. One such example is the ferrocene molecule. Consider the following calculation:

!MR-EOM def2-SVP

%casscf
  nel 6
  norb 5
  mult  1,3
  nroots 5,6
end

%mdci
  DoSingularPT true
  MaxIter 50
end

%mrci
  newblock 1 *
    nroots 18
    refs cas(6,5) end
  end
  newblock 3 *
    nroots 10
    refs cas(6,5) end
  end
end

* xyz 0 1
  Fe 0.000000  0.000000  0.000000
  C  0.000000  1.220080  1.650626
  C -1.160365  0.377025  1.650626
  C -0.717145 -0.987065  1.650626
  C  0.717145 -0.987065  1.650626
  C  1.160365  0.377025  1.650626
  C  0.000000  1.220080 -1.650626
  C  1.160365  0.377025 -1.650626
  C  0.717145 -0.987065 -1.650626
  C -0.717145 -0.987065 -1.650626
  C -1.160365  0.377025 -1.650626
  H  0.000000  2.306051  1.635648
  H -2.193184  0.712609  1.635648
  H -1.355463 -1.865634  1.635648
  H  1.355463 -1.865634  1.635648
  H  2.193184  0.712609  1.635648
  H  0.000000  2.306051 -1.635648
  H  2.193184  0.712609 -1.635648
  H  1.355463 -1.865634 -1.635648
  H -1.355463 -1.865634 -1.635648
  H -2.193184  0.712609 -1.635648
end

The T amplitude iterations do not converge after 50 iterations, even though the singular PT/projection scheme is activated. If we increase SingularPTThresh to 0.05 by adding SingularPTThresh 0.05 to the %mdci block:

%mdci
  DoSingularPT true
  SingularPTThresh 0.05
  MaxIter 50
end

the T amplitude iterations successfully converge in 25 iterations.

In conclusion, it occasionally happens that the T amplitude iterations do not converge. In these cases, the singular PT/projection scheme can be activated (DoSingularPT true) to overcome these convergence difficulties. Sometimes, like in the case of ferrocene, it is necessary to adjust the threshold for the singular PT/projection procedure (SingularPTThresh) to achieve convergence. If the procedure still fails with larger values of the threshold, then it might be necessary to revisit the definition of the state-averaged CAS.

7.40.4. An Orbital Selection Scheme for More Efficient Calculations of Excitation Spectra with MR-EOM¶

As described in Ref. [406], the MR-EOM implementation in ORCA can make use of a sophisticated scheme to discard inactive and virtual orbitals, which are not important for the description of the excited states of interest. The selection of inactive core orbitals is based on the eigenvalues of the core orbital selection density

(7.250)¶\[D_{i'j'}=D^{t}_{i'j'}+\frac{\text{Tr}\left( D^t \right)}{\text{Tr}\left( D^s \right)+\text{Tr}\left( D^u \right)}\left( D^s_{i'j'}+D^u_{i'j'} \right),\]

in which

\[\begin{split}\begin{aligned} D^t_{i'j'}&=\sum_{w,a,b}t_{i'w}^{ab^{(1) }}\left( 2t_{j'w}^{ab^{(1) }}-t_{j'w}^{ba^{(1) }} \right), \\ D^s_{i'j'}&=\sum_{k,w,a}\Big[s_{i'k}^{aw^{(1) }}\left(2_{j'k}^{aw^{(1) }}-s_{kj'}^{aw^{(1) }}\right)+s_{ki'}^{aw^{(1) }}\left(2s_{kj'}^{aw^{(1) }}-s_{j'k}^{aw^{(1) }}\right) \Big], \\ D^u_{i'j'}&=\sum_{k',w,x}u_{i'k'}^{wx^{(1) }}\left( 2u_{j'k'}^{wx^{(1) }}-u_{k'j'}^{wx^{(1) }} \right),\end{aligned}\end{split}\]

are respectively, the contributions from the first order \(t_{i'w}^{ab^{(1) }}\), \(s_{i'k}^{aw^{(1) }}\) and \(u_{i'k'}^{wx^{(1) }}\) amplitudes (i.e. note that all amplitudes have at least one active label). Similarly, the selection of virtual orbitals is based upon the eigenvalues of the virtual orbital selection density

(7.251)¶\[\rho_{ab}=\rho^t_{ab}+\frac{\text{Tr}\left( \rho^t \right)}{\text{Tr}\left( \rho^s \right)} \rho^s_{ab},\]

in which, the contribution \(\rho^t\), from the first order \(T_2\) amplitudes, is given by

(7.252)¶\[\rho^t_{ab}=\sum_{k,w,c}t_{wk}^{ac^{(1) }}\left( 2t_{wk}^{bc^{(1) }}-t_{wk}^{cb^{(1) }}\right) +\sum_{i',w,c}t_{i'w}^{ac^{(1) }}\left( 2t_{i'w}^{bc^{(1) }}-t_{i'w}^{cb^{(1) }} \right),\]

and the contribution \(\rho^s\), from the first order \(S_2\) amplitudes, is given by

(7.253)¶\[\rho^s_{ab}=\sum_{i',k,w}s_{i'k}^{aw^{(1) }}\left( 2s_{i'k}^{bw^{(1) }}-s_{ki'}^{bw^{(1) }} \right)+\sum_{i',w,x}s_{xi'}^{aw^{(1) }}\left( 2s_{xi'}^{bw^{(1) }}-s_{i'x}^{bw^{(1) }} \right).\]

Diagonalization of the core orbital selection density \(D_{i'j'}\) and the virtual orbital selection density \(\rho_{ab}\) then yields two respective sets of eigenvalues \(\{\lambda_{i'}\}\) and \(\{\lambda_a\}\). We have found it useful to compute the ratios,

(7.254)¶\[\mathcal{R}_{\text{core} }=\frac{\sum_{i'=0}^{n_{\text{core} }^{\text{excluded} }} \lambda_{i'} }{\sum_{i'=0}^{n_{\text{core} }} \lambda_{i'} } \times 100\%, \]

(7.255)¶\[\mathcal{R}_{\text{virt} }=\frac{\sum_{a=0}^{n_{\text{virt} }^{\text{excluded} }} \lambda_{a} }{\sum_{a=0}^{n_{\text{virt} }} \lambda_{a} } \times 100\%, \]

of the sum of the excluded eigenvalues to the sum over all eigenvalues. The orbital selection in the core and virtual subspaces is then based upon the values of these ratios, as will be discussed below.

The orbital selection procedure is activated by adding the keyword OrbitalSelection to the first line of input, e.g.

! MR-EOM def2-TZVPP VeryTightSCF OrbitalSelection

There are two threshold parameters CoreThresh and VirtualThresh, which are used to determine which inactive core and virtual orbitals are to be discarded in the orbital selection procedure, respectively. Namely, all inactive core orbitals for which \(\mathcal{R_{\text{core} }}<\) CoreThresh (i.e. \(\mathcal{R_{\text{core} }}\) as defined in Eq. (7.254)) are discarded and all virtual orbitals satisfying the condition \(\mathcal{R_{\text{virt} }}<\) VirtualThresh (i.e. \(\mathcal{R_{\text{virt} }}\) as defined in Eq. (7.255)) are discarded. The default values of these thresholds are CoreThresh = 0.0 (no core orbital selection) and VirtualThresh = 1.0. However, the values of these parameters can easily be changed by redefining them in the %mdci block:

%mdci
  CoreThresh 1.0
  VirtualThresh 1.0
end

Let us consider the calculation of the previous section (A Projection/Singular PT Scheme to Overcome Convergence Issues in the T Amplitude Iterations) on ferrocene, with the orbital selection procedure activated (using the default thresholds):

!MR-EOM def2-SVP OrbitalSelection

%casscf
  nel 6
  norb 5
  mult   1,3
  nroots 5,6
end

%mdci
  DoSingularPT true
  SingularPTThresh 0.05
  MaxIter 50
end

%mrci
  newblock 1 *
    nroots 18
    refs cas(6,5) end
  end
  newblock 3 *
    nroots 10
    refs cas(6,5) end
  end
end

* xyz 0 1
  Fe 0.000000  0.000000  0.000000
  C  0.000000  1.220080  1.650626
  C -1.160365  0.377025  1.650626
  C -0.717145 -0.987065  1.650626
  C  0.717145 -0.987065  1.650626
  C  1.160365  0.377025  1.650626
  C  0.000000  1.220080 -1.650626
  C  1.160365  0.377025 -1.650626
  C  0.717145 -0.987065 -1.650626
  C -0.717145 -0.987065 -1.650626
  C -1.160365  0.377025 -1.650626
  H  0.000000  2.306051  1.635648
  H -2.193184  0.712609  1.635648
  H -1.355463 -1.865634  1.635648
  H  1.355463 -1.865634  1.635648
  H  2.193184  0.712609  1.635648
  H  0.000000  2.306051 -1.635648
  H  2.193184  0.712609 -1.635648
  H  1.355463 -1.865634 -1.635648
  H -1.355463 -1.865634 -1.635648
  H -2.193184  0.712609 -1.635648
end

The details of the orbital selection procedure are printed under the heading ORBITAL SELECTION:

------------------------------------------------
ORBITAL SELECTION
------------------------------------------------

T1 is NOT used in the construction of the orbital selection densities
Factor (in percent) for inactive (core) orbital selection                                    ...      0.000000000
Factor (in percent) for virtual orbital selection                                            ...      1.000000000
Inactive orbitals before selection:    15 ...   44 ( 30 MO's/ 60 electrons)
Virtual orbitals before selection:     50 ...  220 (171 MO's              )
Inactive orbitals after selection:     15 ...   44 ( 30 MO's/ 60 electrons)
Virtual orbitals after selection:      50 ...  126 ( 77 MO's              )

-------------------------------------------
TIMINGS FOR THE ORBITAL SELECTION PROCEDURE
-------------------------------------------

Total Time for Orbital Selection            ...       61.470 sec

First Half Transformation                   ...       58.041 sec ( 94.4%)
Second Half Transformation                  ...        1.789 sec (  2.9%)
Formation of Orbital Selection Densities    ...        1.629 sec (  2.7%)
Core Orbital Selection                      ...        0.000 sec (  0.0%)
Virtual Orbital Selection                   ...        0.003 sec (  0.0%)
Finalization of Orbitals                    ...        0.006 sec (  0.0%)

Comparing the number of virtual orbitals before the orbital selection procedure (171) with the number that are left after orbital selection (77), we see that more than half have been discarded (94). The canonical calculation (without orbital selection) takes 149373 seconds to run and yields the following excitation energies:

-------------------
TRANSITION ENERGIES
-------------------

The lowest energy is  -1648.190045042 Eh

State Mult Irrep Root Block    mEh          eV      1/cm
  1    -1    0    0      0.000      0.000       0.0
  3    -1    0    1     65.110      1.772   14289.9
  3    -1    1    1     65.110      1.772   14289.9
  3    -1    2    1     70.413      1.916   15454.0
  3    -1    3    1     70.413      1.916   15454.0
  3    -1    4    1     95.979      2.612   21065.0
  3    -1    5    1     95.979      2.612   21065.0
  1    -1    1    0    105.302      2.865   23111.1
  1    -1    2    0    105.302      2.865   23111.1
  1    -1    3    0    107.034      2.913   23491.4
  1    -1    4    0    107.034      2.913   23491.4
  1    -1    5    0    160.595      4.370   35246.6
  1    -1    6    0    160.596      4.370   35246.6
  3    -1    6    1    164.694      4.482   36146.1
  3    -1    7    1    165.379      4.500   36296.6
  3    -1    8    1    165.379      4.500   36296.6
  3    -1    9    1    171.464      4.666   37632.1
  1    -1    7    0    208.587      5.676   45779.6
  1    -1    8    0    208.587      5.676   45779.6
  1    -1    9    0    213.093      5.799   46768.6
  1    -1   10    0    213.093      5.799   46768.6
  1    -1   11    0    216.225      5.884   47456.0
  1    -1   12    0    220.230      5.993   48334.9
  1    -1   13    0    220.230      5.993   48334.9
  1    -1   14    0    224.583      6.111   49290.3
  1    -1   15    0    224.583      6.111   49290.3
  1    -1   16    0    237.914      6.474   52216.0
  1    -1   17    0    237.914      6.474   52216.0

In contrast, the calculation with the orbital selection procedure activated runs in 28977 seconds (a factor of 5 speedup) and produces the following excitation energies:

-------------------
TRANSITION ENERGIES
-------------------

The lowest energy is  -1647.788478559 Eh

State Mult Irrep Root Block    mEh          eV      1/cm
  1    -1    0    0      0.000      0.000       0.0
  3    -1    0    1     65.112      1.772   14290.4
  3    -1    1    1     65.134      1.772   14295.3
  3    -1    2    1     70.520      1.919   15477.3
  3    -1    3    1     70.520      1.919   15477.3
  3    -1    4    1     96.105      2.615   21092.7
  3    -1    5    1     96.134      2.616   21099.0
  1    -1    1    0    105.415      2.868   23136.0
  1    -1    2    0    105.450      2.869   23143.5
  1    -1    3    0    107.294      2.920   23548.3
  1    -1    4    0    107.294      2.920   23548.3
  1    -1    5    0    161.082      4.383   35353.4
  1    -1    6    0    161.094      4.384   35356.0
  3    -1    6    1    164.786      4.484   36166.4
  3    -1    7    1    165.465      4.503   36315.4
  3    -1    8    1    165.465      4.503   36315.5
  3    -1    9    1    171.542      4.668   37649.1
  1    -1    7    0    208.853      5.683   45838.0
  1    -1    8    0    208.853      5.683   45838.0
  1    -1    9    0    213.419      5.807   46840.1
  1    -1   10    0    213.419      5.807   46840.1
  1    -1   11    0    216.526      5.892   47521.9
  1    -1   12    0    220.611      6.003   48418.4
  1    -1   13    0    220.611      6.003   48418.5
  1    -1   14    0    225.135      6.126   49411.5
  1    -1   15    0    225.136      6.126   49411.5
  1    -1   16    0    238.388      6.487   52320.1
  1    -1   17    0    238.388      6.487   52320.1

We note that the excitation energies in the orbital selection procedure agree very nicely with those of the canonical calculation. However, the total energies are significantly different, as we currently have not implemented a procedure to correct them. Hence, the following warning is particularly important.

WARNING

The orbital selection procedure should only be used for the calculation of excitation energies. Total energies computed with the orbital selection procedure have not been corrected and can differ greatly from the canonical results.

Before leaving the discussion of the orbital selection procedure, we note that there is also a keyword PrintOrbSelect, which can be added to the %mdci block to print the eigenvalues of the inactive core orbital selection and virtual orbital selection densities and the corresponding values of \(\mathcal{R}_{\text{core} }\) and \(\mathcal{R}_{\text{virt} }\) defined in Eqs. (7.254) and (7.255), respectively. This is useful if one wants to manually select the orbitals to discard in the orbital selection procedure by adjusting the values of CoreThresh and VirtualThresh. We note that the program terminates after printing. In the case of the calculation on ferrocene, if we modify the %mdci block to read

%mdci
  DoSingularPT true
  SingularPTThresh 0.05
  MaxIter 50
  PrintOrbSelect True
end

we find the following information in the ORBITAL SELECTION section of the output (only the first 50 values for the virtual orbital selection density are shown here):

Eigenvalues and corresponding R_core values for the core orbital selection density

Orbital        Eigenvalue        R_core
         0.00026936        0.419318
         0.00027080        0.840883
         0.00038739        1.443947
         0.00038739        2.047011
         0.00040299        2.674355
         0.00040299        3.301700
         0.00077636        4.510285
         0.00086085        5.850394
         0.00086085        7.190503
         0.00091850        8.620358
         0.00091850       10.050213
         0.00112826       11.806598
         0.00115561       13.605563
         0.00137961       15.753236
         0.00137961       17.900908
         0.00139093       20.066210
         0.00139093       22.231512
         0.00143349       24.463072
         0.00143350       26.694633
         0.00148539       29.006985
         0.00148539       31.319338
         0.00173415       34.018940
         0.00224131       37.508054
         0.00224132       40.997171
         0.00533017       49.294785
         0.00533019       57.592429
         0.00658679       67.846267
         0.00662033       78.152314
         0.00701718       89.076149
         0.00701719      100.000000

 Eigenvalues and corresponding R_virt values for the virtual orbital selection density

Orbital        Eigenvalue        R_virt
         0.00000119        0.000450
         0.00000119        0.000899
         0.00000134        0.001404
         0.00000134        0.001909
         0.00000136        0.002423
         0.00000177        0.003091
         0.00000178        0.003764
         0.00000178        0.004437
         0.00000215        0.005248
         0.00000224        0.006096
         0.00000224        0.006944
         0.00000238        0.007844
         0.00000347        0.009154
         0.00000347        0.010465
         0.00000364        0.011841
         0.00000386        0.013299
         0.00000396        0.014793
         0.00000396        0.016287
         0.00000437        0.017937
         0.00000437        0.019587
         0.00000499        0.021472
         0.00000499        0.023357
         0.00000794        0.026354
         0.00000794        0.029352
         0.00000819        0.032447
         0.00000819        0.035543
         0.00000927        0.039044
         0.00000927        0.042546
         0.00001002        0.046332
         0.00001002        0.050119
         0.00001137        0.054415
         0.00001137        0.058711
         0.00001158        0.063086
         0.00001158        0.067461
         0.00001381        0.072678
         0.00001381        0.077894
         0.00001417        0.083249
         0.00001417        0.088604
         0.00001465        0.094137
         0.00001495        0.099785
         0.00001495        0.105432
         0.00001554        0.111302
         0.00001554        0.117172
         0.00001623        0.123303
         0.00001689        0.129685
         0.00001754        0.136310
         0.00001754        0.142934
         0.00001805        0.149752
         0.00001805        0.156570
         0.00002111        0.164546

.
.
.

In conclusion, the orbital selection scheme provides a more efficient way to calculate accurate excitation spectra within the framework of MR-EOM. It can be used to extend the applicability of this approach to larger systems and we expect it to be much more effective in larger systems where the chromophore is localized to a small part of the molecule. We reiterate that it is currently limited to the calculation of excitation energies and should not be used if one is interested in total energies.

7.40.5. Nearly Size Consistent Results with MR-EOM by Employing an MR-CEPA(0) Shift in the Final Diagonalization Procedure¶

One drawback of the MR-EOM methodology is that it is not size-extensive (or size-consistent). The size-extensivity errors arise due to the final uncontracted MR-CI diagonalization step. Namely, they result from the components of the eigenvectors of the transformed Hamiltonian, which lie outside of the CASSCF reference space (e.g. 1h, 1p, etc. configurations). As more of the excitation classes are included through the successive similarity transformations of the Hamiltonian, the size of the final diagonalization manifold is greatly decreased resulting in much smaller size-extensivity errors upon going from MR-EOM-T|T\(^\dagger\)-h-v to MR-EOM-T|T\(^\dagger\)|SXD|U-h-v. To illustrate this, let us consider the O\(_2\)—O\(_2\) dimer where the O\(_2\) molecules are separated by a large distance. For the O\(_2\) monomer, we employ a minimal active space consisting of 2 electrons distributed amongst the two \(\pi^*\) orbitals and we only consider the ground \(^3\Sigma^{-}_g\) state (no state-averaging). In the MR-EOM calculations, we also calculate the higher lying \(^1\Delta_g\) and \(^1\Sigma_g^{+}\) singlet states. For example, the input file for the MR-EOM-T|T\(^\dagger\)|SXD|U-h-v calculation is given by:

!MR-EOM AUG-CC-PVTZ EXTREMESCF

%casscf
  nel 2
  norb 2
  nroots 1
  mult 3
end

%mdci
  MaxIter 300
  STol 1e-12
end

%mrci
  newblock 1 *
    nroots 3
    refs cas(2,2) end
  end
  newblock 3 *
    nroots 1
    refs cas(2,2) end
  end
end

* xyz 0 3
  O  0.00000000    -0.00000000    -0.60500000
  O -0.00000000     0.00000000     0.60500000
*

In the case of the dimer, we take the reference state as the coupled quintet state which is formed as the product \(^3\Sigma_g^+\otimes\,^3\Sigma_g^+\) of the monomer states. We note that at large separation, in the non-interacting limit, the dimer state energies can be decomposed as the sum of monomer state energies. There are various possibilities, taking into account the degeneracies of the various states:

a singlet, a triplet, and a quintet with energy \(E(^3\Sigma_g^{-}+\,^3\Sigma_g^{-})\),
four triplets with energy \(E(^3\Sigma_g^{-}+\,^1\Delta_g)\),
two triplets with energy \(E(^3\Sigma_g^{-}+\,^1\Sigma_g^{+})\),
four singlets with energy \(E(^1\Delta_g+\,^1\Delta_g)\),
four singlets with energy \(E(^1\Delta_g+\,^1\Sigma_g^+)\),
a singlet with energy \(E(^1\Sigma_g^{+}+\,^1\Sigma_g^{+})\).

Hence, in the final diagonalization step of the MR-EOM calculation, we must ask for 10 singlets, 7 triplets and 1 quintet. The input file for the MR-EOM-T|T\(^\dagger|\)SXD|U-h-v calculation on the dimer is given by:

!MR-EOM AUG-CC-PVTZ EXTREMESCF

%casscf
  nel 4
  norb 4
  nroots 1
  mult 5
  etol 1e-13
  gtol 1e-13
end

%mdci
  MaxIter 300
  STol 1e-12
end

%mrci
  newblock 1 *
    nroots 10
    refs cas(4,4) end
  end
  newblock 3 *
    nroots 7
    refs cas(4,4) end
  end
  newblock 5 *
    nroots 1
    refs cas(4,4) end
  end
end

* xyz 0 5
  O                  0.00000000     0.00000000  -500.60500000
  O                  0.00000000    -0.00000000  -499.39500000
  O                 -0.60500000     0.00000000   500.00000000
  O                  0.60500000    -0.00000000   500.00000000
*

In Table Test for size consistency in MR-EOM: Differences in energy (in mE_\text{h}) between the _2—O_2 dimer energies (at large separation) and the sum of the monomer energies for the ground state and various excited states. The results were obtained in an aug-cc-pVTZ basis using minimal active spaces., we have compiled the results of the size consistency test, taking the difference of the dimer state energies (at large separation) and the sum of the monomer state energies (in \(\text{m}E_\text{h}\)). It is evident that as more excitation classes are included in the similarity transformed Hamiltonian and the size of the final diagonalization manifold is decreased, the size-consistency errors decrease. Of particular note are the results for the MR-EOM-T|T\(^\dagger\)|SXD|U-h-v approach (only includes 1h and 1p configurations in the final diagonalization manifold), for which the largest deviation is \(1.25 \times 10^{-2}~\text{m}E_\text{h}\). The much larger deviations for the MR-EOM-T|T\(^\dagger|\)SXD-h-v approach clearly demonstrate the large effect that the 2h excitations have on the size-consistency errors.

Table 7.24 Test for size consistency in MR-EOM: Differences in energy (in mE\(_\text{h}\)) between the \(_2\)—O\(_2\) dimer energies (at large separation) and the sum of the monomer energies for the ground state and various excited states. The results were obtained in an aug-cc-pVTZ basis using minimal active spaces.¶
	T\|T\(^{\dagger}\)-h-v	T\|T\(^{\dagger}\)\|SXD-h-v (with 1h1p)	T\|T\(^{\dagger}\)\|SXD-h-v	T\|T\(^{\dagger}\)\|SXD\|U-h-v
\(\Delta E(^3\Sigma_g^{-}+\,^3\Sigma_g^{-})\)	12.74	2.77	1.11	1.00 \(\times\) \(10^{-5}\)
\(\Delta E(^3\Sigma_g^{-}+\,^1\Delta_g)\)	14.20	3.84	1.85	1.54 \(\times\) \(10^{-4}\)
\(\Delta E(^3\Sigma_g^{-}+\,^1\Sigma_g^{+})\)	17.21	5.52	2.83	4.13 \(\times\) \(10^{-4}\)
\(\Delta E(^1\Delta_g+\,^1\Delta_g)\)	15.69	3.10	2.31	5.26 \(\times\) \(10^{-3}\)
\(\Delta E(^1\Delta_g+\,^1\Sigma_g^{+})\)	18.83	7.52	4.76	5.89 \(\times\) \(10^{-3}\)
\(\Delta E(^1\Sigma_g^{+}+\,^1\Sigma_g^{+})\)	22.34	10.75	7.31	1.25 \(\times\) \(10^{-2}\)

To reduce the size-consistency errors, one can make use of the MR-CEPA(0) shift in the final diagonalization step. This MR-CEPA(0) shift can easily be activated by adding the line

citype mrcepa_0

to the beginning of the %mrci block. The results of the size-consistency test with the use of the MR-CEPA(0) shift are tabulated in Table Test for size consistency in MR-EOM, using the MR-CEPA(0) shift: Differences in energy (in mE_\text{h}) between the O_2—O_2 dimer energies (at large separation) and the sum of the monomer energies for the ground state and various excited states. The results were obtained in an aug-cc-pVTZ basis using minimal active spaces and the MR-CEPA(0) shift was applied in the final diagonalization in each case.. For each of the methods, we see a marked improvement over the results of Table Test for size consistency in MR-EOM: Differences in energy (in mE_\text{h}) between the _2—O_2 dimer energies (at large separation) and the sum of the monomer energies for the ground state and various excited states. The results were obtained in an aug-cc-pVTZ basis using minimal active spaces., which do not make use of the MR-CEPA(0) shift. The greatest improvement occurs in the MR-EOM-T|T\(^\dagger\)|SXD-h-v and the MR-EOM-T|T\(^\dagger\)|SXD|U-h-v results. Namely, the errors in the former case are on the order of nano Hartrees, while the errors in the MR-EOM-T|T\(^\dagger\)|SXD|U-h-v results are not detectable (sub-nano Hartree), as the energy is only printed with nine decimal places. It is interesting to note that upon adding the 1h1p configurations to the diagonalization manifold in the MR-EOM-T|T\(^\dagger\)|SXD-h-v calculations (i.e. with 1h1p), the size-consistency errors increase greatly. Hence, it appears that the use of the MR-CEPA(0) shift is most effective at reducing the size-consistency errors resulting from the presence of the 1h, 1p and 2h configurations in the final diagonalization manifold. In any case, one can easily take advantage of this approach to obtain nearly size-consistent results with both the MR-EOM-T|T\(^\dagger\)|SXD-h-v and MR-EOM-T|T\(^\dagger\)|SXD|U-h-v methods.

Table 7.25 Test for size consistency in MR-EOM, using the MR-CEPA(0) shift: Differences in energy (in mE\(_\text{h}\)) between the O\(_2\)—O\(_2\) dimer energies (at large separation) and the sum of the monomer energies for the ground state and various excited states. The results were obtained in an aug-cc-pVTZ basis using minimal active spaces and the MR-CEPA(0) shift was applied in the final diagonalization in each case.¶
	T\|T\(^{\dagger}\)-h-v	T\|T\(^{\dagger}\)\|SXD-h-v (with 1h1p)	T\|T\(^{\dagger}\)\|SXD-h-v
\(\Delta E(^3\Sigma_g^{-}+\,^3\Sigma_g^{-})\)	2.75 \(\times\) \(10^{-3}\)	0.01	2.00 \(\times\) \(10^{-6}\)
\(\Delta E(^3\Sigma_g^{-}+\,^1\Delta_g)\)	0.06	0.07	0.00
\(\Delta E(^3\Sigma_g^{-}+\,^1\Sigma_g^{+})\)	0.14	0.15	4.00 \(\times\) \(10^{-6}\)
\(\Delta E(^1\Delta_g+\,^1\Delta_g)\)	0.21	0.22	1.00 \(\times\) \(10^{-6}\)
\(\Delta E(^1\Delta_g+\,^1\Sigma_g^{+})\)	0.42	0.44	5.00 \(\times\) \(10^{-6}\)
\(\Delta E(^1\Sigma_g^{+}+\,^1\Sigma_g^{+})\)	0.82	0.87	9.00 \(\times\) \(10^{-6}\)

7.40.6. Perturbative MR-EOM-PT¶

The MR-EOMPT approach was developed for situations where the full accuracy of the iterative MR-EOMCC method is not required. It performs on par with other multireference perturbation theories such as fic-NEVPT2 and does not have the convergence difficulties with the \(\hat T\) and \(\hat S,\hat X,\hat D\) amplitudes like its iterative parent method as these amplitudes are computed in a non-iterative fashion. The only iterative part of the MR-EOMPT method is the calculation of the \(\hat U\) amplitudes since they are quick to converge anyways [501].

The setup procedure for the MR-EOMPT method is the same as for the MR-EOMCC method, and the foregoing also applies to the perturbative variant. Please note that the orbital selection scheme has not been tested with this variant and should be unnecessary anyways since calculations are much faster than with the iterative MR-EOMCC method.

To invoke the new variant, set up the calculation as you would for an MR-EOMCC calculation and then add the keyword DoMREOM_MRPT True to the %mdci block.

The results are interpreted just like results for the iterative MR-EOMCC method. After transforming the Hamiltonian with the perturbatively estimated amplitudes and the final MRCIS diagonalization step, the final state energies are printed along with their reference weights. For reliable results, we again recommend that the reference weight be >90%.