6.8. MR-EOM-CC: Multireference Equation of Motion Coupled-Cluster¶
The Multireference Equation of Motion Coupled-Cluster (MR-EOM-CC) methodology [193, 194, 203, 406, 407, 639] 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 :
a state-averaged CASSCF calculation (CASSCF module),
the solution of amplitude equations and the calculation of the elements of the similarity transformed Hamiltonians (MDCI module),
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. Multireference Equation of Motion Coupled-Cluster (MR-EOM-CC) Theory. We also note that the theoretical details underlying these methods can be found in Ref. [407]. In Sec. Multireference Equation of Motion Coupled-Cluster (MR-EOM-CC) Theory, 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.
6.8.1. 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 Complete Active Space Self-Consistent Field Method
and The Complete Active Space Self-Consistent Field (CASSCF) Module and that a description of the MRCI
module is given in
sections Multireference Configuration Interaction and Pertubation Theory
and The Multireference Correlation Module. 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):
!MR-EOM def2-TZVP VeryTightSCF
%casscf # reference state
nel 6
norb 4
mult 1,3
nroots 3,3
end
%mdci
STol 1e-7
end
%mrci # final roots
newblock 1 *
nroots 6
refs cas(6,4) end
end
newblock 3 *
nroots 6
refs cas(6,4) end
end
end
* xyz 0 1
H 0.000000 0.934473 -0.588078
H 0.000000 -0.934473 -0.588078
C 0.000000 0.000000 0.000000
O 0.000000 0.000000 1.221104
*
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
!MR-EOM-T|Td def2-TZVP VeryTightSCF
runs the MR-EOM-T|T\(^{\dagger}\)-h-v calculation, while
!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:
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
A Projection/Singular PT Scheme to Overcome Convergence Issues in the T Amplitude Iterations.
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
A Tutorial Type Example of a MR Calculation. 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
A Tutorial Type Example of a MR Calculation 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
[194, 203, 406, 407, 529, 530, 639],
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.
6.8.2. 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 (General Description) 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 Multireference Equation of Motion Coupled-Cluster (MR-EOM-CC) Theory.
6.8.3. Perturbative MR-EOM-PT¶
The MR-EOM family of methods now also features an almost fully perturbative approach called MR-EOMPT [501]. 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.