Subscriber access provided by UNIVERSITY OF TOLEDO LIBRARIES
C: Energy Conversion and Storage; Energy and Charge Transport
A Solution for the Trivial Crossing Problem in Surface Hopping Simulations by the Classification on Excited States Zhen Sun, Sheng Li, Shijie Xie, and Zhong An J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.8b00084 • Publication Date (Web): 28 Mar 2018 Downloaded from http://pubs.acs.org on March 28, 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 23 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
A Solution for the Trivial Crossing Problem in Surface Hopping Simulations by the Classification on Excited States Zhen Sun,∗,† Sheng Li,† Shijie Xie,‡ and Zhong An¶ †Department of Physics, Zhejiang Normal University, Jinhua, Zhejiang 321004, China ‡School of Physics, State Key Laboratory of Crystal Materials, Shandong University, Jinan, 250100, China ¶College of Physics, Hebei Normal University, Shijiazhuang 050016, China E-mail:
[email protected] 1
ACS Paragon Plus Environment
The Journal of Physical Chemistry 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
Abstract In this paper, we study excited-state relaxation dynamics in single conjugated polymer chains with surface hopping simulations and propose a solution for the encountered “trivial crossing problem”. By analyzing the transition density matrices, we find that all excited states can be divided into two types. The surface hopping simulations show that the two types of excited states do not interact each other (the nonadiabatic couplings between them are negligible), indicating that the trivial crossings only occur between different types of excited state. The two types of excited states separately span two uncoupled subspaces, and we can let the system evolve only in one of the subspaces. In such a way, the trivial crossings between different types of excited states are completely eliminated. Moreover, this method greatly reduces the computational cost of surface hopping simulations.
Introduction Recent years, the fewest switches surface hopping (FSSH) method, described by Tully, 1 is in rapid developing and becomes one of the most popular methods in the field of nonadiabatic dynamics. 2,3 It has been widely used to simulate various nonadiabatic processes, such as photoexcited dynamics in organic molecules 4–7 and charge transport in organic crystals. 8,9 Moreover, many programs (e.g., Newton-X, 10 MNDO, 11 and NA-ESMD 12 ) have been developed in the past decade, that combine the FSSH method with ab initio electronic structure packages. The popularity and success of the FSSH method is due to its ease of implementation and the fact that it can provide fairly accurate results for many problems. In the FSSH simulations, when the current state approaches an other state which is noninteracting with the current state, the nonadiabatic coupling between them is described as sharp peaks that are strongly localized in time. In such cases, a hop to the new state should occur with unit probability, and hence it is called “trivial crossing”. If a hop does not occur during a trivial crossing, it may lead to unphysical long range charge or energy transfer. 2
ACS Paragon Plus Environment
Page 2 of 23
Page 3 of 23 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
This is the so-called trivial crossing problem. Several algorithms have been developed for detecting trivial crossing in surface hopping simulations. 11,13–17 Tretiak et al. proposed a solution which is based on a Min-Cost algorithm, calculating overlaps between adiabatic states and comparing them with a certain threshold. 13 However, it is hard to eliminate all trivial crossings because there always exist intermediate cases that cannot be distinguished. In addition, when the overlap is lower than the threshold during a time step, it needs to reduce the time step to recalculate the process. This may bring huge extra computational cost. Prezhdo et al. solve Schrödinger equations in a diabatic representation and then transform the solution to the adiabatic representation at each time step. 15 They compare the sum of the surface hopping probabilities calculated from these two representations to detect trivial crossings. As demonstrated in this algorithm, the cleanness of this approach depends on the time step. On the other hand, the diabatic basis are generally not ready-made, but obtained from adiabatic basis transformations. Granucci et al. recommend using a local diabatic basis to propagate the electronic wave function. 16 However, there exist no unique diabatic representation, and the accuracy of simulations performed in the diabatic representation is lower than those based on adiabatic representation. 18 In the current paper, we provide a solution to this problem from a new perspective. In our case, this method can completely eliminate all trivial crossings between non-interacting excited states. Moreover, this method can greatly reduce the computational cost. To do that, we divide all excited states into two uncoupled (non-interacting) subspaces, and the excited states in each subspace can be treated isolatedly. In such a way, the trivial crossings between non-interacting states are completely eliminated.
Methods Over the years, we developed our own code for the FSSH method. Our approach combines the semiempirical Pariser-Parr-Pople (PPP) Hamiltonian 19 with configuration interaction singles
3
ACS Paragon Plus Environment
The Journal of Physical Chemistry 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 4 of 23
(CIS) to calculate the excited-state electronic structure. 20 We consider a one-dimensional model system to mimic a single polymer chain. Each site has a π electron, and all these electrons are described by PPP Hamiltonian. The sites are treated as classical nuclei cores. The total Hamiltonian H = He + Hn , where He and Hn are written as
He = −
X
tµ,µ+1 (c†µ,s cµ+1,s + c†µ+1,s cµ,s ) +
µ,s
X 1 1 (nµ↑ − )(nµ↓ − ) + Vµν (nµ − 1)(nν − 1), 2 2 µ µ,ν 1 X = K (rµ,µ+1 − r0 )2 . 2 µ U
Hn
X
(1) (2)
In eq. (1), c†µ,s and cµ,s are the creation and annihilation operators for a π electron with P spin s at site µ. nµ = s c†µ,s cµ,s is the total number of electrons on site µ. tµ,µ+1 is the nearest hopping integral, tµ,µ+1 = t0 − α(rµ,µ+1 − r0 ) + (−1)µ+1 te . It depends on the hopping integral t0 for equal intersite distances r0 , the distance between neighbor sites rµ,µ+1 ≡ |xµ − xµ+1 |, the electron-phonon coupling constant α, and te the extrinsic transfer term, introduced to lift the ground state degeneracy. U and Vµν are the on-site and inter-site Coulomb interactions, respectively. The Vµν are obtained from a modification of the Ohno 1
2 2 ) , where κ denotes an effective dielectric constant. parameterization Vµν = U/κ(1+0.6117rµν
The value of parameters U and κ are set 8 eV and 2, respectively, derived by Chandross and Mazumdar. 21 In eq. (2), K is the elastic constant. In the present work, the parameters in this model are set to common values: t0 = 2.0 eV, α = 4.0 eV/Å, r0 = 1.50 Å, K = 45.0 2
eV/Å , and te = 0.05 eV. In the FSSH method, the nuclear degrees of freedom are propagated on currently occupied electronic state φK according to the Langevin equation of motion 22
M x¨µ (t) = −
K dEtot − γM x˙ µ (t) + ζ(t), dxµ
(3)
K where Etot is the total energy of K-th excited state, M , x¨µ and x˙ µ represent the mass,
4
ACS Paragon Plus Environment
Page 5 of 23 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
K acceleration and velocity of the µth site, respectively. The analytical calculation of dEtot /dxµ
is outlined in SI. The stochastic force ζ depends on the bath temperature T and the friction coefficient γ. It obeys the fluctuation-dissipation theorem satisfying the condition: hζ(t)ζ(t+ ∆t)i = 2M γkB T δ(∆t). In this work, the bath temperature T is fixed to 300 K, since we do not discuss the temperature effects. A numerical velocity Verlet algorithm is used to integrate Eq. (3). The classical time step ∆t=0.01 fs. To account for quantum transitions among adiabatic excited states, we employ the FSSH algorithm of Tully. The electronic wavefunction of the system is expressed by a linear P combination of the adiabatic CIS states |Φi = I cI |φI i, where {φI } are the eigenvectors of the CIS matrix of He . Using the time-dependent Schrödinger equation, we obtain
c˙K =
X X 1 cI x˙ µ dµIK , EK c K − i¯h µ I6=K
(4)
where dµIK = hφI |dφK /dxµ i, being the nonadiabatic coupling vector, is analytically calculated. The calculation is outlined in SI. Eq. (4) is solved through the standard fourth-order Runge-Kutta algorithm 23 . The electronic evolution requires a smaller time step compared with the nuclear evolution. Here, the quantum time step δt = 0.001 fs. Following Tully’s FSSH algorithm, the switching probability from the currently occupied excited state K to another excited state I is
gKI =
Re[cK c∗I ]
PNq P
j=1 µ ∗ cK cK
x˙ µ dµKI δt
,
(5)
where Nq is the number of quantum steps per classical step. In the simulations, we employ the “instantaneous decoherence” procedure to add the “decoherence effect”. 24 Following this procedure, after each attempted hop, regardless of whether the hop is allowed or forbidden due to energy conservation constrains, the quantum coefficients are reinitialized and a coefficient of one is assigned to the new current state, while the populations of all other states are erased. 5
ACS Paragon Plus Environment
The Journal of Physical Chemistry 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
Results and Discussion The preliminary step for surface hopping simulations is to prepare the initial site coordinates and velocities. To do that, an adiabatic molecular dynamics is run on the ground state at T=300 K. After the polymer is equilibrated, a snapshot is taken every 0.2 ps. The initial site coordinates and velocities are obtained from these snapshots. Excited-state dynamics are started from these initial configurations by populating the initial excited state, i.e., set the initial values of the quantum coefficients.
6
ACS Paragon Plus Environment
Page 6 of 23
Page 7 of 23
20
Site
18
20 (a)
18
18
16
16
14
14
14
12
12
12
10
10
10
8
8
8
6
6
6
4
4
4
2 2
4
6
8 10 12 14 16 18 20
20 18
Site
20 (b)
16
2
18
4
6
8 10 12 14 16 18 20
2 18
16
16
14
14
14
12
12
12
10
10
10
8
8
8
6
6
6
4
4
4
2
2 6
18
4
6
(g)
18
18
16
16
14
14
14
12
12
12
10
10
10
8
8
8
6
6
6
4
4
4
2
2 4
6
8 10 12 14 16 18 20
4
6
8 10 12 14 16 18 20
4
6
8 10 12 14 16 18 20
20 (h)
16
2
8 10 12 14 16 18 20
(f)
2
8 10 12 14 16 18 20
20
20
6
2 2
8 10 12 14 16 18 20
4
20 (e)
16
4
(c)
2 2
20 (d)
2
Site
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
(i)
2 2
4
6
8 10 12 14 16 18 20
Site
Site
2
Site
Figure 1: (a) - (i) correspond to S1 - S9 , respectively, showing the contour maps of transition density matrices of these excited states.
We start with analyzing the excited-state transition density matrix QIµν , which is wildly used in the description of excited-state properties. 25–27 The diagonal elements QIµµ represent the transition density on the µ-th site when the polymer undergoes the transition from ground state to the I-th excited state, whereas the element QIµν (µ 6= ν) represents the joint probability of finding an extra electron on site µ and a hole on site ν. Transition density 7
ACS Paragon Plus Environment
The Journal of Physical Chemistry 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
matrices establish a connection between electron-hole pair distribution in real space and molecular optical properties. The transition density matrices are commonly displayed using contour charge density maps. Fig. 1 shows the contour maps of the transition density matrices for the lowest nine singlet excited states. These transition density matrices are calculated using one of the initial lattice configurations which are prepared above. From these pictures, we find that some of them, including S1 , S3 , S5 , S6 and S9 , have nonzero diagonal elements (QIµµ 6= 0), while the rest of them, including S2 , S4 , S7 and S8 , have zero diagonal elements (QIµµ = 0). We call the first group type I excited state, and the second group type II excited state. Actually, if we examine the transition density matrices of all excited states, we will find that all excited states can be divided into these two types. Moreover, we find that this classification is independent with lattice configuration of the polymer chain. Type I and type II excited states have many different characters. From the pictures in Fig. 1, we can see that the distributions of type I excited states are along the diagonal of the picture, while the distributions of type II excited states along the off-diagonal of the picture. This means that in type II excited states the probabilities of finding the electron and the hole within far distance are higher than in type I excited states, indicating that the energies of type II excited states respond more sensitively to external environment (e.g., geometry changing, inter-chain interactions and external field) relative to the energies of type I excited states. For this reason, type I and type II excited states may cross at energies during excited-state dynamics.
8
ACS Paragon Plus Environment
Page 8 of 23
Page 9 of 23 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
Table 1: CIS-calculated energy levels, types and oscillator strengths of the lowest nine singlet excited states. Excited state Energy(eV) Type S1 2.43 I S2 3.11 II S3 3.31 I S4 3.71 II S5 4.06 I S6 4.15 I S7 4.38 II S8 4.64 II S9 4.75 I
Oscillator strength 9.90 0 0.04 0 2.09 0.13 0 0 0.02
Table I summarizes the CIS-calculated energies relative to the ground state, types, as well as oscillator strengths of the lowest nine singlet excited states. These results are calculated using the lattice configuration as same as in Fig. 1. From this table, we can see that the oscillator strengths of all type I excited states are nonzero, which means that they can be directly generated by absorbing a photon. In contrast, the oscillator strengths of all type II excited states are zero because there is no overlap between electron and hole wavefunctions since QIµµ = 0. Therefore, type II excited states can not be directly generated by absorbing a photon.
9
ACS Paragon Plus Environment
The Journal of Physical Chemistry
0.10 (a)
3.0x10
Nonadiabatic coupling (a.u.)
Nonadiabatic coupling (a.u.)
-7
54
-7
2.0x10
-7
1.0x10
0.0 -7
-1.0x10
-7
-2.0x10
(b)
0.05
0.00
-0.05
-7
-3.0x10
-0.10 0
20
40
60
80
100
0
20
Time (fs)
40
60
80
100
Time (fs)
-198.0 S
(c)
-198.5 Energy (eV)
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 10 of 23
S S
-199.0
S
-199.5
S
1
2
3
4
5
-200.0 -200.5 -201.0 -201.5 0
20
40
60
80
100
Time (fs)
Figure 2: (a), (b) and (c) show time dependence of the nonadiabatic couplings between excited state I and J σIJ and energy levels, respectively. They are obtained after the polymer is initialized in state S5 .
Keep these results in mind, we begin to perform the FSSH simulation. We first focus on the time dependence of the nonadiabatic couplings between different excited states after the system is photoexcited to a certain excited state. The dynamical simulations are restricted to the first five excited states. We assume that the system is vertically excited to state S5 . The time dependence of nonadiabatic couplings between state S5 and other four excited states for a typical trajectory are showed in the top row of Fig. 2. Fig. 2 (a) shows that the nonadiabatic coupling σ54 exhibits sharp peaks extremely localized in time, and they vanish elsewhere. This behavior indicates that state S5 does not interact with state S4 . The sharp
10
ACS Paragon Plus Environment
Page 11 of 23
peaks are caused by the infinite approach of the two states during dynamics. This is testified by the time evolution of energy levels, see Fig. 2 (c). According to our previous discussion, a hop from state S5 to S4 should occur with unit probability when they are approaching each other. However, in the present simulation the “trivial crossing" problem is not treated and the transition probabilities are evaluated according to the normal method. If a hop from state S5 to S4 does not happen during a trivial crossing, the electron/hole distribution on the chain will undergo unphysical drastic changes. As a result, the nonadiabatic coupling between state S5 and other exited states (S3 , S2 and S1 ) will change unnormally and suddenly at those points, as we see in Fig. 2 (b). So the trivial crossings must be properly treated.
0.10
0.10
(b)
Nonadiabatic coupling (a.u.)
(a)
0.05
0.05
0.00
0.00
-0.05
-0.05
-0.10
-0.10 0
20
40
60
80
100
0
0.10
Nonadiabatic coupling (a.u.)
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
20
40
60
80
100
0.10
(d)
(c) 0.05
0.05
0.00
0.00
-0.05
-0.05
-0.10
-0.10 0
20
40
60
80
100
0
20
Time (fs)
40
60
80
100
Time (fs)
Figure 3: Time dependence of the nonadiabatic couplings between excited state I and J σIJ , which are obtained after the polymer is initialized in state S3 .
11
ACS Paragon Plus Environment
The Journal of Physical Chemistry 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
If the polymer is vertically excited to state S3 , the time dependence of all nonadiabatic coupling terms (there are totally ten unknown nonadiabatic coupling terms since σIJ = −σJI ) are shown in Fig. 3. We can see that none of the curves in Fig. 3 have unnormal changes, indicating that there are no trivial crossings happening during the simulation. We also can see that the nonadiabatic coupling terms σ54 , σ52 , σ43 , σ41 , σ32 and σ21 are always zero during the simulation. If we combine the above results with the classification on excited states, we will find that the couplings between type I excited states (S1 , S3 and S5 ) and type II excited states (S2 and S4 ) are always zero during the whole simulation time. This means that type I excited states do not interact with type II excited states. An excited state only interacts with other excited states with the same type. That is to say, trivial crossings only occur between excited states with different types. We try our best to test more excited states, and find that this is a general rule and can be applied to all excited states. We think that the deep reason for that is connected with symmetry of the system. The above results indicate that all type I and all type II excited states span two uncoupled subspaces, and they can be treated separately. Therefore, we can let the system evolve in a subspace which separately contains only type I excited states or type II excited states. By doing that, the trivial crossings are completely eliminated because they only occur between excited states with different type. According to our discussions on the properties of type I and type II excited states, type I or type II excited states can be chosen according to whether the diagonal elements of their transition density matrices being zero or not. Actually, in computer calculations, the diagonal elements of transition density matrices of type II excited states are not zero exactly. To solve this problem, we can set a very small value, for example 10−10 as criteria. If the largest diagonal element of transition density matrix of an excited state is smaller than this value, it is a type II excited state. Otherwise, it is a type I excited state. In general, there are no intermediate cases that the two types of excited states cannot be distinguished because of the huge differences in diagonal elements of transition density matrices between type I and type II excited states. For type I excited states the diagonal
12
ACS Paragon Plus Environment
Page 12 of 23
Page 13 of 23
elements of transition density matrices are typically on the order of 10−1 , while for type II excited states they are typically less than 10−12 . Because of such a big gap, the two types of excited states are easy to be distinguished. It should be stressed that intermediate cases and trivial crossings within the same type of excited states may be possible in general, and other methods discussed in the introduction can be used. However, these cases do not exist in the investigated single conjugated polymers, so the proposed classification technique on excited states is enough to describe the encountered trivial crossing problem. We note that the selection of type I or type II excited states should be done in every time step because an excited state may change its type when the lattice configuration changes. In addition, by using this method, the dimensionality of the problem to be treated is reduced since the system evolves in a subspace. So the computation cost is greatly lowered. This is an other advantage of this method. 1.0
1.0
(a) S
0.8
Population
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
S S
0.6
S S
(b)
S
1
0.8 2
S
3
S 0.6
I 1 I 2 I 3
4
5
0.4
0.4
0.2
0.2
0.0
0.0 0
200
400
600
800
1000
0
200
Time (fs)
400
600
800
1000
Time (fs)
Figure 4: Time dependence of the population (i.e., the occupation probability) of (a) the lowest five singlet excited states without treating the trivial crossing problem, and (b) the lowest three singlet type I excited states for the polymer chain, averaged over 200 trajectories.
We calculate the time dependence of the populations of the lowest five singlet excited states (S1 - S5 ) without treating the trivial crossing problem, as shown in Fig. 4 (a). From Fig. 4 (a), we can see that the system is initialized in state S5 , since the population of state S5 is one at the time of beginning. Then, it decays approximately exponentially, reaching 0.2 at the end of simulation time. With the population of state S5 decreasing, only state S3 and 13
ACS Paragon Plus Environment
The Journal of Physical Chemistry 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
S1 are markedly populated, reaching 0.34 and 0.43, respectively, at the end of simulation time. Since state S5 , S3 and S1 are type I excited states (see Table I), this is consistent with our result that an excited state only interact with other excited states with the same type. However, there are tiny populations (less than 0.02) of type II excited state S4 and S2 involved in the evolution. This is because we do not treat the trivial crossing problem in Fig. 4(a). When state S5 approaches state S4 , the current occupied state is most probably still S5 because the hop probability is assessed by the normal FSSH procedure, instead of setting to unit. In contrast, we calculate the time dependence of the populations of the lowest three type I singlet excited states (SI1 - SI3 ), as shown in Fig. 4 (b). We can see that the evolution tendencies of state SI3 , SI2 and SI1 are similar with state S5 , S3 and S1 , respectively. The differences are that the decay of state SI3 is slower than state S5 , reaching 0.3 at the end of simulation time. Correspondingly, the populations of SI2 and SI1 increase slower than S3 and S1 , respectively, both reaching 0.34. In Fig. 4(a) and 4(b), the similarities between curves SI3 , SI2 , SI1 and curves S5 , S3 , S1 are because states SI3 , SI2 , SI1 are same as states S5 , S3 , S1 , respectively, as shown in Table I. The differences between curves SI3 , SI2 , SI1 and curves S5 , S3 , S1 are because the trivial crossing problem is not treated in Fig. 4(a). As we have known, the miss of trivial crossings may lead to unphysical long range charge or energy transfer. Because of this, the forces from electrons on nuclear may undergo unphysical changes, leading to different nonadiabatic couplings between excited states, and then the populations of them. On the other hand, the miss of trivial crossings may also lead to the transition from state S5 to S2 . This is an other reason of the differences between curves SI3 , SI2 , SI1 and curves S5 , S3 , S1 .
14
ACS Paragon Plus Environment
Page 14 of 23
Page 15 of 23 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
Conclusions In summary, we proposed a solution for the trivial crossing problem in FSSH simulations. Our solution is based on classifying excited states into two types according to whether the diagonal elements of their transition density matrices being zero or not. The FSSH simulations show that the nonadiabatic couplings between the two types of excited states are negligible, and the trivial crossings only occur between excited states with different types. That is to say, the two types of excited states do not interact each other, and they span two uncoupled subspaces. We can let the system evolve in a subspace which contains only one type of excited states. In this way, the trivial crossings between different types of excited states are completely eliminated. It is worth to note that the proposed solution for trivial crossing problem in this paper may not work generally. This approach works well only if the considered excited states can be classified into two types according to whether the diagonal elements of their transition density matrices are zero or not. This may not be the case for all kinds of systems. Our studies show that this method is applicable to most single-chain systems, such as poly-pphenylenevinylene (PPV), polyacetylene, pentacene and so on. While for other systems, one need to check before using this method.
Acknowledgement This work was supported by the National Natural Science Foundation of China (Grant No. 21374105 and No.11574180) and the Zhejiang Provincial Science Foundation of China under Grant LR12B04001). We gratefully acknowledge Prof. Desheng Liu, Dr. Changfeng Fang, Dr. Jianhua Zhao, and Dr. Dawei Kang for providing computational support.
15
ACS Paragon Plus Environment
The Journal of Physical Chemistry 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
References (1) Tully, J. C. Molecular dynamics with electronic transitions. Journal of Chemical Physics 1990, 93, 1061–1071. (2) Wang, L.; Akimov, A. V.; Prezhdo, O. V. Recent Progress in Surface Hopping: 20112015. The journal of physical chemistry letters 2016, (3) Doltsinis, N. L. Nonadiabatic Dynamics: Mean-Field and Surface Hopping. Quantum Simulations of Complex Many 2002, (4) Nelson, T.; Fernandezalberti, S.; Chernyak, V.; Roitberg, A. E.; Tretiak, S. Nonadiabatic excited-state molecular dynamics modeling of photoinduced dynamics in conjugated molecules. Journal of Physical Chemistry B 2011, 115, 5402–14. (5) Zobač, V.; Lewis, J. P.; Abad, E.; Mendietamoreno, J. I.; Hapala, P.; Jelínek, P.; Ortega, J. Photo-induced reactions from efficient molecular dynamics with electronic transitions using the FIREBALL local-orbital density functional theory formalism. Journal of Physics Condensed Matter An Institute of Physics Journal 2015, 27, 175002. (6) Prlj, A.; Curchod, B. F.; Corminboeuf, C. Excited state dynamics of thiophene and bithiophene: new insights into theoretically challenging systems. Physical Chemistry Chemical Physics 2015, 17, 14719–14730. (7) Fazzi, D.; Barbatti, M.; Thiel, W. Modeling ultrafast exciton deactivation in oligothiophenes via nonadiabatic dynamics. Physical Chemistry Chemical Physics 2015, 17, 7787–99. (8) Wang, L.; Beljonne, D. Flexible Surface Hopping Approach to Model the Crossover from Hopping to Band-like Transport in Organic Crystals. Journal of Physical Chemistry Letters 2013, 4, 1888–1894.
16
ACS Paragon Plus Environment
Page 16 of 23
Page 17 of 23 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
(9) Spencer, J.; Gajdos, F.; Blumberger, J. FOB-SH: Fragment orbital-based surface hopping for charge carrier transport in organic and biological molecules and materials. Journal of Chemical Physics 2016, 145, 926–71. (10) Barbatti, M.; Ruckenbauer, M.; Plasser, F.; Pittner, J.; Granucci, G.; Persico, M.; Lischka, H. Newton-X: a surface-hopping program for nonadiabatic molecular dynamics. Wiley Interdisciplinary Reviews Computational Molecular Science 2014, 4, 26–33. (11) Fabiano, E.; Keal, T. W.; Thiel, W. Implementation of surface hopping molecular dynamics using semiempirical methods. Chemical Physics 2008, 349, 334–347. (12) Nelson, T.; Fernandezalberti, S.; Roitberg, A. E.; Tretiak, S. Nonadiabatic ExcitedState Molecular Dynamics: Modeling Photophysics in Organic Conjugated Materials. Accounts of Chemical Research 2014, 47, 1155–64. (13) Fernandezalberti, S.; Roitberg, A. E.; Nelson, T.; Tretiak, S. Identification of unavoided crossings in nonadiabatic photoexcited dynamics involving multiple electronic states in polyatomic conjugated molecules. Journal of Chemical Physics 2012, 137, 460–19249. (14) Nelson, T.; Fernandez-Alberti, S.; Roitberg, A. E.; Tretiak, S. Artifacts due to trivial unavoided crossings in the modeling of photoinduced energy transfer dynamics in extended conjugated molecules. Chemical Physics Letters 2013, 590, 208–213. (15) Wang, L.; Prezhdo, O. V. A Simple Solution to the Trivial Crossing Problem in Surface Hopping. Journal of Physical Chemistry Letters 2014, 5, 713. (16) Granucci, G.; Persico, M.; Toniolo, A. Direct semiclassical simulation of photochemical processes with semiempirical wave functions. Journal of Chemical Physics 2001, 114, 10608–10615. (17) Plasser, F.; Granucci, G.; Pittner, J.; Barbatti, M.; Persico, M.; Lischka, H. Surface hopping dynamics using a locally diabatic formalism: charge transfer in the ethylene 17
ACS Paragon Plus Environment
The Journal of Physical Chemistry 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
dimer cation and excited state dynamics in the 2-pyridone dimer. Journal of Chemical Physics 2012, 137, 22A514. (18) Prezhdo, O. V.; Rossky, P. J. Mean-field molecular dynamics with surface hopping. Journal of Chemical Physics 1997, 107, 825–834. (19) Jug, K. Theoretical basis and design of the PPP model Hamiltonian. International Journal of Quantum Chemistry 1990, 37, 403–414. (20) Gadaczek, I.; Krause, K.; Hintze, K. J.; Bredow, T. MSINDO-sCIS: A new method for the calculation of excited states of large molecules. Journal of chemical theory and computation 2011, 7, 3675–3685. (21) Chandross, M.; Mazumdar, S. Coulomb interactions and linear, nonlinear, and triplet absorption in poly (para-phenylenevinylene). Physical Review B 1997, 55, 1497. (22) Paterlini, M. G.; Ferguson, D. M. Constant temperature simulations using the Langevin equation with velocity Verlet integration. Chemical Physics 1998, 236, 243–252. (23) Hull, T. E.; Enright, W. H.; Jackson, K. R. User’s guide for DVERK-a subroutine for solving non-stiff ODEs. Department of Computer Science. 1976. (24) Nelson, T.; Fernandezalberti, S.; Roitberg, A. E.; Tretiak, S. Nonadiabatic excited-state molecular dynamics: treatment of electronic decoherence. Journal of Chemical Physics 2013, 138, 1061. (25) Tretiak, S.; Mukamel, S. Density matrix analysis and simulation of electronic excitations in conjugated and aggregated molecules. Chemical reviews 2002, 102, 3171–3212. (26) Sun, M.; Kjellberg, P.; Beenken, W. J.; Pullerits, T. Comparison of the electronic structure of PPV and its derivative DIOXA-PPV. Chemical physics 2006, 327, 474– 484.
18
ACS Paragon Plus Environment
Page 18 of 23
Page 19 of 23 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
(27) Bittner, E. R.; Ramon, J. G. S.; Karabunarliev, S. Exciton dissociation dynamics in model donor-acceptor polymer heterojunctions. i. energetics and spectra. The Journal of chemical physics 2005, 122, 214719.
19
ACS Paragon Plus Environment
The Journal of Physical Chemistry
20
Site
18
20 (a)
18
(b)
18
16
16
14
14
14
12
12
12
10
10
10
8
8
8
6
6
6
4
4
4
2
2
2
4
6
8 10 12 14 16 18 20
20 18
Site
20
16
2
2
4
6
8 10 12 14 16 18 20
(d)
18
2
(e)
18
16
16
16
14
14
14
12
12
12
10
10
10
8
8
8
6
6
6
4
4
4
2
2
2
4
6
4
6
8 10 12 14 16 18 20
20
20 18
2
8 10 12 14 16 18 20
(g)
18
18 16
14
14
14
12
12
12
10
10
10
8
8
8
6
6
6
4
4
4
2
2
2
6
8 10 12 14 16 18 20 Site
6
8 10 12 14 16 18 20
4
6
8 10 12 14 16 18 20
4
6
8 10 12 14 16 18 20
20 (h)
16
4
4
(f)
2
16
2
(c)
20
20
2
Site
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 20 of 23
2
4
6
8 10 12 14 16 18 20 Site
ACS Paragon Plus Environment
(i)
2
Site
Page 21 of 23
Nonadiabatic coupling (a.u.)
-7
(a)
3.0x10
54
-7
2.0x10
-7
1.0x10
0.0 -7
-1.0x10
-7
-2.0x10
Nonadiabatic coupling (a.u.)
0.10 (b)
0.05
0.00
-0.05
-7
-3.0x10
-0.10 0
20
40
60
80
100
0
Time (fs)
S
(c)
-198.5
S S
-199.0
S
-199.5
S
-200.0 -200.5 -201.0 -201.5 0
20
40
60
20
40
60
Time (fs)
-198.0
Energy (eV)
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
80
100
Time (fs)
ACS Paragon Plus Environment
1
2
3
4
5
80
100
The Journal of Physical Chemistry
0.10
0.10
(b)
Nonadiabatic coupling (a.u.)
(a)
0.05
0.05
0.00
0.00
-0.05
-0.05
-0.10
-0.10 0
20
40
60
80
100
0.10
Nonadiabatic coupling (a.u.)
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 22 of 23
0
20
40
60
80
100
0.10
(d)
(c) 0.05
0.05
0.00
0.00
-0.05
-0.05
-0.10
-0.10 0
20
40
60
80
100
0
Time (fs)
20
40
60
Time (fs) ACS Paragon Plus Environment
80
100
Page 23 of 23 1.0
Population
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
The Journal of Physical Chemistry 1.0
(b)
(a) S
0.8
S S
0.6
S S
S
1
0.8
S
2
3
S 0.6
0.2
0.2
0.0
0.0 400
600
Time (fs)
800
I 2 I 3
5
0.4
200
1
4
0.4
0
I
1000 0 ACS Paragon Plus Environment
200
400
600
Time (fs)
800
1000