Stochastic and Quasi-Stochastic Hamiltonians for Long-Time

Oct 6, 2017 - This property is utilized to formulate novel stochastic and quasi-stochastic vibronic Hamiltonian methodologies, which open a new route ...
0 downloads 9 Views 3MB Size
Letter Cite This: J. Phys. Chem. Lett. 2017, 8, 5190-5195

pubs.acs.org/JPCL

Stochastic and Quasi-Stochastic Hamiltonians for Long-Time Nonadiabatic Molecular Dynamics Alexey V. Akimov* Department of Chemistry, University at Buffalo, The State University of New York, Buffalo, New York 14260, United States S Supporting Information *

ABSTRACT: In the condensed-matter environments, the vibronic Hamiltonian that describes nonadiabatic dynamics often appears as an erratic entity, and one may assume it can be generated stochastically. This property is utilized to formulate novel stochastic and quasi-stochastic vibronic Hamiltonian methodologies, which open a new route to long-time excited state dynamics in atomistic solid-state systems at negligible computational cost. Using a model mimicking a typical solid-state material in noisy environment, general conclusions regarding the simulation of nonadiabatic dynamics are obtained: (1) including bath is critical to complete excited state relaxation; (2) a totally stochastic modulation of energies and couplings has a net effect of no bath and inhibits relaxation; (3) including a single or several dominant electron−phonon modes may be insufficient to complete the excited state relaxation; (4) only the multiple modes, even those that have negligible weights, can represent both the deterministic modulation of system’s Hamiltonian and stochastic effects of bath. applied, the vibronic Hamiltonian, Hvib ij = Eiδij − iℏdij, behaves quite erratically: the energies of electronic states, Ei, and the scalar (time-derivative) NACs, dij, fluctuate randomly around some average values (e.g., see Figure 6b in ref 38.). It is appears nonoptimal to execute the time-critical electronic structure calculations to get such quasi-random numbers at every point of the MD trajectory. A somewhat more efficient strategy would be to obtain a descriptive statistics of each involved number (energies and couplings) first and use this information to generate the matrix elements of the vibronic Hamiltonian randomly, at the negligible computational costs. The descriptive statistics would consist of the mean values, standard deviations, and, perhaps, higher moments (skewness, curtosis, etc.). In this way, the vibronic Hamiltonian could be constructed stochastically (hence the name) and efficiently despite the length of the NA-MD trajectories. A stochastic potential methodology has been recently proposed by the Chakraborty group [Scher, J. A.; Srihari, A.; Bayne, M. G.; Nangia, S.; Chakraborty, A. “Development of Effective Stochastic Potential Method Using Random Matrix Theory for Embedding Quantum Subystems in Noisy Environment” (under review), 2017], who applied it for fast computations of the absorption spectra in large systems. Stimulated by this idea, I explore and extend it in the context of the NA-MD simulations. In the method of Chakraborty et al., only the extensive sampling of the potential is important, since the ensemble-averaged properties are of interest. In the dynamical simulations, like those undertaken in this study,

G

rowing interests of many researchers in solar energy harvesting and conversion materials have stimulated the development of multiple nonadiabatic molecular dynamics (NA-MD) methodologies for efficient computational modeling of excited states dynamics.1−6 The NA-MD simulations can provide intricate details about the charge carrier and excitation energy transfer mechanisms, nonradiative and radiative excited state relaxation dynamics, coherence loss channels, and the interplay of the above. The computational expenses of the NAMD grow rapidly with the size of the system and the number of excited states to be accounted for, forcing one to develop and use simplified theories, such as the one-electron description of excited states,1,6−11 classical description of nuclei (the so-called classical path approximation, CPA),12 neglect of the electron− phonon back-reaction approximation (NBRA),13,14 and many others.3,15−24 Despite their simplicity, these theories have been proven successful for many classes of systems, especially for solid-state materials,25−30 quantum dots,31−36 and high-temperature states of matter.37 In such environments, the NBR and one-electron approximations are well-justified by the relative rigidity of the system and the presence of many electrons, respectively. Even with the existing approximations, the NA-MD simulations remain limited by the expenses of explicit electronic structure calculations, which restrict the time scale of the processes that can be modeled directly to dozen of picoseconds at most. Simulating the processes that occur concurrently on the hundreds of picoseconds or nanosecond time scales would be technically challenging and would require either an extrapolation of the short-term dynamics or an ad hoc scaling of nonadiabatic couplings (NACs). At the same time, one can notice that in the structurally rigid environments, where the NA-MD with the NBRA can be © 2017 American Chemical Society

Received: August 18, 2017 Accepted: October 6, 2017 Published: October 6, 2017 5190

DOI: 10.1021/acs.jpclett.7b02185 J. Phys. Chem. Lett. 2017, 8, 5190−5195

Letter

The Journal of Physical Chemistry Letters

Figure 1. Comparison of the NBRA-based NA-MD simulations with and without bath. (a) Excited state population dynamics; (b) Trajectories in the phase space (q−coordinate, p - momentum); (c) Evolution of the energies of the ground and excited states; (d) Evolution of the kinetic, potential, total, and extended energy.

in Figure 1 indicate that the thermostat facilitates the excited state relaxation dynamics, especially in the low-dimensional systems as the model used here. The simulations with the thermostat coupled to the nuclear dynamics for the whole duration leads to a complete decay of the excited state population and erratic phase space trajectory. One can conclude that the bath-induced fluctuations assist the system in finding efficient relaxation pathways and their chaotic nature prevents the re-excitations. A dramatic change of the dynamics can be observed if the thermostat is decoupled from the system after the initial 2 ps. The excited state population does not change and remains at the level of 0.5 (Figure 1a). The phase space trajectory becomes regular. The kinetic energy is sufficient for the system to overcome all the ground state PES barriers thermally, so that the trajectory corresponds to the particle moving to the negative direction with the momentum fluctuating regularly (Figure 1b). The PES image along the trajectory becomes regular (Figure 1c). In fact, since the motion is linear and the trajectory coincides with the PES scan coordinate, Figure 1c (after 2 ps) also illustrates the PES profile (scaled in the coordinate direction). Note that the particle crosses the regions with small (large gap) and large NAC magnitudes. However, because of such a quasi-ballistic motion, the system does not spend enough time in the regions of large couplings to lead to significant population transfer probabilities. When thermostat is present, the system can spend longer time at various points due to the effective friction exerted by thermostat. Finally, the regularity of the nuclear dynamics facilitates the reversal of the electronic dynamics and reexcitation. The correctness of the simulation algorithm is confirmed by the conservation of the corresponding extended energy (system + bath) for all times (Figure 1d). The extended energy coincides with the total energy after 2 ps, when thermostat is decoupled from the system. The total energy shows large

the details of the Hamiltonian’s time-dependence is critical, because the latter determines the evolution of the effective electronic coordinates. It will be demonstrated here that generating the Hamiltonian stochastically does not work in NAMD. Instead, one needs a deterministic vibronic Hamiltonian that also mimics the stochasticity introduced by the environment. In the present work, such an approach is termed a quasistochastic Hamiltonian method. Along the way of exploring the two (stochastic and quasi-stochastic) methodologies, some general notes pertaining to the simulations of excited states dynamics in arbitrary systems will be revealed. In this work, a 1D (nuclear degree of freedom, DOF), twolevel system introduced earlier39 is considered (see the SI for further details). The model features the periodic potential energy surfaces (PES) and NACs, mimicking their behavior in realistic solid-state materials. The NA-MD is then performed using the Tully’s surface hopping scheme within the NBRA, that is, the nuclear dynamics is always driven by the groundstate PES (further referred to as the direct simulation) or is fully frozen (in the stochastic or quasi-stochastic simulations). The PES and the initial parameters are chosen such that the classical particle could cross the ground state PES barriers at the temperature of simulation, so that the NBRA is accurate. The simulations are performed using the Libra methodology development library.39 The details of simulations are available in the SI. In this work, the classical DOF can be optionally coupled to a thermostat (Nose-Hoover chain), which together with the periodicity of the PES can mimic the stochasticity of the environment (bath) into which the quantum system is immersed. The unharmonicity and periodicity of the potential are meant to mimic the potential of typical periodic solids, where the NA-MD with the NBRA is often applied. The role of a thermostat in NA-MD calculations has been a subject of interest lately39 and is further investigated in this work using the present model Hamiltonian. The simulations summarized 5191

DOI: 10.1021/acs.jpclett.7b02185 J. Phys. Chem. Lett. 2017, 8, 5190−5195

Letter

The Journal of Physical Chemistry Letters fluctuations during the first 2 ps when the thermostat is in effect, whereas the extended energy is still conserved. The initial 2 ps of the direct simulations with the thermostat are then used to compute the averages ⟨Hij⟩ and variances, ⟨(Hij − ⟨Hij⟩)2⟩ of all matrix elements of the vibronic Hamiltonian. They are used to generate the stochastic vibronic Hamiltonian for the rest of the simulation (e.g., up to 20 ps): Hij(t) = ⟨Hij⟩ + ⟨(Hij − ⟨Hij⟩)2⟩1/2ζt. Here, ζ is a random number that can be chosen in a number of ways: (a) is distributed normally (mean 0, variance 1) or uniformly (on [−1,1]); (b) can be uncorrelated or correlated. The correlation of the random numbers is realized via an algorithm reported by Deserno:40 ζ0

situation when two adiabatic PES come close to each other and allow thermal excitation. Consequently, the excited state repopulation is minimized, as can be seen in Figure 3b. The saturation of the equilibrium population level becomes more apparent. The simulations with all four variants of construction of the stochastic Hamiltonian show similar results. The results in Figure 3 demonstrate that the fully stochastic Hamiltonians are not generally applicable to the dynamical simulations, even when the cropping is applied. The totally stochastic Hamiltonians have worked for the computations of the absorption spectra [Chakraborty et al. under review] as the statistical sampling tools, but only because the details of the dynamics were not important for such type of calculations. This simulation also strongly suggests that the details of the temporal dependence of the vibronic Hamiltonian are important, even if the latter appears to be stochastic. The vibrational modes and their interplay constitute such details. It is generally known that the nuclei drive the electronic transitions, but it may be easy to get confused by the complex internal dynamics of such modes and be tempted to disregard such internal information by using the totally stochastic sampling. The next step is to reintroduce the existing nuclear modes in the dynamics and sample the vibronic Hamiltonians quasistochastically. According to the quasi-stochastic Hamiltonian formulation, the matrix elements of the vibronic Hamiltonian are computed as Hij(t) = ⟨Hij⟩ + 1/s ⟨(Hij − ⟨Hij⟩)2⟩1/2 ∑Nk=1aijk sin(2πωijkt), where s is the standard deviation of the ∑Nk=1aijk sin(2πωijkt) function, ωijk and aijk are the kth frequencies and their relative amplitudes for the matrix element Hij. These properties are computed by taking the Fourier transform of the autocorrelation function of such matrix elements. Finally, N is the number of the most intense frequencies to be included in the calculations. Similar to the stochastic Hamiltonian formulation, the cropping transformation can be applied if needed. The results of the quasi-stochastic Hamiltonian simulations are summarized in Figure 4. One can observe a drastic qualitative difference in the dynamics, in comparison to the stochastic formulation. First of all, it becomes possible to correctly reproduce the direct dynamics of the excited state (Figure 4b). However, the two conditions must be met: (a) a sufficient number of the nuclear modes must be included; (b) the windowing transformation should be applied. Using a single nuclear mode (even though the most dominant) leads to an early saturation or slow decay. This is because the single vibrational mode can not capture the bath modes and leads to a regular electronic dynamics, similar to the one observed in the direct simulations without the bath (e.g., see Figure 1). Using a somewhat larger number of nuclear modes (2−10) can induce coherent oscillations of the population. If the cropping is not applied, the oscillation may be not decaying because the thermal re-excitation has high probability (Figure 4a). Once the cropping is applied, the oscillations have a smaller amplitude and show signs of a slow decay (Figure 4b). The presence of several driving modes provides more opportunities for the system to sample more efficient decay pathways. On the other hand, the limited number of the modes makes the dynamics quite regular. The temporal symmetry of the electronic transitions is not fully broken, leading to the oscillations. Finally, when the number of the nuclear modes is increased over a certain threshold, the complete relaxation of the excited state can be achieved (Figure 4, panels a and b). In the present

= ξ0, ζi + 1 = fζi + 1 − f 2 ξi + 1, where f = e−1/τ, with τ being the correlation time. The correlation time is found from the dominant frequency in the spectra of the corresponding matrix elements. For the present system, the typical dominant frequency is 141 cm−1 for energies and 220 cm−1 for NAC (Figure 2). The details of spectra calculations can be found in the SI.

Figure 2. Fourier transform of the autocorrelation function of the three distinct matrix elements of the vibronic Hamiltonian: energies of the states (real diagonal components) and the nonadiabatic coupling (imaginary off-diagonal component).

The results of the NA-MD simulations with four versions of stochastic Hamiltonian are summarized in Figure 3 (panels a and c). The first 2 ps are obtained directly and illustrate the behavior of the reference Hamiltonian (energies and NAC) and the transitions it induces. One can see that the stochastically constructed Hamiltonian shows much higher fluctuation frequencies, and sometimes the amplitude may exceed the maximal amplitude of the reference Hamiltonian (Figure 3c). The large fluctuations facilitate thermal excitation of the system to the excited state (when the adiabatic PESs get close, so the Boltzmann rejection unlikely). The high oscillation frequency minimizes the bath effects, making the Hamiltonian effectively deterministic (as if without the bath). Correspondingly, the vanishing bath effects make population dynamics similar to the one observed in the reference simulations without the bath (Figure 1a): the excited state population reaches the level somewhat below 0.5 and slowly decays further. However, the population decay shows the signs of saturation and the longtime repopulation of the excited states. To correct the magnitude of the stochastic Hamiltonian fluctuation, the windowing transformation has been applied. The upper and lower bound of each matrix element are computed during the initial 2 ps sampling. When the stochastic Hamiltonian exceeds these boundaries, it is reset to the closest boundary value (Figure 3d). This cropping technique is especially important for the energies, since it helps avoid the 5192

DOI: 10.1021/acs.jpclett.7b02185 J. Phys. Chem. Lett. 2017, 8, 5190−5195

Letter

The Journal of Physical Chemistry Letters

Figure 3. Evolution of the excited state population (a,b) and matrix elements of the vibronic Hamiltonian (c,d) as computed initially with the direct NA-MD method (first 2 ps) and continued with the stochastic Hamiltonian NA-MD (beyond 2 ps). The left column (panels a,c) reports the results without the cropping technique; the right column (panels b, d) reports the results for the analogous simulations, but with inclusion of the cropping technique.

Figure 4. Evolution of the excited state population (a,b) and matrix elements of the vibronic Hamiltonian (c,d) as computed initially with the direct NA-MD method (first 2 ps) and continued with the quasi-stochastic Hamiltonian NA-MD (beyond 2 ps). The left column (a,c) reports the results without the cropping technique; the right column (b,d) reports the results for the analogous simulations, but with inclusion of the cropping technique. The line legends correspond to the number of vibrational modes included in the calculations of quasi-stochastic Hamiltonian.

simulations, this number is estimated to be 25. It is likely that the added extra modes (of the decreasing weights) mimic the bath and lead to longer recurrence times (break temporal symmetry), so that the re-excitation of the system would occur on the time scale much longer than the trajectory length. At this point, there is no general prescription on how many modes should suffice for the complete relaxation. However, having all

the modes (even with small intensities) present in the FT spectrum (Figure 2) should always deliver a desired complete relaxation. At the same time, in the multidimensional systems a limited number of modes may be sufficient since they can act as a bath with respect to all other modes. It is expected that a reasonable relaxation dynamics may be properly described as 5193

DOI: 10.1021/acs.jpclett.7b02185 J. Phys. Chem. Lett. 2017, 8, 5190−5195

Letter

The Journal of Physical Chemistry Letters

based on the modes existing in the reference system (including bath modes) allows one to accurate model the relaxation dynamics. It is found that the inclusion of a small number of the dominant modes alone may not be sufficient, since the recurrences and regular dynamics may appear, inhibiting the electronic transitions. The presence of multiple modes (even nondominant) and energy cropping are essential for accurate simulations. These factors ensure the absence of the short-term recurrences and thermal re-excitation of the system. The modes introduce a quasi-stochasticity in the evolution of the Hamiltonian. The resulting quasi-stochastic Hamiltonian methodology can be useful for accelerating the atomistic ab initio based NA-MD simulations in the condensed-matter systems in which the NBR approximation works well.

long as the retained modes deliver enough complexity for the system’s dynamics to have sufficiently long revival times. One can also note that the cropping transformation helps avoid thermal re-excitation (e.g., N = 25−200 in Figure 4b), but slows the relaxation dynamics a bit (e.g., compare with the analogous line in Figure 4a). Thus, a careful windowing of the quasi-stochastically generated Hamiltonian matrix elements may deliver the best compromise between the accuracy of the computed time scales and the longer-time effects. It is worth mentioning that the time evolution of all matrix elements of the quasi-stochastic vibronic Hamiltonian (Figure 4, panels c and d, after 2 ps) closely resembles that of the originally sampled ones (before 2 ps). One can visually identify several frequencies, which are comparable to the frequencies of the oscillations of the reference matrix elements. Even periodic spikes in the evolution of the NAC are visible. This is in clear contrast to the very dense, high-frequency oscillations observed in the fully stochastic formulations. The explicit consideration of the nuclear modes in this work suggests that the quantization of the nuclear DOFs can be included. For instance, the frequencies obtained in the Fourier transform (Figure 2) of the energy ACF can be used to construct a set of harmonic oscillator energy levels and modulate the dynamics by hopping between these states stochastically. However, the present methodology focuses on the presets of the NBRA (including the classical nuclei approximation) and on the conditions under which the NBRA works. This is usually the case when the nuclear kinetic energy is comparable or larger than the energy fluctuation magnitudes, such that the system is not confined within a single minimum valley, but can visit many of them. Under the hightemperature/small energy barrier approximation, important to NBRA, the dynamics of nuclei is well-described classically. One may envision quantizing the nuclear modes at lower temperatures, but the NBRA may be not work as well as it works in the high temperature regime. The success of the quasi-stochastic Hamiltonian approach emphasizes the importance of including nuclear modes in the NA-MD simulations explicitly. Moreover, the “bath modes” must be included for a proper description of the decay dynamics. These modes, although not dominant and still deterministic, mimic the stochasticity of the bath and introduce the complexity, which helps avoid recurrences (revivals) in the dynamics, and provides the means of effective energy dissipation that in needed for complete population decay. The use of a chain of thermostats coupled to each other with the leading one coupled to the system makes the systems’ dynamics difficult to reverse. In fact, it would be nearly impossible for the system to return to the initial state in the cause of its evolution (Nose−Hoover chain thermostat is known to make the dynamics of nonergodic systems ergodic41). It should be reiterated, that another source of irreversibility in NBRA-based simulations is the constant energy dissipation, which is implicitly encoded in the Boltzmann factor-based hop rejection algorithm (see more in the SI). To recapitulate, in this work, the stochastic and quasistochastic Hamiltonian approaches to NA-MD simulations have been described and studied. It is demonstrated that the fully stochastic Hamiltonians can not properly capture the excited state relaxation, even if the stochastic components are correlated. The random fluctuations of high frequency cancel each other, making the system effectively isolated. On the contrary, the deterministic construction of the Hamiltonian



ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpclett.7b02185. Details of the model Hamiltonian, methodologies, and computational setups (PDF)



AUTHOR INFORMATION

Corresponding Author

*E-mail: alexeyak@buffalo.edu. Twitter: @AkimovLab. ORCID

Alexey V. Akimov: 0000-0002-7815-3731 Notes

The author declares no competing financial interest.



ACKNOWLEDGMENTS The author acknowledges the financial support from the University at Buffalo, The State University of New York startup package. Computational support is provided by the Center for Computational Research at the University at Buffalo.



REFERENCES

(1) Oberhofer, H.; Reuter, K.; Blumberger, J. Charge Transport in Molecular Materials: An Assessment of Computational Methods. Chem. Rev. 2017, 117, 10319−10357. (2) Wang, L.; Akimov, A.; Prezhdo, O. V. Recent Progress in Surface Hopping: 2011−2015. J. Phys. Chem. Lett. 2016, 7, 2100−2112. (3) Wang, L.; Sifain, A. E.; Prezhdo, O. V. Fewest Switches Surface Hopping in Liouville Space. J. Phys. Chem. Lett. 2015, 6, 3827−3833. (4) Kilina, S.; Kilin, D.; Tretiak, S. Light-Driven and Phonon-Assisted Dynamics in Organic and Semiconductor Nanostructures. Chem. Rev. 2015, 115, 5929−5978. (5) Nelson, T.; Fernandez-Alberti, S.; Roitberg, A. E.; Tretiak, S. Nonadiabatic Excited-State Molecular Dynamics: Modeling Photophysics in Organic Conjugated Materials. Acc. Chem. Res. 2014, 47, 1155−1164. (6) Barbatti, M.; Crespo-Otero, R. Surface Hopping Dynamics with DFT Excited States. In Density-Functional Methods for Excited States; Springer International Publishing: Cham, Switzerland, 2014; Vol. 368, pp 415−444. (7) Blumberger, J. Recent Advances in the Theory and Molecular Simulation of Biological Electron Transfer Reactions. Chem. Rev. 2015, 115, 11191−11238. (8) Akimov, A. V.; Prezhdo, O. V. The PYXAID Program for NonAdiabatic Molecular Dynamics in Condensed Matter Systems. J. Chem. Theory Comput. 2013, 9, 4959−4972. (9) Akimov, A. V.; Prezhdo, O. V. Advanced Capabilities of the PYXAID Program: Integration Schemes, Decoherence Effects, Multi5194

DOI: 10.1021/acs.jpclett.7b02185 J. Phys. Chem. Lett. 2017, 8, 5190−5195

Letter

The Journal of Physical Chemistry Letters excitonic States, and Field-Matter Interaction. J. Chem. Theory Comput. 2014, 10, 789−804. (10) Rego, L. G. C.; Hames, B. C.; Mazon, K. T.; Joswig, J.-O. Intramolecular Polarization Induces Electron−Hole Charge Separation in Light-Harvesting Molecular Triads. J. Phys. Chem. C 2014, 118, 126−134. (11) Casalegno, M.; Pastore, R.; Idé, J.; Po, R.; Raos, G. Origin of Charge Separation at Organic Photovoltaic Heterojunctions: A Mesoscale Quantum Mechanical View. J. Phys. Chem. C 2017, 121, 16693−16701. (12) Tully, J. C. Molecular Dynamics with Electronic Transitions. J. Chem. Phys. 1990, 93, 1061−1071. (13) Craig, C.; Duncan, W.; Prezhdo, O. Trajectory Surface Hopping in the Time-Dependent Kohn-Sham Approach for Electron-Nuclear Dynamics. Phys. Rev. Lett. 2005, 95, 163001. (14) Fischer, S. A.; Habenicht, B. F.; Madrid, A. B.; Duncan, W. R.; Prezhdo, O. V.; et al. Regarding the Validity of the Time-Dependent Kohn-Sham Approach for Electron-Nuclear Dynamics via Trajectory Surface Hopping. J. Chem. Phys. 2011, 134, 024102. (15) Subotnik, J. E. Augmented Ehrenfest Dynamics Yields a Rate for Surface Hopping. J. Chem. Phys. 2010, 132, 134112. (16) Shenvi, N.; Subotnik, J. E.; Yang, W. Phase-Corrected Surface Hopping: Correcting the Phase Evolution of the Electronic Wavefunction. J. Chem. Phys. 2011, 135, 024101. (17) Jaeger, H. M.; Fischer, S.; Prezhdo, O. V. Decoherence-Induced Surface Hopping. J. Chem. Phys. 2012, 137, 22A545. (18) Akimov, A. V.; Long, R.; Prezhdo, O. V. Coherence Penalty Functional: A Simple Method for Adding Decoherence in Ehrenfest Dynamics. J. Chem. Phys. 2014, 140, 194107. (19) Akimov, A. V.; Prezhdo, O. V. Second-Quantized Surface Hopping. Phys. Rev. Lett. 2014, 113, 153003. (20) Meek, G. A.; Levine, B. G. Wave Function Continuity and the Diagonal Born-Oppenheimer Correction at Conical Intersections. J. Chem. Phys. 2016, 144, 184109. (21) Meek, G. A.; Levine, B. G. Accurate and Efficient Evaluation of Transition Probabilities at Unavoided Crossings in Ab Initio Multiple Spawning. Chem. Phys. 2015, 460, 117−124. (22) Sifain, A. E.; Wang, L.; Prezhdo, O. V. Communication: Proper Treatment of Classically Forbidden Electronic Transitions Significantly Improves Detailed Balance in Surface Hopping. J. Chem. Phys. 2016, 144, 211102. (23) Akimov, A. V. Nonadiabatic Molecular Dynamics with TightBinding Fragment Molecular Orbitals. J. Chem. Theory Comput. 2016, 12, 5719−5736. (24) da Silva Oliboni, R.; Bortolini, G.; Torres, A.; Rego, L. G. C. A Nonadiabatic Excited State Molecular Mechanics/Extended Hückel Ehrenfest Method. J. Phys. Chem. C 2016, 120, 27688−27698. (25) Jankowska, J.; Prezhdo, O. V. Ferroelectric Alignment of Organic Cations Inhibits Nonradiative Electron−Hole Recombination in Hybrid Perovskites: Ab Initio Nonadiabatic Molecular Dynamics. J. Phys. Chem. Lett. 2017, 8, 812−818. (26) Long, R.; Prezhdo, O. V. Time-Domain Ab Initio Modeling of Electron−Phonon Relaxation in High-Temperature Cuprate Superconductors. J. Phys. Chem. Lett. 2017, 8, 193−198. (27) Long, R.; Fang, W.; Prezhdo, O. V. Moderate Humidity Delays Electron−Hole Recombination in Hybrid Organic−Inorganic Perovskites: Time-Domain Ab Initio Simulations Rationalize Experiments. J. Phys. Chem. Lett. 2016, 7, 3215−3222. (28) Long, R.; Fang, W.; Akimov, A. V. Nonradiative Electron−Hole Recombination Rate Is Greatly Reduced by Defects in Monolayer Black Phosphorus: Ab Initio Time Domain Study. J. Phys. Chem. Lett. 2016, 7, 653−659. (29) Madjet, M. E.-A.; Akimov, A. V.; El-Mellouhi, F.; Berdiyorov, G. R.; Ashhab, S.; Tabet, N.; Kais, S. Enhancing the Carrier Thermalization Time in Organometallic Perovskites by Halide Mixing. Phys. Chem. Chem. Phys. 2016, 18, 5219−5231. (30) Nijamudheen, A.; Akimov, A. V. Excited State Dynamics in 2D Heterostructures: SiR/TiO2 and GeR/TiO2 (R = H, Me) as Promising Photocatalysts. J. Phys. Chem. C 2017, 121, 6520−6532.

(31) Kilina, S. V.; Craig, C. F.; Kilin, D. S.; Prezhdo, O. V. Ab Initio Time-Domain Study of Phonon-Assisted Relaxation of Charge Carriers in a PbSe Quantum Dot. J. Phys. Chem. C 2007, 111, 4871−4878. (32) Kilina, S. V.; Kilin, D. S.; Prezhdo, O. V. Breaking the Phonon Bottleneck in PbSe and CdSe Quantum Dots: Time-Domain Density Functional Theory of Charge Carrier Relaxation. ACS Nano 2009, 3, 93−99. (33) Kilina, S. V.; Neukirch, A. J.; Habenicht, B. F.; Kilin, D. S.; Prezhdo, O. V. Quantum Zeno Effect Rationalizes the Phonon Bottleneck in Semiconductor Quantum Dots. Phys. Rev. Lett. 2013, 110, 180404. (34) Senanayake, R. D.; Akimov, A. V.; Aikens, C. M. Theoretical Investigation of Electron and Nuclear Dynamics in the [Au25(SH)18]−1 Thiolate-Protected Gold Nanocluster. J. Phys. Chem. C 2017, 121, 10653−10662. (35) Dong, S.; Pal, S.; Lian, J.; Chan, Y.; Prezhdo, O. V.; Loh, Z.-H. Sub-Picosecond Auger-Mediated Hole-Trapping Dynamics in Colloidal CdSe/CdS Core/Shell Nanoplatelets. ACS Nano 2016, 10, 9370− 9378. (36) Lin, Y.; Akimov, A. V. Dependence of Nonadiabatic Couplings with Kohn−Sham Orbitals on the Choice of Density Functional: Pure vs Hybrid. J. Phys. Chem. A 2016, 120, 9028−9041. (37) Pradhan, E.; Magyar, R.; Akimov, A. Scaling Relationships for Nonadiabatic Energy Relaxation Times in Warm Dense Matter: Toward Understanding the Equation of State. Phys. Chem. Chem. Phys. 2016, 18, 32466−32476. (38) Akimov, A. V.; Prezhdo, O. V. Nonadiabatic Dynamics of Charge Transfer and Singlet Fission at the Pentacene/C60 Interface. J. Am. Chem. Soc. 2014, 136, 1599−1608. (39) Akimov, A. V. Libra: An Open-Source “methodology Discovery” Library for Quantum and Classical Dynamics Simulations. J. Comput. Chem. 2016, 37, 1626−1649. (40) Deserno, M. How to Generate Exponentially Correlated Gaussian Random Numbers; University of California, Los Angeles, 2015. http:// www.cmu.edu/biolphys/deserno/pdf/corr_gaussian_random.pdf. (41) Tuckerman, M. E.; Liu, Y.; Ciccotti, G.; Martyna, G. J. NonHamiltonian Molecular Dynamics: Generalizing Hamiltonian Phase Space Principles to Non-Hamiltonian Systems. J. Chem. Phys. 2001, 115, 1678−1702.

5195

DOI: 10.1021/acs.jpclett.7b02185 J. Phys. Chem. Lett. 2017, 8, 5190−5195