Exploring the Photophysical Properties of Molecular Systems Using Excited State Accelerated ab Initio Molecular Dynamics.

In the present work, we employ excited state accelerated ab initio molecular dynamics (A-AIMD) to efficiently study the excited state energy landscape and photophysical topology of a variety of molecular systems. In particular, we focus on two important challenges for the modeling of excited electronic states: (i) the identification and characterization of conical intersections and crossing seams, in order to predict different and often competing radiationless decay mechanisms, and (ii) the description of the solvent effect on the absorption and emission spectra of chemical species in solution. In particular, using as examples the Schiff bases formaldimine and salicylidenaniline, we show that A-AIMD can be readily employed to explore the conformational space around crossing seams in molecular systems with very different photochemistry. Using acetone in water as an example, we demonstrate that the enhanced configurational space sampling may be used to accurately and efficiently describe both the prominent features and line-shapes of absorption and emission spectra.

Introduction

Electronically excited states play a prominent role in many different areas of chemical and condensed matter physics and organic chemistry.[1,2] While nonadiabatic processes, such as irradiation induced DNA damage and repair mechanisms,[3,4] are the hallmark of photochemical and photobiological reactions, the fluorescent properties of excited state systems are increasingly being employed as biomolecular probes.[5,6] Despite their obvious importance, a detailed investigation of the excited state remains a challenge to experimentalists and theoreticians alike. Recent developments in laser chemistry have afforded significant advances in the study of light-induced dynamic phenomena.[7] However, due to the short excited state lifetime, the information obtained from such experiments is limited to a very small region of the excited state energy landscape. For those systems exhibiting longer excited state lifetimes, the interpretation of emission spectra is complicated by the fact that a detailed configurational representation of the experimental data requires knowledge of both the excited state and ground state energy surfaces. From a theoretical perspective, the accurate and proper computation of the excited state polyelectronic wave function is a decisive step to obtaining an accurate representation of the photophysical properties of the system at hand. In light of this, a variety of post-Hartree–Fock methods including coupled cluster, configuration interaction,[8] and multiconfigurational self-consistent field theory[9] have been developed. In comparison to the ground state analogue, a single Slater determinant is no longer appropriate to describe the excited state wave function, and a multiconfigurational description is required. Due to the intensive CPU time and large memory requirements, the above methods are often limited to “static” single point energy calculations for small molecular systems. The advent of approximate density-functional theory (DFT)-based methods for the calculation of the excited state wave function, such as restricted open-shell Kohn–Sham (ROKS)[10] theory and, more recently, time-dependent DFT,[11] has facilitated the development of efficient excited-state ab initio molecular dynamics (AIMD) simulations. Such AIMD simulations are readily employed to study the detailed molecular motions of small isolated molecules and complex systems in the excited state. Nevertheless, despite the substantial increase in available computational resources and the development of more efficient simulation algorithms, such as the well-known Car–Parrinello approach,[12] AIMD simulations are generally limited to time-scales of tens to hundreds of picoseconds. This is often too short to fully sample the ground and excited state surfaces.

Recently, we implemented the accelerated molecular dynamics method[13] (AMD) in the framework of ab initio molecular dynamics.[14,15] Accelerated ab initio molecular dynamics, A-AIMD, represents a highly efficient and robust configurational space sampling method that allows for the study of slow molecular motions and rare events. A-AIMD has successfully been applied to study both slow conformational transitions and chemical reactions in the ground state[14,15] occurring on time scales many orders of magnitude longer than those accessible using standard AIMD methods. In the present work, we extend the application of A-AIMD to the study of excited state systems and their photophysical properties. We first demonstrate that A-AIMD can be readily implemented as an efficient tool for studying the excited state energy surface, identifying both stable and metastable local energy minima and quantifying the magnitude of the associated energy barriers that separate these states. We then demonstrate how the excited state A-AIMD method can be used to efficiently map out the photophysical topology of the system of interest (i.e., characterizing the variation in the vertical de-excitation energy gap as a function of the relevant internal degrees of freedom of the system). This information provides a detailed insight into the photophysical properties of the system, allowing the prediction of both excited state reactions and relaxation phenomena.

The work presented in this paper is divided into two sections: In the first section, we focus on systems exhibiting fast radiationless decay processes, focusing on the identification of conical intersections (CIs) and extended crossing seams. In the second section, we demonstrate how A-AIMD can be used to describe the thermally averaged solvent effect on the fluorescence spectra of molecular systems, most notably making use of the enhanced configurational space sampling obtained by A-AIMD to accurately describe both the prominent features and line-shapes of absorption and emission spectra.

Radiationless Decay Process: Identification and Characterization of Extended Crossing Seams and Conical Intersection

Extended crossing seams are sets of molecular geometries where the energy of two electronic states varies, preserving their degeneracy. Among this collection of points, those presenting nonvanishing diabatic couplings[16] are known as conical intersections (CIs), which play an important role in photochemistry[17−19] and photobiology.[20−23] CIs control the deactivation pathway of the electronically excited system back to the ground state. As very little information about CIs can be inferred directly from experimental data, this is clearly one area of research where theory can provide unprecedented information to complement experimental studies, providing answers to important chemical and biochemical questions. The last 20 years has seen considerable progress in the development of new algorithms for localizing CIs algorithms. Traditionally, three methods have been employed: (i) Lagrange–Newton methods that minimize a Lagrangian in the presence of constraints that act to maintain the degeneracy of the electronic states,[24−26] (ii) gradient projection based techniques,[27−29] and (iii) penalty function-based methods.[30] A comparative analysis of these approaches has recently been presented by Thiel et al.[31] However, these procedures are primarily designed for locating single optimized molecular structures, known as minimum energy crossing points (MECPs).[16] However, the location and identification of CIs with these standard methods becomes very challenging when multiple MECPs are required to properly describe the intersection space relevant to the photochemical landscape of interest. More recently, a variety of computational approaches have been developed to address this problem, providing a more rigorous exploration of the conformational subspace around the CI, including the identification of extended conical intersection seams based on the definition of reaction coordinates,[26,32,33] or by direct search without the need of any guess.[34] In the present work, we introduce an alternative direct search approach for the efficient exploration of CIs and extended crossing seams based on an enhanced configurational space sampling method, A-AIMD, which complements previous methods. Using two well-studied small molecular systems, the Schiff bases fomaldimine (FD)[35] and salicyclideneanile (SA),[36] we demonstrate how the efficient and robust enhanced configurational space sampling afforded by the accelerated molecular dynamics approach in the excited state allows for a full characterization of the excited state surface and the photophysical topology of these two systems, including the rapid identification of CIs and extended crossing seams on the (adiabatic) excited state potential energy surface.

Theory and Computational Details

Accelerated ab Intio Molecular Dynamics

A comprehensive account of the accelerated ab initio molecular dynamics (A-AIMD) approach is available in the literature,[14,15] and we only provide a brief outline of the essential details here. In the standard AMD formalism,[13] a continuous non-negative bias potential, ΔV(r), is defined such that when the true (underlying) potential of the system, V(r), lies below a certain, predefined threshold “boost” energy, Eb, the simulation is performed on a modified potential, V*(r) = V(r) + ΔV(r), but when V(r) ≥ Eb, the simulation is performed on the true potential [V*(r) = V(r)]. The modified potential, V*(r), is related to the true potential, V(r), bias potential, ΔV(r), and boost energy, Eb, by[13]and the bias potential, ΔV(), is defined as:In the framework of AIMD, the true potential, V(r), is defined as the density functional energy. The application of the bias potential results in a raising and flattening of the PES, thus enhancing the exchange rate between low energy states. The extent of acceleration (i.e., how aggressively we enhance the configurational space sampling) is determined by the choice of the boost energy, Eb, and the acceleration parameter, α. Configurational space sampling can be enhanced by either increasing the boost energy or decreasing α. The reader is referred to ref (13) for more details.

Obtaining Accurate Free Energy Surfaces Using A-AIMD

As the bias potential destabilizes low energy regions of conformational space on the potential energy landscape, the population statistics on the modified potential are inherently incorrect. However, as it is known at each step in the simulation exactly how much the system is destabilized, one can readily retrieve the correct population statistics by reweighting each point in the configuration space on the modified potential by the strength of the Boltzmann factor of the bias energy, exp[βΔV(r)], at that particular point. As a result, A-AIMD yields a canonical average of an observable such that thermodynamic and other equilibrium properties can be accurately determined.

Using the canonical Boltzmann reweighting formulism described above, representative free energy surfaces for the system of interest can be obtained as follows: Having defined an appropriate set of internal degrees of freedom (see the Results and Discussions section for specific details), the configurational subspace is divided into bins, and the structures collected across the A-AIMD trajectories are allocated to their respective bin (j). The effective population statistic for each structure, i, is given by exp[βΔV(i)]. For each bin, the population statistics are summed across all A-AIMD trajectories to give an effective total population in that bin, pop(j), and the free energy surface, ΔG(j), is obtained aswhere pop(j)max is the effective total population of the most populated bin.

In the case of FD, we show that the ground and excited state free energy surfaces obtained from the A-AIMD simulations are directly compared to a metadynamics[37] reference (see the Computational Details section).

Mapping out the Photophysical Topology

Having efficiently and exhaustively sampled the configurational space of the molecular system of interest in the excited state, the photophysical topology of the system was mapped out by identifying the appropriate internal degrees of freedom and dividing this configurational subspace into bins. Each structure in the A-AIMD simulations was then allocated to a given bin. Given the large number of structures obtained from the A-AIMD trajectories, 10 structures in each bin were randomly chosen, and a ground state energy calculation was performed. The vertical de-excitation energy gap associated with each bin was then averaged over the 10 representative structures. The specific internal degrees of freedom used for each system are defined explicitly in the Results and Discussions section.

Computational Details

All molecular dynamics simulations were performed using an in-house modified version of the AIMD 3.13 package.[38] In the case of formaldimine (FD), the system was placed at the center of an orthorhombic box with side lengths (17.6 × 18.3 × 14.0 au3). Similarly, for salicylidenaniline (SA), the system was placed at the center of a box with side lengths (22.9 × 35.6 × 14.0 au3). In all ground state (S0) simulations, the electronic structure problem was solved with density functional theory (DFT), and the restricted open-shell Kohn–Sham (ROKS) method was implemented for simulations in the singlet (S1) excited state. For all systems, the Becke (B) exchange[39] and Lee, Yang, Parr correlation functional[40] were employed. Core electrons were treated using the norm-conserving pseudopotentials of Troullier and Martins, and the valence orbitals were expanded in a plane-wave basis up to an energy cutoff of 70 Ry. In all ground and excited state (accelerated and nonaccelerated) simulations, a fictitious mass of 400 au was ascribed to the electronic degrees of freedom, and the coupled equations of motion were solved using the velocity Verlet algorithm with a time step of 4 au. All AIMD simulations (accelerated and nonaccelerated) were carried out at T = 300 K. Standard AIMD simulations were performed using a Nosé–Hoover chain thermostat[41] on the ions, while in the case of the accelerated AIMD simulations, a Nosé–Hoover thermostat was applied to both the electronic and nuclear degrees of freedom.

The ground and excited state metadynamics[37] simulations used for comparative analysis in the FD study (see Results and Discussions) were performed under exactly the same physical conditions as the standard AIMD simulations described above. The angular and torsional coordinates (θ,ϕ) were employed as the collective variables. Metadynamics trajectories were run for 1 000 000 MD steps (∼100 ps) and Gaussian hills with a half-width of 2.5° and a depth of 0.2 kcal/mol were added every 200 MD steps.

Before proceeding to the Results and Discussions section, it should be noted that there are several limitations to the method presented in this work: First, it is clear that the accuracy of the relevant excited state free energy surfaces and photophysical topologies, including the identified CIs and extended conical intersection seams, is strongly dependent on the integrity of the employed polyelectronic wave function representation. In the present examples, we have used the restricted open-shell Kohn–Sham (ROKS) method to describe the excited state. The enhanced configurational space sampling afforded by A-AIMD obviously does not have any effect on the accuracy of the polyelectronic wave function representation, and, as in all cases when using DFT-based methods, it is important to validate the results using more sophisticated post-Hartree–Fock methods (see Introduction). Indeed, we would like to point out that for the two systems investigated here, FD and SA, the accuracy of the ROKS approach has already been validated by comparison to more sophisticated multireference configuration interaction methods.[36,42,43] Second, the excited state A-AIMD trajectories presented here are strictly performed under adiabatic conditions, and therefore explore the adiabatic excited state energy surface. As a result, the excited state populations and the corresponding free energies in the absence of relaxation are only approximations that do not have a direct physical meaning. In addition, no direct information is obtained concerning the specific dynamic processes that occur after vertical excitation that drive the system toward the conical intersection. Nevertheless, A-AIMD simulations provide important information about the topology of the excited and ground state surfaces and allow converged sampling of both surfaces. Finally, it is also well-known that radiationless decay processes occur not only close to the crossing seam, where the energy gap between the intersecting states is small, but also in regions of configurational space away from it.[3,44,45] For these cases, a variety of theoretical approaches, such as the method of fewest switches (MFS),[46] which involve the calculation of a transition probability between the relevant electronic states have been developed to address such relaxation mechanisms. With specific regard to the MSF approach, the transition probability is strongly dependent on the temporal velocity of the wave function, and as the evolution of the system in A-AIMD occurs on a nonlinear time scale, the application of such methods is beyond the scope of the present study.

Results and Discussions

Conical Intersections and Photoisomerization in Formaldimine

Formaldimine (FD), the smallest unprotonated Schiff base, is the prototypical system for studying cis–trans photoisomerization around a double bond. The photophysical properties and associated nonadiabatic photoisomerization mechanism for this system have been rigorously studied using both nonadiabatic ab initio MD[35] and high-level “static” multireference configuration interaction (MRCI) calculations.[42,43] In the ground state, S0, FD adopts a planar geometry as depicted in Figure 1a. After vertical excitation to the lowest excited singlet state, S1, the system rapidly relaxes toward the local energy minimum on the excited state PES and the N–H bond vector twists out of the molecular plane. On approaching the orthogonal twist geometry, the system enters the CI region resulting in strong nonadiabatic coupling between the S1 and S0 states. After a surface hop from the S1 to the S0 state, FD returns to a planar geometry, leading either to the photoisomerized product or to the regeneration of the reactant state (see Figure 1b).

Figure 1

(a) Molecular structures of the two isomeric forms of formaldimine. The degrees of freedom ϕ and θ discussed in this work are indicated. (b) Scheme representing the photoreaction pathways of formaldimine. (c) Molecular structures of the two tautomeric forms of salicylideneaniline. The degrees of freedom (ω1,r1) and (ω2,r2) discussed in this work are indicated (dihedral angles highlighted in red). (d) Scheme representing the photochemical landscape of salicylideneaniline, including the ESIPT and the nonadiabatic pathways.

In order to study the configurational space sampling properties of FD, a multiple-copy ab initio MD simulation strategy was employed: After performing a standard geometry optimization and equilibration procedure, five independent AIMD simulations were performed for 500 000 MD steps (the equivalent of ∼50 ps) at 300 K in both the ground and excited state. Analysis of the resulting trajectories revealed that the most significant variations in the configurational geometry involved the H–C=N–H dihedral angle, ϕ, and the C=N–H angle, θ. Using these two internal degrees of freedom, the configurational space sampling obtained from all five independent AIMD simulations in the ground and excited state are depicted in green in Figure 2a and b, respectively. Notably, in both electronic states, the configurational space sampling afforded by the standard AIMD simulations is rather limited: In the ground state, the H–C=N–H dihedral angle varies between −25 and +25° and the C=N–H angle varies between 100 and 125°. These variations are the product of local thermal fluctuations at 300 K about the geometry-optimized structure with coordinates [ϕ,θ] = [0.0,111.36]. The observed thermal fluctuations in the excited state are slightly more pronounced with the ϕ angle varying by ±40° and the θ angle varying by ±20° about the geometry-optimized structure located at [ϕ,θ] = [124.69,101.71].

Figure 2

(a,b) Conformational space sampled in formaldimine by standard AIMD (green) and A-AIMD (red) simulations, on the ground (left) and first singlet excited (right) electronic states. (c,d) Free energy surfaces derived from the conformational space sampled by A-AIMD simulations in a and b, respectively. The most energetically stable state has been defined as the most sampled state. (e,f) Free energy surfaces obtained from metadynamics simulations. The ϕ and θ coordinates have been used as the two collective coordinates. Energies are expressed in kcal/mol.

The standard AIMD trajectories described above served as the initiation point for the accelerated ab initio MD simulations. Using a fixed value for the acceleration parameter, α = 0.01 au, a series of short [100 000 MD step] A-AIMD simulations were performed across a broad range of acceleration levels in both the ground and excited state. The acceleration level was controlled by varying the boost energy, [Eb – V0], where V0 is the average density functional energy obtained from the standard AIMD simulations. The choice of the acceleration level is exceedingly important: If the acceleration level is too low, ostensibly no enhanced configurational space sampling is observed. In contrast, applying too high an acceleration level results in a molecular explosion or dissociation of the system. In the present case of FD, [Eb – V0] was varied between 0.050 au and 0.100 au in steps of 0.005 au. The optimal acceleration levels for enhanced configurational space sampling were found to be [Eb – V0] = 0.085 au and [Eb – V0] = 0.060 au for the ground and excited state, respectively. These “optimal” acceleration levels afforded the most efficient configurational space sampling, while not rendering dissociation of the system. Once the optimal acceleration levels in the ground and excited state had been determined, five independent A-AIMD simulations were performed for 500 000 MD steps.

The configurational space sampling obtained from the five independent A-AIMD simulations in the ground and excited state is depicted in red in Figure 2a and b, respectively. Clearly, the accelerated MD simulations provide a much more complete and robust description of the conformational space sampling properties of FD compared to the standard AIMD trajectories. Most importantly, in the A-AIMD simulations, we observe multiple isomerization events in both the ground and excited state. Interestingly, the A-AIMD simulations reveal significant differences in the specific adiabatic isomerization pathways in the two electronic states: In the ground state, as the ϕ angle flips from 0° to ±180°, the C=N–H angle, θ, is seen to significantly increase. At the orthogonal twist geometry (ϕ = ±90°), the C=N–H angle adopts values between +140° and +170°. By contrast, in the excited state adiabatic isomerization process, the θ angle adopts values between 90° and 140° at ϕ = 0.

The associated free energy surfaces in the (ϕ,θ) configurational subspace (see Theory and Computational Details) for the ground and excited state of FD are depicted in Figure 2c and d, respectively. The lowest free energy barrier for adiabatic isomerization on the ground state surface was found to be 0.052 au, and the associated transition state is located at [ϕ,θ] = [±100,150]. Similarly, in the S1 excited state, the lowest free energy barrier for isomerization is 0.037 au, and the associated transition state is located at [ϕ,θ] = [0,120]. The free energy barriers obtained from the A-AIMD simulations presented here are in excellent agreement with previous theoretical and experimental studies, where the free energy of isomerization in the S0 and S1 states was found to be 0.048 au and 0.040 au, respectively.[10,42]

In order to provide a reference to validate the A-AIMD free energy surfaces shown in Figures 2c and d, we performed two metadynamics[37] simulations using the ϕ and θ internal degrees of freedom as the collective variables (see Theory and Computational Details). The free energy surfaces for the ground and excited state obtained from the metadynamics simulations are depicted in Figure 2e and f, respectively. Clearly, the adiabatic free energy surfaces provided by both methods are very similar, and the extent of configurational space sampling obtained from the A-AIMD and metadynamics simulations is also similar. This is an encouraging result when one considers that in the accelerated molecular dynamics approach there is no a priori definition of a reaction coordinate or a set of collective variables.

The comparative analysis of the present A-AMD approach and the more well-known metadynamics method clearly demonstrates that one can obtain an accurate energetic description of the excited state energy surface, which is both an important and necessary prerequisite to acquiring a detailed picture of the photophysical topology of the system. The excited state free energy surface and associated population statistics acquire much more significance when studying a system with considerably longer excited state lifetimes, such as in the case of the acetone-in-water system presented later in this paper.

The efficiency of the configurational space sampling obtained from the A-AIMD simulations can be readily inferred from the magnitude of the free energy barriers that are traversed during the simulations of FD in the ground and excited state. According to transition state theory (TST), the predicted rate constant, k, for a transition between two states separated by a free energy barrier, ΔG, is given bywhere κ is the transmission coefficient and kB is the Boltzmann constant. As discussed above, the free energy barrier for adiabatic isomerization of FD in the ground state is 0.052 au. In the fast limit (κ = 1), the TST predicted rate constant for this process is therefore 1.62 × 10–11 s–1. Given that a 50-ps ground state A-AIMD trajectory produces ∼4–5 isomerization events, on average, in the A-AIMD simulations, an isomerization event is observed every 100 000 MD steps (∼10 ps). In light of this observation, the application of the bias potential enhances the configurational space sampling rate by a factor of 1021.

In order to map out the photophysical topology of FD, the configurational subspace (ϕ,θ) shown in Figure 1b was divided into bins; the vertical de-excitation energy gap across the entire excited state configurational space was calculated (see Theory and Computational Details). Two large crossing seam regions are clearly visible, located at [±50 < ϕ < ±150, 70 < θ < 120]. Notably, outside of these two regions, the vertical de-excitation energy gap significantly increases to values on the order of 20–30 kcal/mol. Closer inspection of Figure 3a reveals that there exists a considerable amount of fine-structure on the vertical deactivation energy surface. This observation raises the important question as to whether this fine structure is really a physical phenomenon of the system or whether it is merely statistical noise invoked by the application of the bias potential in the A-AIMD simulations. In order to address this issue, 200 excited state structures along the ϕ coordinate were selected from the A-AIMD simulations with a fixed C=N–H angle, θ = 90. A partial geometry optimization procedure was performed for each structure, in which the internal (ϕ,θ) coordinates were constrained and the vertical de-excitation gap was then recalculated. Figure 3b depicts the S1–S0 energy profile along the ϕ coordinate obtained from the A-AIMD structures collected “on-the-fly” and their geometry-optimized counterparts. The two resulting energy profiles are remarkably similar, confirming the fact that small local molecular distortions arising from the application of the bias potential in the A-AIMD simulations do not compromise the accuracy of the vertical deactivation energy statistics. A very interesting feature in the energy–gap profile depicted in Figure 3b is observed deep in the crossing seam region: Between coordinates [−110 < ϕ < −75] and again at [75 < θ < 110], the vertical deactivation energy gap is seen to increase by ∼2 kcal/mol. This small variation in the energy gap is due to a change in the pyramidalization angle about the C atom, which reaches its maximum at a pyramidalization angle of ∼27°. This result is in agreement with a previous study by Doltsinis et al., who observed a loss of degeneracy of the crossing seam at a pyramidalization angle of ∼30°.[35] This result provides additional validation to the degree of detail that can be achieved in both the conformational space and free energy sampling, otherwise not possible using standard MD or single point electronic calculations. A second interesting feature shown in Figure 3a and b is the observation of negative [S1–S0] energy gaps for a limited region of the A-AMD-sampled configurational space for FD. Such small negative energy gaps have been observed previously[35] and are invariably the result of an error arising from the fact that the ground and excited electronic state calculations are performed independently. As such, we would like to point out that these small artifacts are not a result of the application of the A-AMD biasing potential but rather a product of the polyelectronic wave function representation employed in this work.

Figure 3

(a) Variation in the energy gap [S1–S0] for formaldimine as a function of the internal coordinates ϕ and θ, obtained from single point calculations in randomly chosen structures sampled across the A-AIMD trajectory (see Theory and Computational Details). The solid line indicates the cut of the surface for which the direct relaxed difference has been computed. Dark blue indicates regions of energetic degeneracy between S0 and S1. Evolution toward red indicates a loss of energetic degeneracy. (b) Direct (dashed line) and relaxed (solid line) S1–S0 energy gap profiles along the surface cut at ϕ = 90°. Energies are expressed in kcal/mol.

Conical Intersections and Proton Transfer in Salicylideneaniline

Using our A-AIMD approach on the FD system, we were able to exhaustively investigate the photophysical properties of the system and identify crossing seams from the data collected in S1 at a higher sampling rate compared to standard AIMD approaches. In light of this successful initial study, we moved on to a more complex system. The aromatic Schiff base salicylideneaniline (SA) represents a classic example of an excited state intramolecular proton transfer (ESIPT) reaction. The reactant and product forms of SA, corresponding to the enol and cis-keto isomeric forms, respectively, are connected by an ESIPT event after photoexcitation to S1, as shown in Figure 1c. At this point, the electronically excited cis-keto tautomer can either undergo a cis–trans isomerization to form a trans-keto tautomer or directly deactivate to S0. Despite this simple chemical scheme, much controversy has arisen in the past 40 years concerning the experimental determination of the time scale for the ESIPT process in SA, with measurements ranging from 750 fs to less than 50 fs. The reader is referred to ref (36) for an updated list of references on the SA studies. Sekiya and collaborators proposed the presence of two nonadiabatic processes competing with the proton transfer event:[47] The first process involves a highly deformed enol conformation, while the second relaxation mechanism is associated with the cis–trans isomerization path of the cis-keto tautomer, as depicted in Figure 1d. Such competing processes explained the discrepancy in the experimental observations. Very recently, “static” TDDFT calculations (validated with high-level multiconfigurational ab initio calculations) and quantum dynamical simulations confirmed the presence of two CI regions associated with both the enol and cis-keto states.[36]

In order to study the competition between the two relaxation pathways on the excited state surface, we performed a series of AIMD and A-AIMD simulations similar to those described previously for FD. A standard geometry optimization and equilibration procedure was performed on the enol tautomer in the excited state. The geometry optimization led directly to the cis-keto tautomer of SA, indicating that the enol state is not a stable species in S1. Five independent AIMD simulations were performed for a total of 1 500 000 MD steps (∼150 ps) at 300 K on the cis-keto tautomer in S1. The system remained in the cis-keto state during the simulations, and we were not able to observe a reverse ESIPT event. The O–H atomic distance was seen to fluctuate around the cis-keto geometry-optimized value (∼2 Å) in all the ∼150 ps trajectories, as depicted in green in Figure 4a. The principal geometrical fluctuations observed in the simulations can be described by two independent variables: the C–C=C–N dihedral angle (highlighted in red in Figure 1c), ω2, which varies between −50 and +50 degrees, and the C=C bond distance, r2, which varies between 1.32 and 1.53 Å. Notably, although large amplitude vibrations were seen along the ω2 coordinate (Figure 4c), cis–trans isomerization was not observed in the standard AIMD simulations.

Figure 4

(a) Evolution of the O–H distance of salicylideneaniline during ∼150 ps of standard AIMD (green) and A-AIMD (red) simulations in the S1 state. The initial structure of the simulations corresponds to a cis-keto tautomer. (b,c) Conformational space sampled in salicylideneaniline using standard AIMD (green) and A-AIMD (red) simulations in the S1 state in the enol and cis-keto regions, respectively. The intersection of the solid lines indicates the location of the two MECPs described in ref (36).

In order to improve the conformational space sampling, a set of short A-AIMD simulations were initiated. Using a fixed value of α = 0.01 au for the acceleration parameter, a series of 500 000 MD-step A-AIMD simulations were performed across a broad range of acceleration levels for the cis-keto state in S1, whereby [Eb – V0] was varied between 0.050 au and 0.200 au in steps of 0.05 au. The optimal acceleration level for enhanced configurational space sampling was found to be [Eb – V0] = 0.150 au for the excited state. This optimal acceleration level was employed for five independent accelerated ab initio MD simulations, each performed for a total of 1 500 000 MD steps (∼150 ps).

The analysis of our A-AIMD simulations indicated that the configurational space sampled within ∼150 ps covered both the cis-keto and the enol states. The O–H bond distance presented large fluctuations during the first 100 ps, ranging between 1.2 and 5.5 Å. After ∼110 ps, a reverse ESIPT event was observed, which led to a decrease of the O–H distance to ∼1.0 Å, as depicted in red in Figure 4a. A difference of ∼10 kcal/mol between the enol and cis-keto regions in the excited state was estimated, which explains why we did not observe a reverse ESIPT event during the standard AIMD simulations. By contrast, however, the application of the bias potential in the A-AIMD simulations facilitated the observation of a reverse ESIPT event after only 110 ps (see Figure 4a). The newly formed enol tautomer remained stable for the last ∼40 ps of the simulation, affording both ample conformational and energy statistics for the enol state. As such, two chemically different but relevant sections of the PES were successfully sampled from a single molecular geometry.

The configurational space sampled in the enol region revealed significant variations in both the C–C=N–C dihedral angle (highlighted in red in Figure 1c), ω1, and the C=N bond distance, r1, which may be compared directly to the prominent internal degrees of freedom identified for the cis-keto state [ω2, r2]. Interestingly, the rotation and stretching along double bonds in excited electronic states has been previously proposed as a common feature in nonadiabatic deactivation pathways for complex organic compounds.[48] The variation of these degrees of freedom sampled in the five independent A-AIMD simulations in the excited state are depicted in red in Figure 4b and c. As expected, the accelerated MD simulations provide considerably more exhaustive configurational space sampling compared to their AIMD counterparts. In the A-AIMD simulations, we observed multiple rotation events around the dihedral angles ω1 and ω2 in the enol and cis-keto regions, respectively, and complete rotation events through 360° were observed in both cases. Simultaneously, rather large fluctuations in the bond distances (1.22–1.67 Å for r1 and 1.23–1.75 Å for r2) were also observed. The intersection of the two straight lines depicted in Figure 4b and c defines the coordinates for the two MECPs described in a previous study of SA.[36] These two MECPs are located at the coordinates (92, 1.37) and (89, 1.47) for the enol and cis-keto regions, respectively. Notably, the conformational space sampled in the standard AIMD simulations for the cis-keto state (Figure 4c) does not include these coordinates, while the A-AIMD simulations readily encompass a considerably larger region of the configurational subspace including the associated MECP coordinates. The variation in the underlying PES explored in our A-AIMD simulations indicates that the rotational barriers associated with the dihedral angles ω1 and ω2 are 16.75 and 18.03 kcal/mol, respectively. As such, the twisted conformations for both the enol and cis-keto tautomers are clearly inaccessible in the standard AIMD simulations but can be readily observed in the A-AIMD trajectories.

The vertical de-excitation energy gaps mapped across the entire excited state configurational subspaces (ω1, r1) and (ω2, r2) for the enol and cis-keto tautomeric states are depicted in Figures 5a and b. Two large conical intersection regions are clearly visible on each plot, located at [±75 < ω1 < ±110, 1.40 < r1 < 1.55] and [±60 < ω2 < ±120, 1.40 < r2 < 1.65] for the enol and cis-keto states, respectively. These results show that a complete characterization of both crossing seams can be obtained from the A-AIMD simulations, even though they belong to very different regions of the PES.

Figure 5

(a,b) Variation in the energy gap [S1 – S0] of salicylideneaniline, in the enol and cis-keto regions, respectively. The surfaces are projected on the (ω1,r1) and (ω2,r2) subspaces and have been obtained from single point calculations in randomly chosen structures sampled in the A-AIMD trajectories (see Theory and Computational Details). Energies are expressed in kcal/mol. Dark blue indicates regions of energetic degeneracy between S0 and S1. Evolution toward red indicates a loss of energetic degeneracy.

The extended crossing seams present in the de-excitation energy maps (Figures 4 and 5) obtained from the A-AIMD simulations clearly suggest that there exist two competing radiationless relaxation decay processes which completely describe the photophysical properties of SA after vertical excitation: In the first of these relaxation decay mechanisms, the vertical excitation energy initiates a spontaneous proton transfer event which leads directly to the formation of the (excited) cis-keto state (Figure 1d), and a subsequent rotation in ω2 brings the system to the extended conical intersection seam (CI2) depicted in Figure 5b, resulting in either a cis-keto or trans-keto ground state photoisomerized product. In a second competing relaxation mechanism, the vertical excitation energy induces a rotation about the ω1 dihedral angle which breaks the hydrogen bond and results in the formation of a highly twisted enol conformation. In the twisted enol state, the C=N bond is weakened, and the r1 bond length increases up to 1.67 (see Figure 4b), bringing the system into the conical intersection region (CI1). Notably, the correlation between the ω1 and r1 coordinates assists the evolution of the enol tautomer toward the extended CI seam, an observation that cannot be ascertained from the MECP result but is readily confirmed by the exhaustive configurational space sampling provided by the A-AIMD simulations. The existence of two competing radiationaless relaxation decay processes predicted from this study readily explains the experimental observation that the quantum yield of the cis/trans-keto products is lower than one would expect when only considering a single nonadiabatic relaxation mechanism associated with the enol-to-cis-keto proton transfer event.

Using A-AIMD Simulations to Recapitulate and Interpret Absorption and Emission Spectra

In section I, we focused on the use of excited state A-AIMD to efficiently identify and explore CIs and conical intersection seams allowing for the prediction of diverse and often competing radiationless decay mechanisms in systems exhibiting extremely short excited state lifetimes. We now extend our study to the direct recapitulation and interpretation of absorption and emission spectra, whereby we make use of both the enhanced configurational space sampling properties and the approximate free energy statistics afforded by the accelerated molecular dynamics method. The ab initio prediction of absorption and emission spectra in the condensed phase represents a considerable challenge, since it requires both the use of an accurate quantum mechanical method and a correct description of the solvent effect and temperature. At room temperature, our ability to simulate all the relevant conformational states of a molecule and the phase space statistics of the solvent is often limited. Using a proto-typical system, acetone in bulk water, we now demonstrate how the application of A-AIMD leads to an improved description of both absorption and emission spectra, when compared to conventional AIMD simulations. Acetone in water was chosen for this test as it has been studied extensively both theoretically and experimentally.

AIMD and A-AIMD simulations were performed on a periodically repeating cubic unit cell of edge length 22.7 au containing a single acetone molecule and 34 H2O molecules, which had already been brought to equilibrium under the correct density conditions using classical molecular dynamics. Initially, the simulation system was further equilibrated with AIMD at 300 K for 15 ps. Production run AIMD and A-AIMD trajectories were carried out for 10 ps in both the ground and first-excited state. All simulation details including the choice of density functional (BLYP) and employed pseudopotentials, the fictitious electron mass (400 au), and the time step (4 au) were identical to those described previously for the FD and SA systems. For the A-AIMD simulations, the acceleration parameters were set to α = 0.1 au, and [Eb – V0] = 0.15 au, based on a previous A-AIMD study of bulk water.[15] It should be noted that, while the application of a larger bias potential is possible, using lower acceleration levels guarantees that the numerical error remains minimal during the free energy reweighting process. For both ground and excited state trajectories, the atomic coordinates were saved each 10 MD steps, and the respective vertical excitation and vertical de-excitation energy gap for each simulation frame was calculated in order to recapitulate the absorption and emission spectra. In the case of the accelerated MD simulations, the relative intensity of the absorption/emission spectra was obtained to a first approximation using the population statistics obtained from the canonical Boltzmann reweighting procedure described previously.

In Figure 6b, the absorption and emission spectra calculated from the A-AIMD trajectories are shown. A modest, but nevertheless not insignificant, improvement in both the theoretical recapitulation of both these spectra was obtained from the A-AIMD simulations: Average values for the emission spectrum [2.58(0.26) and 2.50(0.25) eV for A-AIMD and AIMD, respectively] should be compared to a range of values, between 2.99[49] and 3.06 eV,[50,51] measured experimentally The deviation between ROKS and experiments has been explained in a previous QM/MM study[52] where ROKS results for acetone in water were found to be systematically red-shifted with respect to TDDFT calculations by a constant value, estimated around −0.28 eV.[52] Average values for the absorption spectrum [4.11(0.12) and 3.95(0.11) eV for A-AIMD and AIMD, respectively] are also in reasonable agreement with the experimental observation, 4.69 eV.[53] In addition to improving the average values of the absorption and emission distributions, A-AIMD also improves the quality of the population statistics and hence calculated spectral lineshapes. In particular, irregularity in the spectral histograms suggests that 10 ps of AIMD are too short to obtain converged sampling, whereas clear improvement can be seen with A-AIMD simulations. As an example of the enhanced configurational space sampling observed in the A-AIMD simulations, Figure 6a shows a comparison for the methyl group rotation observed in both the standard AIMD and A-AIMD simulations in the excited state. Notably enhanced phase space sampling for the solvent water molecules around the acetone system was also observed in the A-AIMD simulations.

Figure 6

(a) The rotation of the methyl group (shown here for a representative excited state simulation on S1), one of the motions that occurs more rapidly in A-AIMD (blue) when compared to AIMD simulations (black). (b) Absorption (full lines) and fluorescence (broken lines) spectra of acetone, as calculated with A-AIMD (blue) and a conventional AIMD (black) using the same number of MD steps.

Conclusions

In this paper, we have performed excited state accelerated ab initio molecular dynamics (A-AIMD) to explore the excited state energy landscape and photophysical properties of a variety of systems. In an initial case study on formaldamine, we have found that A-AIMD can be readily employed to enhance configurational space sampling for the study of molecular systems in the excited state, crossing energy barriers on the excited state energy landscape of up to 30 kcal/mol and affording a quantitative description of the excited state free energy surface. Notably, in comparison to other popular enhanced sampling methods that involve the specific definition of a reaction coordinate, A-AIMD does not require any a priori understanding of the underlying free energy surface. The method was subsequently used to study the photophysical topologies of salicylideneaniline (SA). It allowed the identification and characterization of conical intersections (CIs) and extended conical intersection seams and the prediction of two different, competing radiationless decay processes. For acetone in water, we demonstrated that the enhanced configurational space sampling obtained by A-AIMD may be used to compute converged absorption and emission spectra. The A-AIMD method is general and may be implemented in combination with any methodological representation of the polyelectronic wave function. Given the rapid and sustained increase in available computational resources, it may soon be possible to combine A-AIMD with more sophisticated post-Hartree–Fock representations of the polyelectronic wave function in order to efficiently obtain an accurate and complete description of both ground and excited state energy surfaces.