Xiaofang Wanga,
Jirui Guoa,
Tanping Li*a and
Zhiyi Weib
aSchool of Physics and Optoelectronic Engineering, Xidian University, Xi'an, 710071, People's Republic of China. E-mail: tpli@xidian.edu.cn
bBeijing National Laboratory for Condensed Matter Physics, Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China
First published on 13th May 2020
The connections between the non-equilibrium solvation dynamics upon optical transitions and the system's equilibrium fluctuations are explored in aqueous liquid. Linear response theory correlates time-dependent fluorescence with the equilibrium time correlation functions. In the previous work [T. Li, J. Chem. Theory Comput., 2017, 13, 1867], Stokes shift was explicitly decomposed into the contributions of various order time correlation functions on the excited state surface. Gaussian fluctuations of the solute–solvent interactions validate linear response theory. Correspondingly, the deviation of the Gaussian statistics causes the inefficiency of linear response evaluation. The above mechanism is thoroughly tested in this study. By employing molecular simulations, multiple non-equilibrium processes, not necessarily initiated from the ground state equilibrium minimum, were examined for tryptophan. Both the success and failure of linear response theory are found for this simple system and the mechanism is analyzed. These observations, assisted by the width dynamics, the initial state linear response approach, and the variation of the solvation structures, integrally verify the virtue of the excited state Gaussian statistics on the dynamics of Stokes shift.
The application of linear response theory on Stokes shift has achieved remarkable success.9,13–15,17–22 Nevertheless, the validity of linear response theory is interpreted differently in literatures.23–29 Using the equilibrium properties of the unperturbed surface, the ground state linear response approach is employed in numerous systems.13,15,16,30 This treatment, originating from the perturbation theory,12 assumes the perturbation to be much less than the Boltzmann energy kBT and is often inapplicable for the Stokes shift. Gaussian statistics on the excited state surface ensures the rigorous validity of linear response theory, even though the initial configurations of the non-equilibrium process are far from the equilibrium minimum.25,31 This argument in favour of the excited-state linear response theory has been well-supported by a number of studies. For example, Li's theoretical investigations show the consistent dynamics between the Stokes shift and excited state equilibrium fluctuations for the protein Staphylococcus nuclease, which are clearly different from the ground state equilibrium ones.26,32 Heid et al. have thoroughly tested the ground and excited state linear response theory for different benzene-like solutes and concluded that the ground state linear response approach is not recommended.27,33
Laird and Thompson have previously shown that linear response theory is validated due to the factorizing of the high order time correlation functions into the lower ones.25,31,34 Li and co-workers further explicitly decomposed the total energy shift into the summation of all the time correlation functions.26,35 Whether those orders with the dominant contributions obey the same dynamics as the linear one, as characterized by Gaussian statistics, sets the validity of linear response theory. This mechanism is supported by investigations in the protein.26 In this context, we further examined the linear response for a single tryptophan in aqueous solution using molecular simulations (MD). The fluorescence Stokes shift experiments of tryptophan have been performed extensively.8,36–40 The internal conversion between the two excited state singlets of the tryptophan indole occurs within sub-picoseconds. The lifetime of the chromophore is longer than 500 picoseconds.36 These properties enable tryptophan to be an ideal optical probe to study the hydration dynamics in bulk water and proteins. Molecular dynamics simulations of tryptophan in bulk water were performed owing to the following merits. Firstly, the simulation results can be compared with the experimental results. Furthermore, it is feasible to acquire sufficient sampling data on the configurational space for such a simple system, which is always a bottleneck in the complex proteins. The multiple non-equilibrium MD simulations were designed to propagate on the excited state surface with the initial configurations not necessarily sampled from the ground state equilibrium minimum. The connections with the system's equilibrium statistics were examined to analyze the origin of linear response theory. The remainder of this paper is organized as follows. The descriptions of the theory and the simulation details are firstly presented. The simulation results are shown in the discussions. Finally, the results are summarized in the conclusion session.
In literature, Stokes shift is also evaluated by .5,25 By applying Taylor expansion on eδΔE, ΔES(t) is correlated with the various order equilibrium time correlation 〈δΔE(t)δΔE(0)n〉e on the excited state surface. Our previous work further rearranges the time correlation function into26
![]() | (1) |
The accumulation 〈δΔE(t)δΔE(0)n〉e is thus converted into an integration over δΔE. δΔE(0,t) is interpreted as follows. For a given value of δΔE, the related configurations are sampled over the excited state equilibrium fluctuations. Upon the evolution of the time interval t, each configuration results in a solute–solvent interaction. δΔE(0,t) is referred to as the average of the above propagated energies. Therefore, δΔE(0,t) relies on both δΔE and time t. The value of δΔE(0,0) at time zero is naturally δΔE itself. Pn(δΔE) = δΔEnPe(δΔE) is the weight function and Pe(δΔE) is the excited state equilibrium distribution. Sn(t) is thus pictured as the ensemble summation over a branch of δΔE(0,t) weighted with the weight function Pn(δΔE). Sn(t) relaxes and decays to zero in the long time limit. The decay rate is characterized by . Obviously, those δΔE, like the stationary points on the weight function Pn(δΔE), take the most significant contributions over the summation and dominantly control the dynamics of Sn(t). It is noted that Sn(0) = ∫d(δΔE)δΔEPn(δΔE) can be calculated for a given distribution Pe(δΔE) and the associated weight function Pn(δΔE).
Dynamic Stokes shift is correlated with all the time correlation functions Sn(t).25,31 The explicit decomposition can be found in our previous work,26 and here is a concise presentation. For the non-equilibrium relaxation propagating on the excited state surface, the initial distribution of the solute–solvent interactions is represented as the summation of all the weight functions by
PI(δΔE) = ∑fnPn(δΔE) | (2) |
ΔES(t) = ∑fnSn(t) = ∑fnSn(0)cn(t) | (3) |
The n-th order contributes fnSn(0) to the total Stokes shift with the dynamics characterized by cn(t). Giving Gaussian statistics of the solute–solvent interaction on the excited state surface, Wick's theorem enables the even order S2n(t) to disappear, and all the odd S2n+1(t) to exhibit the same decay rates by25,41
c2n+1(t) = c1(t) | (4) |
Therefore, the Stokes shift merely arises from the odd order time correlation function with the contribution f2n+1S2n+1(0) from the order 2n + 1. Both f2n+1 and S2n+1(0) can be evaluated by the Gaussian type weight function , in which
is the excited state Gaussian distribution with the width σ. According to eqn (4), the solvation time correlation function remains consistent with the odd order time correlation functions, namely the linear one by
S(t) = c1(t) | (5) |
The above derivations validate the linear response theory of Stokes shift with the excited state Gaussian fluctuations.
We further examined the decay rates of the odd order S2n+1(t). By equivalently accumulating 〈δΔE(t)δΔE(0)2n+1〉e over the excited state equilibrium, S2n+1(t) was derived for the orders 1, 3, 5, 7 and 9. Ideally, one would like to acquire various order time correlation functions, as many as possible, to examine the dynamics of the spontaneous fluctuations. However, the accurate evaluation of S2n+1(t) requests abundant samplings on those δΔE at and around the stationary points of the weight function, which is obviously difficult for the higher orders. Hence, we merely calculated the time correlation functions up to the 9-th order. The normalized form c2n+1(t) was calculated and is displayed in Fig. 2. All the curves surprisingly exhibit almost the same decay rates. The relaxation processes are characterized by an inertial decay on the femtosecond time scale, and the dynamics on a few picoseconds correlated to the rotational/diffusional motion of the water molecules. The identical relaxation processes of these orders verify the Gaussian statistics of the energy fluctuations on the excited state surface. Its critical role in the non-equilibrium relaxations was further examined in the following context.
![]() | ||
Fig. 2 Normalized time correlation functions c2n+1(t) of the odd orders. The decay rates are the same up to the 9-th order. |
Non-equilibrium simulations were performed by propagating trajectories on the excited-state surface with the initial configurations sampled from the equilibrium fluctuations of the states g, gT1, gq1 and gq2. We firstly inspected the total energy shifts by applying ΔES(0) = ∫d(δΔE)δΔEPI(δΔE), which are 16.0, 16.0, 10.8 and 46 for the above non-equilibrium processes, respectively. By further projecting the initial distribution PI(δΔE) into the Gaussian type weight functions (eqn (2)), ΔES(0) was decomposed into the contributions of the odd order time correlation functions by f2n+1S2n+1(0) (Fig. 4). The peak energy of PI(δΔE), namely the center of the profile, notifies the order of the weight function making the maximum contribution. For example, the center of PI(δΔE) is 10.8 for the state gq1, residing closely to the stationary point (10.1) of the 7-th weight function. Correspondingly, f7S7(0) of the 7th order contributes maximally among all the time correlation functions. The same correspondences are observed for the states g, gT1, and gq2. For both the states g and gT1, the maximum contribution to the energy decomposition arises from the 17-th order, which is attributed to the close residence between the center (16.0) of their PI(δΔE) and the stationary point (15.7) of the 17-th weight function. However, about 80% of the contributions to ΔES(0) arise from the order 9–25 for the state g, and from the order 7–29 for the state gT1. The broader distribution PI(δΔE) of the latter enables it to correlate with more weight functions than the ground state ones. For the state gq2, the orders of 121–163 contribute 82% to the total energy shift. The maximum contribution is from the 145-th order. The associated stationary point (45.8) is close to the center (46.0) of the initial distribution PI(δΔE).
The time relaxations of ΔES(t) were obtained for the non-equilibrium processes, in which the inital configurations were sampled from the states g, gq1, gT1 and gq2 respectively (Fig. 5). The normalized solvation time correlation functions S(t) were compared with the excited state time correlation functions c1(t) to check the validity of linear response theory. As a result of Gaussian statistics, Stokes shift is represented as the summation of all the odd order time correlation functions by
ΔES(t) = f2n+1S2n+1(0)c2n+1(t) | (6) |
The energy allocations of the total decay ΔES(0) are shown in Fig. 4. For the states g, gq1, gT1, those orders taking the most dominant contributions are the lower ones. Their equilibrium time correlation functions c2n+1(t) obey the same dynamics as the linear one c1(t) owing to Gaussian statistics of the energy variable δΔE. Therefore, the non-equilibrium solvation time correlation S(t) remains consistent with c1(t) as expected by linear response theory. In contrast, the non-equilibrium dynamics of the state gq2 are evidently different from c1(t). In this case, most of the energy contributions are from the higher orders. For example, the maximum contribution is from the 145-th order. For those very high orders, the equilibrium time correlation functions of c2n+1(t) can hardly be accumulated due to sampling difficulties. We speculate that their dynamics are different from the linear order c1(t) due to the deviation of the Gaussian statistics, causing discrepancy between the non-equilibrium dynamics and the linear response expectation. Therefore, both the success and failure of linear response theory were observed and the mechanisms were interpreted. We also compared the simulation with the experimental results, in which the relaxation ΔES(t) of the state g corresponds to the experimental measurements of Stokes shift. The simulated solvation time correlation function of the state g is fitted with the function , leading to the inertial decay τ1 of 170 femtosecond and τ2 of 1.0 picosecond. The simulation results are in good agreement with the experimental measurements, which shows the dynamics of 1.6 picoseconds upon the initial ultrafast decay for the tryptophan solution.45 The molecular mechanism for the biphasic relaxation is shown in Fig. 5 (bottom panel). The time scale τ1 is mostly attributed to the reorientation motion of the neighboring water molecules of indole, and τ2 is linked with the rotational and diffusional motions of water.
For the non-equilibrium solvation time correlation function S(t) of the state gq2, we further probed the molecular origin for the failure of linear response theory. The solvation structures in the proximity of indole were examined for gq2 and the excited state, which provided the most relevant information on the initial and final configuration over the non-equilibrium process. The artificial charge modulation on the state gq2 makes a significant change in its dipole moment compared with the excited state, which is 9.7 D and 5.1 D, respectively, with an angle flip of 137° (Fig. S2 in ESI†). As a result, the solvation structure in proximity to indole was significantly rearranged. Fig. 6 displays the typical snapshots for the two states. The average number of water molecules within 3 Å of the nitrogen atom is 2.3 and 1.1 for the state gq2 and the excited state, respectively. Therefore, the H-bond structure between the solute and the solvents vary significantly over non-equilibrium relaxation. This gross change in the solvation structure, as reported in literatures,24,46,47 correlates with the discrepancy between the non-equilibrium and equilibrium dynamics and accounts for the breakdown of linear response theory.
The mechanism of linear response theory has been carefully examined in literature. Heid and co-workers found that a large intermediate broadening for the width of the energy distribution correlates with a failure of the time correlation function to describe the non-equilibrium event.27 Our observations coincide with their conclusion. Among the non-equilibrium processes, the only case with the failure of linear response theory shows the strong peculiarity of width dynamics, namely intermediate broadening, and all the others with the success of linear response theory show either relatively stable or gradual narrowing of width. Therefore, the abrupt change, and not the absolute change, of the width links to the failure of linear response theory. How this specificity of the width dynamics inherently correlates with the non-Gaussian statistics, as assigned to be the origin of nonlinear behavior in this study, needs further investigations and illumination.
Footnote |
† Electronic supplementary information (ESI) available. See DOI: 10.1039/d0ra01227k |
This journal is © The Royal Society of Chemistry 2020 |