# Nuclear Ensemble Approach with Importance Sampling

Fábris Kossoski, Mario Barbatti

► To cite this version:

Fábris Kossoski, Mario Barbatti. Nuclear Ensemble Approach with Importance Sampling. Journal of Chemical Theory and Computation, 2018, 14 (6), pp.3173-3183. hal-01965459

HAL Id: hal-01965459

https://amu.hal.science/hal-01965459

Submitted on 26 Dec 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Nuclear Ensemble Approach with Importance Sampling

F. Kossoski*,†,‡ and M. Barbatti*,‡

†Instituto de Física “Gleb Wataghin”, Universidade Estadual de Campinas, 13083-859, Campinas, São Paulo, Brazil

‡Aix Marseille Univ, CNRS, ICR, Marseille, France

E-mail: fabris@if.usp.br; mario.barbatti@univ-amu.fr

## Abstract

We show that the importance sampling technique can effectively augment the range of problems where the nuclear ensemble approach can be applied. A sampling probability distribution function initially determines the collection of initial conditions for which calculations are performed, as usual. Then, results for a distinct target distribution are computed by introducing compensating importance sampling weights for each sampled point. This mapping between the two probability distributions can be performed whenever they are both explicitly constructed. Perhaps most notably, this procedure allows for the computation of temperature dependent observables. As a test case, we investigated the UV absorption spectra of phenol, which has been shown to have a marked temperature dependence. Application of the proposed technique to a range that covers 500 K provides results that converge to those obtained with conventional sampling. We further show that an overall improved rate of convergence is obtained when sampling is performed at intermediate temperatures. The comparison between calculated and the available measured cross sections is very satisfactory, as the main features of the spectra are correctly reproduced. As a second test case, one of Tully's classical models
was revisited, and we show that the computation of dynamical observables also profits
from the importance sampling technique. In summary, the strategy developed here can
be employed to assess the role of temperature for any property calculated within the
nuclear ensemble method, with the same computational cost as doing so for a single
temperature.

# 1 Introduction

In photophysics and photochemistry, the nuclear ensemble approach consists in representing
the molecular vibrational wave functions by a set of classical nuclear configurations. One can
say it is equivalent to the quasi-classical trajectory method 1-6 when applied to unimolecular
problems. For calculations of static properties (when there is no need for time propagation),
the term “trajectories” is somewhat misleading, and the former term might be more appro-
priate. Here we aim exclusively at unimolecular processes, namely, absorption spectroscopy
and excited-state dynamics, and since both static and dynamical properties are involved, we
sticked with the “nuclear ensemble approach” expression, rather than “quasi-classical trajec-
tory method”. For spectra calculations, there exists a construction that leads to the nuclear
ensemble method.⁷ On top of that, the method counts with a quite intuitive interpretation,
can be easily implemented, and can offer a very successful description of various types of
spectra. 7-11 Adiabatic and non-adiabatic dynamics that employ a classical nuclei approxi-
mation can also be framed as a nuclear ensemble approach. In this case, however, a set of
initial coordinates and momenta have to be sampled, while for spectra calculations, only a set of
nuclear configurations is required. For either static or dynamical simulations, one computes
a given physical observable of interest for a set of independent phase space points. These
follow a certain probability distribution function (PDF), which can be constructed based
on a quantum or classical description of the nuclei.¹² Once the calculations are done, the
individual computed values are simply averaged, which is taken as the estimate for its real
value.
While being a very straightforward approach, converged results require an extensive number of sampled nuclear configurations (or initial conditions for the dynamics), depending on the type of calculation and the desired degree of precision. Here we have considered the situation in which the same observable should be computed for different PDFs. In principle one could perform several sets of calculations, but the cost of doing so for each PDF clearly makes this attempt unaffordable. We propose a strategy that makes use of the importance sampling technique, allowing for the consideration of different target PDFs, while the sampling and the calculations are performed for only one PDF.

The first situation that comes to mind where the underlying PDFs change is when temperature is varied. An increase in temperature allows the nuclei to visit phase space regions further away from the equilibrium geometry, which should impact on the sampling probabilities. At a fundamental research level, a temperature dependent observable can reveal details about the physico-chemical process of interest, and examples of that are endless. We briefly review some specific motivation related to the role of temperature on molecular photoabsorption, since this topic will be discussed in more depth here. Modeling the dynamics and interactions of atmospheric pollutants requires the understanding of many different processes. Of special interest are the rates for photolysis, which is very often the dominant decaying mechanism.<sup>13-15</sup> Therefore, it is critical to understand how the molecule absorbs a photon at the UV-visible spectral range and to have reliable absorption cross sections. Its dependence on wavelength and temperature are equally important pieces of information,<sup>13-15</sup> despite being usually overlooked. Accounting for the temperature effects on the cross sections can lead to considerable modifications on the lifetimes of atmospherically relevant compounds.<sup>15,16</sup> As a related motivation, absorption spectroscopy is one of the main tools used for combustion and gasification analysis, and for detecting and measuring concentration of pollutants in flue gases.<sup>17</sup> However, the lack of data on the temperature dependent spectra of important industrial compounds limits its application to lower temperatures. Since elevated temperatures is a commonplace for these applications, more expensive and complicated techniques must be
devised instead.¹⁸ More complete databases on higher-temperature cross sections should help promoting the use of absorption spectroscopy in different industrial applications. Furthermore, we mention that absorption of UV radiation is of central importance in astrochemistry, and once more, temperature effects can play key roles.¹⁹⁻²¹ Despite several experimental results that have one of the above mentioned points as a motivation,¹³⁻¹⁵,¹⁷,¹⁹,²²,²³ there are very few results from the side of theory.²⁴,²⁵

As a first test case for the proposed importance sampling strategy, we investigated the temperature dependence of the photoabsorption spectra of phenol, in the 4.0–6.1 eV interval. This particular choice was motivated by the available experimental data¹⁷ that covered a wide temperature range (296–773 K), showing a quite noticeable temperature effect. The electronic excited states of phenol have been extensively studied, by means of vertical excitations calculations,²⁶⁻³¹ absorption spectroscopy,¹⁷,³¹⁻³⁴ and electron energy loss spectroscopy.³⁰,³⁵ From the theoretical side, results on the absorption spectra are still lacking. More generally, we notice a huge gap from theory regarding the role of temperature on UV-vis absorption cross sections. While this first test case involves the calculation of a static property, the importance sampling strategy should also be valid for dynamics simulations. Thus, in order to show the general applicability of the method, we have investigated one of Tully’s classical models³⁶ as a second test case. Our main goal here concerns the presentation and validation of the idea that simulations based on nuclear ensembles can largely profit from the importance sampling technique. A second goal is to describe the temperature dependent absorption spectra of phenol, by comparing it with the existing experimental data.

# 2 Theory and calculations

## 2.1 Nuclear ensemble approach

The central concept behind the nuclear ensemble approach is that a physical observable $\mu$ can be computed by integrating its value over a set of representative nuclear configurations.⁷
This statement can be cast into the following expression:

$$ \mu = \int f(\mathbf{x}) P(\mathbf{x}) d\mathbf{x} = E_P[f(\mathbf{X})], \qquad (1) $$

where $f$ is the function to be integrated in phase space $\mathbf{x} = (\mathbf{q}, \mathbf{p})$, $\mathbf{q}$ and $\mathbf{p}$ representing the set of nuclear coordinates and momenta, respectively. These, in turn, follow a normalized PDF $P(\mathbf{x})$, while $E_P(\cdot)$ indicates the expectation value when its argument follows this PDF. The capital letter is used here to represent that the variable should be sampled according to a given distribution, i.e., $(\mathbf{X} \sim P)$.

In practice, one actually computes an estimate $\hat{\mu}$ by performing an average over a set of $N_p$ sampled points:

$$ \hat{\mu} = \frac{1}{N_p} \sum_{l=1}^{N_p} f(\mathbf{X}_l), \quad \mathbf{X} \sim P. \qquad (2) $$

For the calculation of static observables, where $f$ depends only on the set of initial coordinates $\mathbf{q}$, the momenta dependence of the PDF can be integrated out. Whenever the PDF is separable into coordinate and momenta, which is often the case, this simplifies to:

$$ \mu = \iint f(\mathbf{q}) P(\mathbf{q}, \mathbf{p}) d\mathbf{q} d\mathbf{p} = \int f(\mathbf{q}) P(\mathbf{q}) d\mathbf{q}, \qquad (3) $$

where it is implicit that the second $P$ represents a function of $\mathbf{q}$ only.

The PDF $P(\mathbf{x})$ can be built based on a quantum or classical description of the nuclei. The formally correct one is to use a quantum quasi-probability distribution function, and the Wigner distribution tends to be the natural choice.³⁷⁻³⁹ There are also other possibilities to represent the quantum distribution, based on the Husimi distribution, for example.⁴⁰,⁴¹ An alternative to the use of these quasi-distributions is to ignore the interdependency of the conjugated coordinate and momenta, and to use instead the coordinate and momenta wavefunctions as the PDFs.⁴²,⁴³ Once the distribution is constructed, sampling of each coordinate and momenta pair can be performed in a correlated on uncorrelated way.³⁹ In the
former case, only coordinate is sampled, while momenta is selected (with a random direction) so that the energy of the mode matches the exact quantum value. In the second option, sampling is performed independently for each variable, and a distribution of vibrational energies is obtained. The molecular vibrational states are usually described as a set of normal modes obtained in the harmonic approximation, for which the Wigner distribution is known analytically. For a system at a finite temperature $T$, the molecule is in a mixture of vibrationally excited states, weighted by the Boltzmann factors. It can be shown that the harmonic oscillator, normalized Wigner distribution for finite temperatures is given by:³⁷

$$P(\mathbf{q}) = \prod_{i=1}^{M} \left( \frac{\alpha_i}{\pi \hbar} \right) \exp \left( - \frac{\alpha_i \mu_i \omega_i q_i^2}{\hbar} - \frac{pi^2}{\mu_i \omega_i \hbar} \right) \quad (4)$$

where the product runs over $M$ vibrational modes, $q_i$ is the coordinate, $p_i$ is the momentum, $\mu_i$ is the reduced mass and $\omega_i$ is the harmonic frequency of the $i$-th mode. The temperature $T$ appears in the factor $\alpha_i = \tanh(\hbar\omega_i/2k_B T)$, where $k_B$ is Boltzmann constant. The quantum sampling method is usually preferred for small to medium sized molecules, when the calculation of frequencies and normal modes is cheap in comparison to the set of calculations for the ensemble of configurations. One drawback is that it does not account for anharmonic effects, which are particularly important whenever one has many low-frequency modes, as is the case for larger molecular systems. We note, however, that this limitation stems from the harmonic description of the vibrational states, and not from the use of a distribution that tries to mimic the quantum nature of the nuclei.

An alternative way to perform the sampling does not require an explicit functional form for the PDF. Instead, one can perform molecular dynamics for the electronic source state, and configurations are sampled as a series of snapshots from the simulation. The dynamics can be performed in the canonical or microcanonical ensemble, and with the average energy corresponding to the classical thermal energy or to the zero-point energy (plus thermal contribution if the case). Such classically-based sampling has the advantage that anharmonicities are naturally accounted for, and it is usually the preferred choice when dealing
with larger compounds or simulations in condensed phase.⁴⁴ However, since this method
does not sample from classically forbidden regions, the molecular motion is unrealistically
restricted, which might lead to artifacts in the calculations. Besides, it completely ignores
zero-point energy effects. A related way to perform the sampling is to rely on the classical
distribution of the harmonic potential. Results computed with different sampling methods
for the initial conditions⁸,¹¹,³⁹,⁴⁵⁻⁴⁷ show that overall better agreement with experiment is
achieved when a quantum-based distribution is used, instead of the classically-based ther-
mal sampling alternative.

## 2.2 Importance sampling

Suppose we have already computed the observable $\hat{\mu}$ for a given PDF $Q(\mathbf{x})$ and now we would like to evaluate it for a distinct PDF $P(\mathbf{x})$. This would be the case if we are interested in evaluating how $\hat{\mu}$ varies with temperature. There are other situations which can also be framed in terms of a change in the PDF, which will be discussed later. The most immediate way to proceed would be to sample a separate set of points from the target distribution $P(\mathbf{x})$, and evaluate the function $f(\mathbf{x})$ for these new configurations. This can rapidly become very expensive, even unfeasible if many distinct PDFs are to be sampled (a wide range of temperatures, for instance). However, in most cases of interest the PDFs are not expected to differ significantly, such that the probabilities for sampling any particular point should not vary much from one PDF to another. At this point, intuition tells that we would be overshooting in performing many rounds of samplings and calculations for each distinct PDF. There should be a way to better exploit the invaluable calculations already performed for the first PDF, and somehow reuse it for other distributions of interest.

We propose a strategy that does exactly that, which is based on the variance reduction
importance sampling (IS) technique.⁴⁸ Basically, after the function $f$ is evaluated for a set
of points $\mathbf{X}_t$ obtained from a given sampling distribution $Q(\mathbf{x})$, the expected value $\hat{\mu}$ can be
readily computed for a distinct target distribution $P(\mathbf{x})$ by introducing appropriate weights
for each sampled point. The following expression synthesizes the idea:

$$ \mu = \int f(\mathbf{x}) P(\mathbf{x}) d\mathbf{x} = \int f(\mathbf{x}) \frac{P(\mathbf{x})}{Q(\mathbf{x})} Q(\mathbf{x}) d\mathbf{x} = \\ = \int f(\mathbf{x}) w(\mathbf{x}) Q(\mathbf{x}) d\mathbf{x} = E_Q [f(\mathbf{X}) w(\mathbf{X})]. \qquad (5) $$

The sampling distribution $Q(\mathbf{x})$ is introduced in both numerator and denominator. In the numerator it plays the role of the actual PDF from which $\mathbf{X}$ has already been sampled, while in the denominator it is combined with $P(\mathbf{x})$ and gives rise to the IS weighting factors $w(\mathbf{x}) = P(\mathbf{x})/Q(\mathbf{x})$. Therefore, the sampling and the function evaluations are performed only once, according to the sampling distribution $Q(\mathbf{x})$, and for each target distribution $P(\mathbf{x})$ of interest one simply needs to apply the weights $w(\mathbf{x})$. We stress that no approximation is involved here, as the introduced weights exactly compensate for the different sampling and target distributions. For an actual collection of sampled points, the expected value is thus given by:

$$ \hat{\mu} = \frac{1}{N_p} \sum_{l=1}^{N_p} f(\mathbf{X}_l) w(\mathbf{X}_l), \mathbf{X} \sim Q. \qquad (6) $$

Even though our strategy is based on the IS technique, there is a key distinction between our motivation to use it and the customary reason to do so.⁴⁸ Usually, IS is employed when the PDF $P$ provides a poor sampling for the function to be evaluated $f$, the common example being when only a small portion of the sampling space contribute decisively to the integral. Therefore, one needs to search for a better suited PDF $Q$, that concentrates the sampling on the important regions of the integrand, hence the name of the technique. In some sense our case is the opposite of that, since for each target distribution $P$ we actually use a worse sampling distribution $Q$. The great advantage is that we sample and perform the expensive calculations of $f(\mathbf{X})$ only once, allowing the IS-based strategy to evaluate the integral for any other target distribution. In order to further illustrate this point, we
comment on previous applications of the IS technique for bimolecular collisions<sup>4,49,50</sup> as an illustrative example. In these problems, one is interested in computing rate constants, which require performing dynamics simulations for various sampled collision energies. However, reactions might not take place at all for low collision energies, and therefore it would be wasteful to perform costly trajectory simulations in this regime. With this motivation, it has been proposed to sample the collision energies from a biased distribution (avoiding lower energies), and then using the IS technique accordingly.<sup>4,49,50</sup> This has the great advantage of improving the statistics of the final results, but normally one would still have to repeat such calculations for each new sampling distribution. Our proposed usage of the IS tool allows for a straightforward mapping between the set of results from any pair of distributions of interest. It should be possible, for example, to investigate the temperature dependence of the rate constant while sampling at one single temperature, which is fundamentally different from the previous applications of IS to the problem of bimolecular collisions.<sup>4,49,50</sup>

We will discuss in greater detail results for the temperature dependence of a static observable (the absorption cross section). In this case, the difference between sampling and target PDFs lies only on the temperature, and working with the Wigner distribution (eq. 4), the IS weights can be computed as:

$$w(\mathbf{q}) = \prod_{i=1}^{M} \sqrt{\frac{\alpha_i^t}{\alpha_i^s}} \exp \left[ -(\alpha_i^t - \alpha_i^s) \frac{\mu_i \omega_i q_{il}^2}{\hbar} \right] \quad (7)$$

where $T^s$ is the sampling temperature, $T^t$ is the target temperature, $\alpha_i^s = \tanh(\hbar\omega_i/2k_B T^s)$ and $\alpha_i^t = \tanh(\hbar\omega_i/2k_B T^t)$. We note that the momentum dependence was integrated out, as in eq. 3. If a dynamical observable should be computed, then the IS weights would be given by:

$$w(\mathbf{q}, \mathbf{p}) = \prod_{i=1}^{M} \frac{\alpha_i^t}{\alpha_i^s} \exp \left[ -(\alpha_i^t - \alpha_i^s) \left( \frac{\mu_i \omega_i q_i^2}{\hbar} + \frac{p_i^2}{\mu_i \omega_i \hbar} \right) \right] \quad (8)$$
## 2.3 Photoabsorption spectra calculation

Within the nuclear ensemble approach, the photoabsorption cross section is computed according to:⁷

$$ \sigma(E) = \frac{\pi e^2 \hbar}{2mc\varepsilon_0 E} \sum_n^{N_{\text{fs}}} \frac{1}{N_p} \sum_l^{N_p} \Delta E_{0n}(\mathbf{q}_l) f_{0n}(\mathbf{q}_l) \times g(E - \Delta E_{0n}(\mathbf{q}_l, \delta)) \quad (9) $$

where $\Delta E_{0n}$ is the excitation energy and $f_{0n}$ is the oscillator strength of the n-th state, $N_{\text{fs}}$ is the number of states for which excitations are considered, $N_p$ is the number of nuclear configurations and $g$ is a phenomenological normalized line profile, usually taken as a Gaussian or a Lorentzian. Here we employed a Lorentzian line profile with $\delta = 0.05$ eV.

Excited states were described with the algebraic diagrammatic construction to the second order (ADC(2)),⁵¹,⁵² with the resolution-of-the-identity approximation.⁵³ Møller-Plesset Perturbation Theory (MP2) was used for the geometry optimization of the ground state. Electronic vertical excitations and oscillator strengths were computed for 12 states of singlet multiplicity. We employed the aug-cc-pVDZ basis set supplemented with doubly augmented functions for the hydrogen of the hydroxyl group, and for the carbon atoms closest and furthest from the it. For the energy range we are interested (up to 6.1 eV) this intermediate (d)-aug-cc-pVDZ set provided vertical excitation energies identical to the ones obtained with the full d-aug-cc-pVDZ to within 0.01 eV, while counting with only 78% the number of basis functions. The Turbomole software⁵⁴ was used for all electronic structure calculations.

4400 nuclear configurations were sampled according to the finite-temperature uncorrelated Wigner distribution (eq. 4), for three different temperatures: 296 K, 573 K and 773 K. The normal modes were obtained at the same level as the geometry optimization. Both the sampling and the absorption cross sections calculations were performed with the Newton-X program.⁵⁵,⁵⁶ The cross sections were also computed for target temperatures different from the sampling temperatures according to the proposed IS strategy. This allowed us to eval-
uate whether and how efficiently the IS results converge towards the conventional sampling
results.

## 2.4 Dynamics calculations

As already pointed out, the IS technique can be adopted to any method that is based on
an ensemble of independent points in phase space, each point representing a set of nuclei
positions and momenta. This is the case for the popular fewest-switches trajectory surface-
hopping method,³⁶ which we have employed for our dynamics case study. We have investi-
gated the first of Tully's classical models,³⁶ with the original set of parameters. It consists in
two states with parallel potential energy curves, except at the region where they interact and
form a single avoided crossing (with a corresponding peaked non-adiabatic coupling). The
system is initially in the ground state in the asymptotic region, moving towards the avoided
crossing. The initial momentum controls the hopping probabilities in the interaction region,
which in turn determine the final probabilities for transmission and reflection on the lower
state and for transmission on the upper state.³⁶

Since we are interested in evaluating how a set of results can be mapped between different
sampling and target distributions, the natural choice is to assess scenarios with varying
temperatures. While for a harmonic potential the finite-temperature Wigner distribution is
given by a Gaussian, the free-particle case is of course governed by the Maxwell-Boltzmann's
PDF. Any physical observable can be computed with the aid of eq. 6, where the IS weights
are given by the ratio of two Maxwell-Boltzmann's PDFs, in analogy to eq. 7.

Three temperatures were considered (4,000 K, 10,000 K and 15,000 K), which correspond
to average momenta around 8.0, 12.7 and 15.6 (in atomic units). That is the same range of
momenta as in Tully's tests.³⁶ For each temperature, 30,000 trajectories were propagated,
all of them starting far away from the interacting region (at -6.0a₀), and with velocities
sampled according to the Maxwell-Boltzmann's distribution. The dynamics run for 100 fs,
with a time step of 0.5 fs, which was small enough to ensure good energy conservation even
for the most energetic trajectory. The time-dependent Schroedinger equation was integrated with a 0.025 fs time step (with cubic interpolation from the classical dynamics), while a decoherence correction ($\alpha$ = 0.1 hartree) was employed.⁵⁷ The calculations were performed with the Newton-X program.⁵⁵,⁵⁶

# 3 Results and discussion

## 3.1 Importance sampling vs. conventional sampling

The absorption cross sections computed for each sampling temperature (296 K, 573 K and 773 K) were taken as a reference result ($\sigma^0$) for comparison to the IS results. For each sampling temperature, the IS strategy was applied to compute the cross sections $\sigma$ for the other two target temperatures. We then compared both set of results by computing the normalized root mean square deviation (NRMSD):

$$ \text{nrms} = \frac{1}{\sigma_{\text{max}}^0} \sqrt{\frac{\sum_{k=1}^{N_{\text{pts}}} [\sigma^0(E_k) - \sigma(E_k)]^2}{N_{\text{pts}}}} \quad (10) $$

where $\sigma_{\text{max}}^0$ is the maximum value of the reference cross section and the sum runs over the grid of $N_{\text{pts}}$ energy points $E_k$. Since we are interested in describing the two first absorption bands of phenol, which are located in the 4–6.1 eV range, only points lying inside this interval entered in the calculation of the NRMSDs.

The cross sections for the nine combinations of sampling and target temperatures are displayed in Fig. 1. By visually inspecting the curves, there is already a good indication that both conventional and importance sampling approaches converge to the same results. This is confirmed in Fig. 2, which shows for a sampling temperature of 773 K that the NRMSDs consistently decreases as the number of sampling points increase. The same behavior was observed for the other sampling temperatures. The comparison among cross sections com-
puted for the same target temperature and different sampling temperatures (see Fig. 1) also reveals that the curves become noisier as these temperatures get apart. This can be seen more directly in Fig. 2, which show that the NRMSDs are systematically smaller when target and sampling temperatures are the same, increasing as they deviate from each other. We also notice that a polynomial fit (in the 20–3500 interval) to the NRMSD of the conventional sampling curve gives a $N_p^{-0.53}$ behavior, very close to the theoretically expected $N_p^{-1/2}$ scaling.

<figure><img src="image_2.png" /><figcaption>Figure 1: Absorption cross sections of phenol, computed for combinations of 296 K, 573 K and 773 K as sampling ($T^s$) and target ($T^t$) temperatures.</figcaption></figure>

When using the IS strategy, the set of sampled configuration is obviously sub-optimal if the target PDF differs from the sampling PDF. Therefore, the price one pays when extrapolating the set of results from the later to the former distributions is a loss in the rate of convergence of the results, as clearly seen in Fig. 2. Following the same definition usually employed in IS, we introduce the effective sample size as:

$$N_p^{\text{eff}} = \frac{\left( \sum_{l=1}^{N_p} w(\mathbf{X}_l) \right)^2}{\sum_{l=1}^{N_p} w(\mathbf{X}_l)^2} \qquad (11)$$
<figure><img src="image_3.png" /><figcaption>Figure 2: NRMSDs computed for three target temperatures, where the sampling was performed at 773 K. On the left, the x-axis represents the actual number of configurations, and on the right it denotes the effective number of configurations, as defined in eq. 11.</figcaption></figure>

This quantity is therefore related to the distribution of IS weights $w(\mathbf{X}_l)$. The closer sampling and target PDFs are, the sharper will be the distribution of $w(\mathbf{X}_l)$ around 1, such that $N_p^{\text{eff}}$ approaches $N_p$. On the other hand, if IS is used for two quite distinct PDFs, most of the configurations will be of very little importance, while a few will contribute with significant weights, and in this case $N_p^{\text{eff}}$ would be much smaller than $N_p$. In Fig. 2, we also show how the NRMSDs converge as a function of the effective sample size, which was done by simply rescaling the x-axis by the ratio $r = N_p^{\text{eff}}/N_p$. The rate of convergence becomes the same, regardless of the combination of sampling and target PDFs, which supports the interpretation that $N_p^{\text{eff}}$ represents an effective number of configurations. In other words, this is the number of points that would have to be sampled from the target PDF in order to achieve the same level of precision (in a statistical sense) as that obtained from the actual sampling PDF.

All our IS results do converge towards the reference results, but it is also important to understand how sampling and target temperatures influence the convergence. Using the above expressions for $N_p^{\text{eff}}$ and $w(\mathbf{X}_l)$ for the case of phenol, we evaluated how the ratio $r = N_p^{\text{eff}}/N_p$ behaves as a function of the target temperature, for the three sampling temperatures,
and the results are displayed in Fig. 3. When target and samplings temperatures coincide, $r = 1$, and as they depart from each other $r$ gets smaller and smaller. At the lowest temperatures, it stops dropping and converges to a fixed value at 0 K, which reflects the fact that the occupation numbers of the vibrational levels change very slowly at this regime. Despite the apparent Gaussian profile of the curves, there is in fact an asymmetry that makes the drop faster towards higher temperatures. For a target temperature of 534.5 K, midway the lowest and highest sampling temperatures considered, $r = 0.38$ when sampling is performed at 296 K, a quite smaller ratio than when it is done at 773 K ($r = 0.60$). Yet, when sampling and target temperatures differ considerably (again, 296 K and 773 K), $r$ is rather small nevertheless, whether sampling is performed at the low ($r = 0.07$) or at the high ($r = 0.19$) temperature. A much better way to proceed is to perform the sampling at some intermediate temperature. We did so for 573 K (results also shown in Fig. 3), simply because the spectra was also measured at this temperature, and the ratios for both the lowest ($r = 0.50$) and the highest temperatures ($r = 0.62$) increase a lot and also become closer to each other. Finally, we also compare in Fig. 3 the effective sample size whether a static (thus only coordinate sampling) or a dynamical (both coordinate and momenta) observable is being computed, when sampling at 573 K. The ratio drops more sharply for dynamical calculations, since there is double the number of variables for which IS weights have to be computed (compare eqs. 7 and 8).

When sampling and target PDFs differ more considerably, the effective sample size might not be as large as one needs. It is always possible to perform an additional set of calculations by sampling at the target temperature, thus supplementing the results obtained with the IS strategy. Yet, there is no reason not to make use of previously performed calculations, as IS can always map results between two known distributions, with virtually no cost.
<figure><img src="image_4.png" /><figcaption>Figure 3: Ratio between effective and total sample size, as a function of target temperatures, for three sampling temperatures. The IS weights were computed by considering only coordinates (full lines), or both coordinates and momenta (dashed line).</figcaption></figure>

## 3.2 Temperature dependence of the absorption spectra

The computed excitation energies, vertical (f) and ensemble averaged (⟨f⟩) length-gauge oscillator strengths and characters of the four first singlet states are shown in Table 1. Among the many available theoretical reports,²⁶⁻³¹ we also present in Tab. 1 data from the most recent one,³¹ which employed the EOM-CCSD method. Overall, our results are in good agreement with the previous ones. However, ADC(2) energies are systematically red-shifted (~0.3 eV in average), and the order of the 2A'' and 3A' states is inverted, relative to more accurate calculations.²⁶,³¹

In Fig. 4, the contribution from each adiabatic state to the absorption cross sections is displayed. The first band arises exclusively from excitation to the S₁ state (4π–1π*). The S₂ state (4π–3s/σ*<sub>OH</sub>) contributes with a low-intensity, very broad and asymmetrical structure. Our computed spectra reveals that both S₃ and S₄ states account for the second band, with a somewhat larger component of the later. At the equilibrium geometry, S₃ and S₄ states have well-defined characters, 4π–3p/σ*<sub>OH</sub> and 4π–2π*, respectively, with the later presenting a much larger oscillator strength. We checked the dominant natural transition orbitals for a couple of distorted geometries, and found very mixed 3p/σ*<sub>OH</sub>/2π* characters,
<table><caption>Table 1: Vertical excitation energies (ΔE, in units of eV), vertical (f) and ensemble averaged (⟨f⟩) length-gauge oscillator strengths, for the four lowest-lying singlet excited states of phenol. The dominant configuration for these excitations involve a transition from the highest-lying occupied orbital (4π).</caption><thead><tr><th></th><th>ΔE<sup>a</sup></th><th>f<sup>a</sup></th><th>⟨f⟩<sup>a</sup></th><th>ΔE<sup>b</sup></th><th>f<sup>b</sup></th></tr></thead><tbody><tr><td>2A' (1π<sup>*</sup>)</td><td>4.85</td><td>0.038</td><td>0.029</td><td>4.95</td><td>0.021</td></tr><tr><td>1A'' (3s/σ<sub>OH</sub><sup>*</sup>)</td><td>5.33</td><td>0.000</td><td>0.005</td><td>5.70</td><td>0.000</td></tr><tr><td>2A'' (3p/σ<sub>OH</sub><sup>*</sup>)</td><td>5.80</td><td>0.004</td><td>0.036</td><td>6.18</td><td>0.004</td></tr><tr><td>3A' (2π<sup>*</sup>)</td><td>5.92</td><td>0.059</td><td>0.047</td><td>6.16</td><td>0.034</td></tr></tbody></table>

<sup>a</sup> ADC(2)/(d)-aug-cc-pVDZ level (this work).

<sup>b</sup> EOM-CCSD/aug-cc-pVTZ+R level (from Ref. 31).

<figure><img src="image_5.png" /><figcaption>Figure 4: Photoabsorption cross sections of phenol computed at 296 K, highlighting the contribution from each adiabatic state.</figcaption></figure>

which indicates pronounced intensity borrowing between these two states. However, this aspect alone cannot explain the intensity of the band, since ⟨f⟩ is around 30% larger than f, when summing both S₃ and S₄ contributions. This difference must come from the non-equilibrium configurations, which induce an overall increase in the oscillator strengths. Our interpretation of the spectra is consistent with previous reports.30,31 The second band should arise mostly from excitations counting with a 2π* component, but we call attention to the fact that due to the strong mixing observed for the second band, excitations with a 3p/σ<sub>OH</sub><sup>*</sup> character could also play a role. Higher lying excited also have a non-negligible contribution.
Theoretical and experimental¹⁷ absorption cross sections are compared in Fig. 5, for 296 K, 573 K and 773 K. The experimental data was convoluted with a 0.05 eV line width Gaussian profile. The individual comparisons are very favorable. Positions, shapes and relative intensities of the two bands are reproduced at a semi-quantitative level. The agreement on the absolute cross sections is also quite satisfactory, considering that it is much harder to compute accurate oscillator strengths than excitation energies. As shown in Tab. 1, ADC(2) results seems to underestimate the excitation energies and overestimate the vertical oscillator strengths, when compared to the more sophisticated EOM-CCSD results. These factors should account in great part for the ~0.25 eV difference on the position of the second band and for the difference between theoretical and experimental cross section magnitudes.

<figure><img src="image_6.png" /><figcaption>Figure 5: Photoabsorption cross sections of phenol at different temperatures, according to our calculations and to available experimental results.¹⁷</figcaption></figure>

Most importantly for our main discussion, we are able to capture the main trends of the temperature dependent spectra. An increase in temperature displaces the first band to lower energies, while its maximum intensity is slightly reduced. This, in turn, induces an increase of the cross section in the 4–4.4 eV region, and a displacement of the minimum between the two bands to lower energies. The stronger band gets broader while its maximum intensity does not change much, thus following a larger cross section in the 5.1–5.65 eV and
the 5.9–6.1 eV ranges. All these features have been observed in the measured spectra, with the only exception on the intensity of the second band, which do increase with temperature in the experiment. However, this is a reasonably modest disagreement, as the measured magnitudes increase by around 15% (from 296 K to 773 K) at 5.9 eV, while we found a small decrease of around 5%.

With a 0.05 eV line width, the cross section curves still present structures that do not allow for the identification of a single maximum. In order to better quantify the energy shifts in the spectra, both theoretical and experimental cross sections were convoluted with a larger line width of 0.1 eV. Back to the first band, we found a 0.11 eV redshift from the computed vertical excitation (4.85 eV) to the band maximum (4.74 eV), for the zero temperature cross sections (not shown). Redshifts of the same order have been reported for other compounds,¹¹,⁵⁸ and this behavior seems to be very general.⁵⁹ The difference between calculated band maximum and vertical excitation can be used to estimate an “experimental vertical excitation”, by displacing the experimental band maximum by the same amount.⁷,¹¹ The band maximum is further red-shifted by 0.02 eV when moving from 0 K to room temperature (296 K), and by 0.09 eV when going up to 773 K, whose difference (0.07 eV) compares very favorably to the experimentally convoluted 0.06 eV shift. Despite being a small effect, this finite temperature correction should also be included when comparing vertical and band maximum energies, since measurements are usually performed at room temperature. Considering the convoluted experimental band center of 4.62 eV, and the 0.13 eV theoretical room temperature redshift, the “experimental vertical excitation” would lie at 4.75 eV, in very good agreement with both ADC(2) (4.85 eV) and EOM-CCSD (4.95 eV) results. Furthermore, the theoretical ratio for the intensities of the first band (1:0.89:0.82) is in quite good agreement with the convoluted experimental result (1:0.93:0.86). We also compared vertical excitations and band maximum for the S₃ and S₄ states, and found zero-temperature redshifts of 0.14 eV and 0.12 eV, and a small room temperature correction of 0.01 eV for the former. The theoretical redshift of 0.13 eV (weighting both states) thus displaces the measured maximum (6.00
eV) to an “experimental vertical excitation” of 6.13 eV. Again, agreement is very good with the vertical excitation energy (weighted by the oscillator-strengths) computed with ADC(2) (5.91 eV) and EOM-CCSD (6.16 eV). The case of S₂ is atypical, given its skewed shape towards higher energies. While vertical and band maximum position are virtually the same, there is a 0.19 eV redshift between the vertical excitation energy (5.33 eV) and the median of the cross section curve (5.14 eV). This value is equivalent to the average excitation energy of the ensemble (weighted by the oscillator-strengths), and therefore it more closely represents the vibrationally averaged excitation energy than the position of the band maximum.⁵⁹ It is interesting to notice that, despite having the smallest redshift (in the conventional sense), the S₂ state is the one that stabilizes the most when considering an ensemble of nuclear configurations.

Our implementation of the IS algorithm also allows us to select individual temperatures for each normal mode. This procedure can point out which vibrations promote the observed changes in the absorption spectra. In Fig. 6, we demonstrate exactly that, by showing the spectra computed at 296 K, at 773 K and for two combinations of hot and cold vibrations. By heating only the *a''* modes (10 out of 33), the difference between the 773 K and 296 K cross sections drops to about a third, and by further including the 4 first *a'* modes (14 in total), only a sixth of the original difference remains. This demonstrates that the band shifts and broadenings arise mostly when the molecular vibrations induce out-of-plane distortions.

### 3.3 Tully's single avoided crossing model

The fraction of trajectories lying on the ground state are depicted in Fig. 7, as computed for dynamics simulation where momenta sampling was performed at 4,000 K, 10,000 K and 15,000 K. The overall Gaussian decay of the population is associated with the stochastic hopping to the excite state, while the small bump around 40 fs can be traced back to the chance of hopping back to the ground state. Still, the ground state is repopulated in a much
<figure><img src="image_7.png" /></figure>

Figure 6: Photoabsorption cross sections of phenol, computed at 296 K, 773 K, and at two intermediate cases: when $a''$ modes are at 773 K and $a'$ modes are at 296 K, and when the 14 lower frequency modes (which include all the $a''$) are at 773 K and the remaining 19 are at 296 K.

<figure><img src="image_8.png" /></figure>

Figure 7: Fraction of trajectories on the ground state, for coincident sampling and target temperatures of 4,000 K (blue), 10,000 K (orange) and 15,000 K (red).

longer timescale, which is associated to trajectories that remain momentarily trapped in the excited state. As the temperature increases, there are more trajectories with larger initial velocities (and fewer with smaller), such that the avoided crossing region is reached faster. Also, there is an increase in the hopping probabilities and therefore a decrease on the final population of the ground state.
<figure><img src="image_9.png" /><figcaption>Figure 8: Fraction of trajectories on the ground state, for sampling temperatures of 4,000 K (blue), 10,000 K (orange) and 15,000 K (red), and for target temperatures of 10,000 K (left) and 4,000 K (right).</figcaption></figure>

Concerning our main discussion, the most relevant result from the dynamics simulations are shown in Fig. 8, which is a zoomed in version of Fig. 7 with the addition of the IS results. There is an excellent level of agreement between the curves obtained with IS sampling and with the standard conventional sampling. In the present case, $r = 0.84$ for the sampling and target temperatures of 15,000 K and 10,000 K, while $r = 0.31$ for the 15,000/4,000 K pair and $r = 0.51$ for the 10,000/4,000 K combination. The larger effective number of trajectories of the later case explains its better convergence (see Fig. 8), which is a consequence of the closer sampling and target temperatures. We computed the NRMSD associated with the curves in Fig. 8 (not shown) and we found that it decreases consistently when the effective number of trajectory grows, which demonstrates that the proposed strategy is also valid for the calculation of dynamical properties.

## 3.4 Other potential applications

We have presented a detailed discussion on how IS can be applied to investigate the influence of temperature on photoabsorption cross sections. Naturally, the strategy can be extended to study other static process, such as fluorescence, photoionization and phosphorescence, which can present marked temperature dependent effects.⁶⁰⁻⁶² Additionally, the method can
be used to inquire the role of temperature for any dynamical observable framed within the
nuclear ensemble approach, such as photo-induced dynamics and all its decaying mechanisms
(internal conversion, dissociation, intersystem crossing). Since we can also control the oc-
cupation numbers of each mode, investigations of vibrationally mediated processes<sup>63</sup> should
also profit from the proposed IS strategy. Furthermore, we recall that bimolecular colli-
sions are very often described with trajectory quasiclassical simulations, which also require
sampling initial conditions from a given PDF.⁵ There are, of course, important differences
on how they are generated for isolated molecules and for a collision problem,<sup>5,64,65</sup> but in
principle the presented method should likewise be of interest for the later class of problems.
Besides all these direct potential applications, we elaborate on the sequence about other
situations in which the technique might be valuable.

A very interesting possibility concerns the level of description of the molecular vibrations,
and how to more easily assess its impact on nuclear ensemble calculations. The harmonic
approximation is usually employed, for which the PDFs are Gaussian-shaped around the
equilibrium geometry. By including higher-order terms or using more intricate potentials,
the PDFs also acquire a more complicated profile.⁴²,⁴³ That would lead to more accurate
probabilities for sampling at the most distorted geometries, which are under or oversampled
in the harmonic approximation, depending on how the potentials deviate from harmonicity.
This aspect should be specially important for larger molecules, where many low-frequency
vibrations are present. If one is interested in evaluating how two levels of description for the
vibrations affect the computation of an observable, the conventional way to proceed would
be to perform two separate rounds of calculations, by sampling from the PDFs associated
to each set of potentials. Instead, the sampling could be done only once (say, for the har-
monic potentials), and then the IS strategy could be applied to compute the results for any
other potentials. That could allow for systematic and inexpensive studies about the role of
anharmonicities in the selection of initial condition, ultimately leading to a more accurate
description of the process being investigated.
A related and simpler possibility is to evaluate the impact of a different set of vibrational frequencies. That would be the case if one introduces a scaling factor to each frequency in an attempt to recover some anharmonic effect. Then, target $\omega_i^t$ and sampling $\omega_i^s$ frequencies differ, and by also considering variations in temperature, the IS weights for a static property would be given by:

$$w(\mathbf{q}) = \prod_{i=1}^{M} \sqrt{\frac{\alpha_i^t \omega_i^t}{\alpha_i^s \omega_i^s}} \exp \left[ - \left( \alpha_i^t \omega_i^t - \alpha_i^s \omega_i^s \right) \frac{\mu_i q_i^2}{\hbar} \right], \quad (12)$$

where in this case $\alpha_i^t = \tanh(\hbar\omega_i^t/2k_B T^t)$. We did some experimentation for phenol, and found very marginal changes in the spectra (not shown). In any case, the low frequency modes were the most sensitive to the frequency scaling.

The effect of different isotopologues on the selection of initial conditions could also be investigated with the IS technique. A different isotope cause a small shift on the normal modes frequencies, as well as an overall scaling of the PDFs, due to the $\sqrt{\mu_i}$ factor. One could readily investigate this isotope effect on the absorption spectra, for example, which might be of relevance for the isotope fractionation of atmospheric compounds, 25,66,67 and for photodesorption rates of molecular ices. 19

Finally, we also considered whether IS could be used to map results obtained from a quantum based PDF to a classical one. For one single mode, the portion of the zero-temperature quantum PDF that lies inside the classically allowed region is given by $\int_{-q_{cl}}^{+q_{cl}} P(q)dq$, where $P(q)$ is the quantum (Wigner) distribution and $q_{cl}$ are the classical turning points. In the harmonic approximation this corresponds to $\text{erf}(1) \sim 0.84$, where $\text{erf}$ is the error function. Hence, there is a reasonable overlap between the two distributions. When we consider $M$ normal modes, the probability that a point sampled with a quantum PDF also lies within the classical volume is given by $\text{erf}(1)^M$, thus decreasing exponentially with the number of degrees of freedom. Therefore, there is an increasingly large chance that at least for one mode the sampled point will be in a classically forbidden region. Only one component is all
that is needed to turn the phase space point inaccessible by a classical PDF. This is quite a
revealing result. It implies that quantum and classical distribution effectively sample from
different regions of the phase space. In the case of phenol, for example, $\operatorname{erf}(1)^{33} \sim 4 \times 10^{-4}$,
and if we were to perform dynamics, than the classical-quantum overlap would be an im-
pressively small $\operatorname{erf}(1)^{2 \times 33} \sim 10^{-5}$. Yet, there is an extensive number of papers$^{44,68-72}$ where
good results can be achieved even when the nuclear configuration/initial conditions phase
space is sampled classically. This indicates that very often the computed observables should
be considerably homogeneous throughout the classical, inner part of the phase space and the
quantum exclusive, outer portion of it.

# 4 Conclusions

We proposed the adoption of the importance sampling technique into the nuclear ensemble
approach, which allows for a mapping between results computed for one sampling probability
distribution function to any target distribution of interest. The only requisite being the
relative probabilities for sampling from the two distributions, which can be readily computed
when a quantum based distribution is used. The theoretical framework is very general, and it
can be employed for the calculation of any static or dynamical observable. It should also be
useful for assessing models for the nuclei motion that go beyond the harmonic approximation.
Incorporation of isotope effects should also be possible with the proposed strategy. We
also discussed how quantum and classical distributions for the initial conditions can have a
surprisingly small overlap.

We studied how varying temperatures affect the photoabsorption cross sections of phenol,
and our calculations were able to capture the main features of the measured spectra. By using
the IS technique the cross sections did converge towards those obtained with conventional
sampling. Convergence for each individual temperature is slower, as expected, but when
considering the whole range of temperatures the best practice is to adopt IS. For the same
level of precision, applying IS leads to considerable savings in total computational time in comparison to performing a series of separate samplings. We further showed that a better overall rate of convergence is attained when the sampling temperature lies at the middle of the range of temperatures of interest. The IS technique was further put under test for trajectory surface hopping dynamics, which was used to investigate the role of temperature in Tully's single avoided crossing model. The same set of dynamical observables were obtained, whether they were computed with conventional sampling or with importance sampling.

The IS algorithm was coded as a modular part and integrated into the Newton-X package.⁵⁵,⁵⁶ Both photoabsorption spectra and dynamics simulations can make use of it, as the user is now able to set up target temperatures, frequencies and occupation numbers for each normal mode. We believe our proposed strategy might serve as a valuable tool for the study of several physico-chemical processes, as it effectively widens the applicability of the nuclear ensemble approach.

## Acknowledgement

F. K. acknowledges financial support from São Paulo Research Foundation (FAPESP), under grant No. 2016/21880-7. M. B. work was supported by Excellence Initiative of Aix-Marseille University (A*MIDEX) and the project Equip@Meso (ANR-10-EQPX-29-01), both funded by the French Government "Investissements d'Avenir" program. MB also acknowledges funding from the WSPLIT project (ANR-17-CE05-0005-01).

## References

(1) Bunker, D. L. *J. Chem. Phys.* **1962**, *37*, 393-403.

(2) Karplus, M.; Porter, R. N.; Sharma, R. D. *J. Chem. Phys.* **1964**, *40*, 2033-2034.

(3) Schatz, G. C.; Bowman, J. M.; Kuppermann, A. *J. Chem. Phys.* **1975**, *63*, 685-696.
(4) Bernstein, R. *Atom-molecule collision theory: a guide for the experimentalist*; Physics of atoms and molecules; Plenum Press, 1979.

(5) Paranjothy, M.; Sun, R.; Zhuang, Y.; Hase, W. L. *WIREs: Comput. Mol. Sci.* **2013**, *3*, 296–316.

(6) Pratihar, S.; Ma, X.; Homayoon, Z.; Barnes, G. L.; Hase, W. L. *J. Am. Chem. Soc.* **2017**, *139*, 3570–3590.

(7) Crespo-Otero, R.; Barbatti, M. *Theor. Chem. Acc.* **2012**, *131*, 1237.

(8) Svoboda, O.; Ončák, M.; Slavíček, P. *J. Chem. Phys.* **2011**, *135*, 154301.

(9) Barbatti, M. *Phys. Chem. Chem. Phys.* **2011**, *13*, 4686.

(10) Li, Q.; Mennucci, B.; Robb, M. A.; Blancafort, L.; Curutchet, C. *J. Chem. Theory Comput.* **2015**, *11*, 1674–1682.

(11) Barbatti, M.; Sen, K. *Int. J. Quantum Chem.* **2016**, *116*, 762–771.

(12) Persico, M.; Granucci, G. *Theor. Chem. Acc.* **2014**, *133*, 1526.

(13) McGillen, M. R.; Fleming, E. L.; Jackman, C. H.; Burkholder, J. B. *Geophys. Res. Lett.* **2013**, *40*, 4772–4776.

(14) Bernard, F.; McGillen, M. R.; Fleming, E. L.; Jackman, C. H.; Burkholder, J. B. *J. Photochem. Photobiol.*, **A** **2015**, *306*, 13–20.

(15) Papadimitriou, V. C.; McGillen, M. R.; Fleming, E. L.; Jackman, C. H.; Burkholder, J. B. *Geophys. Res. Lett.* **2013**, *40*, 440–445.

(16) Blitz, M. A.; Heard, D. E.; Pilling, M. J. *J. Phys. Chem. A* **2006**, *110*, 6742–6756.

(17) Grosch, H.; Sárossy, Z.; Egsgaard, H.; Fateev, A. *J. Quant. Spectrosc. Radiat. Transfer* **2015**, *156*, 17–23.
(18) Dufour, A.; Girods, P.; Masson, E.; Normand, S.; Rogaume, Y.; Zoulalian, A. *J. Chromatogr. A* **2007**, *1164*, 240–247.

(19) Cruz-Diaz, G. A.; Caro, G. M. M.; Chen, Y.-J. *Mon. Notices Royal Astron. Soc.* **2014**, *439*, 2370–2376.

(20) Blasberger, A.; Behar, E.; Perets, H. B.; Brosch, N.; Tielens, A. G. G. M. *Astrophys. J.* **2016**, *836*, 1–8.

(21) Öberg, K. I. *Chem. Rev.* **2016**, *116*, 9631–9663.

(22) Grosch, H.; Fateev, A.; Clausen, S. *J. Quant. Spectrosc. Radiat. Transfer* **2015**, *154*, 28–34.

(23) Chen, F.; Wu, C. *J. Quant. Spectrosc. Radiat. Transfer* **2004**, *85*, 195–209.

(24) Schmidt, J. A.; Johnson, M. S.; McBane, G. C.; Schinke, R. *J. Chem. Phys.* **2012**, *137*, 054313.

(25) Schmidt, J. A.; Johnson, M. S.; Schinke, R. *Proc. Natl. Acad. Sci. U.S.A.* **2013**, *110*, 17691–17696.

(26) Lorentzon, J.; Malmqvist, P. Å.; Fülscher, M.; Roos, B. O. *Theor. Chim. Acta* **1995**, *91*, 91–108.

(27) Daigoku, K.; ichi Ishiuchi, S.; Sakai, M.; Fujii, M.; Hashimoto, K. *J. Chem. Phys.* **2003**, *119*, 5149–5158.

(28) Lan, Z.; Domcke, W.; Vallet, V.; Sobolewski, A. L.; Mahapatra, S. *J. Chem. Phys.* **2005**, *122*, 224315.

(29) Pino, G. A.; Oldani, A. N.; Marceca, E.; Fujii, M.; Ishiuchi, S.-I.; Miyazaki, M.; Broquier, M.; Dedonder, C.; Jouvet, C. *J. Chem. Phys.* **2010**, *133*, 124313.
(30) Jones, D. B.; da Silva, G. B.; Neves, R. F. C.; Duque, H. V.; Chiari, L.; de Oliveira, E. M.; Lopes, M. C. A.; da Costa, R. F.; Varella, M. T. d. N.; Bettega, M. H. F.; Lima, M. A. P.; Brunger, M. J. *J. Chem. Phys.* **2014**, *141*, 074314.

(31) Limão-Vieira, P.; Duflot, D.; Ferreira da Silva, F.; Lange, E.; Jones, N. C.; Hoffmann, S. V.; Śmiałek, M. A.; Jones, D. B.; Brunger, M. J. *J. Chem. Phys.* **2016**, *145*, 034302.

(32) Kimura, K.; Nagakura, S. *Mol. Phys.* **1965**, *9*, 117-135.

(33) Nix, M. G. D.; Devine, A. L.; Cronin, B.; Dixon, R. N.; Ashfold, M. N. R. *J. Chem. Phys.* **2006**, *125*, 133318.

(34) Miura, M.; Aoki, Y.; Champagne, B. *J. Chem. Phys.* **2007**, *127*, 084103.

(35) Ari, T.; Güven, H.; Ecevit, N. *J. Electron. Spectrosc. Relat. Phenom.* **1995**, *73*, 13-23.

(36) Tully, J. C. *J. Chem. Phys.* **1990**, *93*, 1061.

(37) Wigner, E. *Phys. Rev.* **1932**, *40*, 749-759.

(38) Hillery, M.; O'Connell, R. F.; Scully, M. O.; Wigner, E. P. *Phys. Rep.* **1984**, *106*, 121-167.

(39) Sun, L.; Hase, W. L. *J. Chem. Phys.* **2010**, *133*, 044313.

(40) Christoffel, K. M.; Bowman, J. M. *J. Chem. Phys.* **1996**, *104*, 8348-8356.

(41) Kube, S.; Lasser, C.; Weber, M. *J. Comput. Phys.* **2009**, *228*, 1947-1962.

(42) García-Vela, A. *J. Chem. Phys.* **2000**, *112*, 8302-8311.

(43) García-Vela, A. *Chem. Phys.* **2002**, *285*, 245-253.

(44) Cusati, T.; Granucci, G.; Persico, M. *J. Am. Chem. Soc.* **2011**, *133*, 5109-5123.
(45) Pattengill, M. *Chem. Phys.* **1983**, *75*, 59-66.

(46) Pattengill, M. *Chem. Phys. Lett.* **1984**, *105*, 651-654.

(47) Han, Y.-C.; Bowman, J. M. *Chem. Phys. Lett.* **2013**, *556*, 39-43.

(48) Tokdar, S. T.; Kass, R. E. *WIREs: Comput. Stat.* **2010**, *2*, 54-60.

(49) Faist, M. B.; Muckerman, J. T.; Schubert, F. E. *J. Chem. Phys.* **1978**, *69*, 4087-4096.

(50) Faist, J. T. M.; B., M. *J. Phys. Chemistry* **1979**, *83*.

(51) Trofimov, A. B.; Schirmer, J. *J. Phys. B: At., Mol. Opt. Phys.* **1995**, *28*, 2299.

(52) Schirmer, J. *Phys. Rev. A* **1982**, *26*, 2395-2416.

(53) Hättig, C.; Köhn, A. *J. Chem. Phys.* **2002**, *117*, 6939-6951.

(54) Ahlrichs, R.; Bär, M.; Häser, M.; Horn, H.; Kölmel, C. *Chem. Phys. Lett.* **1989**, *162*, 165-169.

(55) Barbatti, M.; Ruckenbauer, M.; Plasser, F.; Pittner, J.; Granucci, G.; Persico, M.; Lischka, H. *WIREs: Comput. Mol. Sci.* **2014**, *4*, 26-33.

(56) Barbatti, M.; Granucci, G.; Ruckenbauer, M.; Plasser, F.; Crespo-Otero, R.; Pittner, J.; Persico, M.; Lischka, H., NEWTON - X: A package for Newtonian dynamics close to the crossing seam. Version 2. 2016; http://www.newtonx.org.

(57) Granucci, G.; Persico, M. *J. Chem. Phys.* **2007**, *126*.

(58) Patrick, C. E.; Giustino, F. *Nat. Commun.* **2013**, *4*, 2006.

(59) Ferrer, F. J.; Cerezo, J.; Stendardo, E.; Improta, R.; Santoro, F. *J. Chem. Theory Comput.* **2013**, *9*, 2072-2082.

(60) Sajoto, T.; Djurovich, P. I.; Tamayo, A. B.; Oxgaard, J.; Goddard, W. A.; Thompson, M. E. *J. Am. Chem. Soc.* **2009**, *131*, 9813-9822, PMID: 19537720.
(61) Wang, C.-R.; Nguyen, J.; Lu, Q.-B. *J. Am. Chem. Soc.* **2009**, *131*, 11320–2.

(62) Lambrecht, D. S.; Clark, G. N. I.; Head-Gordon, T.; Head-Gordon, M. *J. Phys. Chem. A* **2011**, *115*, 5928–5935.

(63) Rosenwaks, S. *Vibrationally Mediated Photodissociation*; The Royal Society of Chemistry, 2009; pp P001–P204.

(64) Bowman, J. M.; Kuppermann, A.; Schatz, G. C. *Chem. Phys. Lett.* **1973**, *19*, 21–25.

(65) Nagy, T.; Vikár, A.; Lendvay, G. *J. Chem. Phys.* **2016**, *144*.

(66) Danielache, S. O.; Nanbu, S.; Eskebjerg, C.; Johnson, M. S.; Yoshida, N. *J. Chem. Phys.* **2009**, *131*.

(67) Hattori, S.; Danielache, S. O.; Johnson, M. S.; Schmidt, J. A.; Kjaergaard, H. G.; Toyoda, S.; Ueno, Y.; Yoshida, N. *Atmospheric Chem. Phys.* **2011**, *11*, 10293–10303.

(68) Cattaneo, P.; Persico, M. *J. Am. Chem. Soc.* **2001**, *123*, 7638–7645.

(69) Wohlgemuth, M.; Bonai-Kouteck, V.; Mitri, R. *J. Chem. Phys.* **2011**, *135*.

(70) Long, R.; Prezhdo, O. V. *J. Am. Chem. Soc.* **2011**, *133*, 19240–19249.

(71) Spörkel, L.; Cui, G.; Thiel, W. *J. Phys. Chem. A* **2013**, *117*, 4574–4583.

(72) Mai, S.; Gattuso, H.; Fumanal, M.; Muñoz-Losa, A.; Monari, A.; Daniel, C.; González, L. *Phys. Chem. Chem. Phys.* **2017**, *3*, 21–25.
Graphical TOC Entry

<figure><img src="image_10.png" /></figure>