2.3.4. Solvent Effects

Finally, a brief discussion of the environmental effects is appropriate, as ESIPT systems are usually investigated in solution or in other complex condensed-phase environments. In particular, the polarity of the surrounding medium has been observed to have a strong impact on the ESIPT reaction efficiency [73,76].

Thus far, several different approaches have been employed to tackle environmental effects on ESIPT, including microsolvation [43,69], the conductor-like screening model (COSMO) [48,94,95], the polarizable continuum model (PCM) [76,96,97], the solvent model density (SMD) method [46,98], and the integral equation formalism version of PCM (IEF-PCM) [99–101]. The latter approach, particularly popular recently [33,37,73,102], has been applied e.g. by Wang et al. to the BTS system [39] in methylene chloride, yielding very high accuracy predictions for excitation and emission wavelengths, with divergence from the experimental values measured in just a few nm. In addition to the general purpose methods listed above, state-specific PCM treatments of correlated linear response (cLR) [103] and the vertical excitation model within the unrelaxed density approximation (VEM-UD) [104,105] have been successfully applied to study ESIPT by Vérité et al. [40], who pointed out the advantages that these approaches bring for the description of charge-transfer states in ESIPT reactions. Nevertheless, the explicit inclusion of (typically few) solvent molecules is necessary in certain cases, especially for protic solvents and solvents exhibiting proton-accepting properties, since resulting competition between intra- and intermolecular hydrogen bond formation may drastically affect the ESIPT reaction yield [88,90,106].

#### *2.4. Summary of the Static ESIPT Investigation Methods*

To summarize the section dedicated to the static ESIPT investigation protocol, we again underline its strengths as being a relatively affordable and ye<sup>t</sup> informative approach, designed to provide a fundamental characterization of ESIPT, including the system's absorption and emission properties, as well as information on the topography of GS and ES PESs over pre-selected reaction coordinates. Due to its inherent compatibility with a grea<sup>t</sup> variety of electronic structure methods, this protocol allows researchers to take advantage of new developments in electronic structure theory and, thus, constantly provides opportunities for cutting-edge studies of ESIPT in all types of molecular systems.

At the same time, it should be noted that, under certain circumstances, the investigation of static ESIPT paths may not be sufficient. In particular, systems undergoing multiple PT reactions are typically challenging to be accurately studied with this protocol due to the large computational cost of multi-dimensional PES scans, on the one hand, and the critical role of the sequence of individual processes missed at this level, on the other hand. Another situation, in which special precautions should be taken, is when ESIPT occurs within a dense manifold of electronic states, such as in situations in which a competition between various photochemical transformations is to be expected; in these cases, one may need to explicitly determine the relative efficiency of each channel, which usually requires the inclusion of nuclear-dynamic effects.

#### **3. Nonadiabatic Molecular Dynamic Approaches**

New opportunities of delving deeper into the course of the ESIPT reaction open when one turns toward dynamic approaches. In the most general view, this refers to a large (and growing) number of methods allowing for real-time simulations of molecular systems' evolution in terms of their electronic and nuclear structure, beyond the static picture.

The time-dependent Schrödinger equation is the core and starting point for the dynamic methods, yielding two broad families of approaches, differed by the level of the applied approximations. The first family consists of fully quantum dynamic (QD) methods, in which electronic and nuclear degrees of freedom are treated both at the quantummechanical level. The second family is built on nonadiabatic mixed quantum–classical (NA-MQC) dynamics methods, in which the nuclei, which are much slower than the electrons, are treated at the classical or semi-classical level, propagating under the Newton equation of motion. The quantum–electronic and classical–nuclear subsystems are in this picture connected by the nonadiabatic coupling, ensuring the self-consistency of the description of the total molecular system. We start our discussion of ESIPT nonadiabatic molecular dynamics studies with the NA-MQC methods, which will be subsequently followed by the analysis of the QD performance, in line with the increasing level of the method exactness.

#### *3.1. Mixed Quantum–Classical Dynamic Calculations in ESIPT Studies*

The NA-MQC methods, recently summarized in an excellent review by Crespo-Otero and Barbatti [107], can be divided into several groups, out of which the trajectory surface hopping (TSH) [108–111], ab initio multiple spawning (AIMS) [112,113], and, most recently, the nuclear–electronic orbital Ehrenfest (NEO-Ehrenfest) [114] methods are, to the best of our knowledge, the ones that have been successfully employed in dynamic ESIPT investigations thus far. Below, a brief characterization of these approaches is provided, along with an illustration of their performance for the description of the ESIPT process. Additionally, a schematic representation of their underlying mechanisms is presented in Figure 3.

**Figure 3.** Schematic illustration of the mechanisms of the discussed NA-MQC dynamic approaches: trajectory surface hopping (TSH), ab inito multiple spawning (AIMS), and nuclear–electronic orbital Ehrenfest (NEO-Ehrenfest).

#### 3.1.1. Trajectory Surface Hopping Approach

The trajectory surface hopping method, especially under Tully's fewest-switches (FSSH) [108] algorithm and under other algorithms based on the Landau–Zener (LZ) model [115,116], is the most widely used NA-MQC approach in ESIPT studies thus far. TSH relies on the modeling of the real-time evolution of the molecular system by a set of independent classical trajectories, which together are assumed to represent a nuclear wave packet in an approximate (statistical) way. While the trajectories are propagated on individual Born–Oppenheimer adiabatic PESs, nonadiabatic transitions between these surfaces are possible in regions characterized by large interstate nonadiabatic couplings (NACs). The interstate transitions are controlled by a stochastic algorithm (FSSH) or induced in minimum-energy-gap regions (LZ methods), with the "hopping" probability proportional to the NAC. Importantly, the TSH method can be implemented as an "on-thefly" approach [117], which means that the actual PES, on which the system is propagating, does not have to be known in advance, and electronic properties, such as energies or gradients, are calculated along the NA-MQC path as needed. It should be noted, however, that usually, many trajectories are required for reliable and converged TSH results [107].

As of today, the TSH method has been implemented in a number of dedicated software packages, including Newton X [118,119], Shark [120], Jade [121], etc. [107]. In all cases, the NA-MQC dynamics protocol has to be paired with an electronic structure method, and its choice needs to be made with grea<sup>t</sup> care, as it directly impacts the quality of the results, as well as the simulation cost. One needs to be aware, however, that due to the inherent mixed-classical nature of the TSH approach, certain effects that may play an important role in the ESIPT reaction cannot be reproduced at the TSH level of theory. This includes all phenomena stemming from the nonlocality of the true nuclear wave function, which is reduced to a single point on adiabatic PES within the TSH picture. In particular, proton tunneling, wave-packet interference, and decoherence effects are not included, the latter being partially restored in the TSH simulations via the introduction of various kinds of decoherence corrections [122].

In terms of recent applications of the TSH methodology to particular ESIPT studies, Li et al. reported interesting results explaining 3-hydroxyflavone dual fluorescence in solvents containing protic contamination with a competition between intra- and intermolecular excited-state PT reactions [90], discussing an effect of the number of explicitly

included water molecules on the simulation outcomes. In this case, the applied FSSH/TD-DFT methodology allowed for high-quality predictions of the electronic excitation energies (typical error below 0.2 eV) and also yielded ESIPT timescale in very good agreemen<sup>t</sup> with available experimental data, with a deviation of less than 10 fs. Another challenging aspect of dealing with a large number of possible photo-reaction products has been tackled by Tuna et al., who employed a robust multiconfiguration interaction variant of the orthogonalization-corrected semi-empirical OM2 approach to model ESIPT-driven photochemistry of urocanic acid [123]. The same method has also been applied by Xia et al. to study relaxation mechanisms in the isolated benzodiazepinone molecule, in which several interconnected relaxation channels come to play [124]. Furthermore, we recently performed a TSH study at the TDA-DFT level to analyze the impact of the character of the lowest excited state on the ESIPT process efficiency [84], eventually confirming the important role of the *ππ*\* states.

#### 3.1.2. Ab Initio Multiple Spawning Approach

Another NA-MQC approach that has found applications in time-resolved ESIPT simulations is the ab initio multiple spawning method. AIMS originates from the formally exact full multiple spawning methodology [125,126]. Its core concept relies on representing the nuclear wave function with partially coupled traveling Gaussian functions, having a finite width both in position and momentum coordinates and interacting during the dynamics. Importantly, the total number of the "on-the-fly" propagated Gaussian functions changes in time since on each passage through a PES region characterized with strong NAC, a new Gaussian is spawned (hence, the S in AIMS).

Similar to the TSH case, AIMS simulations require combining the particular MS protocol with a suitable electronic structure method. As for the AIMS code itself, as of today, it is available within several software packages, including GAMESS [127,128], MOL-PRO [129,130], and MOPAC [131,132]. Technically, AIMS involves a higher computational cost than TSH, ye<sup>t</sup> it should be considered a superior approach, inherently including decoherence effects, and yielding a correct description of some non-local phenomena. At the same time, due to certain intrinsic limitations of the AIMS approach, the tunneling effect, although theoretically possible to be covered through the intrastate spawning procedure [112,125], is not reproduced at this level of theory [113].

Turning to recent interesting applications of the AIMS in ESIPT studies, Pijeau et al. investigated the photophysics of the paradigmatic salicylideneaniline (SA) system [133], focusing on the effect of nonplanarity on ESIPT and on the total deactivation mechanism. In this study, the AIMS protocol has been connected with the floating occupation molecular orbital complete active space configuration interaction (FOMO-CASCI) method, with further wave function-in-DFT embedding. The same group also tackled the hydroxyphenyl benzothiazole (HBT) system at this level of theory, obtaining very good agreemen<sup>t</sup> with the experimental results [134].

## 3.1.3. Nuclear–Electronic Orbital Ehrenfest Approach

Recently, a new NA-MQC dynamic approach, NEO-Ehrenfest, aiming toward the further enhanced recovery of nonlocal effects, has been developed. Within this method, protons are treated quantum mechanically on an equal footing with the electrons, yielding automatic inclusion of the ZPE, quantized vibrational levels, and tunneling effects associated with these species [114]. The NEO-Ehrenfest approach, specifically tailored to provide a high-level description of the ESIPT and PCET processes [135], is built on the concept of semi-classical traveling proton basis functions, which, on the one hand, provide means for the quantum-mechanical representation of protons, as has been demonstrated before for the time-independent case [136], and, on the other hand, enable the description of its long-range displacements.

In a recent pioneering NEO-Ehrenfest study by Zhao et al. the ESIPT process in ohydroxy-benzaldehyde has been investigated [135]. Upon comparison of results obtained

using the NEO-Ehrenfest and the traditional Ehrenfest approach [137] with all-classical nuclei, the proton transfer reaction acceleration in the quantum case has been observed, which has been ascribed to the delocalization of the proton wave function, resulting in a smaller necessary displacement of the proton-accepting and proton-donating centers. Moreover, the kinetic isotope effect upon deuterium substitution has been reproduced at this level of theory.

#### 3.1.4. Summary of the NA-MQC Dynamic ESIPT Simulations

To summarize the section dedicated to mixed quantum–classical ESIPT studies, we again highlight the grea<sup>t</sup> contributions of the NA-MQD dynamic methods for the field. By allowing real-time picturing of the proton transfer process, characteristic timescales and unforeseen reaction mechanisms can be modeled at this level of theoretical description. While the methods share the mixed quantum–classical nature, they also still bear important differences, making them possible methods of choice for different conditions. In particular, TSH is a robust and probably most universal tool, reliable for the modeling of barrierless ESIPT, including complex situations, in which multiple PTs or competition from other photoreaction channels needs to be taken into account. The AIMS approach, formally more exact, may also be generally applied to this class of processes, as long as it does not become prohibitively expensive due to the extended molecular system size. At the same time, when nuclear quantum effects of protons are expected to play a role, such as in the barrier-restricted ESIPT case, the NEO-Ehrenfest method may be considered a good choice.

#### *3.2. Quantum Dynamics Methods for ESIPT Simulations*

Despite many useful conclusions on the ESIPT reaction course that may be taken from the NA-MQC dynamics, there are situations in which one needs to advance even further with the level of the system's dynamic description, up to the point of full quantum treatment of all the species, including nuclei. As has been already pointed out, the most typical reason of adopting this approach is when tunneling through an energy barrier along the ESIPT path needs to be included, i.e., when highly accurate rates or proton-transfer equilibrium have to be characterized. Another situation calling for the QD treatment is when a strongly nonadiabatic ESIPT mechanism is expected, e.g., when trivial interstate crossings are present, potentially threatening the correct NA-MQC dynamics performance. The latter problem, however, has been, in recent years, partially resolved by the successful design of correction strategies to several NA-MQC protocols [138–140].
