Nuclear Ensemble Approach with Importance Sampling

We show that the importance sampling technique can eﬀectively 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 proﬁts 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.


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 appropriate. 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 trajectory method". For spectra calculations, there exists a construction that leads to the nuclear ensemble method. 7 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][8][9][10][11] Adiabatic and non-adiabatic dynamics that employ a classical nuclei approximation can also be framed as a nuclear ensemble approach. In this case, however, a set of initial coordinates and momenta have to 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. 12 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. [13][14][15] 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, [13][14][15] 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. 15,16 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. 17 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 3 devised instead. 18 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. [19][20][21] Despite several experimental results that have one of the above mentioned points as a motivation, [13][14][15]17,19,22,23 there are very few results from the side of theory. 24,25 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 17 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, 26-31 absorption spectroscopy, 17,31-34 and electron energy loss spectroscopy. 30,35 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 36 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.

Nuclear ensemble approach
The central concept behind the nuclear ensemble approach is that a physical observable µ can be computed by integrating its value over a set of representative nuclear configurations. 7 4 This statement can be cast into the following expression: where f is the function to be integrated in phase space x = (q, p), q and p representing the set of nuclear coordinates and momenta, respectively. These, in turn, follow a normalized PDF P (x), while E P (·) 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., (X ∼ P ).
In practice, one actually computes an estimateμ by performing an average over a set of N p sampled points:μ For the calculation of static observables, where f depends only on the set of initial coordinates 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: where it is implicit that the second P represents a function of q only.
The PDF P (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. [37][38][39] There are also other possibilities to represent the quantum distribution, based on the Husimi distribution, for example. 40,41 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. 42,43 Once the distribution is constructed, sampling of each coordinate and momenta pair can be performed in a correlated on uncorrelated way. 39 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: 37 where the product runs over M vibrational modes, q i is the coordinate, p i is the momentum, µ i is the reduced mass and ω i is the harmonic frequency of the i-th mode. The temperature T appears in the factor α i = tanh(hω 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 6 with larger compounds or simulations in condensed phase. 44 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 8,11,39,[45][46][47] show that overall better agreement with experiment is achieved when a quantum-based distribution is used, instead of the classically-based thermal sampling alternative.

Importance sampling
Suppose we have already computed the observableμ for a given PDF Q(x) and now we would like to evaluate it for a distinct PDF P (x). This would be the case if we are interested in evaluating howμ 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 (x), and evaluate the function f (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. 48 Basically, after the function f is evaluated for a set of points X l obtained from a given sampling distribution Q(x), the expected valueμ can be readily computed for a distinct target distribution P (x) by introducing appropriate weights for each sampled point. The following expression synthesizes the idea: The sampling distribution Q(x) is introduced in both numerator and denominator. In the numerator it plays the role of the actual PDF from which X has already been sampled, while in the denominator it is combined with P (x) and gives rise to the IS weighting factors . Therefore, the sampling and the function evaluations are performed only once, according to the sampling distribution Q(x), and for each target distribution P (x) of interest one simply needs to apply the weights w(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 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. 48 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 (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 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: where T s is the sampling temperature, T t is the target temperature, 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:

Photoabsorption spectra calculation
Within the nuclear ensemble approach, the photoabsorption cross section is computed according to: 7 where ∆E 0n is the excitation energy and f 0n is the oscillator strength of the n-th state, N 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 δ =0.05 eV.
Excited states were described with the algebraic diagrammatic construction to the second order (ADC (2) 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. 55,56 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.

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 surfacehopping method, 36 which we have employed for our dynamics case study. We have investigated the first of Tully's classical models, 36 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. 36 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).
where σ 0 max is the maximum value of the reference cross section and the sum runs over the grid of N 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 con- 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: This quantity is therefore related to the distribution of IS weights w(X l ). The closer sampling and target PDFs are, the sharper will be the distribution of w(X l ) around 1, such that N eff p 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 eff p 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 eff p /N p . The rate of convergence becomes the same, regardless of the combination of sampling and target PDFs, which supports the interpretation that N 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 eff p and w(X l ) for the case of phenol, we evaluated how the ratio r = N eff p /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

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, [26][27][28][29][30][31] we also present in Tab. 1 data from the most recent one, 31 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. 26,31 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 1 state (4π−1π * ).
The S 2 state (4π−3s/σ * OH ) contributes with a low-intensity, very broad and asymmetrical structure. Our computed spectra reveals that both S 3 and S 4 states account for the second band, with a somewhat larger component of the later. At the equilibrium geometry, S 3 and S 4 states have well-defined characters, 4π−3p/σ * OH 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/σ * OH /2π * characters, .  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 3 and S 4 contributions. This difference must come from the nonequilibrium 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/σ * OH character could also play a role. Higher lying excited also have a non-negligible contribution.  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 18 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 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, 11,58 and this behavior seems to be very general. 59 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. 7,11 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 We also compared vertical excitations and band maximum for the S 3 and S 4 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 2 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. 59 It is interesting to notice that, despite having the smallest redshift (in the conventional sense), the S 2 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.

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 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. 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.

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. [60][61][62] 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 occupation numbers of each mode, investigations of vibrationally mediated processes 63 should also profit from the proposed IS strategy. Furthermore, we recall that bimolecular collisions are very often described with trajectory quasiclassical simulations, which also require sampling initial conditions from a given PDF. 5 There are, of course, important differences on how they are generated for isolated molecules and for a collision problem, 5,64,65 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. 42,43 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 harmonic 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 ω t i and sampling ω s i frequencies differ, and by also considering variations in temperature, the IS weights for a static property would be given by: where in this case α t i = tanh(hω t i /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 √ µ 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 +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 erf(1) ∼ 0.84, where 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 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, erf(1) 33 ∼ 4 × 10 −4 , and if we were to perform dynamics, than the classical-quantum overlap would be an impressively small erf(1) 2×33 ∼ 10 −5 . Yet, there is an extensive number of papers 44,[68][69][70][71][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.

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. 55,56 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.