A simple and powerful method for comparing many-electron wavefunctions constructed at different levels of theory is presented. By using wavefunction overlaps, it is possible to analyze the effects of varying wavefunction models, molecular orbitals, and one-electron basis sets. The computation of wavefunction overlaps eliminates the inherent ambiguity connected to more rudimentary wavefunction analysis protocols, such as visualization of orbitals or comparing selected physical observables. Instead, wavefunction overlaps allow processing the many-electron wavefunctions in their full inherent complexity. The presented method is particularly effective for excited state calculations as it allows for automatic monitoring of changes in the ordering of the excited states. A numerical demonstration based on multireference computations of two test systems, the selenoacrolein molecule and an iridium complex, is presented.

A common task in quantum chemistry is the comparison of electronic structure computations performed at different levels of theory. This is unavoidable when comparing a cheaper method against a high-level benchmark or when testing a newly developed computational protocol. Specifically for excited state calculations, it is important to have a quick and reliable method to monitor whether changes in the wavefunction characters occur between calculations carried out with different methods. A common practice is to compare the energies and selected other physical observables derived from the particular wavefunctions. While this is a straightforward and physically sound approach, it requires a careful choice of the properties to be considered. Alternatively, qualitative insight into the wavefunction character can be obtained through visualization of the canonical orbitals or by application of some specific visualization protocols.1–3 However, the assignment of state characters based on orbitals alone is challenging in many cases, e.g., when the orbital character is not clearly visible, when the orbitals are delocalized, when many configurations are involved, or when partial double excitation character is present. To overcome some of these problems, different phenomenological descriptors have been introduced that monitor diverse specific properties, such as the number of unpaired electrons,4–6 double excitation character,3,7 charge transfer,8–12 excitonic effects,13–15 collectivity,16 entanglement,17,18 and orbital relaxation.19 While these certainly provide important insight, the sheer number of such descriptors shows how challenging it is to grasp many-electron wavefunctions in their full complexity. In this communication we want to promote a more fundamental and unbiased mode of comparing two many-electron wavefunctions: Computing the scalar product between the two wavefunctions, i.e., their overlap. Despite its simplicity, this idea does not seem to have received appreciable attention so far, aside from some formal discussions20,21 and a partially related investigation of transition density matrices between arbitrary state functions.22 

The starting point for the discussion is the electronic Schrödinger equation

H ̂ Ψ = E Ψ ,
(1)

where H ̂ is the electronic Hamilton operator, E is an energy eigenvalue, and Ψ is the corresponding antisymmetric many-electron eigenfunction. Given two sets of individually orthonormal approximate eigenfunctions of H ̂ , denoted { Ψ I } and Ψ J , we propose to compute the overlap between two such functions,

S I J = Ψ I Ψ J ,
(2)

to answer the question of how much they resemble each other.

Computing wavefunction overlaps is a rather involved task in the general case of varying orbitals that requires either an explicit construction of the Slater determinants23–25 or a transformation of the orbital bases.26,27 We have recently introduced an algorithm for the computation of such overlaps, using a very general formalism based on Slater determinants.28 Aside from its original intention of computing overlaps between varying geometries along a molecular dynamics trajectory, this formalism allows to compute the overlap between any pair of wavefunctions given in a Slater determinant expansion. The code is used here to compute overlaps between wavefunctions constructed using varying one-electron basis sets, molecular orbitals (MOs), and wavefunction expansions.

A detailed discussion of the properties of the overlap matrix S is presented elsewhere.28 Here, only one crucial property is highlighted. For this purpose, the wavefunction Ψ I is expanded in the form

Ψ I = J = 1 N Ψ J Ψ J Ψ I + 1 J = 1 N Ψ J Ψ J Ψ I = J = 1 N Ψ J S I J + Ψ I .
(3)

Here, N denotes the dimension of the space spanned by the Ψ J , and Ψ I is the component of Ψ I that belongs to the orthogonal complement of this space. The wavefunction Ψ I is thus decomposed into the individual projection components and the “missing part” Ψ I . As a next step, the above equation is projected onto Ψ I ,

1 = Ψ I Ψ I = J = 1 N S I J 2 + Ψ I 2 .
(4)

This shows that the combined weight of the projection components is normalized to unity and consequently that the sum of the squared overlap values along a column (or a row) of the overlap matrix is less or equal to one. To visualize this decomposition, we plot the SIJ2 values in pie charts with a missing piece corresponding to the Ψ I 2 term (see below).

Two examples29 are provided to illustrate the power of this wavefunction analysis protocol (cf. Fig. 1): selenoacrolein, a molecule with rather simple electronic structure,28,30 and an iridium complex whose theoretical description is highly challenging.19,31 For selenoacrolein, a geometry close to an avoided crossing between the lowest two triplet states is chosen by setting the C = C torsion angle to 55°. It is well known32–34 that the relative energies of and ππ states strongly depend on the treatment of dynamic correlation and therefore we expect the state-ordering to be sensitive to the method. To evaluate this hypothesis, the T1 and T2 states are computed at five different levels of theory using various complete active space self-consistent field (CASSCF) and multireference configuration interaction (MR-CI) expansions.35–38 First, an active space of 6 electrons in 5 orbitals and the ANO-RCC-VDZP basis set (abbreviated D)39 are chosen and the maximum excitation level is varied using the CASSCF(6,5), MR-CIS(6,5), and MR-CISD(6,5) levels of theory. Then the active space is increased to the MR-CISD(8,7) level, and the basis set is enlarged to the triple-ζ level (MR-CISD(6,5)/T). Corresponding energies, dipole moments, and qualitative wavefunction characters are presented in Table I. The energy difference between the T1 and T2 states is below 0.1 eV in all cases, in agreement with the assumption of a nearby crossing. However, discrepancies are found with respect to the characters of the states. At the CASSCF(6,5) level the T1 state is of ππ character while MR-CIS predicts it to be . In contrast, all MR-CISD calculations show a mixture with predominant ππ character. The dipole moments are consistent with respect to the qualitative state characters: they are somewhat above 2 D for the predominant ππ states and well below 1 D for the states.

FIG. 1.

Structural formulas of (a) selenoacrolein and (b) Ir(C3H3N)3 as studied in this work.

FIG. 1.

Structural formulas of (a) selenoacrolein and (b) Ir(C3H3N)3 as studied in this work.

Close modal
TABLE I.

Energies (E, eV) relative to the ground state minimum, dipole moments (μ, D), and wavefunction characters for the lowest two triplet states of selenoacrolein at a torsion angle of 55° computed at different levels of theory.

T1 T2
Methoda E μ Char. E μ Char.
CASSCF(6,5)/D  2.399  2.24  ππ  2.491  0.86   
MR-CIS(6,5)/D  2.670  0.37    2.758  2.42  ππ 
MR-CISD(6,5)/D  2.505  2.17  ππ/  2.561  0.67  /ππ 
MR-CISD(8,7)/D  2.547  2.24  ππ/  2.612  0.59  /ππ 
MR-CISD(6,5)/T  2.484  2.05  ππ/  2.529  0.80  /ππ 
T1 T2
Methoda E μ Char. E μ Char.
CASSCF(6,5)/D  2.399  2.24  ππ  2.491  0.86   
MR-CIS(6,5)/D  2.670  0.37    2.758  2.42  ππ 
MR-CISD(6,5)/D  2.505  2.17  ππ/  2.561  0.67  /ππ 
MR-CISD(8,7)/D  2.547  2.24  ππ/  2.612  0.59  /ππ 
MR-CISD(6,5)/T  2.484  2.05  ππ/  2.529  0.80  /ππ 
a

Basis set: D = ANO-RCC-VDZP, T = ANO-RCC-VTZP.

Despite the fact that the analysis of Table I is probably sufficient here, there are inherent limitations with respect to the strategies used. First, a comparison of orbital transitions requires that the orbitals produced by the different methods resemble each other, and if the wavefunctions are mixed, e.g., having ππ/ character, additional care is necessary. Then, the application of operator expectation values requires a priori knowledge of the relevant properties to analyze. It is, thus, clearly beneficial to have an automated strategy for comparing two wavefunctions that does not require any problem-specific adjustment or a priori knowledge of the results. We argue that wavefunction overlaps provide just that. As an illustration, the overlaps of the T1 and T2 MR-CISD(6,5)/D wavefunctions with respect to the other computational methods are shown in the pie charts in Fig. 2. These overlaps evidence the effects of lowering the excitation level, on the one hand, and of increasing the active space or basis set, on the other hand. The two columns of Fig. 2 correspond to the wavefunctions, Ψ 1 and Ψ 2 , computed with the indicated respective methods. The pie charts represent the square of the overlap values of these states with respect to the Ψ 1 (blue) and Ψ 2 (red) states computed at the MR-CISD(6,5)/D level. Each pie chart thus corresponds to one column of the overlap matrix. The weight of the orthogonal complement, i.e., the Ψ I 2 term in Eq. (4), is shown as the missing piece of the pie. At the CASSCF(6,5)/D level (Fig. 2(a)), one immediately sees that the state ordering is the same as provided by MR-CISD(6,5)/D, in agreement with Table I. There is a squared overlap of 82% between the T1 wavefunctions and 4% mixing with T2. However, the orthogonal complement (14%) plays an important role, highlighting that dynamic correlation alters the wavefunction. The T2 state is basically the mirror image of the T1 state, due to the simple nature of the excited states in this case. The MR-CIS(6,5)/D level of theory (Fig. 2(b)) induces a reversal in state ordering and the squared Ψ 2 Ψ 1 term amounts to 58% while Ψ 1 Ψ 1 2 = 31 % . The weight of the orthogonal complement is 11%, slightly smaller than CASSCF (14%). The T2 state is again the mirror image of T1. In the case of MR-CISD(8,7)/D (Fig. 2(c)) a 99% overlap is observed for both the T1 and T2 states and it is clear that the increase of the active space is irrelevant here. A somewhat more pronounced effect is observed upon increasing the basis set to the triple-ζ level (MR-CISD(6,5)/T, Fig. 2(d)) where only a 96% squared overlap is obtained. While the previous analysis of Table I is consistent with these results, the overlap computation provides the added benefit of a well-defined quantitative analysis.

FIG. 2.

Squared overlap terms Ψ I Ψ J 2 between wavefunctions computed at different levels of theory using the example of the two lowest triplet states of selenoacrolein at a torsion angle of 55°. The pie charts correspond to individual ket states computed by the method indicated above them. The overlap terms of the Ψ 1 and Ψ 2 wavefunctions computed by the MR-CISD(6,5)/D method are shown in blue and red, respectively, while the orthogonal complement corresponds to the sector missing in the chart.

FIG. 2.

Squared overlap terms Ψ I Ψ J 2 between wavefunctions computed at different levels of theory using the example of the two lowest triplet states of selenoacrolein at a torsion angle of 55°. The pie charts correspond to individual ket states computed by the method indicated above them. The overlap terms of the Ψ 1 and Ψ 2 wavefunctions computed by the MR-CISD(6,5)/D method are shown in blue and red, respectively, while the orthogonal complement corresponds to the sector missing in the chart.

Close modal

A second, more demanding example is provided by analyzing the excited state wavefunctions of the fac-tris(3-iminoprop-1-en-1-ido)iridum complex [Ir(C3H3N)3] shown in Fig. 1(b). The wavefunctions of transition metal complexes are well-known to strongly depend on the computational method.40,41 Therefore, two active spaces (using 12 electrons in either 9 or 12 orbitals) and two different excitation levels are considered. Overall, we have performed CASSCF(12,9), CASSCF(12,12), MR-CIS(12,9), and MR-CIS(12,12) calculations, in connection with the LANL2DZ pseudopotential.42,43 In all cases state-averaging over the lowest ten singlet states, corresponding to four states of A and six of E symmetry, is carried out. For simplicity, the analysis is restricted to the four non-degenerate A states. As shown in Table II, already the excitation energies of these states exhibit strong deviations depending on the method used. The 21A energies vary between 4.17 eV for CASSCF(12,9) and 3.56 eV for CASSCF(12,12), and in the case of the 41A state the difference between CASSCF(12,9) and MR-CIS(12,12) exceeds 1 eV. These excited states are combinations of ligand-centered (LC) and metal-to-ligand charge transfer (MLCT) states. However, the precise character of each state is difficult to pinpoint because several excited state configurations are involved in most cases. The relevant orbitals, in turn, possess mixed Ir-d/ligand-π character. A better comparison of the results can be achieved by inspecting the Mulliken charges on the iridium atom qIr. For the ground state, a variation from 1.20 to 1.30 is observed when going from CASSCF(12,9) to MR-CIS(12,12). For the first and second excited states the opposite trend is observed and the CASSCF method gives higher charges. Both trends combined mean that the proportion of MLCT is significantly enhanced by CASSCF (see also Ref. 31). For example, in the case of the 11A → 21A transition the charge transferred from iridium to the ligands amounts to 0.46 e in the case of CASSCF(12,9) while only half of this value is obtained for MR-CIS(12,9). The 41A states generally possess reduced MLCT character and again some discrepancies between the methods are observed.

TABLE II.

Excitation energies (ΔE, eV) and Mulliken charges on the iridium atom (qIr, e) for the first four totally symmetric singlet states of Ir(C3H3N)3 computed at the CASSCF and MR-CIS levels of theory using different active spaces.

CAS(12,9) CAS(12,12) MR-CI(12,9) MR-CI(12,12)
ΔE qIr ΔE qIr ΔE qIr ΔE qIr
11A  …  1.20  …  1.23  …  1.27  …  1.30 
21A  4.17  1.66  3.56  1.63  3.97  1.50  3.76  1.54 
31A  4.49  1.60  3.93  1.58  4.66  1.47  4.46  1.50 
41A  6.19  1.35  5.36  1.37  5.01  1.42  4.96  1.44 
CAS(12,9) CAS(12,12) MR-CI(12,9) MR-CI(12,12)
ΔE qIr ΔE qIr ΔE qIr ΔE qIr
11A  …  1.20  …  1.23  …  1.27  …  1.30 
21A  4.17  1.66  3.56  1.63  3.97  1.50  3.76  1.54 
31A  4.49  1.60  3.93  1.58  4.66  1.47  4.46  1.50 
41A  6.19  1.35  5.36  1.37  5.01  1.42  4.96  1.44 

The analysis of the excitation energies and the Mulliken charges certainly provides some insight into the states of Ir(C3H3N)3; however, some important questions remain unanswered. Of particular interest is the problem of whether the overall state-ordering is preserved in all calculations despite the large shifts in energies. Equally important is to know whether there is a convergence of the wavefunctions when improving the level of theory. These two questions can be easily addressed by computing wavefunction overlaps. In this particular case, the MR-CIS(12,12) level of theory is used as a reference. The results are presented in Figure 3. As opposed to Figure 2, it is observed that the overall state ordering agrees in all cases. However, large pieces of the pies are missing illustrating that no convergence is reached. The lowest level used here, CASSCF(12,9), provides a reasonable description of the ground state with a projection weight of 84% while the agreement for the excited states is significantly lower with a minimal value of 57% in the case of the 41A state. Some mixing among the excited states further suggests that the CASSCF(12,9) description is insufficient. The agreement is somewhat better with CASSCF(12,12). Here, the projection weights to MR-CIS(12,12) vary between 87% and 70% and no significant mixing between the states is observed. Better agreement is observed between MR-CIS(12,9) and MR-CIS(12,12), but even in this case the overlaps are not higher than about 90%. This analysis reveals that discrepancies are present for the different wavefunction models used here to describe Ir(C3H3N)3. They are related to variations in the amount of MLCT character, and a more detailed analysis shows that the MLCT character is in turn affected by orbital relaxation effects.19,31

FIG. 3.

Squared overlap terms Ψ I Ψ J 2 between wavefunctions computed at different levels of theory using the example of the four lowest singlet states of the Ir(C3H3N)3 complex. The pie charts correspond to individual ket states computed by the method indicated above them. The overlap terms of the Ψ 0 , Ψ 1 , Ψ 2 , and Ψ 3 wavefunction computed by the MR-CIS(12,12) method are shown in blue, red, green, and yellow while the orthogonal complement corresponds to the sector missing in the chart.

FIG. 3.

Squared overlap terms Ψ I Ψ J 2 between wavefunctions computed at different levels of theory using the example of the four lowest singlet states of the Ir(C3H3N)3 complex. The pie charts correspond to individual ket states computed by the method indicated above them. The overlap terms of the Ψ 0 , Ψ 1 , Ψ 2 , and Ψ 3 wavefunction computed by the MR-CIS(12,12) method are shown in blue, red, green, and yellow while the orthogonal complement corresponds to the sector missing in the chart.

Close modal

In conclusion, it has been shown that the computation of wavefunction overlaps provides a simple and powerful method for comparing wavefunctions constructed at different levels of theory, that is especially useful in the case of excited state computations. Wavefunction overlaps are suitable to examine changes in wavefunction character not only at a qualitative, but also at a quantitative level. This method is directly applicable to models producing explicit wavefunctions, i.e., the CI and multiconfigurational SCF methods, as used above. Extensions to time-dependent density functional theory, coupled cluster, and other response theory methods are straightforward using approximations that have been described previously23,24,44 allowing for at least semi-quantitative results in these cases. There are no formal limitations on the wavefunctions compared, as long as they possess the same number of α and β electrons, and can be expanded in Slater determinants. Comparisons between different electronic structure packages are accessible if proper care is taken of the ordering and normalization conventions of the basis functions. The tools described here will be released as part of the Sharc molecular dynamics package45,46 and as a standalone module. Interfaces to the Columbus,47Molcas,48ADF,49 and Turbomole50 electronic structure packages are available.

This paper is based on work supported by the VSC Research Center funded by the Austrian Federal Ministry of Science, Research, and Economy (bmwfw). The computational results presented have been achieved in part using the Vienna Scientific Cluster (VSC), Project Nos. 70719 and 70726. We thank COST Action No. CM1405 (MOLIM) for providing a basis for early discussions on the topic of this manuscript.

1.
M.
Head-Gordon
,
A. M.
Grana
,
D.
Maurice
, and
C. A.
White
,
J. Chem. Phys.
99
,
14261
(
1995
).
2.
R. L.
Martin
,
J. Chem. Phys.
118
,
4775
(
2003
).
3.
F.
Plasser
,
M.
Wormit
, and
A.
Dreuw
,
J. Chem. Phys.
141
,
024106
(
2014
).
4.
K.
Takatsuka
,
T.
Fueno
, and
K.
Yamaguchi
,
Theor. Chim. Acta
183
,
175
(
1978
).
5.
M.
Head-Gordon
,
Chem. Phys. Lett.
372
,
508
(
2003
).
6.
S.
Grimme
and
A.
Hansen
,
Angew. Chem., Int. Ed.
54
,
12308
(
2015
).
7.
S.
Matsika
,
X.
Feng
,
A. V.
Luzanov
, and
A. I.
Krylov
,
J. Phys. Chem. A
118
,
11943
(
2014
).
8.
A. V.
Luzanov
and
O. A.
Zhikol
,
Int. J. Quantum Chem.
110
,
902
(
2010
).
9.
F.
Plasser
and
H.
Lischka
,
J. Chem. Theory Comput.
8
,
2777
(
2012
).
10.
T.
Etienne
,
X.
Assfeld
, and
A.
Monari
,
J. Chem. Theory Comput.
10
,
3906
(
2014
).
11.
F.
Plasser
,
B.
Thomitzni
,
S. A.
Bäppler
,
J.
Wenzel
,
D. R.
Rehn
,
M.
Wormit
, and
A.
Dreuw
,
J. Comput. Chem.
36
,
1609
(
2015
).
12.
C.
Adamo
,
T.
Le Bahers
,
M.
Savarese
,
L.
Wilbraham
,
G.
Garcia
,
R.
Fukuda
,
M.
Ehara
,
N.
Rega
, and
I.
Ciofini
,
Coord. Chem. Rev.
304-305
,
166
(
2015
).
13.
S.
Tretiak
and
S.
Mukamel
,
Chem. Rev.
102
,
3171
(
2002
).
14.
S. A.
Bäppler
,
F.
Plasser
,
M.
Wormit
, and
A.
Dreuw
,
Phys. Rev. A
90
,
052521
(
2014
).
15.
D.
Casanova
and
A. I.
Krylov
,
J. Chem. Phys.
144
,
014102
(
2016
).
16.
A. V.
Luzanov
,
A. A.
Sukhorukov
, and
V. E.
Umanskii
,
Theor. Exp. Chem.
10
,
354
(
1976
).
17.
A. V.
Luzanov
and
O. V.
Prezhdo
,
Mol. Phys.
105
,
2879
(
2007
).
18.
F.
Plasser
,
J. Chem. Phys.
144
,
194107
(
2016
).
19.
F.
Plasser
,
S. A.
Bäppler
,
M.
Wormit
, and
A.
Dreuw
,
J. Chem. Phys.
141
,
024107
(
2014
).
20.
W.
Kohn
,
Rev. Mod. Phys.
71
,
1253
(
1999
).
21.
G. K.-L.
Chan
,
WIREs: Comput. Mol. Sci.
2
,
907
(
2012
).
22.
P.-Å.
Malmqvist
and
V.
Veryazov
,
Mol. Phys.
110
,
2455
(
2012
).
23.
E.
Tapavicza
,
I.
Tavernelli
, and
U.
Rothlisberger
,
Phys. Rev. Lett.
98
,
23001
(
2007
).
24.
R.
Mitrić
,
U.
Werner
, and
V.
Bonačić-Koutecký
,
J. Chem. Phys.
129
,
164118
(
2008
).
25.
J.
Pittner
,
H.
Lischka
, and
M.
Barbatti
,
Chem. Phys.
356
,
147
(
2009
).
26.
P.-Å.
Malmqvist
and
B. O.
Roos
,
Chem. Phys. Lett.
155
,
189
(
1989
).
27.
A.
Mitrushchenkov
and
H.-J.
Werner
,
Mol. Phys.
105
,
1239
(
2007
).
28.
F.
Plasser
,
M.
Ruckenbauer
,
S.
Mai
,
M.
Oppel
,
P.
Marquetand
, and
L.
González
,
J. Chem. Theory Comput.
12
,
1207
(
2016
).
29.
See supplementary material at http://dx.doi.org/10.1063/1.4958462 for computational details and for Cartesian geometries of the systems studied.
30.
S.
Mai
,
T.
Müller
,
F.
Plasser
,
P.
Marquetand
,
H.
Lischka
, and
L.
González
,
J. Chem. Phys.
141
,
074105
(
2014
).
31.
F.
Plasser
and
A.
Dreuw
,
J. Phys. Chem. A
119
,
1023
(
2015
).
32.
S.
Perun
,
A. L.
Sobolewski
, and
W.
Domcke
,
J. Am. Chem. Soc.
127
,
6257
(
2005
).
33.
P. G.
Szalay
,
T.
Watson
,
A.
Perera
,
V. F.
Lotrich
, and
R. J.
Bartlett
,
J. Phys. Chem. A
116
,
6702
(
2012
).
34.
J.-M.
Mewes
,
V.
Jovanović
,
C. M.
Marian
, and
A.
Dreuw
,
Phys. Chem. Chem. Phys.
16
,
12393
(
2014
).
35.
H.
Dachsel
,
H.
Lischka
,
R.
Shepard
,
J.
Nieplocha
, and
R. J.
Harrison
,
J. Comput. Chem.
18
,
430
(
1997
).
36.
H.
Lischka
,
R.
Shepard
,
R. M.
Pitzer
,
I.
Shavitt
,
M.
Dallos
,
T.
Müller
,
P. G.
Szalay
,
M.
Seth
,
G. S.
Kedziora
,
S.
Yabushita
, and
Z. Y.
Zhang
,
Phys. Chem. Chem. Phys.
3
,
664
(
2001
).
37.
T.
Müller
,
J. Phys. Chem. A
113
,
12729
(
2009
).
38.
H.
Lischka
,
T.
Müller
,
P. G.
Szalay
,
I.
Shavitt
,
R. M.
Pitzer
, and
R.
Shepard
,
WIREs: Comput. Mol. Sci.
1
,
191
(
2011
).
39.
B. O.
Roos
,
R.
Lindh
,
P.-Å.
Malmqvist
,
V.
Veryazov
, and
P.-O.
Widmark
,
J. Phys. Chem. A
108
,
2851
(
2004
).
40.
S.
Vancoillie
,
H.
Zhao
,
V. T.
Tran
,
M. F. A.
Hendrickx
, and
K.
Pierloot
,
J. Chem. Theory Comput.
7
,
3961
(
2011
).
41.
D.
Escudero
and
L.
González
,
J. Chem. Theory Comput.
8
,
203
(
2012
).
42.
P. J.
Hay
and
W. R.
Wadt
,
J. Chem. Phys.
82
,
299
(
1985
).
43.
P.
Hariharan
and
J.
Pople
,
Theor. Chim. Acta
28
,
213
(
1973
).
44.
F.
Plasser
,
R.
Crespo-Otero
,
M.
Pederzoli
,
J.
Pittner
,
H.
Lischka
, and
M.
Barbatti
,
J. Chem. Theory Comput.
10
,
1395
(
2014
).
45.
S.
Mai
,
M.
Richter
,
M.
Ruckenbauer
,
M.
Oppel
,
P.
Marquetand
, and
L.
González
, SHARC: Surface hopping including arbitrary couplings—Program package for non-adiabatic dynamics, 2014, http://sharc-md.org.
46.
S.
Mai
,
P.
Marquetand
, and
L.
González
,
Int. J. Quant. Chem.
115
,
1215
(
2015
).
47.
H.
Lischka
,
R.
Shepard
,
I.
Shavitt
,
R. M.
Pitzer
,
M.
Dallos
,
T.
Muller
,
P. G.
Szalay
,
F. B.
Brown
,
R.
Ahlrichs
,
H. J.
Boehm
,
A.
Chang
,
D. C.
Comeau
,
R.
Gdanitz
,
H.
Dachsel
,
C.
Ehrhardt
,
M.
Ernzerhof
,
P.
Hoechtl
,
S.
Irle
,
G.
Kedziora
,
T.
Kovar
,
V.
Parasuk
,
M. J. M.
Pepper
,
P.
Scharf
,
H.
Schiffer
,
M.
Schindler
,
M.
Schueler
,
M.
Seth
,
E. A.
Stahlberg
,
J.-G.
Zhao
,
S.
Yabushita
,
Z.
Zhang
,
M.
Barbatti
,
S.
Matsika
,
M.
Schuurmann
,
D. R.
Yarkony
,
S. R.
Brozell
,
E. V.
Beck
,
J.-P.
Blaudeau
,
M.
Ruckenbauer
,
B.
Sellner
,
F.
Plasser
, and
J. J.
Szymczak
, Columbus: Anab initioelectronic structure program, release 7.0, 2015, www.univie.ac.at/columbus.
48.
F.
Aquilante
,
J.
Autschbach
,
R. K.
Carlson
,
L. F.
Chibotaru
,
M. G.
Delcey
,
L.
De Vico
,
I.
Fdez Galvan
,
N.
Ferre
,
L. M.
Frutos
,
L.
Gagliardi
,
M.
Garavelli
,
A.
Giussani
,
C. E.
Hoyer
,
G.
Li Manni
,
H.
Lischka
,
D.
Ma
,
P.-Å.
Malmqvist
,
T.
Müller
,
A.
Nenov
,
M.
Olivucci
,
T. B.
Pedersen
,
D.
Peng
,
F.
Plasser
,
B.
Pritchard
,
M.
Reiher
,
I.
Rivalta
,
I.
Schapiro
,
J.
Segarra-Marti
,
M.
Stenrup
,
D. G.
Truhlar
,
L.
Ungur
,
A.
Valentini
,
S.
Vancoillie
,
V.
Veryazov
,
V. P.
Vysotskiy
,
O.
Weingart
,
F.
Zapata
, and
R.
Lindh
,
J. Comput. Chem.
37
,
506
(
2016
).
49.
G.
Te Velde
,
F. M.
Bickelhaupt
,
E. J.
Baerends
,
C.
Fonseca Guerra
,
S. J. A.
van Gisbergen
,
J. G.
Snijders
, and
T.
Ziegler
,
J. Comput. Chem.
22
,
931
(
2001
).
50.
TURBOMOLE V7.0.1, a development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH (1989-2007), TURBOMOLE GmbH, 2015, http://www.turbomole.com.

Supplementary Material