Quantum Dynamics of Photoactive Transition Metal

Complexes. A Case Study of Model Reduction

arXiv:2305.07682v1 [physics.chem-ph] 12 May 2023

Olga S. Bokarevaa,b, Oliver Kühnc,∗

a Leibniz Institute for Catalysis (LIKAT) Albert-Einstein-Str. 29A 18059 Rostock Germany
b Institute of Physics University of Kassel Heinrich-Plett-Str. 40 34132 Kassel Germany
c Institute of Physics University of Rostock Albert-Einstein-Str. 23-24 D-18059 Rostock

Germany

Abstract
Transition metal complexes for photochemical applications often feature a high
density of electron-vibrational states characterized by nonadiabatic and spin-
orbit couplings. Overall, the dynamics after photoexcitation is shaped by rapid
transitions between states of different character and multiplicity. Even though
transient absorption experiments enable characterization in terms of kinetic
rates, the complexity of the systems usually prevents a more detailed analysis.
Quantum dynamics simulations using quantum chemically determined model
Hamiltonians may provide such details. In particular, one is tempted to pursue
a model reduction, such as to identify couplings or vibrational modes most rele-
vant for the dynamics. Here, we address how such an endeavor is challenged by
the particular nature of transition metal complexes. For that purpose, we per-
formed quantum dynamics simulations for a recently studied iron(II) homoleptic
complex.
Keywords: quantum dynamics, transition metals, linear vibronic coupling,
correlation effects

1. Introduction

Transition metal (TM) complexes are important for a broad range of catalyti-
cal, medical, biological, and material science applications such as photocatalysis,
molecular memory, cell imaging, etc., see e.g., Refs. [1, 2, 3, 4]. The typical size
of many natural and artificial photocatalysts and their high density of electronic
and vibrational states, together with spin-orbit, relativistic, and environmental
effects pose a challenge to quantum chemical and quantum dynamical simula-
tions [5, 6, 7, 8, 9]. Therefore, most of the reported simulations of TM excited
state dynamics employ the trajectory surface hopping (TSH) method, see. e.g.,
Refs. [10, 11, 12, 13, 14, 15]. This quasi-classical approach enjoys popularity

∗ Corresponding author
Email address: oliver.kuehn@uni-rostock.de (Oliver Kühn )

Preprint submitted to Elsevier May 16, 2023

due to its clear concept and the possibility of studying complex systems with
numerous electronic and vibrational degrees of freedom (DOFs). Of course, one
might ask whether all states and modes are equally important for the dynamics.
The answer to this question often will provide a better mechanistic understand-
ing in addition to a means for model reduction towards a setup more feasible
for quantum dynamics. Methods for identifying essential nuclear motions are
well established in the field of molecular dynamics [16, 17, 18, 19, 20]. They all
rely solely on the distribution of nuclear geometries, including minimum energy
seams of crossing and minimum energy conical intersections [21]. However, as
pointed out by Mai and González [22], displacements alone are not sufficient to
judge on the importance of a particular mode, as the impact of the amplitude
of a given nuclear motion on the wave packet evolution may be very complex.
Instead, information about electronic energy gaps and state couplings has to
be taken into account. A workflow for model reduction starting from a TSH
simulation has been proposed by Gómez et al. [12], which additionally included
feedback on the reduced dimensional quantum simulation, such as to improve
TSH parameters. Trajectory surface hopping was also used by other authors to
identify relevant modes for quantum dynamics simulations of TM compounds
[13].
As excited dynamics takes place on the potential energy surfaces (PESs) pre-
computed with the chosen quantum-chemical approach or on-the-fly, the out-
comes of such computations may strongly vary with the applied computational
scheme. An exemplary homoleptic iron(II) complex FePS with two push-pull
ligands [23] has been recently studied by us using the combination of TSH
excited states dynamics with optimally-tuned long-range separated functionals
[24]. Although the concept of optimal tuning was introduced to reliably describe
charge-transfer properties such as metal-to-ligand charge transfer (MLCT) ex-
citations, the appearance of multiple low-lying metal-centered states makes the
choice of an appropriate range-separation scheme more complicated. Being for-
mally local, those MC states may suffer from uncompensated self-interaction
error and, thus, require special attention beyond the classical optimal tuning.
Following this, additional criteria may be desired to adjust the computational
scheme. Assisted by comparison with the experimental absorption spectra and
CASPT2 calculations, two suitable sets of range-separation parameters were
chosen for excited state surface hopping dynamics simulations. The suggested
dynamics mechanisms were notably different for both sets of parameters. Using
the set of optimal parameters from the classical ∆SCF method (set A) resulted
in the population of 3MLCT states that remained stable within the 2 ps sim-
ulation time, which is at variance with transient absorption experiments [23].
Applying multi-reference CASPT2-derived parameter set, TSH calculations led
to relaxation processes into hot singlets states. In the course of the deactiva-
tion, the system stayed in regions of high density of states, reachable for further
transfer into quintet states in accordance with experimental findings [23].
In what follows, we will use the CASPT2-based optimized functional to il-
lustrate some aspects of quantum dynamics simulations of TM compounds alike
FePS. It should be stressed that the purpose is not to obtain converged results

2
for the total system, although, in principle, this should be feasible. Instead,
we will highlight some issue which one may encounter when trying to design a
reduced model for a system with a large density of states and a large number
of couplings of comparable magnitude. In what follows, we first outline the lin-
ear vibronic coupling (LVC) model used for the quantum dynamics simulations.
Next, we address the Multi-layer Multiconfiguration Time-dependent Hartree
(ML-MCTDH) method. In the applications section, we first characterize the
LVC model for the FePS system. Possible model reduction is investigated us-
ing the Time-dependent Hartree (TDH) limit of MCTDH. For the thus identified
reduced model, ML-MCTDH simulations are performed to address correlation
effects. Finally, conclusions are given.

2. Theory

2.1. Linear Vibronic Coupling

In contrast to trajectory-based methods, quantum dynamics simulations re-
quire a priori knowledge of the PES, including possible electronic state couplings.
For TM complexes, this inevitably means that one has to resort to model PES,
the most prominent one being the vibronic coupling model [25, 26, 27]. It rests
on a harmonic oscillator description of nuclear motion in terms of normal modes,
combined with a low-order expansion for the displacement of the PESs and the
state couplings. Specifically, the linear LVC Hamiltonian reads, denoting the
set of f normal mode coordinates as Q,
X
H= [δmn Hm (Q) + Vmn (Q)]|mihn| (1)
mn

with the Hamiltonian for electronic state m

1X X
Hm (Q) = Em + h̄ωξ (−∂ξ2 + Q2ξ ) + κm,ξ Qξ (2)
2
ξ ξ

and the interstate couplings

X
(0)
Vmn (Q) = Vmn + λmn,ξ Qξ . (3)
ξ

Here, κm,ξ and λmn,ξ are the intra- and interstate coupling constants, respec-
(0) (SOC)
tively, related to mode Qξ . Further, Vmn = Vmn denotes the spin-orbital
coupling (SOC), assuming that its coordinate dependence is negligible. The
Hamiltonian, Eqs. (1 - 3), can be supplemented by
X
Hf (t) = −E(t) dmn |mihn| (4)
mn

describing the interaction with an external time-dependent electric field E(t) in

dipole approximation (dmn is the vector of the dipole matrix elements).

3
 Note that the LVC Hamiltonian is not only used for quantum dynamics
simulations (see, e.g. Refs. [28, 29]), if applicable, it also provides an ideal
framework for TSH. Fortunately, this often holds true for TM complexes where
on-the-fly generation of PES is computationally very expensive [14, 24]. This
model is only valid if no large amplitude motions in the excited state dynamics
are present, for diagnostics of accuracy, see the exemplary study in Ref. [30].
In cases where bond breaking is simulated, the LVC model for MCTDH may
be complemented by the corresponding dissociative mode, see e.g. studies of
CO-release from TM-complexes [31, 32].

2.2. ML-MCTDH Method

The method of choice for high-dimensional quantum dynamics simulations
of LVC Hamiltonians is the MCTDH approach [33, 34, 35, 36, 37], which is
briefly sketched in what follows. Given the expansion of the state vector in the
diabatic electronic basis
X
|Ψ(Q; t)i = ψm (Q; t) |mi . (5)
m

The state-specific vibrational wave packets are represented as follows (skipping

the state index m for simplicity)

X
n1 X
nf
ψ(Q; t) = ... A1j1 ...jf (t) · φ1;1 1;f
j1 (Q1 , t) . . . φjf (Qf , t) . . (6)
j1 =1 jf =1

Here, φ1,κ
λ (Qκ , t) are the so-called single-particle functions (SPFs) which form
a time-dependent basis for expanding the wave packet in terms of different
Hartree products of SPFs. Eq. (6) provides a compact representation such that
the number of SPFs per DOF, nκ , usually will be well below that required for a
static basis. This ansatz can, in principle, be converged towards the numerically
exact limit. On the other hand side, there is the TDH limit, where only a single
SPF per DOF is taken into account. In practice, the SPFs are expanded into a
(κ)
time-independent (primitive) basis set, χj (Qκ ), according to

X
Nκ
(κ)
φ1;κ
λ (Qκ , t) = A2;κ
λ;j (t)χj (Qκ ) . (7)
j=1

MCTDH can be seen as a two-layer approach. This view becomes clear upon
interpretation of A2;κ
λ;j (t) as an additional set of expansion coefficients to repre-
sent the time-dependent SPFs of the upper (first) layer in a time-independent
primitive basis in the lower (second) layer. In A2;κ λ;j (t) the superscripts 2 and
κ correspond to the second layer and the κth mode, respectively, whereas the
subscript λ refers to the indices of the respective SPF of the first layer. The ex-
pansion coefficients in Eqs. (6) and (7) can be obtained using the time-dependent
Dirac-Frenkel variational principle [34].

4
 ML-MCTDH builds on the idea of combining strongly correlated DOFs into
p logical coordinates (combined modes or MCTDH particles) yielding the set
{q11 , . . . , qp1 } with qκ1 = {Q1κ , . . . , Qdκ }, where the superscript refers to the first
layer, and dκ is the dimension of the κth particle. The wave packet within the
combined modes’ picture for the uppermost (first) layer can be written as

X
n1 X
np
Ψ(q11 , . . . , qp1 , t) = ... A1j1 ...jp (t) φ1;1 1 1;p 1
j1 (q1 , t) . . . φjp (qp , t) . (8)
j1 =1 jp =1

In ML-MCTDH, the multi-dimensional SPFs for the particles, φ1;κ 1

λ (qκ , t), are
considered as multi-dimensional wave packets of the second layer and expanded
according to

φ1;κ 1
λ (qκ , t) = φ2;κ 2;κ 2;κ
λ (Q1 , . . . , Qdκ , t)
X
nκ,1
X
nκ,dκ
= ... A2;κ 2;κ;1
λ;j1 ...jd (t) φj1 (Q2;κ 2;κ;dκ
1 , t) . . . φjdκ (Q2;κ
dκ , t) .(9)
κ
j1 jdκ

In the third layer, these wave packets are represented in the primitive basis, i.e.

X
Nα
(α)
X
κ−1
φ2;κ;σ
λ (Q2;κ
σ , t) = A3;κ;σ 2;κ
λ;j (t) χj (Qσ ), α=σ+ di . (10)
j=1 i

While this example included three layers, the general idea can easily be applied
to a case with more layers containing particles with decreasing dimensionality
when stepping towards the bottom layer. Within this scheme, electronic states
are incorporated by introducing an electronic coordinate, Qel , in the uppermost
layer. By construction, ML-MCTDH allows for a numerical solution of the time-
dependent Schrödinger equation with controlled convergence. Finally we note
that the effect of an external laser field, Eq. (4), can be incorporated into the
MCTDH scheme a straightforward manner [38]. For the simulations present
below the Heidelberg program package was used [39].

3. Application

3.1. LVC Hamiltonian for FePS

In what follows, we will apply MCTDH to study the quantum dynamics of
an exemplary iron(II) homoleptic complex [Fe(cpmp)2 ]2+ , where cpmp = 6,2”-
carboxypyridyl-2,2’-methylamine-pyridyl-pyridine, shown as an inset in Figure 3
and abbreviated as FePS [23]. The choice is motivated by our previous inves-
tigation of the influence of the parametrization of the range-separated hybrid
functional on the photoinduced dynamics as described by TSH. However, the
present focus will not be on a quantitative comparison with the TSH dynamics
reported in Ref. [24], but on coping with the high dimensionality of the LVC
Hamiltonian for FePS in general. As an example, we use the parameter set of

5
the functional derived based on a comparison with CASPT2 results, which gave
a better agreement between calculated and experimental dynamics. In num-
bers, taking into account the lowest 9 singlet and 14 triplet states (energies up
to about 2.8 eV), there are in total 51 electronic states. Further, the molecule
has 212 vibrational DOFs.
Let us start with an analysis of the LVC Hamiltonian reported in Ref. [24].
Assuming an energy threshold of 0.0003 eV, the LVC Hamiltonian contains
4722 κm,ξ terms, 20119 λmn,ξ terms, and 2464 SOC terms. Figure 1 shows the
distribution of LVC parameters, which indicates that the majority of couplings
is rather small. In Figure 2a, we have plotted the distribution of zeroth-order
states, defined as spin-free states without nonadiabatic state coupling (λmn,ξ =
0). Apparently, the number of states increases rather rapidly with energy.

103
a 104 b c
103 102
102
count

count

count
102
101 101 101

100 100
0.00 0.08 0.16 0.24 0.32 0.00 0.02 0.04 0.06 0.08 0.00 0.01 0.02 0.03 0.04 0.05
κ, eV λ, eV |SOC|, eV

Figure 1: Distribution of LVC and SOC coupling parameters: (a) κ = |κm,ξ |, (b) λ = |λmn,ξ |,
(SOC)
and (c) |SOC| = |Vmn |.

Since it is known that TM compounds often feature a rather rapid ISC

dynamics, it is interesting to inspect the crossing points between singlet and
triplet states. According to the LVC Hamiltonian, for a pair of zeroth-order
states, they are given by qξ∗ = (En − Em )/(κm,ξ − κn,ξ ). The distribution of
crossing points for FePS is shown in Figure 2b. To put these numbers into
perspective, one should note that the mode elongation at the zero point energy
ranges between ±1 (in dimensionless units). The actual values for FePS are
mostly much larger. The reason is to be found in the rather small displacements,
cf. Figure 1. This observation is no contradiction to a rapid ISC since the
potential curves run almost parallel until they eventually cross. Since SOC
is assumed to be not coordinate dependent, this facilitates an efficient ISC if
energies of spin-free states are energetically not too far apart.

3.2. TDH Dynamics

Obtaining converged MCTDH dynamics for the given model on a few pi-
cosecond time scale should be possible with the multilayer approach but is not
the focus of this contribution. Instead, we will use the TDH limit of the MCTDH
ansatz, which allows for full-dimensional quantum simulations [40]. The TDH
limit facilitates the exploration of overall aspects of the dynamics of different
reduced models on equal footing. Note that for ML-MCTDH calculations, each

6
 1000
(a)

800

600
counts

400

200

0
1.6 1.8 2.0 2.2 2.4 2.6 2.8
energy (eV)

400 (b)

300
counts

200

100

0
-40 -20 0 20 40
q*

Figure 2: Analysis of the LVC model for FePS [24]. (a) Distribution of diabatic electron-
vibrational spin-free singlet and triplet states (bin width 5 meV). (b) Distribution of state
crossing for all possible transitions between spin-free singlet and triplet states (bin width 0.5).

model should be converged separately. The analysis in the previous section

showed that the LVC model is characterized by a high density of states and a
large number of couplings terms having a comparable magnitude. In this sec-
tion, we will address the question of whether this situation allows for a model
reduction by defining certain cut-offs for the coupling parameters.
Figure 3 compares the experimental absorption spectrum with the one cal-
culated with the correlation function approach using TDH dynamics of the full
model. Further, an assignment of spin-free electronic transitions is given. The
overall lineshape of the spectrum is reasonably reproduced; in addition to the
inherent broadening due to the large density of states, we have added a phe-
nomenological broadening of γ −1 = 10 fs to account for solvation effects. As
compared with experiments, there is a shift of about 0.3 eV. Since this range of
the experimental spectrum does not convey many details, model reduction will
be investigated for laser-driven dynamics.
For the following simulations, we will consider the x-component of the dipole

7
Figure 3: Upper panel: Absorption spectrum of FePS as calculated by the TDH approach
to the dipole autocorrelation function as well as experimental results [23]. Also shown is the
laser pulse spectrum (gray shaded area). Lower panels: Decomposition of electronic transitions
based on density matrix analysis [41].

operator only (which is of higher magnitude than the other components). Fur-
ther we assume a Gaussian pulse shape, i.e.

E(t) = E0 cos(Ωt) exp{−(t − t0 )2 /2σ 2 } . (11)

Here, Ω, t0 , and σ are the carrier frequency, the pulse center, and the pulse
width, respectively. The transition energy at the peak of absorption due to
the x-component of the dipole is at 2.27 eV. For the width parameter, we use
24 fs yielding a full-width a half maximum of the pulse intensity of 40 fs. The
field amplitude has been chosen as 0.001 Hartree/Bohr, which gave a population
transfer of about 25-30% to the excited singlet states for the TDH dynamics.
Figure 4 compares the population of all singlet states (except the ground
state) for different models. The reference case (full model) is characterized by
an initial rise due to the laser pulse and an immediate rapid decay. (Note that
this decay is not due to stimulated emission.) Subsequently, one observes a
slower decay with a partial revival around 600-700 fs as well as superimposed
rapid oscillations. The dynamics of the other three models shown in the figure
essentially agree during the first 200 fs. Putting the threshold for the nonadia-
batic couplings to |λmn,ξ | = 1 meV (I, blue line), the decay after 200 fs is a bit
slower, and the partial revival doesn’t occur. Using, in addition, the threshold
|κm,ξ | = 5 meV (II, orange line) brings back the revival. Putting in addition a
(SOC)
threshold for SOC of |Vmn | = 0.5 meV (III, red line), there is again no pro-
nounced revival. It should be stressed that these parameters have been taken

8
 0.3
ref
I
II
0.25 III

singlet population
0.2

0.15

0.1

0.05

0
0 200 400 600 800 1000
time (fs)

Figure 4: TDH dynamics of summed excited singlet state populations after excitation with
a x-polarized Gaussian laser pulse centered at t0 = 100 fs (h̄Ω = 2.27 eV, σ = 24 fs, E0 =
0.001 Hartree/Bohr). The different curves correspond to the full model (denoted ’ref’) and
various reduced modes as explained in the text.

after consideration of certain larger cut-off values such as |λmn,ξ | = 5 meV,

which gave a pronounced difference in the dynamics. (Note that for the pur-
pose of illustration, it was not attempted to pinpoint these parameters in more
detail.) At first glance, the small magnitude of these cut-off parameters is a
bit surprising. In fact, inspecting Figure 1 one could have been tempted to
select just a few modes having with the largest couplings. Nevertheless, in the
present case of FePS, the dynamics after laser excitation appears to be shaped
by the combined effect of many small-magnitude couplings. We note that the
absorption spectrum is essentially indistinguishable for these four models. The
resulting LVC Hamiltonian with all three coupling parameters reduced contains
2941 κm,ξ terms, 13243 λmn,ξ terms, and 2166 SOC terms.
Next, we explore a possible reduction of the number of vibrational DOFs.
Intuitively one would be tempted to compare the coordinate expectation val-
ues with the values of the crossing points qξ∗ between singlet and triplet states.
However, according to Figure 2b, the majority of qξ∗ takes rather large values
which are not reached during the dynamics. As a technical note, we would like
to mention that for the simulation, the single-set method has been used where
a given SPF is optimized for all electronic states and coordinate expectation
values are calculated accordingly. Figure 5 compares the reduced model iden-
tified above with models having different cut-offs for the average normal mode
elongation. Namely, the cut-off has been set to 0.5 (24 modes), 0.2 (70 modes),
and 0.1 (131 modes). First, we note that in contrast to Figure 4, the initial
dynamics during the first 200 fs is influenced by the cut-off. In fact, inspecting
the absorption spectrum shown in the inset, one notices marked changes upon
reduction of the number of modes. The change in cut-off modifies the overlap

9
 0.3
ref (212 modes)
55 modes
70 modes
131 modes
0.25

absorption (arb. un.)

5
singlet population
0.2 4

0.15 1

0
1.6 1.8 2 2.2 2.4 2.6 2.8 3

energy (eV)
0.1

0.05

0
0 200 400 600 800 1000
time (fs)

Figure 5: Comparison of the population dynamics for model III of Fig. 4 (212 modes) with
various models having a reduced number of vibrational modes as indicated; for selection
details see text. The inset shows the absorption spectrum for the four modes according to the
x-component of the dipole.

with the pulse spectrum and, thus, the wave packet created on the excited state
manifold. To have reasonably accurate initial excitation and, at the same time,
a considerable reduction of the number of modes, we identified the model with
70 modes as a compromise. Again, we emphasize that the thresholds are rather
small and not that much different, i.e. it is the net effect of a large number of
small-amplitude vibrations, which is responsible for the observed dynamics.

3.3. ML-MCTDH Dynamics

TDH is the mean-field approximation to MCTDH and, as such, neglects
correlations in the dynamics of the different DOFs. One could argue that in
cases of a large number of small couplings combined with a large density of
states, the mean-field approximation should work well. However, since FePS
already gave some surprises, it will be better to compare TDH with MCTDH
dynamics, in particular, since ML-MCTDH provides a means to obtain numer-
ically converged results. For the purpose of illustration, we will use the final
70 modes model obtained from the TDH simulations in the previous section.
ML-MCTDH requires to identify a ML-tree to structure the wave packet. This
procedure is a nontrivial task, and different trees might lead to rather different
CPU usage for obtaining the same accuracy [42]. Presently, we did not attempt
to optimize the tree structure. Instead, we started with two branches, one in-
cluding all modes with an average elongation beyond 0.5 (21 modes) and one
including the remaining 49 modes. The first branch was further subdivided into
three branches, each containing two two-dimensional and one three-dimensional

10
particle. The second branch was first subdivided into four subbranches, each of
them containing three further branches. On the lowest level, there were 23 two-
dimensional and one three-dimensional particle. The number of SPFs in each
layer has been chosen such that the largest population of the least occupied
natural orbital was smaller than 2%.
0.25 0.14
(a) (b)

0.12
0.2

0.1
singlet population

1
MLCT
0.15

population
0.08

ML-MCTDH
3
0.06 MLCT
0.1

0.04

0.05 TDH 3
MC
0.02
1
MC
0 0
0 200 400 600 800 1000 0 200 400 600 800 1000
time (fs) time (fs)

Figure 6: Comparison between TDH and ML-MCTDH dynamics after laser pulse excitation
for the reduced model introduced in Section 3.2. (a) Total population of excited singlet states.
(b) Contributions of individual types of states to the ML-MCTDH dynamics as indicated.

Figure 6a compares TDH and ML-MCTDH dynamics for the total excited
singlet state population. Even though the linear absorption spectrum is es-
sentially identical, the difference in population dynamics is dramatic. First,
ML-MCTDH does not show the pronounced initial peak during and immedi-
ately after the laser pulse as is observed for TDH. More notably, however, the
decay of the total singlet state population is much slower in the ML-MCTDH
case. In fact, this brings the present simulation in closer agreement with the
TSH results reported in Ref. [24]. Of course, one might argue that at this point,
one should go back to the model reduction performed in Section 3.2 and to re-
peat it using ML-MCTDH. However, we don’t expect a qualitative difference,
i.e., the illustration of the challenges to model reduction performed with TDH
will still be valid.
In Figure 6b, ML-MCTDH populations of individual sets of states are shown.
The rapid initial population of triple states is typical for TM compounds and es-
sentially occurs during the laser pulse action. Here, we notice that the most/least
population is in the 1 MLCT/1 MC states. The population of 3 MLCT/3 MC
states is about equal and in between the 1 MLCT and 1 MC cases. This general
behavior is also in accord with the TSH results of Ref. [24]. Note that one
would not expect a quantitative agreement, not at least because the description
of the laser pulse excitation differs in the two simulations (cf. systematic study
of this aspect in Ref. [43]).

11
4. Conclusions

We have studied the laser-driven dynamics of the exemplary FePS com-

plex whose photophysical and quantum chemical properties were previously re-
ported along with TSH simulations. FePS shows the high density of electron-
vibrational states typical for such TM complexes. Using an LVC Hamiltonian,
the reduction of the number of electronic coupling parameters was explored,
followed by a reduction of the number of vibrational modes. The structureless
linear absorption spectrum was not useful in this respect, and, therefore, the
TDH dynamics was analyzed. It was found that the dynamics is influenced by
a large number of small couplings, which prevents the identification of a rep-
resentative model with just a few DOFs. For the thus established model, still
containing 70 vibrational modes and thousands of coupling terms, the effect of
correlations on the dynamics was investigated. Here one would have expected
that given the nature of the model, effects beyond the mean-field description
play a minor role only. Surprisingly, only the inclusion of correlations brings
the dynamics in good agreement with previously reported TSH simulations and
experiments.
We expect that the challenges for model reduction discussed here for an
exemplary case are typical for similar TM complexes.

Acknowledgement

The authors are grateful to Dr. J.-P. Zobel for providing the LVC Hamil-
tonian and L. Kleindienst for his support with preparing MCTDH input files
from LVC data. This work has been financially supported by the Deutsche
Forschungsgemeinschaft [Priority Program SPP 2102 “Light controlled reactiv-
ity of metal complexes” (grant KU952/12-1)].

References

[1] A. J. Esswein, D. G. Nocera, Hydrogen Production by Molecular Photo-

catalysis, Chem. Rev. 107 (2007) 4022–4047.
[2] R. E. Blankenship, D. M. Tiede, J. Barber, G. W. Brudvig, G. Fleming,
M. Ghirardi, M. R. Gunner, W. Junge, D. M. Kramer, A. Melis, T. A.
Moore, C. C. Moser, D. G. Nocera, A. J. Nozik, D. R. Ort, W. W. Parson,
R. C. Prince, R. T. Sayre, Comparing Photosynthetic and Photovoltaic
Efficiencies and Recognizing the Potential for Improvement, Science 332
(2011) 805–809.
[3] J. Jiao, J. He, M. Li, J. Yang, H. Yang, X. Wang, S. Yang, A porphyrin-
based metallacage for enhanced photodynamic therapy, Nanoscale 14
(2022).
[4] M. A. Halcrow (Ed.), Spin-Crossover Materials: Properties and Applica-
tions, Wiley, Chichester, West Sussex, United Kingdom, 2013.

12
 [5] L. González, D. Escudero, L. Serrano-Andrés, Progress and Challenges in
the Calculation of Electronic Excited States, ChemPhysChem 13 (2012)
28–51.
[6] T. J. Penfold, E. Gindensperger, C. Daniel, C. M. Marian, Spin-Vibronic
Mechanism for Intersystem Crossing, Chemical Reviews 118 (2018) 6975–
7025.
[7] C. Daniel, Ultrafast processes: Coordination chemistry and quantum the-
ory, Phys. Chem. Chem. Phys. 23 (2020) 43–58.
[8] P. S. Wagenknecht, P. C. Ford, Metal centered ligand field excited states:
Their roles in the design and performance of transition metal based photo-
chemical molecular devices, Coord. Chem. Rev. 255 (2011) 591–616.
[9] J. P. Zobel, L. González, The Quest to Simulate Excited-State Dynamics
of Transition Metal Complexes, JACS Au 1 (2021) 1116–1140.
[10] J. C. Tully, Molecular dynamics with electronic transitions, The Journal of
Chemical Physics 93 (1990) 1061–1071.
[11] S. Mai, P. Marquetand, L. González, Nonadiabatic dynamics: The SHARC
approach, WIREs Comput. Mol. Sci. 8 (2018) e1370.
[12] S. Gómez, M. Heindl, A. Szabadi, L. González, From Surface Hopping to
Quantum Dynamics and Back. Finding Essential Electronic and Nuclear
Degrees of Freedom and Optimal Surface Hopping Parameters, J. Phys.
Chem. A 123 (2019) 8321–8332.

[13] M. Pápai, Toward Simulation of Fe(II) Low-Spin → High-Spin Photo-

switching by Synergistic Spin-Vibronic Dynamics, J. Chem. Theory Com-
put. 18 (2022) 1329–1339.
[14] J. P. Zobel, M. Heindl, F. Plasser, S. Mai, L. González, Surface Hopping
Dynamics on Vibronic Coupling Models, Acc. Chem. Res. 54 (2021) 3760–
3771.
[15] F. Plasser, S. Mai, M. Fumanal, E. Gindensperger, C. Daniel, L. González,
Strong Influence of Decoherence Corrections and Momentum Rescaling
in Surface Hopping Dynamics of Transition Metal Complexes, J. Chem.
Theor. Comp. 15 (2019) 5031–5045.
[16] A. M. Virshup, J. Chen, T. J. Martı́nez, Nonlinear dimensionality reduction
for nonadiabatic dynamics: The influence of conical intersection topogra-
phy on population transfer rates, J. Chem. Phys. 137 (2012) 22A519.
[17] J. P. P. Zauleck, S. Thallmair, M. Loipersberger, R. de Vivie-Riedle, Two
New Methods To Generate Internal Coordinates for Molecular Wave Packet
Dynamics in Reduced Dimensions, J. Chem. Theor. Comp. 12 (2016) 5698–
5708.

13
[18] G. Capano, T. J. Penfold, M. Chergui, I. Tavernelli, Photophysics of a cop-
per phenanthroline elucidated by trajectory and wavepacket-based quan-
tum dynamics: A synergetic approach, Phys. Chem. Chem. Phys. 19 (2017)
19590–19600.
[19] X. Li, Y. Xie, D. Hu, Z. Lan, Analysis of the Geometrical Evolution in On-
the-Fly Surface-Hopping Nonadiabatic Dynamics with Machine Learning
Dimensionality Reduction Approaches: Classical Multidimensional Scaling
and Isometric Feature Mapping, J. Chem. Theor. Comp. 13 (2017) 4611–
4623.
[20] S. R. Hare, L. A. Bratholm, D. R. Glowacki, B. K. Carpenter, Low dimen-
sional representations along intrinsic reaction coordinates and molecular
dynamics trajectories using interatomic distance matrices, Chem. Sci. 10
(2019) 9954–9968.
[21] Y. Harabuchi, J. Eng, E. Gindensperger, T. Taketsugu, S. Maeda,
C. Daniel, Exploring the Mechanism of Ultrafast Intersystem Crossing
in Rhenium(I) Carbonyl Bipyridine Halide Complexes: Key Vibrational
Modes and Spin–Vibronic Quantum Dynamics, J. Chem. Theor. Comp. 12
(2016) 2335–2345.
[22] S. Mai, L. González, Identification of important normal modes in nonadia-
batic dynamics simulations by coherence, correlation, and frequency anal-
yses, J. Chem. Phys. 151 (2019) 244115.
[23] J. Moll, R. Naumann, L. Sorge, C. Förster, N. Gessner, L. Burkhardt,
N. Ugur, W. Nürnberger, P. Seidel, C. Ramanan, M. Bauer, H. K., Pseudo-
ocathedral iron(II) complexes with near-degenerate charge transfer and lig-
and field states at the Franck-Condon geometry, Chem. Eur. J. 28 (2022)
e2022ß1858.
[24] J. P. Zobel, A. Kruse, O. Baig, S. Lochbrunner, S. I. Bokarev, O. Kühn,
L. González, O. S. Bokareva, Can range-separated functionals be optimally
tuned to predict spectra and excited state dynamics in photoactive iron
complexes?, Chem. Sci. 14 (2023) 1491–1502.
[25] H. Köppel, W. Domcke, L. S. Cederbaum, Multimode molecular dynamics
beyond the Born-Oppenheimer approximation, Adv Chem Phys 57 (1984)
59.
[26] T. J. Penfold, E. Gindensperger, C. Daniel, C. M. Marian, Spin-Vibronic
Mechanism for Intersystem Crossing, Chem. Rev. 118 (2018) 6975–7025.
[27] V. May, O. Kühn, Charge and Energy Transfer Dynamics in Molecular
Systems, Fourth Edition, Wiley-VCH, Weinheim, 2023.
[28] M. Fumanal, E. Gindensperger, C. Daniel, Ultrafast Intersystem Crossing
vs Internal Conversion in α-Diimine Transition Metal Complexes: Quan-
tum Evidence, J. Phys. Chem. Lett. 9 (2018) 5189–5195.

14
[29] M. Pápai, T. Rozgonyi, T. J. Penfold, M. M. Nielsen, K. B. Møller, Sim-
ulation of ultrafast excited-state dynamics and elastic x-ray scattering by
quantum wavepacket dynamics, J. Chem. Phys. 151 (2019) 104307.
[30] T. J. Penfold, J. Eng, Mind the GAP: Quantifying the breakdown of the
linear vibronic coupling Hamiltonian, Physical Chemistry Chemical Physics
25 (2023) 7195–7204.
[31] K. Falahati, H. Tamura, I. Burghardt, M. Huix-Rotllant, Ultrafast carbon
monoxide photolysis and heme spin-crossover in myoglobin via nonadia-
batic quantum dynamics, Nat. Comm. 9 (2018) 4502.
[32] M. Fumanal, C. Daniel, E. Gindensperger, Excited-state dynamics of
[Mn(im)(CO) 3 (phen)] + : PhotoCORM, catalyst, luminescent probe?,
J. Chem. Phys. 154 (2021) 154102.
[33] H. D. Meyer, U. Manthe, L. S. Cederbaum, The multi-configurational time-
dependent Hartree approach, Chem. Phys. Lett. 165 (1990) 73–78.
[34] M. H. Beck, A. Jäckle, G. A. Worth, H.-D. Meyer, The multiconfiguration
time-dependent Hartree method: A highly efficient algorithm for propagat-
ing wavepackets, Phys. Rep. 324 (2000) 1–105.
[35] H. Wang, M. Thoss, Multilayer formulation of the multiconfiguration time-
dependent Hartree theory, J. Chem. Phys. 119 (2003) 1289–1299.
[36] H.-D. Meyer, Studying molecular quantum dynamics with the multicon-
figuration time-dependent Hartree method, WIREs Comput. Mol. Sci. 2
(2011) 351–374.
[37] O. Vendrell, H.-D. Meyer, Multilayer multiconfiguration time-dependent
Hartree method: Implementation and applications to a Henon-Heiles
Hamiltonian and to pyrazine, J. Chem. Phys. 134 (2011) 044135.
[38] H. Naundorf, G. A. Worth, H.-D. Meyer, O. Kühn, Multiconfiguration
Time-Dependent Hartree Dynamics on an ab Initio Reaction Surface: Ul-
trafast Laser-Driven Proton Motion in Phthalic Acid Monomethylester, J.
Phys. Chem. A 106 (2002) 719–724.
[39] G. A. Worth, M. H. Beck, A. Jäckle, H.-D. Meyer, The MCTDH Package,
Version 8.2, (2000). H.-D. Meyer, Version 8.3 (2002), Version 8.4 (2007). O.
Vendrell and H.-D. Meyer, Version 8.5 (2013), used Version 8.6.3 (2023),
See http://mctdh.uni-hd.de/.
[40] GK. Paramonov, H. Naundorf, O. Kühn, Ultrafast multidimensional dy-
namics of strong hydrogen bonds - A Cartesian reaction surface approach,
Eur. Phys. J. D 14 (2001) 205–215.
[41] F. Plasser, TheoDORE: A toolbox for a detailed and automated analysis of
electronic excited state computations, J. Chem. Phys. 152 (2020) 084108.

15
[42] M. Schröter, S. D. Ivanov, J. Schulze, S. P. Polyutov, Y. Yan, T. Pullerits,
O. Kühn, Exciton-vibrational coupling in the dynamics and spectroscopy
of Frenkel excitons in molecular aggregates, Phys. Rep. 567 (2015) 1–78.
[43] M. Heindl, L. González, Validating fewest-switches surface hopping in the
presence of laser fields, The Journal of Chemical Physics 154 (2021) 144102.

16