A Simple Phase Correction Makes a Big Difference in Nonadiabatic

Subscriber access provided by University of Sunderland

Spectroscopy and Photochemistry; General Theory

A Simple Phase Correction Makes a Big Difference in Nonadiabatic Molecular Dynamics Alexey V. Akimov J. Phys. Chem. Lett., Just Accepted Manuscript • DOI: 10.1021/acs.jpclett.8b02826 • Publication Date (Web): 04 Oct 2018 Downloaded from http://pubs.acs.org on October 7, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry Letters

A Simple Phase Correction Makes a Big Difference in Nonadiabatic Molecular Dynamics Alexey V. Akimov*

Department of Chemistry, University at Buffalo, The State University of New York, Buffalo, NY, 14260 USA.

ACS Paragon Plus Environment

1

The Journal of Physical Chemistry Letters 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 30

Abstract The outcomes of nonadiabatic molecular dynamics (NA-MD) calculations are modulated by the parameters entering the time-dependent Schrodinger equation (TDSE). The adiabatic states are commonly used as the basis in which the TD-SE is integrated. However, the phase inconsistencies of such states along the nuclear trajectories obtained in NA-MD simulations may render the wavefunction and other relevant properties ill-behaving, adversely affecting the dynamics. This work illustrates the consequence of adiabatic state phase inconsistencies in nonadiabatic Ehrenfest dynamics. A simple phase correction approach is proposed and is demonstrated to alter the dynamics to make it consistent with the reference calculations done in the phaseconsistent diabatic representation.


2

Page 3 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


TOC GRAPHICS

Nonadiabatic molecular dynamics (NA-MD) has become a popular tool to study charge and energy transfer processes in many solar energy materials,1–7 biological systems,8,9 and even exotic states of matter.10 The quantum-classical trajectory surface hopping (TSH) algorithm of Tully11 is among the most commonly used approximations in a wide spectrum of available computational techniques, thanks to its computational efficiency, satisfactory accuracy, and simplicity of its implementation. The computational scheme


3


Page 4 of 30

involves solving the time-dependent Schrodinger equation (TD-SE) in the basis of electronic states, |𝜓⟩ = (𝜓1,𝜓2,…, 𝜓𝑁), that are parametrically dependent on time via the motion of nuclei (classical path approximation, CPA). The evolution of nuclei is described by the classical Hamiltonian equations of motion (classical nuclei) with either statespecific (as in TSH) or state-averaged (Ehrenfest) forces. The electronic degrees of freedom (DOFs) are represented by the amplitudes of electronic basis states in a selected representation “rep”, 𝐶 = (𝑐1,𝑐2,…,𝑐𝑁)𝑇. The evolution of thes DOFs is determined by a simplified TD-SE: ∂𝐶𝑟𝑒𝑝

𝑖ℏ𝑆𝑟𝑒𝑝

∂𝑡

𝑝𝑛

= 𝐻𝑟𝑒𝑝𝐶𝑟𝑒𝑝 ―𝑖ℏ∑𝑛𝐷𝑟𝑒𝑝,𝑛𝑚𝑛𝐶𝑟𝑒𝑝.

(1)

Here, 𝐶𝑟𝑒𝑝 is the vector of the electronic basis state amplitudes in a selected representation, “rep”. 𝑆𝑟𝑒𝑝 and 𝐻𝑟𝑒𝑝 are the overlap and electronic Hamiltonian matrices in a given representation “rep”. 𝐷𝑟𝑒𝑝,𝑛 are the matrices of derivative nonadiabatic couplings (DNACs) for all nuclear DOFs, 𝑛 = 1,…𝑁𝑛𝑢𝑐𝑙. Commonly, the adiabatic (“rep=adi”) or diabatic (“rep=dia”) representations are used.


4

Page 5 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The solution of TD-SE, Eq. 1, requires the computation of the DNACs. In many atomistic implementation, the DNACs are not readily available. Moreover, they may be 𝑝𝑛

unnecessary, considering the fact that only the sum ∑𝑛𝐷𝑟𝑒𝑝,𝑛𝑚𝑛, which represents the time-

⟨

│∂𝑡∂ │𝜓𝑟𝑒𝑝⟩, is sufficient. It is common approach to

derivative NACs (TNACs), 𝑑 = 𝜓𝑟𝑒𝑝

compute NACs using the overlaps between the electronic wavefunctions at the adjacent time steps, 𝑆(𝑡,𝑡′) = ⟨𝜓𝑟𝑒𝑝(𝑡)│𝜓𝑟𝑒𝑝(𝑡′)⟩, following the numerical prescription of HammessSchiffer and Tully (HST)12:

(

𝑑 𝑡+

) = 𝑆(𝑡,𝑡 + Δ𝑡)2Δ𝑡― 𝑆(𝑡 + 𝛥𝑡,𝑡).

Δ𝑡 2

(2)

This is the point where a potential problem may arise that has a capacity to significantly affect the computed dynamics. Namely, the adiabatic wavefunction at different points (e.g. |𝜓𝑎𝑑𝑖(𝑡)⟩ and |𝜓𝑎𝑑𝑖(𝑡′)⟩) are not generally consistent in their phases. This problem may originate from the details of the eigensolver implementations and post-processing of the eigenvectors. With different mathematical packages adopting different standards and conventions, the phase behavior of the computed adiabatic wavefunction may vary from one computational system to another, with little to no expectation of consistency. Most of


5


Page 6 of 30

the time, the numerical solvers do not have any means of “unifying” the phases of wavefunction obtained for distinct nuclear geometries. Fundamentally, nothing prevents the solutions of a stationary SE to have distinct phases at distinct geometries. Consequently, the adiabatic wavefunctions that serve in NA-MD as the basis states to solve the TD-SE are defined only up to a complex phase. Ensuring the phase consistency of the wavefunctions used in Eq. 2 and in computing other properties entering Eq. 1 is necessary as far as the electron-nuclear dynamics is concerned. The straightforward application of the HST scheme with the “raw”, unregularized adiabatic basis states can make the NACs ill-defined and potentially affect the electron-nuclear dynamics of the system to the point of unreliability of the results. The problem of the phase consistency in numerical NA-MD calculations has been in the air for a while, but is still hasn’t been addressed in some implementation.13 Typically, the phase inconsistency is discussed in context of spurious state reordering that takes place in the course of the dynamics. A number of schemes has been proposed to resolve the state reordering problem, but not necessarily to correct the phases. Fernandez-Alberti et al.14 introduced the “mincost” algorithm to track the identities of adiabatic states in the


6

Page 7 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


TSH context, although the underlying coherent evolution stays unchanged. Similar “maximal overlap” procedures to resolve spurious reordering of states in TSH implementations are discussed by a number of other authors.15–18 Ryabinkin et al.19 suggested an elegant phase matching and reordering procedure. However, it mainly forces the wavefunction to choose a positive or negative overall sign, and may not be suitable when the wavefunctions must smoothly change the complex phase. This may be a potential issue when complex-values wavefunctions are used. This is a common situation with the periodic DFT calculations, which utilize complex (planewave) wavefunctions. The situation may be especially critical if overlaps of wavefunctions with distinct crystal momenta (k-points) are considered. It should be noted that many atomicorbital-based quantum chemistry codes do use real-valued basis functions, so the wavefunction correction by a factor of ± 1 is suitable in many cases. Mai et al.20 developed a more elaborate 3-step integrator to smoothly integrate TD-SE in the presence of trivial crossings. Meek and Levine21,22 developed a unitary interpolation scheme to facilitate an accurate integration of the TD-SE in the adiabatic representation (when the magnitude of the NACs may be large). Although the technique itself does not


7


Page 8 of 30

provide the way to regularize the phases of the wavefunction at the adjacent points in time, it relies on the nature of the starting and final wavefunctions (at the two points in time). If these functions are inconsistent in phase, the overall integration of TD-SE may be affected. The phase inconsistencies may be avoided by integrating TD-SE in the diabatic representation, e.g. when the diabatic Hamiltonian of a system is known18,23,24 or by using local diabatization.25–31 Nonetheless, the widespread use of the adiabatic representation demands simple and robust strategies for ensuring the phase consistencies. Although different degree of awareness regarding the phase corrections has been expressed, the effects of not addressing the phase inconsistencies in the NAMD simulations have not been under the focus. Thus, a further insight into this topic together with a simplified correction are needed. In this work, I introduce a simple phase correction method to ensure the phase consistency of the wavefunction used in the NA-MD schemes. This work focuses on the Ehrenfest dynamics, which in a simplified form serves as a backbone of the THS simulations and has a well-controlled mathematical foundation. Using the controlled model problems, I illustrate how the neglect of the phase consistency may change the


8

Page 9 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


nonadiabatic Ehrenfest dynamics integrated in the adiabatic representation. When the adiabatic states (and properties) are used as they are (directly obtained from the stationary SE eigenvalue problem), the dynamics obtained significantly deviates from the one in the diabatic representation (after the projection to the adiabatic one), which is always phase-consistent. The proposed phase correction scheme makes the two types dynamics consistent with each other. To describe the phase correction algorithm, I rely on equations and definitions already introduced above. More detailed definitions and conventions are summarized in the Supporting Information, Sections S1-S6. Assume that the adiabatic wavefunction for a state 𝑖 at times 𝑡 and 𝑡′ are |𝜓𝑖(𝑡)⟩ and |𝜓𝑖(𝑡′)⟩, respectively. Each function can be represented as a product of a phase factor, 𝑒𝑖𝜙𝑖(𝑡), and a phase-independent component, |𝜒𝑖(𝑡)⟩, such that |𝜓𝑖(𝑡)⟩ = 𝑒𝑖𝜙𝑖(𝑡)|𝜒𝑖(𝑡)⟩ and |𝜓𝑖(𝑡′)⟩ = 𝑒𝑖𝜙𝑖(𝑡′)|𝜒𝑖(𝑡′)⟩. To ensure the wavefunctions at times 𝑡 and 𝑡′ are consistent in phase, we introduce the phase correction factor: ⟨𝜓𝑖(𝑡)│𝜓𝑖(𝑡′)⟩

𝑓𝑖 = |⟨𝜓𝑖(𝑡)│𝜓𝑖(𝑡′)⟩|.

(3)


9


Page 10 of 30

This phase factor shows how much the phase of the wavefunction 𝜓𝑖 changes when it evolves from time 𝑡 to time 𝑡′. Indeed, considering our definitions above, we can see that 𝑓𝑖 = 𝑒𝑖[𝜙(𝑡′) ― 𝜙(𝑡)]. The phase-corrected functions at time 𝑡′, |𝜓𝑖(𝑡′)⟩, are then obtained by multiplying the “raw” functions |𝜓𝑖(𝑡′)⟩ by the complex conjugate of 𝑓𝑖: |𝜓𝑖(𝑡′)⟩ = |𝜓𝑖(𝑡′)⟩𝑓𝑖∗ .

(4)

Such a corrected function is consistent in phase with the function |𝜓𝑖(𝑡)⟩, which had been phase-corrected already to make it consistent with its predecessor. Indeed, according to the above definition, we obtain: |𝜓𝑖(𝑡′)⟩ = 𝑓𝑖∗ |𝜓𝑖(𝑡′)⟩ = 𝑒 ―𝑖[𝜙(𝑡′) ― 𝜙(𝑡)]𝑒𝑖𝜙(𝑡′) |𝜒𝑖(𝑡′)⟩ = 𝑒𝑖𝜙(𝑡)|𝜒𝑖(𝑡′)⟩. That is the overall phase factor of the function at time 𝑡′ is the same as that of the function at the time 𝑡: 𝑒𝑖𝜙(𝑡). The only change of the wavefunction is due to phase-independent part: |𝜒𝑖(𝑡)⟩→|𝜒𝑖(𝑡′)⟩, as presumed to be the case, e.g. when Eq. 2 is applied. Assume the adiabatic states are given as linear combination of the diabatic states as encoded by the diabatic-to-adiabatic transformation matrix, 𝑈: |𝜓𝑎𝑑𝑖⟩ = |𝜓𝑑𝑖𝑎⟩𝑈. Then, the phase correction described above consists of scaling the columns of the matrix 𝑈(𝑡′) by


10

Page 11 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the computed phase factors 𝑓𝑖 such that each column 𝑖 of the matrix is multiplied by the factor 𝑓𝑖∗ . Furthermore, the correction Eq. 4 should be reflected in the properties dependent on the wavefunctions. Namely, the NACs (DNACs in model calculations such as those considered the present work or TNACs such as those that appear in atomistic calculations), electronic couplings, etc. should be updated according to the indicated phase corrections of the wavefunctions. For instance, the DNACs are transformed as 𝑑𝛼𝑖𝑗 = ⟨𝜓𝑖│∇𝛼│𝜓𝑗⟩→𝑓𝑖⟨𝜓𝑖│∇𝛼│𝜓𝑗⟩𝑓𝑗∗ . That is the element 𝑑𝛼𝑖𝑗 is rescaled by the factor 𝑓𝑖𝑓𝑗∗ . Consequently, the diagonal properties, such as excited state forces (e.g. as in TSH) are not changed (for 𝑓𝑖𝑓𝑖∗ = 1). Thus, the phase correction may be less dramatic in the TSH nuclear dynamics. On the contrary, the nuclear forces in Ehrenfest dynamics involve offdiagonal contributions from DNACs, in which case the phase inconsistency would play a non-negligible role. Finally, to maintain the continuity of the overall wavefunction, Eq. S1, the amplitudes of the adiabatic states should be corrected as 𝐶𝑎𝑑𝑖,𝑖→𝐶𝑎𝑑𝑖,𝑖 = 𝑓𝑖𝐶𝑎𝑑𝑖,𝑖, to compensate for the phase correction of the adiabatic basis states. In practical implementation, the cumulative and relative phase corrections should be distinguished


11


Page 12 of 30

and used to transform the adiabatic states and their amplitudes, respectively. The details are elaborated in the Supporting Information, section S7. Two points should be kept in mind when considering the phase correction. First, the problem is associated with the adiabatic representation. The diabatic states have fixed phases by definition, so they are always phase-consistent. This means that the dynamics in the diabatic representation can serve as a reference to verify any phase correction algorithm. It is important that such a reference calculation relies only on the properties directly available (or evolved) in the diabatic representation, not those converted from the properties given in other (e.g. adiabatic) representations. Second, the phase inconsistencies usually develop when adiabatic states are nearly degenerate and may change their nature/ordering. This situation corresponds to the case when diabatic PESs cross each other. To demonstrate the utility of the phase correction algorithm, two 2-level and one 3-level model problems are considered. In the diabatic representation, the model Hamiltonians are given by: 𝐻𝑖𝑖 = 𝐴𝑖cos (𝜔𝑖𝑥 + 𝛿𝑖) + 𝐵𝑖,

(5a)


12

Page 13 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


𝐻𝑖𝑗 = 𝑉𝑖𝑗 = 𝑐𝑜𝑛𝑠𝑡.

(5b)

The models are distinct from each other by the parameters. The Model I is defined by: 𝐴0 = 0.1, 𝐴1 = 0.1, 𝐵0 = 0.0, 𝐵1 = ―0.1, 𝜔0 = 0.25, 𝜔1 = 0.25, 𝛿0 = 0.0, 𝛿1 = 𝜋/2, 𝑉01 = 0.01 . All parameters are given in the atomic units. These parameters make the diabatic surfaces cross in multiple spots (Figure 1a) and the adiabatic surfaces to be rather close to each other, delivering a notable NA coupling (Figure 1d). The Model II has the same parameters as Model I, except for 𝐵1 = 0.2, which eliminates the possibility of the diabatic PES crossings, even though there are still non-negligible nonadiabatic couplings (Figure 1, panels b and e). This model is designed to highlight the observation that the phase correction is only relevant when there are multiple crossings of diabatic PESs. Finally, the Model III considers 3 diabatic PESs with an intricate crossing pattern (Figure 1, panels c and f). The Hamiltonian matrix elements are defined by 𝐴0 = 0.1, 𝐴1 = ―0.1, 𝐴2 = ―0.2, 𝐵0 = 0.0, 𝐵1 = 0.1, 𝐵2 = 0.15, 𝜔0 = 0.25, 𝜔1 = 0.25, 𝜔2 = 0.20, 𝛿0 = 0.0, 𝛿1 = 0.0, 𝛿2 = 0.0, 𝑉01 = 0.05, 𝑉02 = 0.05, 𝑉12 = 0.03.


13


Page 14 of 30

(a)

(b)

(c)

(d)

(e)

(f)

Figure 1. The three model problems considered: (a, d) Model I; (b, e) Model II; (c, f) Model III. (a-c) PES, 𝐻𝑖𝑖, is the diabatic representation; (d-f) PES, 𝐸𝑖, in the adiabatic representation.

The nuclear subsystem is represented by a single degree of freedom (“a particle in 1D”, although the interpretation of this nuclear DOF can be more general than just a Cartesian coordinate) with an effective mass of 2,000 a.u., which roughly corresponds to the mass of the hydrogen atom. In all simulations, the particle is initially positioned at 𝑞0 = 0.1 𝑎.𝑢.,


14

Page 15 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


although the specific value is not critical for the purpose of the present study. Together with the chosen model parameters, this mass corresponds to the range of typical frequencies of organic molecules, e.g. 𝜔0 𝐴0/𝑚 ≈ 7.1 ∗ 10 ―4𝑎.𝑢. ―1, that is approximately 35 fs or 953 𝑐𝑚 ―1. The coupled electron-nuclear dynamics is integrated for 10000 steps with the integration timestep of 1 a.u. to avoid possible complications of the analysis due to the numerical integration inaccuracies. The electronic DOFs (amplitudes) are obtained by integrating TD-SE, Eq. 1 with all the quantities expressed in the selected representation. The nuclei follow the Ehrenfest forces, which can be computed using the representationspecific properties. The dynamics in two representations is executed: adiabatic and diabatic. The quantities relevant to other representation (non-dynamical) can be computed using the representation transformation relationship. For instance, if the dynamics is integrated in the diabatic representation, the diabatic populations (diagonal elements of the corresponding density matrix) are directly available, whereas the adiabatic populations are obtained by first performing the amplitudes transformation and then computing the needed density matrix. The Ehrenfest method with the detailed


15


Page 16 of 30

equations of motion, expressions for all key properties in both representations as well as transformations of various properties between the two representations are summarized in the Supporting Information (Sections S1 – S6). The methods and models are implemented in the Libra package (e.g. commit 196f181),32 which is used to performed all the computations. The Ehrenfest dynamics obtained in the diabatic representation is considered a reference, since there is no phase inconsistencies between the diabatic states at various geometries. The Ehrenfest dynamics in the adiabatic representation is computed with and without the proposed phase correction for all models. In addition, various initial conditions are considered: a) the initial wavefunction is selected either as one of the pure adiabatic states or as a coherent superposition of several; b) the initial momentum of the classical DOF is varied from 0 to 25 a.u. The combination of various models, initial condition, and simulation protocols (diabatic/adiabatic with and without phase correction) resulted in large number of computations, enough to ensure the robustness of the developed correction. The implementation of the computational and plotting scripts are available via


16

Page 17 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


GitHub.33 The typical results of the calculations are illustrated in Figure 2. A number of other representative figures is provided in the Supporting Information (Figures S2 – S10).

(a)

(b)

(c)

(d)


17


(e)

Page 18 of 30

(f)

Figure 2. Comparison of the results computed for the Model I with initial conditions 𝑝0 = 0 𝑎.𝑢. and Ψ(𝑡 = 0) = 𝜓0. (a, b) adiabatic population of the lowest state; (c, d) diabatic population of the lowest state; (e, f) phase space portrait. Left panels (a, c, e): no phase correction in the adiabatic dynamics; right panels (b, d, f): the phase correction is included in the adiabatic dynamics.

As it is apparent from Figure 2, the lack of phase consistency of the adiabatic states in the “raw” Ehrenfest dynamics (which utilizes the quantities directly available from the generalized eigenvalue solver, as implemented in Eigen3 package34) leads to a pronounced difference in the evolution of the adiabatic populations (Figure 2a). The integration in the diabatic representation (green line) shows only a slight variation of the


18

Page 19 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


population at certain points, whereas the overall population remains very close to the unity over the integration interval considered. At the same time, the integration in the adiabatic representation (red line) leads to a notable drop of this population by some 0.15 units. The application of the phase correction makes the adiabatic representation dynamics fully agree with the diabatic representation (Figure 2b). The remarkable difference in the coherent dynamics observed in Figure 2a could be further exaggerated in the TSH calculations, as our earlier experience shows.35 This effect originates from the asymmetry in the transition rates in the course of the TSH process: the detailed balance is present in the TSH and not in the Ehrenfest36 method, makes the re-population of higher-energy states less likely in the TSH than in the Ehrenfest. Therefore, any slight differences in the down-hopping probabilities produced during a coherent Ehrenfest-type dynamics are multiplied by the TSH algorithm, that effectively lets the states evolve down in energy, and inhibits transitions up in energy. The accumulation of such slight differences at every timestep can significantly change the TSH-based dynamics, leading to a notable change in the overall timescales. Thus, the phase correction may have even greater implications in TSH calculations than it does for


19


Page 20 of 30

the Ehrenfest dynamics. This complication can become even more cumbersome considering the complexity of the adiabatic PESs crossings and the need for additional tracking of the state identities.23 Analogous to the adiabatic populations, the diabatic populations behave very differently in the two types of simulations (Figure 2c): the population changes smoothly when the integration is done in the diabatic representation, but develops a highly oscillatory fine structure if the integration is done in the adiabatic representation. Moreover, the overall shapes of the population curves may or may not coincide. The oscillatory structure observed can be attributed to spurious phase changes of the adiabatic states. Because of the high frequency of such oscillations, the mean values of the computed quantities (e.g. forces, NACs) may be close to each other when the dynamics is run in each of the two representations. This is likely the reason for the general coincidence of the evolution of the diabatic populations for some parts of the trajectory. Yet, these slight deviations can lead to an overall deviation of the two evolutions in the longer runs or under more unfavorable conditions. The phase correction removes all the oscillatory features and makes the two evolutions fully coincide (Figure 2d).


20

Page 21 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Not only the electronic DOFs are affected by the phase inconsistencies of the adiabatic states. The evolution of nuclear DOFs is different in diabatic vs. adiabatic representation dynamics, if no phase corrections is applied (Figure 2e). The difference of the phase space portraits shown in Figure 2e can be much more dramatic, depending on the particular condition (e.g. see Figure S4 or S9). The phase correction helps remove the differences once again (Figure 2f). As was already noted above, the phase correction would not be necessary in the adiabatic dynamics, because the forces are insensitive to phase of the adiabatic wavefunctions. The phase inconsistency appears when the diabatic PESs intersect, which leaves the eigensolver confused about the sign (more generally, phase) assignment of the resulting adiabatic states. Thus, the problem is not anticipated when the diabatic surfaces do not cross. The calculations on the Model II (Figure S8) indeed confirm the consistency of the two types of dynamics even when the adiabatic-representation dynamics is done without the phase correction. To further validate the proposed phase correction approach, a 3level model (Model III) is considered. Within this model, all pairs of the diabatic states do intersect at certain geometries, making the model representative of a typical situation


21


Page 22 of 30

observed in modeling NA-MD in the condensed-matter systems, where multiple crossings are possible (both in adiabatic and diabatic representations). The results of the calculations (Figures S9-S10) do indeed follow the expectations already discussed above. For certain models and conditions, the phase-corrected adiabatic-representation based dynamics may show differences from the diabatic-representation dynamics. The observations that: a) this occurs at relatively large nuclear momenta (energy) and b) starting only after many integration steps (with the perfect agreement before the point of initial deviation) suggest that the origin of the discrepancy for such cases may be in the accumulation of the integration errors or the lack of the integration accuracy (when NACs become large). This effect can likely be fixed by the Meek-Levine interpolation method,21,22 but using phase-corrected adiabatic states as the input to the procedure. Generally, the PESs topology may be quite elaborate, with the identity (character) of electronic states changing after the point of the adiabatic surfaces’ intersection.23 A correct tracking of the individual states and their identities is of key importance in nonadiabatic dynamics.14,23 The neglect of this effect may lead to erroneously accelerated charge carrier diffusion (e.g. via non-physical long-range charge carrier hops)


22

Page 23 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


and biased dynamics of excited states population. Under such circumstances, the phasecorrection algorithm needs to be augmented with a proper state identity tracking algorithm, such as the recently developed crossing-corrected FSSH.23 In the present work, an alternative state reordering algorithm has been implemented, although its discussion is beyond the scope of this work. The simulations undertaken with and without this state reordering algorithm lead to the identical results, indicating that in this work we either have not run into the complex crossing problems, such as those discussed by Qiu et al.,23 or the state reordering had only a minor effect on the Ehrenfest dynamics. The latter expectation is based on the fact that energies and forces are determined by all the states, and the equal-energy states (e.g. at crossings) may have similar adiabatic amplitudes associated with them. Thus, the spurious reordering might be automatically accounted for by the mean-field nature of the Ehrenfest dynamics. However, this topic still needs a closer attention and it is likely that under certain situation, a more advanced approach to account for spurious state reordering within Ehrenfest or TSH dynamics may be needed.


23


Page 24 of 30

To summarize, it is illustrated that the phase inconsistencies of the “raw” adiabatic states may lead to a significant disturbance of the coherent dynamics in comparison to the reference phase-consistent dynamics performed in the diabatic representation. It is noted that such discrepancy may be magnified in the subsequent TSH calculations, although the investigation of this topic is beyond the scope of the present work and may need to account for a spurious state reordering. In this work, it is emphasized that the adiabatic wavefunctions must be corrected to ensure their phase consistency, before they are used in NA-MD calculations, e.g. to compute the NACs using the HST scheme, Eq. 2. A simple phase correction is introduced and justified. The application of the correction to nonadiabatic Ehrenfest dynamics shows excellent agreement with reference phaseinsensitive calculations on a number of model problems. The proposed phase correction is simple to implement which may help utilize it with the existing atomistic implementations of NA-MD for the latter to produce more accurate and reliable dynamics in systems with dense manifolds of excited states.

Associated Content


24

Page 25 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Supporting Information Available: Derivations and definitions relevant to Ehrenfest dynamics and representation transformations, results of Model I-III calculations with different initial conditions. (PDF)

Author Information Corresponding Author * Email: [email protected] Twitter: @AkimovLab

Notes The author declares no competing financial interests.

Acknowledgements

The author acknowledges the financial support from the University at Buffalo, The State University of New York startup package.

References

(1)

Grimaldi, G.; Crisp, R. W.; Brinck, S. ten; Zapata, F.; Ouwendorp, M. van; Renaud, N.; Kirkwood, N.; Evers, W. H.; Kinge, S.; Infante, I.; et al. Hot-Electron Transfer in Quantum-Dot Heterojunction Films. Nat. Commun. 2018, 9, 2310.


25


(2)

Page 26 of 30

He, J.; Vasenko, A. S.; Long, R.; Prezhdo, O. V. Halide Composition Controls Electron–Hole Recombination in Cesium–Lead Halide Perovskite Quantum Dots: A Time Domain Ab Initio Study. J. Phys. Chem. Lett. 2018, 9, 1872–1879.

(3)

Nijamudheen, A.; Akimov, A. V. Criticality of Symmetry in Rational Design of Chalcogenide Perovskites. J. Phys. Chem. Lett. 2018, 9, 248–257.

(4)

Sifain, A. E.; Bjorgaard, J. A.; Nelson, T. R.; Nebgen, B. T.; White, A. J.; Gifford, B. J.; Gao, D. W.; Prezhdo, O. V.; Fernandez-Alberti, S.; Roitberg, A. E.; et al. Photoexcited Nonadiabatic Dynamics of Solvated Push-Pull π-Conjugated Oligomers with the NEXMD Software. J. Chem. Theory Comput. 2018, 14, 3955– 3966.

(5)

Li, W.; Zhou, L.; Prezhdo, O. V.; Akimov, A. V. Spin–Orbit Interactions Greatly Accelerate Nonradiative Dynamics in Lead Halide Perovskites. ACS Energy Lett. 2018, 3, 2159–2166.

(6)

Yazdani, N.; Bozyigit, D.; Vuttivorakulchai, K.; Luisier, M.; Infante, I.; Wood, V. Tuning Electron–Phonon Interactions in Nanocrystals through Surface Termination.

Nano Lett. 2018, 18, 2233–2242. (7)

Li, L.; Kanai, Y. Dependence of Hot Electron Transfer on Surface Coverage and Adsorbate Species at Semiconductor–molecule Interfaces. Phys. Chem. Chem.

Phys. 2018, 20, 12986–12991. (8)

Crespo-Otero, R.; Barbatti, M. Recent Advances and Perspectives on Nonadiabatic Mixed Quantum–Classical Dynamics. Chem. Rev. 2018, 118, 7026–7068.

(9)

Barbatti, M.; Crespo-Otero, R. Surface Hopping Dynamics with DFT Excited States.

Top. Curr. Chem. 2016, 368, 415–444. (10) 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. (11) Tully, J. C. Molecular Dynamics with Electronic Transitions. J. Chem. Phys. 1990,

93, 1061–1071. (12) Hammes-Schiffer, S.; Tully, J. C. Proton Transfer in Solution: Molecular Dynamics with Quantum Transitions. J. Chem. Phys. 1994, 101, 4657–4667.


26

Page 27 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(13) Akimov, A. V.; Prezhdo, O. V. The PYXAID Program for Non-Adiabatic Molecular Dynamics in Condensed Matter Systems. J. Chem. Theory Comput. 2013, 9, 4959– 4972. (14) Fernandez-Alberti, S.; Roitberg, A. E.; Nelson, T.; Tretiak, S. Identification of Unavoided Crossings in Nonadiabatic Photoexcited Dynamics Involving Multiple Electronic States in Polyatomic Conjugated Molecules. J. Chem. Phys. 2012, 137, 014512. (15) Fabiano, E.; Keal, T. W.; Thiel, W. Implementation of Surface Hopping Molecular Dynamics Using Semiempirical Methods. Chem. Phys. 2008, 349, 334–347. (16) Bastida, A.; Soler, M. A.; Zúñiga, J.; Requena, A.; Kalstein, A.; Fernández-Alberti, S. Instantaneous Normal Modes, Resonances, and Decay Channels in the Vibrational Relaxation of the Amide I Mode of N-Methylacetamide-D in Liquid Deuterated Water. J. Chem. Phys. 2010, 132, 224501. (17) Kalstein, A.; Fernández-Alberti, S.; Bastida, A.; Soler, M. A.; Farag, M. H.; Zúñiga, J.; Requena, A. Vibrational Dynamics of Polyatomic Molecules in Solution: Assignment, Time Evolution and Mixing of Instantaneous Normal Modes. Theor.

Chem. Acc. 2011, 128, 769–782. (18) Wang, L.; Prezhdo, O. V. A Simple Solution to the Trivial Crossing Problem in Surface Hopping. J. Phys. Chem. Lett. 2014, 5, 713–719. (19) Ryabinkin, I. G.; Nagesh, J.; Izmaylov, A. F. Fast Numerical Evaluation of TimeDerivative Non-Adiabatic Couplings for Mixed Quantum-Classical Methods. J.

Phys. Chem. Lett. 2015, 6, 4200–4203. (20) Mai, S.; Marquetand, P.; González, L. A General Method to Describe Intersystem Crossing Dynamics in Trajectory Surface Hopping. Int. J. Quantum Chem. 2015,

115, 1215–1231. (21) Meek, G. A.; Levine, B. G. Evaluation of the Time-Derivative Coupling for Accurate Electronic State Transition Probabilities from Numerical Simulations. J. Phys.

Chem. Lett. 2014, 5, 2351–2356. (22) 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.


27


Page 28 of 30

(23) Qiu, J.; Bai, X.; Wang, L. Crossing Classified and Corrected Fewest Switches Surface Hopping. J. Phys. Chem. Lett. 2018, 9, 4319–4325. (24) Wang, L.; Beljonne, D. Flexible Surface Hopping Approach to Model the Crossover from Hopping to Band-like Transport in Organic Crystals. J. Phys. Chem. Lett. 2013,

4, 1888–1894. (25) Granucci, G.; Persico, M.; Toniolo, A. Direct Semiclassical Simulation of Photochemical Processes with Semiempirical Wave Functions. J. Chem. Phys. 2001, 114, 10608–10615. (26) Voityuk, A. A. Estimation of Electronic Coupling for Singlet Excitation Energy Transfer. J. Phys. Chem. C 2014, 118, 1478–1483. (27) Voityuk, A. A. Electronic Coupling for Charge Transfer in Donor–bridge–acceptor Systems. Performance of the Two-State FCD Model. Phys. Chem. Chem. Phys. 2012, 14, 13789-13793. (28) Voityuk, A. A. Electronic Couplings for Photoinduced Electron Transfer and Excitation Energy Transfer Computed Using Excited States of Noninteracting Molecules. J. Phys. Chem. A 2017, 121, 5414–5419. (29) Hoyer, C. E.; Xu, X.; Ma, D.; Gagliardi, L.; Truhlar, D. G. Diabatization Based on the Dipole and Quadrupole: The DQ Method. J. Chem. Phys. 2014, 141, 114104. (30) Hoyer, C. E.; Parker, K.; Gagliardi, L.; Truhlar, D. G. The DQ and DQΦ Electronic Structure Diabatization Methods: Validation for General Applications. J. Chem.

Phys. 2016, 144, 194101. (31) Li, S. L.; Truhlar, D. G.; Schmidt, M. W.; Gordon, M. S. Model Space Diabatization for Quantum Photochemistry. J. Chem. Phys. 2015, 142, 064106. (32) Akimov, A. V. Libra: An Open-Source “Methodology Discovery” Library for Quantum and Classical Dynamics Simulations. J. Comput. Chem. 2016, 37, 1626–1649. (33) Akimov, A. V. https://github.com/AkimovLab/Project_Ehrenfest_phase_corrections (accessed Oct 3, 2018). (34) Eigen http://eigen.tuxfamily.org/index.php?title=Main_Page (accessed Mar 29, 2016). (35) Akimov, A. V. Nonadiabatic Molecular Dynamics with Tight-Binding Fragment Molecular Orbitals. J. Chem. Theory Comput. 2016, 12, 5719–5736.


28

Page 29 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(36) Parandekar, P. V.; Tully, J. C. Detailed Balance in Ehrenfest Mixed QuantumClassical Dynamics. J. Chem. Theory Comput. 2006, 2, 229–235.


29


50x50mm (300 x 300 DPI)


Page 30 of 30

A Simple Phase Correction Makes a Big Difference in Nonadiabatic

Recommend Documents