Clarification of Nonadiabatic Chemical Dynamics
by the Zhu-Nakamura Theory of Nonadiabatic
Transition: From Tri-atomic Systems
to Reactions in Solutions
Toshimasa Ishida1) , Shinkoh Nanbu2) ∗, and Hiroki Nakamura3) †1)
Department of Chemistry, University of Malaya, 50603,
Kuala Lumpur, Malaysia.
e-mail: ishida@fukui.kyoto-u.ac.jp
2)
Department of Materials and Life Sciences,
Faculty of Science and Technology, Sophia University,
7-1 Kioi-cho, Chiyoda-ku, Tokyo, 102-8554, Japan.
e-mail: shinkoh.nanbu@sophia.ac.jp
3)Institute for Molecular Science,
National Institutes of Natural Sciences,
Myodaiji, Okazaki 444-8585, Japan.
e-mail: nakamura@kba.biglobe.ne.jp
23/Feb./2017

Abstract
It is now confirmed that the Zhu-Nakamura (ZN) theory of nonadiabatic
transition is useful to investigate various nonadiabatic chemical dynamics.
The theory, being one-dimensional, presents a whole set of analytical formulas that enables us to treat the dynamics eﬃciently. It is also quite significant
that classically forbidden transitions can be dealt with analytically. The theory can be combined with the trajectory surface hopping (TSH) method (ZNTSH) and is demonstrated to be useful to clarify the dynamics of not only
simple tri-atomic reactions but also large chemical systems. The whole setof analytical formulas directly applicable to practical systems is summarized∗

Corresponding author. E-mail: shinkoh.nanbu@sophia.ac.jp
Present address:Nakamura Institute of Chemical Dynamics, 3-10-20 Tatsumi Higashi, Okazaki,
Japan
†

1

and the applications to polyatomic systems are illustrated. Examples of polyatomic molecules are H2 SO4 , NH3 , indolylmaleimide, cyclohexadiene (CHD),
and retinal. The Fortran code for the whole set of Zhu-Nakamura formulas isprovided in Appendix for the convenience of a reader who is interested in using
them. The ZN-TSH method can be combined with the QM/MM method toclarify reaction dynamics in the surrounding environment. This is named as
ZN-QM/MM-TSH. The particle-mesh Ewald (PME) method can also be combined with ZN-TSH to clarify reaction dynamics in solutions. This is named
as ZN-PME-TSH. Formulations of these methods are presented together withpractical applications. Examples treated by ZN-QM/MM-TSH are photoisomerization dynamics of retinal chromophore embedded in the protein environment. The diﬀerences in the isomerization mechanisms between rhodopsin
and isorhodopsin are made clear. The faster and more eﬃcient isomerizationof rhodopsin compared to isorhodopsin is nicely reproduced. Examples ofreactions in solutions are photoisomerizations of retinal and CHD. The experimentally observed long life time of the excited state of retinal is reproduced
and is found to be due to the long-range solvation eﬀect. The solvent dependent branching ratios of CHD:hexatriene(HT) are clarified for the ethanol
and hexane solvents by the ZN-PME-TSH method. Both ZN-QM/MM-TSH
and ZN-PME-TSH are thus demonstrated to be promising methods to deal
with a wide range of nonadiabatic dynamics in large chemical and biologicalsystems.

1. Introduction
It is well known now that a variety of chemical dynamic processes, for instancephoto-chemical processes such as photo-isomerization, proceed through conical intersections (CI) of potential energy surfaces[1]. The transition at CI is nothing
but nonadiabatic transition between adiabatic potential energy surfaces[2]. The
significance of nonadiabatic chemical dynamics is well recognized and the field is
expected to be further developed[3]. For instance, various molecular functions suchas photochromism[4], molecular switch[2, 5] and molecular machine[6, 7] may becontrolled by the nonadiabatic transition. The concept of ”nonadiabatic transition” is actually very general and plays crucial roles not only in chemistry but alsoin physics, biology and other fields[2]. Examples are nuclear reactions[8, 9], surface physics [10], solid state physics[11], tunnel junctions[12, 13, 14], quantum dots
[15, 16, 17, 18], solar cell[19], neutrino conversion [20, 21], photosynthesis[22, 23]
and various biological processes[24, 25].
Because of the importance of the transition, there have been developed many
analytical theories [2, 26, 27, 28, 29] since the pioneering works done by Landau[30],
Zener[31] and Stückelberg[32]. There are time-dependent and time-independent
versions in which the latter one is more diﬃcult to formulate and relevant to chemical
processes[2]. Zhu and Nakamura formulated a theory (ZN theory) which provides awhole set of analytical formulas not only for the transition probabilities but also for2the phases induced by the transition[2, 33]. There are two types of potential curve
crossings: one is Landau-Zener type (see Fig.1) in which two diabatic potentials
cross with the same sign of slopes and the other is nonadiabatic tunnel type (seeFig.2) in which two diabatic potentials cross with opposite signs of slopes and apotential barrier is created. The ZN theory not only largely improves the famous
Landau-Zener formula in the classically allowed region (energy ≥ crossing energy),
but also provides the analytical formulas for the classically forbidden transitions,covering the whole range of energy. It should be noted that the classically forbidden
transitions are not negligible at all. The whole expressions can be evaluated withuse of the knowledge of adiabatic potentials and no diabatization is needed. This isvery convenient, because quantum chemical ab initio computations provide adiabatic
potentials uniquely and diabatization procedure cannot be unique.
Since real chemical processes proceed, naturally, in a multi-dimensional space of
nuclear coordinates, the classical or the semiclassical mechanics is usually employed
to describe nuclear motions, while quantum mechanical, i.e., quantum chemical,
methods are used for electronic structure computations. The full quantum mechanical treatment such as the wave packet propagation method or multi-channel timedependent Hartree (MCTDH) method has been devised and successfully applied tovarious processes[34], but it would be a formidable task for real high dimensional
chemical and biological systems, including accurate ab initio computations of global
potential energy surfaces. The idea of surface hopping based on classical trajectories of nuclear motion is usually employed[35, 36]. The fewest switches surface
hopping (FSSH) method proposed by Tully[37] has been widely used with various
modifications[38, 39, 40, 41]. Basically, the time-dependent coupled equations are
solved numerically along classical trajectories with the nonadiabatic coupling employed as the cause of the nonadiabatic transition. The ZN theory can be combinedwith this idea of surface hopping. The nonadiabatic transition is assumed to occur at
the avoided crossing of adiabatic potential energy curves in the direction of nonadiabatic coupling vector and thus the decoherence problem[41] does not arise. Besides,
classically forbidden transitions can be treated, as mentioned above, which is notpossible by other available semiclassical methods based on classical trajectories. The
quantum chemical ab initio computations of potential energy surfaces of high quality
can be carried out on-the-fly. Furthermore, the method can be combined with the
QM/MM method or the particle-mesh Ewald summation method to treat large scale
systems such as photo-isomerization of retinal in protein environment and reactions
in solutions.
This review article is organized as follows: the next section summarizes the whole
set of ZN formulas directly applicable to practical systems. One numerical example
is presented for one-dimensional model potentials to demonstrate how the theory
works well even for classically forbidden transitions. The Fortran code for the formulas is provided in Appendix for the convenience of the reader who is interested inapplying the ZN formulas. In Section 3 various methods to treat chemical dynamics with use of the ZN formulas are explained. The first one is the Zhu-Nakamura
3trajectory surface hopping (ZN-TSH) method by which the over-all reaction probability is calculated from the ratio between the number of reactive trajectories andthe total number of trajectories. The eﬀects of phases are not taken into account inthis treatment. The second one is the Herman-Kluk type semiclassical initial value
representation method[42, 43, 44] with use of the ZN formulas (ZN-HKSCIVR), in
which the phases are incorporated. This is appropriate for evaluating transition
probability amplitude and absorption spectrum. In the third subsection a method
to evaluate reaction rate constant of nonadiabatic dynamics is presented. Especially,
the simple rate constant expression that improves the famous Marcus formula forelectron transfer is provided. The last subsection describes how to treat the effects of environment. Namely, the combinations of the ZN theory with QM/MM
and particle-mesh Ewald summation methods are discussed. In Section 4 numericalapplications of the above mentioned methods are presented for various practical processes such as chemical reactions and absorption spectrum in tri- and tetra-atomic
systems, reaction rate constant of electron transfer, photochemical dynamics of polyatomic molecules, photo-isomerization of retinal in vacuo and in protein environmentand reactions in solutions. Section 5 concludes the paper by discussing future perspectives of the field of nonadiabatic chemical dynamics. Semiclassical molecular
dynamics (MD) simulation method can be developed by taking into account not
only nonadiabatic transitions but also quantum mechanical tunneling eﬀects[45].

2. Zhu-Nakamura formulas[2]A whole set of formulas of the Zhu–Nakamura theory are summarized here. It should
be noted that the theory has the following nice features and can be easily utilized
for various applications:(1) All the probabilities are expressed in simple analytical forms.
(2) All the phases induced by nonadiabatic transitions are provided in compact
analytical forms.
(3) All the basic parameters can be directly estimated from adiabatic potentials on
the real axis.
This means that (i) no non-unique diabatization procedure is needed, (ii) no complex
calculus is necessary, (iii) no nonadiabatic coupling information is required, and (iv)the theory works for whole range of energy and coupling strength.
The formulas presented here contain some empirical corrections so that they cancover practically the whole range of parameters a2 and b2 (see below). Some complex
integrals up to the complex crossing point are simplified. These are explained ineach relevant section.
When the diabatic potentials are available, the two basic dimensionless parameters a2 (≥ 0) and b2 (−∞ < b2 < ∞) are defined by
a2 =h̄2 F (F1 − F2 )
2µ
8VX3
4

(2. 1)

andb2 = (E − EX )F1 − F2,
2F VX(2. 2)

where√µ is the mass, Fj (j = 1, 2) is the slope of the j−th diabatic potential with
F = |F1 F2 |, VX is the diabatic coupling and EX is the energy at the avoided
crossing point (see Eq.(2. 6)). Without loss of generality F1 − F2 is assumed to bepositive. These parameters eﬀectively represent the coupling strength and energy:
a2 ≫ 1, a2 ∼ 1 and a2 ≪ 1 correspond to weak, intermediate, and strong coupling
regime, respectively.
As mentioned above, there are two types of crossings: Landau-Zener type andnonadiabatic tunneling type.
(A) Landau-Zener case (see Fig.1)
In terms of adiabatic potentials Ej (R)(j = 1, 2) these parameters a2 and b2 can be
estimated directly from them asa2 =√d2 − 1

b2 =
where

√

h̄2
(0)

(0)

µ(T2 − T1 )2 (E2 (R0 ) − E1 (R0 ))d2 − 1(0),E − (E2 (R0 ) + E1 (R0 ))/2,
(E2 (R0 ) − E1 (R0 ))/2
(0)

(0)

(2. 3)(2. 4)(0)[E2 (T1 ) − E1 (T1 )][E2 (T2 ) − E1 (T2 )]d =,[E2 (R0 ) − E1 (R0 )]22(2. 5)

where the position R0 is defined as the minimum separation of the two adiabatic
potentials and EX is defined by
EX = [E1 (R0 ) + E2 (R0 )]/2 = E1 (T10 ) = E2 (T20 ).

(2. 6)It should be noted that even when the diabatic potentials are available, it is better to
use above expressions defined in terms of adiabatic potentials, since it is confirmedthat they are more accurate.
(B) Nonadiabatic tunneling case (see Fig.2)
The parameters a2 and b2 originally defined by Eqs.(2. 1) and (2. 2) can be used
when diabatic potentials are available; but these can also be estimated directly fromadiabatic potentials by(1 − γ 2 )h̄2
(2. 7)a2 =µ(Rb − Rt )2 (Eb − Et )
and
b2 =

E − (Eb + Et )/2,(Eb − Et )/25(2. 8)

Potential energy

T2

T1

E
T2(0)

EX

T1(0)

E2(R)
E1(R)
R0

R

Figure 1: Landau-Zener type potential curve crossing.
where Eb (Et ) is the bottom (top) of the upper (lower) adiabatic potential andγ=

Eb − Et.
E2 ([Rb + Rt ]/2) − E1 ([Rb + Rt ]/2)

(2. 9)

When Rb = Rt , γ = 1 and
2

a2 =



2



2

h̄
∂ E2 (R)∂ E1 (R).
−
2
24µ(Eb − Et )
∂R
∂R
R=Rb
R=Rt

(2. 10)2.1. Classically allowed transition
2.1.1. Landau-Zener case (E ≥ E2 (R0 ))(see Fig.1)
The nonadiabatic transition probability for one passage of the crossing point is given
by





pZN = exp −

1/2

2
π 

√
4a|b| 1 + 1 + b−4 (0.4a2 + 0.7)

.

(2. 11)It should be noted that the Landau-Zener formula is simply given by
pLZ = exp[−

6

π
],
4a|b|

(2. 12)E2 (f)Potential energy

E2(R)
T1

T2

T 1l

T 2r

E

Eb
Et
C
D

E
A
B

E1(R)
Rt Rb

R

Figure 2: Nonadiabatic-Tunneling type potential curve crossing.
which goes to zero as |b| → 0. The factor 0.4a2 + 0.7 is originally unity and isconfirmed to be slightly better than unity (see (a) of [33]).
When the phases induced by nonadiabatic transition are needed, the following
transition matrix IX at the avoided crossing point R0 should be used:
(

)

( )

C
A
= IX
,
D
B

(2. 13)where A and B (C and D) are the coeﬃcients of the wave functions at R = R0 + 0(R = R0 − 0). This matrix gives the transition amplitude at the crossing point andexplicitly defined as
)
(√
√
ZN
1 − pZN ei(ψZN −σZN )− pZN eiσ0√,
(2. 14)IX =√ZN
pZN e−iσ01 − pZN e−i(ψZN −σZN )where
ψZN = σZN + ϕS
with

(

δψ
δψ δψ
ln
ϕS = − +πππ)(2. 15)(δψ− arg Γ i
π

)

−

π
.
4

(2. 16)The various parameters are defined as follows:δψ = δZN (1 +5a1/2
10−σZN )a1/2 + 0.8
7(2. 17)

and

∫ R∗

σ0ZN + iδ0ZN ≡

R0√[K1 (R) − K2 (R)]dR ≃2π F−C + iF+C,
4a F+2 + F−2(2. 18)where the final approximate expression is derived by using the linear potential modeland 3-point quadrature (see (i) of [33]),√

Kj (R) =
√√

F± =

2µ
(E − Ej (R)),
h̄2(2. 19)(b2 + γ1 )2 + γ2 ± (b2 + γ1 )√√+(b2 − γ1 )2 + γ2 ± (b2 − γ1 ) ,

(

])

[

0.16bx,b −→ b − √ 4
b +1
√
(
)
0.45 d2C
F− = F− γ2 −→1 + 1.5e2.2bx |bx |0.57
F+C = F+

and

2

2

√
γ1 = 0.9 d2 − 1,bx = b2 − 0.9553,(2. 20)√γ2 = 7 d2 /16.

(2. 21)(2. 22)(2. 23)The parameters σZN and δZN are given as follows:
∫ R0

σZN =

T1

K1 (R)dR −

∫ R0
T2

K2 (R)dR + σ0ZN

(2. 24)andδZN = δ0ZN .

(2. 25)The correction factor in Eq.(2. 17) is introduced empirically (see (c) of [33]).
2.1.2. Nonadiabatic tunneling case (E ≥ Eb )(see Fig.2)
The transition probability for one passage of the crossing point is given as

pZN = exp 
−



2
π 
√
4ab 1 + 1 − b−4 (0.72 − 0.62a1.43 )1/2 
 
.

(2. 26)

The factor 0.72 − 0.6a1.43 is originally unity and is confirmed to be slightly better
than unity (see (b) of [33]). When E ≤ E2 (∞), the overall transmission probability
from left to right or vice versa is the physically meaningful quantity and is given by
4 cos2 (σZN − ϕ̄S )P12 =,4 cos2 (σZN − ϕ̄S ) + (pZN )2 /(1 − pZN )8(2. 27)

where

∫ T2

σZN =

K2 (R)dR

T1(2. 28)with K2 (R) defined by Eq.(2. 19),ϕ̄S = ϕS + h1 ,(2. 29)0.23a1/2
40−σZN ,(2. 30)a1/2 + 0.75
δZN δZN
δZN
δZN
πϕS = −+ln() − arg Γ(i)− ,(2. 31)
π
π √π
π
4
√
6 + 10 1 − b−4
π
√
δZN =
.
(2. 32)16ab 1 + 1 − b−4
The correction factor h1 is introduced empirically (see (c) of [33]). The expressionof δZN is derived by using the linear potential model (see (d) of [33]). It should
be noted that when ψZN ≡ σZN − ϕ̄S = (n + 1/2)π(n = 0, 1, 2, , , , ), P12 becomeszero. Namely, the intriguing phenomenon of complete reflection occurs, when thiscondition is satisfied.
When the phases are needed, the transition matrix IX at E ≥ Eb defined below
should be used,
)
(√
√
ZN
1 − pZN eiϕSpZN eiσ0√,(2. 33)IX =√ZN− pZN e−iσ01 − pZN e−iϕS
where ϕS is the same as Eq.(2. 31). At E > E2 (∞), this IX is used in the same way
as IX in the Landau-Zener case. When the trapping by the upper adiabatic potential
occurs (E ≤ E2 (∞)), we have to use the following reduced scattering matrix (S R )that describes the overall transmission/reflection amplitudes from the entrance tothe exit:( i∆11)1e
U2 ei∆12SR =(2. 34)ei∆22
1 + U1 U2 U2 ei∆12with
U2 = 2i Im(U1 )/(|U1 |2 − 1),(2. 35)
h1 =

where
∫ Rb

∆11 = 2

T1
∫ T2

∆22 = 2(σ0ZN =

Rb − Rt2){

Rb∆12 = σZN ,(2. 36)

K2 (R)dR − 2σ0ZN ,(2. 37)

K2 (R)dR + 2σ0ZN ,(2. 38)}1 [K1 (Rt ) − K2 (Rb )]2
K1 (Rt ) + K2 (Rb ) +,3 K1 (Rt ) + K2 (Rb )(2. 39)√U1 = i 1 − pZN exp[i(σZN − ϕ̄S )](2. 40)and σZN [ϕS ] is given by Eq.(2. 28) [Eq.(2. 31)]. The index of the S R -matrix
corresponds to left or right side of the barrier.
92.2. Classically forbidden transition
2.2.1. Landau-Zener case (E ≤ E2 (R0 )) (see Fig.1)
The overall transition probability from the state 1 to state 2 is given by

where

P12 = 4pZN (1 − pZN ) sin2 (ϕS + σ0ZN ),(2. 41)pZN = [1 + B(σZN /π) exp(2δZN ) − g1 sin2 (σZN )]−1 ,(2. 42)

ϕS and σ0ZN are given by Eq.(2. 16) and Eq.(2. 18),
2πx2x e−2xxΓ2 (x)(2. 43)3σZN
1ln(1.2 + a2 ) − 2 .
πδZNa(2. 44)B(x) =andg1 =The parameters σZN and δZN are given below (Eq.(2. 50)-Eq.(2. 53)).
When the phases induced by the transition are needed, the reduced scatteringmatrix (S R ) which provides the transition amplitude from the turning point on theentrance adiabatic surface to that on the exit adiabatic surface should be used:T
S R = IX
· IX ,(2. 45)where the superscript T means ”transpose” and the transition matrix IX is formallygiven by Eq.(2. 14) with the parameters ψZN and σ0ZN replaced by the following
expressions:ψZN = arg(U ),(2. 46)
where

√

e−δZN




,
Re U = cos(σZN )  B(σZN /π)eδZN − g2 sin2 (σZN ) √B(σZN /π) {

Im U = sin(σZN ) B(σZN /π)e2δZN − g22 sin2 (σZN ) cos2 (σZN )(2. 47)e−2δZN
B(σZN /π)

+ 2g1 cos2 (σZN ) − g2

}1/2

(2. 48)with σ0ZN is given by Eq.(2. 18) and
g2 = 1.8(a2 )0.23 e−δZN .

10(2. 49)

The factors g1 and g2 are originally unity. The first (second) factor of g1 is responsible
for a2 ≫ 1(a2 ≪ 1) (see (a) of [33]). The parameters σZN and δZN are given as follows:(1) At E2 (R0 ) > E > E1 (R0 ),
∫ R0

σZN =
and

T1

∫ T2

δZN =

R0

(2) At E ≤ E1 (R0 ),K1 (R)dR + σ0ZN

(2. 50)K2 (R)dR + δ0ZN .

(2. 51)σZN = σ0ZN

and
δZN = −∫ T1
R0

|K1 (R)|dR +

∫ T2
R0

(2. 52)|K2 (R)|dR + δ0ZN .

(2. 53)2.2.2. Nonadiabatic tunneling caseFirst, the overall transmission probability (P12 ) from left to right or vice versa is thephysically meaningful probability.
(1) At Eb ≥ E ≥ Et
W2,
1 + W2

P12 =
where

(2. 54)[]t3
h2 ∫ ∞b2h3tcos
W = 2/3− 2/3 t − 2/3
dt,
a
3
a
a h4 + a1/3 t
00.38
2h2 = 1 + 2 (1 + b2 )1.2−0.4b ,a√a2 − 3b2 √1.23 + b2
h3 = √
a2 + 3and√h4 = 0.61 2 + b2 .

(2. 55)(2. 56)(2. 57)(2. 58)The factors h2 , h3 and h4 are originally 1/2, 1, and 1, respectively, and are confirmed
to cover the whole range of a2 (see (c) of [33]).
(2) At E ≤ Et
P12 =
with

B(σc /π)e−2δZN√√[1 + (0.5 a2 /[ a2 + 1])B(σc /π)e−2δZN ]2 + B(σc /π)e−2δZN(2. 59)σc = σZN (1 − 0.32 × 10−2/a e−δZN ),(2. 60)

2

∫ Tr

δZN =

1

T1l

|K1 (R)|dR,

11(2. 61)

and

√

σZN =

√

6 + 10 1 − 1/b4

π 18a|b| 2√1+

1 − 1/b4,(2. 62)where the function B(x) is given by Eq.(2. 43). It should be noted that when
a2 → 0, namely, the upper adiabatic potential goes away, we haveP12 =e−2δZN
,
1 + e−2δZN(2. 63)

which agrees with the ordinary single potential barrier penetration probability. The
2factor 0.32 × 10−2/a in Eq.(2. 60) is originally 0 and confirmed to be better (see(c) of [33]). The simplified expression of σZN is derived by using the linear potential
model (see (d) of [33]).
When the phases are needed, what we have to use is the reduced scatteringmatrix (S R ) defined by Eq.(2. 34) with U2 in terms of U1 is given by Eq.(2. 35).
The parameters ∆12 , ∆11 , ∆22 and U1 are defined as follows:(1) At Eb ≥ E ≥ Et
∆11 = σZN − 2σ0ZN∆12 = σZN ,and ∆22 = σZN + 2σ0ZN ,1
σ0ZN = − (Rt − Rb )K1 (Rt )(1 + b2 ),3√
U1 = i[ 1 + W 2 eiϕ − 1]/W,(

ϕ = σZN + arg Γ

1
δZN
+i
2
π

)

(

−

δZN
δZN
ln
π
π
(

)

+

(2. 64)(2. 65)(2. 66)

δZN− h5 ,
π(2. 67))1

4
a0.7 (a0.7 + 0.35)
100b2
2h5 = 0.34
(0.42
+
b
)
2
+
,
(2. 68)a2.1 + 0.73
100 + a2
[
]
√
1
1
2 41
σZN = − √0.057(1 + b ) +(1 − b2 ) 5 + 3b2 ,(2. 69)
3
a2
[
]
√
1
1
1
0.057(1 − b2 ) 4 +δZN = √
(1 + b2 ) 5 − 3b2 ,(2. 70)3
a2
The empirical corrections used for h5 , σZN and δZN are introduced by using the linear
potential model (see (c) of [33]).
(2) At E ≤ Et
∆12 = ∆11 = ∆22 = −2σZN ,(2. 71)√ √ ( )
 0.5 a2
σc −δZNeδZN√
B
Re U1 = sin(2σc )e+√,(2. 72)
 1 + a2πB(σc /π) 

Im U1

v
u
u (Re U1 )2= cos(2σ )t+c

2sin (2σc )1
cos2 (2σc )
12

−

Re U1
1
,
2 sin(σc ) cos(σc )(2. 73)

The parameters
σc and√√ δZN are given by Eqs.(2. 60) and (2. 61), respectively. The
2factor 0.5 a /(1 + a2 ) is originally 0.25 and confirmed to be better (see (c) of[33]).

2.3. Numerical example
One example is shown here (see (g) of [33]). The following 4-state problem definedin atomic units is employed.
V1 (R)
V2 (R)
V3 (R)
V4 (R)====

0.037 exp[−1.3(R − 3.25)] − 0.034,
0.037 exp[−1.3(R − 3.25)] − 0.012,
0.4057[1 − exp[−0.344(R − 3)]]2 − 0.03,
0.4057[1 − exp[−0.344(R − 3)]]2 .

(2. 74)The couplings are given by
2V0, V12 (R) = V34 (R) = 0.
1 + exp(R − 3)
(2. 75)This model system represents some states of O2 (reduced mass = 29377.3 me ) [46].
Fig.3 shows the potential curves in the case of V0 = 0.002 Eh . Fig.4 indicates the

Potential energy / Eh

V13 (R) = V14 (R) = V23 (R) = V24 (R) =

3
4

1
2

R / a0Figure 3: Four-state model potentials of Eq.(2. 74) with V0 = 0.002 Eh . Reproduced
from Ref.(g) in [33] by permission of American Institue of Physics (AIP).
transition probability from state 1 to 2 in the case of V0 = 0.004 Eh in the energy13Transition probability P12

Et(3)

Eb(4)

E / EhFigure 4: Transition probability vs. energy for the potential system of Fig.3 with(4)
(3)
V0 = 0.004 Eh . Etop = 0.01009 Eh , Ebottom = 0.0158 Eh and the starting energy
Estart = 0.00568 Eh . Reproduced from Ref.(g) in [33] by permission of AIP.
range covering the classically forbidden region of the avoided crossing between thehighest and second highest adiabatic potentials. The parameters a2j (j = 1 − 4)corresponding to V0 = 0.004 Eh are given as follows, where j represents the crossing
position (see Fig.3):(a21 , a22 , a23 , a24 ) = (0.98, 0.55, 0.43, 0.24).

(2. 76)It should be noted that the energy of the top (bottom) of the third (fourth) adiabatic
potential is equal to
(3)

(4)

Etop = 0.01009 Eh , Ebottom = 0.0158 Eh(2. 77)and the starting energy of Fig.4 is 0.00568 Eh . It is clearly seen that the present
semiclassical theory works very well, reproducing the detailed resonances. There aremany more examples to demonstrate the accuracy of the theory[33].

3. How to treat dynamics
3.1. Zhu-Nakamura trajectory surface hopping(ZN-TSH) method
The ZN-TSH is the simplest method to estimate reaction probabilities by runningclassical trajectories on adiabatic potential energy surfaces. The probability is given
14

by the ratio of the number of relevant reactive trajectories to the total number oftrajectories. As is well known, if the total number of trajectories is large enough,the ratio gives the proper estimate of the transition probability. The phases along
trajectories are not taken into account in this method.
Nonadiabatic transitions in the vicinity of conical intersections are treated asfollows: the adiabatic potential energy diﬀerence (∆E) between the relevant surface
and the nearest potential energy surface is computed along each classical trajectory
and the minimum ∆E position is detected. That provides the position of transition.
There the nonadiabatic coupling vector is evaluated and the direction of the vector
provides the direction of transition. Two potential energy curves along this direction
are used to evaluate the corresponding nonadiabatic transition probability by using
the relevant Zhu-Nakamura formula. As was explained in the previous section, thenecessary parameters can be estimated from the adiabatic potential curves. Whenthe translational energy along the transition direction is larger than the energy at
the avoided crossing point, the transition is the ordinary classically allowed one. The
relevant transition probability is pZN . If the former energy turns out to be lower
than the latter, the transition is classically forbidden and the corresponding ZhuNakamura probability P12 should be used. The judgment whether the transitionis carried out or not is made by generating a random number prandom . There are
two versions: one is ant-eater procedure and the other is ants procedure [47]. In theformer case the transition is carried out when the Zhu-Nakamura probability is larger
than prandom and in the latter case the trajectory is branched into two in which eachtrajectory is associated with the branching ratio given by the corresponding ZhuNakamura probability. Since the number of branched trajectories can be large in theants procedure, the ant-eater procedure is simpler and is usually employed. Whenthe transition is made, the kinetic energy along the transition direction changes
because of the electronic energy diﬀerence ∆E and the appropriate modificationof kinetic energy of the trajectory should be made. As can be understood in theabove explanation, the nonadiabatic transition is considered only at the minimum
energy separation and the nonadiabatic coupling vector is needed only to decide
the transition direction there. Besides the transition probabilities can be estimatedfrom the analytical formulas. Thus the method is much simpler than solving thetime-dependent coupled diﬀerential equations with use of the nonadiabatic couplings
all the way. It should be noted that the classically forbidden transitions cannot beproperly treated by solving the time-dependent coupled equations based on classicaltrajectories.
If it is too much cpu-time consuming to compute nonadiabatic couplings, thefollowing approximation can be employed[48]. In the vicinity of conical intersection
we can assume the following diabatic Hamiltonian,
(

H=

V1
V

V
V2

)

(∑
Ax
= ∑j j j

j Bj xj

15

∑
)
Bj xj
j
∑
.

−

j Aj xj(3. 1)

Then the nonadiabatic coupling vector ⃗e = {ej } is equal to
ej =

∂V
1
2
(V1 − V2 ) − V ( ∂V
− ∂V
)
∂xj
∂xj
∂xj
2
(∆E)

∑

k (Ak Bj − Aj Bk )xk
,
(∆E)2

=
where
∆E =

√ ∑

(

Ak xk )2 + (k∑

Bk xk )2 .

(3. 2)(3. 3)k

On the other hand, the second derivative of ∆E is given by]∑∑∑∂ 2 ∆E−1 [ ∑
2
2
=
(
A
x
)
B
B
+(
B
x
)
A
A
−(
A
x
)(
B
x
)(A
B
+B
A
)
.
k
k
i
j
k
k
i
j
k
k
k
k
i
j
i
j
∂xi ∂xj(∆E)3 kkkk
(3. 4)
Thus, up to an irrelevant scale factor we have∂ 2 ∆E∝ ei ej .
∂xi ∂xj(3. 5)

This gives a rank-1 matrix and the eigenvector of its only one non-zero eigenvalue
provides the direction of the nonadiabatic coupling vector. This approximation
works well, as will be demonstrated in Section 4.1.

3.2. Herman-Kluk type semiclassical initial value
representation method with use of the Zhu-Nakamura
formulas (ZN-HKSCIVR)
When coherence plays crucial roles in chemical dynamics, the transition amplitude
should be evaluated with various phases taken into account. The best way to do
this is to employ the initial value representation (IVR) devised by Miller[42]. Thiswas combined with the idea of frozen Gaussian propagation approximation[43, 44].
Here we use this method.
The initial wave function ψ(r0 , t = 0) is expanded in terms of frozen Gaussian
wave packets and the latter are propagated along classical trajectories. The total
wave function at time t, ψ(r, t), is expressed as∫
dq0 dp0g(r; qt , pt )Cq0 ,p0 ,t exp[iSq0 ,p0 ,t ] dr0 g ∗ (r0 ; q0 , p0 )ψ(r0 , t = 0),
ψ(r, t) =traj (2π)N(3. 6)
where N is the dimensionality of configuration space, Sq0 ,p0 ,t is the action integral
along the classical trajectory propagated from (q0 , p0 , t = 0) to (qt , pt , t), and∫16

Cq0 ,p0 ,t is the Herman-Kluk pre-exponential factor along the trajectory. The frozen
Gaussian wave packets g(r; q, p) are defined byg(r; q, p) = (2γ N/4) exp[−γ(r − q)2 /h̄ + ip · (r − q)/h̄],πh̄

(3. 7)

where γ is a constant parameter common to all the wave packets. The pre-exponential
factor Cq,p,t is defined by[

Cq,p,t = Det(

∂pt
∂qt
∂qt
i ∂pt ]1/2
+
− 2iγ
+
) .
∂p0 ∂q0∂p0 2γ ∂q0(3. 8)

As is easily understood, the above formulation is for the dynamics on adiabatic
potential energy surface. In order to treat nonadiabatic dynamics the corresponding transition amplitude from Zhu-Nakamura theory is incorporated into the above
formulation. In the case of classically allowed transition from the state j to i, theamplitude Iij from the transition matrix IX should be incorporated (see Eqs.(2. 14)
and (2. 33)). In the case of classically forbidden transition, the reduced scatteringmatrix element SijR should be used (see Eqs.(2. 34) and (2. 45)). Once the total
wave function ψ(t) is known, the transition amplitude for any transition can be simply calculated by taking the projection of the above total wave function ψ(r, t) onto
the final state ψf (r). It should be noted that if a classical trajectory goes around a
conical intersection, the geometrical phase, or Berry phase[49] should be added.
In the case of photo-absorption the total absorption cross section is given by
2παa20 ω ∫ ∞dt exp[iEt/h̄] < ψ(0)|ψ(t) >,σ(ω) =h̄−∞(3. 9)

where α, a0 and ω are fine structure constant, Bohr radius, and light frequency,
respectively. The wave function ψ(0) represents the initial state on the excitedpotential energy surface and is given by|ψ(0) >= µ|ϕi >,(3. 10)where ϕi and µ represent the initial state on the ground potential energy surface
and the transition dipole moment. The total wave function ψ(t) is calculated by
the above mentioned Herman-Kluk type IVR method. When the photo-excitedstate is nonadiabatically coupled with some other excited states through conicalintersections, this ψ(t) can be calculated by the ZN-HKSCIVR method. Photoabsorption of NH3 described in the next Section 4.2 corresponds to this case. Eq.(3.
9) indicates that the total absorption cross section is given by the Fourier transformof the auto-correlation function < ψ(0)|ψ(t) >[50].
In the above treatment of nonadiabatic transition, it is assumed that the energy
at the transition point is given by that of classical trajectory of the center of eachfrozen Gaussian wave packet. More accurate treatment to take into account the factthat each frozen Gaussian has energy distribution is to use the energy normalized17eigenfunction expansion method[51]. The frozen Gaussian wave packet just before
the transition is expanded in terms of the energy normalized eigenfunction on thecorresponding potential energy surface. The expansion coeﬃcient is transformed
into the new one by multiplying the Zhu-Nakamura transition amplitude. The wave
function right after the transition is obtained from the new coeﬃcients and theenergy normalized eigenfunctions on the new potential energy surface. In actual
applications the simplified version with use of the energy of the center of each frozen
Gaussian wave packet would be good enough.

3.3. Reaction rate constant of nonadiabatic dynamics
The nonadiabatic transition state theory (NA-TST) has been developed by many
authors by using the perturbative treatment of diabatic coupling or the LandauZener formula (for instance, [52, 53, 54, 55, 56]). Since electron transfer process is
one of the most typical nonadiabatic dynamics, the theory has also been developed
for the process[57, 58, 59] basically by using the perturbation theory. The ZhuNakamura theory can improve the situation by covering the whole range of couplingstrength from diabatic to adiabatic regime[60]. The semiclassical instanton theory
has also been extended to nonadiabatic dynamics[61].
Here, we employ the ordinary formulation of thermal rate constant, which is
given by1 1 ∫∞
dE exp[−βE]N (E),
(3. 11)k(T ) =
Zr 2πh̄ 0
where T is the temperature with β = 1/κT , Zr is the partition function of reactants and N (E) is the cumulative reaction probability. In the case of nonadiabatic
dynamics, the latter is defined by
N (E) =∫
1
⃗ · p⃗)PZN (E⊥ (⃗q)),d⃗pd⃗qδ[E − H(⃗p, ⃗q)]δ[f (⃗q)](∇f(2πh̄)F −1(3. 12)where F is the number of dimension, f (⃗q) = 0 defines the dividing surface, whichis actually taken to be the crossing seam surface, and E⊥ (⃗q) represents the kinetic
energy in the direction perpendicular to the crossing seam surface. The momentum vector is now expressed as p⃗ = (p∥(1) , , , , p∥(F −1) , p⊥ ), where p⃗∥ and p⊥ are themomenta parallel and perpendicular to the crossing seam surface, respectively. By
introducing the energies E∥(j) = p2∥(j) /2 and E⊥ = p2⊥ /2, where the mass is assumed⃗ · p⃗ = p⊥ and the energy
to be unity or the momenta are mass scaled, and using ∇fconservation,[F −2]
E = H = E⊥ + E∥
+ E∥(F −1) + V (⃗q† )(3. 13)
with
[K]
E∥ ≡

F∑
−1
j=K+1

18

E∥(j) ,

(3. 14)the integrals over momenta can be rewritten as
Int(E) ≡∫dp⊥ p⊥

F∏−1 ∫dp∥(j) δ(E − H)

j=1F∏−2 ∫ E
1 F −1 ∫ E−V∥
dE⊥= (√ )
−V
2
j=1 0

[j−1]

[F −3]

with E∥(F −1) = E∥dE∥(j) √

1
E∥(j)

√

1

E∥(F −1)(3. 15)− E∥(F −2) . Using the beta function,∫ 1

B(p, q) =0

dttp−1 (1 − t)q−1 ,(3. 16)

we obtain∫ E−V

Int(E) =

−V

where

dE⊥ (E − E⊥ − V )(F −3)/2

F∏
−2

1 m
B( , ),
2 2
m=1

F∏
−2

1 m
[Γ(1/2)]F −1
B( , ) =.
2 2
Γ[(F − 1)/2]m=1(3. 17)(3. 18)Then the cumulative reaction probability, i.e., the micro-canonical rate constant, isfinally given by∫
π (F −1)/2d⃗qδ[f (⃗q)]
N (E) = (2πh̄)−F +1 (F −1)/22
Γ[(F − 1)/2]

×

∫ E−V
−V

dE⊥ (E − E⊥ − V )(F −3)/2 PZN (E⊥ ).

(3. 19)Since the order of integration over E and E⊥ can be changed as,
∫ ∞
0

∫ E−V

dE

−V

dE⊥ =

∫ ∞
−V

dE⊥

∫ ∞

dE,

(3. 20)E⊥ +V

the thermal rate constant is finally obtained ask(T ) =

∫
∫ ∞
1 1
1
d⃗
q
dEs exp[−βE]PZN (Es ),†Zr 2πh̄ (8πβh̄2 )(F −1)/2
0(3. 21)

where Es = E⊥ + V† , ⃗q† represents the (F − 1)−dimensional coordinate space and V†
is the potential along the crossing seam. The rate constant can be simply rewritten
as
∫
1
1 1
d⃗q† PZN (β; ⃗q† ),(3. 22)k(T ) =
Zr 2πβh̄ (8πβh̄2 )(F −1)/2∫

where
PZN (β; ⃗q† ) = β

dEs exp[−β(Es − V† )]PZN (Es ; ⃗q† ).
19(3. 23)(0)The minimum energy crossing point (⃗q† ) (MECP) approximation is often employed,(0)in which the transition is assumed to occur only at (⃗q† )[53, 60, 55, 56]. In this casethe rate constant is approximated by
∫
1 1
1
(0)
PZN (β; ⃗q† = ⃗q† ) d⃗q† exp[−β(V† )],k(T ) =
Zr 2πβh̄ (8πβh̄2 )(F −1)/2(3. 24)

As is well known, electron transfer is one of the most famous nonadiabatic chemical processes and the Marcus theory is most celebrated[57, 58]. The Marcus theory
is based on the first order perturbation of electronic coupling. By using the ZN theory, we can improve this theory (see (b) and (e) of [60]). Since the electron transfer
is formulated in the free energy space, we introduce the free energy profile by (seeEq.(3. 12))∫exp[−βFr (ξ)] =⃗ (⃗q)|,d⃗q exp[−βVr (⃗q)]δ[ξ − f (⃗q)]|∇f(3. 25)where the index r represents ”reactant”, ξ is the free energy and ξ = ξ0 = f (⃗q)
defines the crossing seam surface. As in the case of Marcus theory, we assume theparabolic approximation,
1
1
Fr (ξ) = ω 2 (ξ − ξ01 )2 , Fp (ξ) = (ξ − ξ02 )2 + ∆G,(3. 26)2
2
where Fp is the product potential, ξ01 (ξ02 ) is the position of the donor (acceptor)
free energy minimum, and ∆G represents the exothermicity. The crossing point ξ0
is given by1∆Gξ0 = (ξ01 + ξ02 ) + 2(3. 27)2ω δξ0with δξ0 = ξ01 − ξ02 . The Marcus formula is given byH2
kM arcus (T ) = DA
h̄

√

πβ
(λ + ∆G)2
exp[−β],λ4λ(3. 28)where1
λ = ω 2 (δξ0 )2 .
2Then the improved formula is given byk(T ) = κkM arcus (T ),
where
h̄
κ=
2
2πHDA

√

λ
P̄ (β; ξ0 ).
πβ(3. 29)(3. 30)(3. 31)The energy averaged transition probability P̄ (β; ξ0 ) is given byP̄ (β; ξ0 ) = β∫ ∞0dE exp[−βE]PZN (E; ξ0 ).

(3. 32)

Both of Marcus’s normal (NT type) and inverted (LZ type) cases can be treatedwith use of the corresponding ZN probability formulas.
20

3.4. How to treat the eﬀects of environment
The environmental eﬀects are quite crucial to comprehend chemical dynamics in wet
chemistry. The so called dielectric medium model or the polarized continuum model(PCM)[62, 63, 64] has been widely employed to take into account the solvation effects. The solvent in this model is treated as a continuum medium with a constantpermittivity and the solute molecule is put in a certain empirical cavity. The dynamical solvation eﬀect cannot be treated by the PCM model, unfortunately. In thismodel, the liquid structure is discussed in terms of the radial distribution function(RDF). Chandler and coworkers have developed the reference interaction site model(RISM) theory[65, 66] based on the Ornstein-Zernike equation. Hirata and coworkers combined this RISM theory with quantum chemistry. This is called RISM-selfconsistent field (RISM-SCF) approach[67, 68, 69]. Since the solution structure andelectronic states are iteratively optimized in RISM-SCF, the solvation eﬀects arereasonablly well taken into account by using the quantum chemical calculations.
The RDF obtained from the RISM-SCF method is, however, one-dimensional averaged distribution in spite of the fact that the dynamical solvation eﬀect occurs inthree-dimensional space. One can employ the three-dimensional RISM-SCF method,which has been extended to the simulation of dynamics[70]. But this is very muchtime-consuming, unfortunately.
On the other hand, the molecular mechanical (MM) simulation is quite popularin molecular dynamics of condensed matter. Analytical potential functions with
empirical parameters are used in the MM method. One of the most famous force
fields is AMBER[71, 72, 73]. The computational cost of the MM method is proportional to the square of number of atoms, and it could be reduced by the neighbor-list
method[74]. This is quite a reasonable approach for obtaining the statistical information. The target of the conventional MM approach is therefore focused on theelectronic ground state of the system. The conventional formulas and parameters inthis MM approach are designed to reproduce equilibrium molecular conformation in
the electronic ground state. Furthermore, the periodic boundary condition (PBC)
is used to take into account the solution or crystal environment, and the solvent
is represented as a finite system surrounded by its replicas. The minimum-image(MI) convention and the cutoﬀ-radius are also introduced to calculate the totalenergy and force in spite of the fact that the Coulombic interaction goes to zerovery slowly. In order to treat the long-range electrostatic interaction, the Ewald
summation method has been implemented. The MD simulation together with the
Ewald summation method can achieve explicit treatment of the long range Coulombic interaction compared to the implicit way of treatment by the RISM theory (seeFig.5)[75, 76, 77].
3.4.1. QM/MM method
Morokuma et al. have proposed a sophisticated approach for treating large chemicalsystems by using quantum chemical and molecular mechanical (QM/MM) hybrid21

(a)

(b)

M
R (real) \M

M

PRGHO

R\M

I (image) \R

I \R

Minimum image convention

Ewald summation method

Figure 5: The diﬀerence between (a) MI convention and (b) Ewald summation
method. The green and red circles represent particles in the region M and in its
outside, respectively. The lines connecting circles indicate interactions taken into
account.
model known as ONIOM (our Own N-layer Integrated molecular Orbitals and molecular Mechanics) theory[78, 79]. The potential energy of the chemical reaction part
(QM layer) is computed by ab initio quantum chemical calculation, and the remaining nonreactive surrounding part (MM layer) is given by an empirical molecularmechanics (MM) potential. The ONIOM potential energy is defined by the extrapolation form,
QM
MM
MM
EONIOM = Ereal+ Emodel− Emodel,(3. 33)where “real” and “model”, respectively, refer to the full system and the QM layer.
QM
MM
Ereal
is the MM potential energy of the full system, Emodelis the QM energy of theMMmodel layer, and Emodel is the MM energy of the same model layer. In this approach
QM
MM
and Emodel; one is based on the mechanicalthere are two methods to compute Emodel
embedding (ME) method, while the other is based on the electrostatic embedding(EE) method[79].
The use of molecular dynamics (MD) simulation to comprehend chemical dynamics is quite popular nowadays. In the on-the-fly MD (direct dynamics) simulation,
the potential energy and its force are provided by ab initio quantum chemical calculations at each time-step. The ONIOM-MD simulation can be performed in thesimilar way[80, 81]. If the boundary between reaction center and its surrounding
environment is clearly known a priori, the ONIOM-MD simulation method can beapplied in a straightforward way. If this is not the case and the boundary is not pre22

dictable, the exchange-algorithm of solvent molecules (ONIOM-XS) and the adaptive multiscale QM/MM-MD method can be used[82, 83, 84]. Namely, the ONIOM
method enables us to deal with large complex chemical and biological systems.
It is known now that various quantum mechanical eﬀects play key roles even incomplex systems [85]-[94] and that the polarized continuum model (PCM) cannot
necessarily reproduce the eﬀects. York and coworkers have extended the particlemesh Ewald summation method (PME) so as to be combined with QM/MM[95].
The ONIOM scheme[78, 79] is superior to the simple QM/MM scheme with respect
to the calculation of the interaction between QM and MM layers; thus the PME
method should be incorporated into the ONIOM scheme[96, 97].
3.4.2. Particle-mesh Ewald summation method
A.Formulas of force field
The AMBER force field is composed of five terms, which are energies of (i) bond,(ii) angle, (iii) dihedral angle, (iv) van der Waals interaction, and (v) Coulombic
interaction:
V =

∑
bonds

kr (r − req )2 +∑kθ (θ − θeq )2 +angles+∑Vn[1 + cos(nφ − γ)]dihedrals 2

∑
i<j

[

]

[

]

∑ q i qj
Bij
Aij
− 6 +
, (3. 34)
12
Rij
Rij
i<j εRij

where r, θ, ϕ and Rij are the bond length, bond angle, dihedral angle and the
distance between atom i and j, respectively. The parameters req and θeq refer tobond distance and bond angle at equilibrium molecular geometry, kr and kθ are therespective force constants. The angle γ is the oﬀset of rotation, qi and qj are atomic
charges of atoms i and j, and ϵ is 4π times of permittivity of vacuum.
If we use the periodic boundary condition, the Coulombic interaction is rewrittenas1 ∑∑1 ∑ ∑∑
qi q jqi qj
U=+,(3. 35)2 i j̸=i |Ri − Rj | 2 N̸=0 i j̸=i |Ri − Rj − N|
where Ri and Rj are the coordinates of atoms i and j, respectively, and N is theorigin of replica, namely, N = ξa + ηb + ζc, where ξ, η, ζ are integers and a, b, c
are the lattice vectors.
B.Definition of model, real and image layers
The three-layer model is defined in the ONIOM approach as shown in Fig.6. The
layers M (model) and R (real) represent the reaction center and the environment,
respectively. The unit cell in the present treatment is R. The surrounding image
layer denoted by I is provided by the periodic replicas of the real layer (R). We havetherefore the relationship M ⊂ R ⊂ I . The subspace R\M (= R − M ) represents
the outer space of M within R and is the region where the electrostatic embedding23

scheme is employed.

Figure 6: Three-layer model based on the ONIOM scheme. The unit cell is composedof M (model) and R (real), and the image layer I is composed of the periodic replicas
of R.
C.The original ONIOM scheme(i) Coulombic interaction
In the original ONIOM scheme, the PBC is not employed and the total Coulombic
energy is given by a finite sum
∑
1∑
qB
qA
,
(3. 36)U=
2 Ar
AB
B(A̸=B)
where qA and qB are eﬀective charges of particles A and B (∈ R), and rAB = |rB −rA |.
If both A and B belong to M , the interaction is taken into account by quantum
chemical calculations (QM). In the case that A ∈ M and B ∈ R\M , on the other
hand, the pairwise interaction is directly estimated as given above, which is called
the mechanical embedding (ME) method. In the more sophisticated procedure called
electrostatic embedding (EE) method, the following modified one-electron operator
for the i-th electron is used in the electronic Hamiltonian of the subspace M :∑ qB,(3. 37)hQM−
ir
B∈R\M iB
where hQM
is the original one-electron operator. This EE method is employed inithe following sections unless otherwise noted.

24

(ii) van der Waals (vdW) interaction
The vdW interaction is expressed as
UV DW =

∑ [ AAB 12∑(

A B(A̸=B)

rAB

) −(

BAB 6 ]
) .
rAB

(3. 38)As mentioned above, the cut oﬀ procedure is employed and the summation in Eq.(3.
38) is finite in the ONIOM scheme. The parameters AAB and BAB are given by
AMBER force field[71].
The ONIOM energy is finally given by Eq.(3. 33) which is pictorially shown in
Fig.7(a). The ONIOM energy gradient that is needed for classical trajectory calculation is straightforwardly obtained from Eq.(3. 33) by term-wise diﬀerentiation
as
QM
MM
MM
∂EONIOM∂Ereal∂Emodel∂Emodel=+−.
(3. 39)
∂R
∂R
∂R
∂R

D. Particle-mesh Ewald summation approach with use of the ONIOM
scheme(i) Coulombic interaction
The Coulombic interaction defined in the scheme of periodic boundary condition(PBC) given by Eq.(3. 35) is rewritten as
UEwald =

∑
qB
1 ∑∑
qA
,
2 N A
r
ABN
B(̸=A)(3. 40)where the lattice vector of the cell is N, and rABN = |rB − rA − N| is the distance
between particle A ∈ R and the particle B ∈ I\R. This can be divided into
the following two terms according to the particle-mesh Ewald summation (PME)approach[76, 77],
UEwald = Uint + Uext .
(3. 41)
The first term Uint represents the interactions in the unit cell,
Uint =

∑ qB
1∑
,
qA
2 A
r
B(̸=A) AB(3. 42)where A and B ∈ R. This term is the same as the corresponding energy in the
ONIOM approach and is estimated by the quantum chemical calculation within theEE scheme.
On the other hand, Uext represents the interaction between the unit cell (R) andthe replicas(I\R),
1 ∑ ∑ ∑ qB,
(3. 43)qA
Uext =
2 N(̸=0) A
B rABN25

where A ∈ R and B ∈ I\R.
When A ∈ R and B ∈ R, the computational procedure of Coulombic interaction is same as that in the ONIOM scheme. The Coulombic potential energies for
B ∈ I\R can only be taken into account in the PME-ONIOM scheme. In the caseof A ∈ R\M , the Coulombic potential energies are calculated by PME wherever theparticle B is located. The Coulombic potential energy for A ∈ I\R is not calculated
wherever B is located.
(ii) vdW interaction
The vdW interaction energy in PBC is given by
UV DW =∑∑ ∑ [(N

A B(̸=A)AAB 12
BAB 6 ]
) −(
) .
rABN
rABN

(3. 44)Since UV DW decays rapidly with the distance, the summation over N is taken only
over the nearest neighbor replicas, namely |N| = ±1. This is the so called minimum
image (MI) convention. The parameters AAB and BAB are obtained also from thegeneral AMBER force field[71, 72, 73]. When A ∈ R and B ∈ R, the computational
method is same as that in the original ONIOM scheme. The vdW potential energy
for B ∈ I\R is calculated by the MM method. But this interaction could usually
be cut oﬀ, since the subspace I\R is quite far from M.
In the case of A ∈ R\M and B ∈ R, the interaction is calculated by the MM
method. The interaction for B ∈ I\R is also calculated by MM, but the summation
is taken over the minimum image (MI) (|N| = ±1).
The PME-ONIOM procedure is summarized as follows[96]:• E MM
real by the PME approach.
• E MM
model by the direct summation of the interaction energy.
• E QM
model by the quantum chemical calculation with use of the electrostatic embedding scheme.
Since the periodic boundary condition (PBC) provides the additional interactionenergy ∆E between the model (M ) and its outside, the PME-ONIOM energy isfinally given by
QM
MMMM
EPME−ONIOM = Ereal+ Emodel− Emodel+ ∆E(3. 45)which is pictorially shown in Fig.7. The orange-colored part circumscribing the green
QM part is of electrostatic embedding (EE) scheme given by Eq.(3. 37); the quantum
chemical calculation is performed by directly taking into account the Coulombic
interaction between the electron in the subspace M and the charge qB in R\M andI \R, where the polar solvent molecules are described by charge distribution. On
26

(a)

QM with EE
MM with Coulomb int.
MM
Ereal
 Emodel EmodelEONIOM

MM

MM
QM

QM

MM
Coulomb int. in realEE

(by direct sum )(b)
Coulomb int.
QM with EE
MM with Coulomb int.
PME
MM
Eimage Erea Emodewith real  Ereal  Emodelll

EPME ONIOM
MM
QM

PME

MM

MM

QM
EE

PME

Coulomb int. in real(by direct sum)

Figure 7: The pictorial scheme of (a) original ONIOM scheme with electrostatic
embedding (EE) and (b) PME-ONIOM approach with EE. These figures are basedon the practical scheme rather than the above mentioned theoretical concept.
the other hand, the QM calculation without the EE scheme is carried out by themechanical embedding (ME) scheme. Thus, the PME method makes it possible for
the Coulombic interaction to directly aﬀect the QM part.
The diﬀerence between the original ONIOM scheme and the PME-ONIOM scheme
is the existence of the term of ∆E. The polarity of solution can be properly taken
into account by the present EE scheme. In the present PME-ONIOM scheme, thesubspace I\R may be divided into two layers: one is described by the electrostatic
embedding method (R\M ) and the other is described by the mechanical embedding
method (I\R). In this sense the present PME-ONIOM scheme can be combinedwith the three-layered ONIOM. In the same way as in Section 3.4.2 C, the energy
gradient is straightforwardly calculated fromQM
MM
MM
∂Ereal∂Emodel∂Emodel∂∆E∂EPME−ONIOM
=
+
−
+
.
∂R
∂R
∂R
∂R
∂R

27

(3. 46)

4. Clarification of dynamics
4.1. Reactions in tri-atomic systems
Applications of the ZN-TSH method to the reactions in tri-atomic systems are re+ported here by taking DH+
2 and CH2 systems as examples. The DH2 is the wellknown most basic system that has a conical intersection at infinity with the crossing
seam almost parallel to the reaction coordinate. In this kind of simple systems we
can use analytical expressions of potential energy surfaces. The potential energy surfaces used for DH+
2 are the diabatic representation developed by Ushakov et al.[98].
In the actual TSH calculations the adiabatic potentials are used by diagonalizing
the diabatic ones and the crossing seam surface is expressed by analytical functions.
Fig.8 shows numerical results of the cumulative reaction probability for the specified
initial rovibrational state for H2 (v = 1, j = 2) + D+ → HD+ + H[99]. QM is theresult of quantum mechanical numerical solutions of coupled Schrödinger equations.
The new (old) TSH means the TSH calculations with use of the ZN (LZ) formulas.
Since resonances are not the target of interest here, we can say that the ZN-TSH
works well to reproduce the overall reaction probability. It should be noted thatthe nonadiabatic transitions are classically forbidden when the vibrational quantum
number v is less than or equal to 3. This means that the ZN-TSH works well even forclassically forbidden transitions. It should also be noted that the LZ-TSH does notwork even at high energies. This is because the transitions may become classicallyforbidden due to energy transfer in the multi-dimensional space. The cross section
calculations for H+ + D2 and D+ + H2 are carried out with use of the Zhu-Nakamura
theory by Li and Han[100].
The second example is the reaction C + H2 (v = 1, j = 0) → CH + H. This has
a conical intersection in the reaction zone and both LZ and NT types of transitions
occur. The model potential energy surfaces are constructed by using the DIM (diatomics in molecule) method[48]. Here the generalized TSH method explained in
Section 3.1 is employed without using analytical expression of crossing seam surface.
The approximate way of estimating nonadiabatic coupling vector with use of Hessian
is demonstrated to work well. One of the numerical results is shown in Fig.9. Thisshows the cumulative reaction probability for the above mentioned process. There
are many resonances because of the deep potential well in the ground state potential energy surface. Again resonances are not the target of interest. The ZN-TSH
method works acceptably well. The open symbol represents the result with use of
the approximate nonadiabatic coupling vector.

4.2. Absorption spectrum of NH3
The photoabsorption dynamics of ammonia (NH3 ) is investigated by the ZN-HKSCIVR
approach with the full dimensional degrees of freedom taken into account[101]. It is
well-known that the progressions in the umbrella mode (ν2 ) in the S1 ← S0 absorp28

Cumulative probability

Collision energy /eV
Figure 8: Cumulative reaction probability for the specified initial rovibrational state(v = 1, j = 2). The new (old) TSH means the present results with use of the ZhuNakamura (Landau-Zener) formulas. Reproduced from Ref.[99] by permission ofAmerican Institue of Physics (AIP).
tion band are observed in the Franck-Condon region [102, 103, 104]. The scenario
of photodissociation of ammonia molecule is as follows:(1) The excited molecule is trapped by the shallow well just after excitation to S1 .
(2) The ν2 -progression features the broad peaks partly due to tunneling eﬀect by
the predissociation in the shallow well.
(3) The excited molecule passes through the conical intersection between S1 andS0 along the dissociation pathway to NH2 + H.
Since this is a problem of spectrum, ZN-HKSCIVR is a suitable approach. The
Table 1: The vertical excitation energy for S1 − S0 of NH3 with use of two diﬀerent
methods
ThoeryVertical excitation energy
MRCI/aug-cc-pVDZ
6.66
SA-CASSCF/aug-cc-pVDZ
5.93

state-averaged 8-electron and 7-orbital complete active space self-consistent field
29

Cumulative probability

C+H2(v=1, j=0, sur=1) H + CH(sur=2)

Quantum
TSH: Nonadiabatic Coupling Vector
TSH: Difference of Hessian

Total energy / eV
Figure 9: Cumulative reaction probability for the specified initial rovibrational state(v = 1, j = 0). The solid (open) square is the result with use of the proper nonadiabatic coupling vector (approximate vector obtained from Hessian. Reproduced
from Ref.[48] by permission of AIP.
(CAS(8e,7o)-SCF) scheme was used in the ab initio part of the on-the-fly dynamics.
The Hessian matrices are required to propagate the frozen Gaussian wavepackets,as explained in Section 3.2. To estimate the error due to SA-CASSCF in the S1 -S0
excitation energy, the MRCI calculation was performed at the same level of basis
set. The result is listed in Table 1. Since the estimated error is 0.73 eV, the resultfor ν2 -progression is shifted by this amount.
Figure 10 shows the full-dimensional spectrum comprised by forty trajectories.
It should be noted that the nonadiabatic transition along the dissociation path
provides the background broad peak that is missing in the one-surface treatment
done by Guo et al. [104]. Table 2 shows the peak energies of the ZN-HKSCIVR
and experimental spectra. The experimental measurement was performed at 16 K
[102]. The peak energies of the present theoretical results reasonably well agree withthe experimental data. The intervals of the peaks of ND3 are clearly smaller thanthose of NH3 . This indicates that the quantization of vibrational motions are well
represented by the semiclassical method. The energy diﬀerence between the theory
and experiment is approximately 300cm−1 . Another significant isotope eﬀect is thenarrowing of peak widths observed experimentally[102, 103]. In this sense it would
be interesting to improve the present treatment by taking into account the quantum
mechanical tunneling eﬀect. This can be done by applying the newly developed
method[45].
30Figure 10: Photoabsorption spectra of (a) NH3 and (b) ND3

4.3. Reaction rate of electron transfer
The nonadiabatic transition state theory (NA-TST) explained in Section 3.3 hasbeen applied to electron transfer processes. The first example is to demonstratethat the NA-TST based on the Zhu-Nakamura theory works well over the whole
range of electronic coupling strength (see (c) of [60]). Fig.11 displays the calculated
results of five diﬀerent approaches as a function of electronic coupling (HAB ) in theunit of ω at low temperature T = 0.3ω, where ω = 500 cm−1 is the frequency of
harmonic oscillators. The potential system used is a simple one-dimensional modelof two shifted harmonic oscillators, the parameters of which are taken from Ref.[105].
The potential energy at crossing V0 , the eﬀective mass m, and the exothermicity ∆G
are V0 =1500 cm−1 , m = 1836 me . and ∆G = 0, respectively. It should be notedthat the simplicity of the model does not destroy the generality. As is seen, thepresent NA-TST theory based on the ZN theory (solid line) agrees very well withthe exact quantum mechanical results (dots) in the whole range of the couplingstrength. The Bixon-Jortner theory[58] (dot line) works well in the weak coupling
regime (HAB < 0.2ω), but fails at larger couplings. The worse results by the Marcus
theory (dash line) and by the Marcus formula with use of the LZ (dash-dot line)
probability indicate the importance of nonadiabatic tunneling. Namely, not only
the nuclear tunneling but also the coupling between electronic and nuclear motions
play crucial roles.
The second example is the reproduction of the experiment done by Nelsen et.
al[106]. They measured the electron transfer rate of 2,7-dinitronaphthalene in three31

Table 2: The peak positions (ν2 -progression) of theoretical and experimental data
for NH3 and ND3 .

Band
00
21
22
23
24
25
26
27
28
29
210
211
212
213

NH3
Exp./cm−1 Theo.∗ /cm−1
46222
46495
47057
47374
47964
48273
48869
49172
49783
50112
50730
51012
51656
51952
52543
52892
53469
53872
54454
54853
55380
55834
56342
56836
57300
57857
58285
58879

ND3
Band Exp./cm−1 Theo.∗ /cm−1
00
46701
46823
21
47369
47538
22
48052
48233
3
2
48697
48928
24
49357
49602
5
2
50062
50277
26
50748
50971
7
2
51451
51666
8
2
52168
52381
29
52874
53097
210
53582
53832
211
54306
54568
12
2
55054
55304
13
2
55751
56019

*The peaks shifted by 0.73 eV.
kinds of solvents by using ESR. They used the Marcus-Hush theory[57] and the
Bixon-Jortner theory[58], but could not reproduce the experimental results, since
the electronic coupling strength is in the intermediate range. Fig.12 shows thenumerical results. Calculations are made based on the one-dimensional model withuse of the parameters taken from Ref.[106] except for the eﬀective frequency ω,which is 1200 cm−1 , 950 cm−1 , 800 cm−1 for the three solvents, CH3 CN, DMF(dimethylformamide) and PrCN (butyronitrile), respectively. As can be seen, theexperiment is nicely reproduced.
Applications of the present theory to charge transfer in organic semiconductors
and organic molecules for solar cells are made by Nan et al.[107] and by Zhao andLiang[108] with use of quantum chemical electronic structure calculations. The
numerical methods of Monte Carlo and umbrella sampling are employed to attack
complex systems by Zhao et al. (see (d) of [60]). The determination of electronic
structure parameters such as the reorganization energy λ and the coupling strength
by quantum chemical calculations is not easy, unfortunately. Use of spectroscopic
data is sometimes useful as is discussed by Nelsen et al.[109].

32

Rate / a.u.

HAB
Figure 11: Electron transfer rate as a function of electronic coupling strength HAB
in the units of ω at low temperature T = 0.3ω and ω = 500 cm−1 . Dot:quantum mechanical exact result, solid line:present NA-TST based on ZN theory, dot line:BixonJortner theory, Dash line: Marcus theory, dash-dot:LZ formula.

4.4. Photodissociation of H2 SO4
The isotopic imprint of sulfur bearing compounds is an important tool to understand the sulfur cycle of planetary atmosphere. Sulfur isotope analysis has been
employed for many years to study the rise of oxygen on the Earth during the Archean
period[110, 111, 112] and to understand the size, strength and emission features ofstratospheric volcanic eruptions[113, 114]. Sulfur isotopes and non-mass dependent(NMD) fractionation eﬀects have been applied to the study of meteorites[115] and to
the study which links the deep sulfur cycle and the Martian atmosphere[116]. Recent
experimental and theoretical studies[117] have presented the strong evidence of theinvolvement of intersystem crossings, i.e., potential surface crossings, as a potential
origin of NMD fractionation for photoexcited (SO∗2 ) species and for their significancefor explaining the isotopic imprint of the geological record. Furthermore, a recent
study[118] has reported the importance of potential surface crossings for the NMD
eﬀect during the photodissociation of SO.
Although the cyclic reactions have been proposed in several studies, the propertyof each reaction still has some ambiguity even in the reaction channel; for example,
what kind of electronic states (excited/ground state or singlet/triplet state . . . ) of
product are correlated with the parent molecule in photodissociation process, etc.
Concerning the sulfuric acid (H2 SO4 ), the information of electronically excited statesis still an open question, albeit the sulfuric acid is a representative one of the most33

CH3CNDMF
PrCN

Rate(s-1)

1.0 u1011

1.0 u1010

1.0 u109

1.0 u108

1000/T(K)
Figure 12: Arrhenius plot of electron transfer rate of 2,7-dinitronaphthalene in threekinds of solvents. Dots are experimental results of Nelsen et al.[106]. Lines are the results of present improved Marcus formula, Eq.(3. 30). Red line and symbol:CH3 CN,
Green line and symbol:DMF, Purple line and small symbol:PrCN. The eﬀective frequencies used are 1200 cm−1 , 950 cm−1 , 800 cm−1 for CH3 CN, DMF and PrCN,
respectively. Reproduced from Ref.(c) of [60] by permission of World Scientific

34

(a)

(b)

Figure 13: The time-evolution of the potential energy surfaces of a typical trajectory
after photoexcitations of (a) S1 ← S0 and (b) S2 ← S0 . The solid and dash lines
represent the occupied state or unoccupied state, respectively.
oxidized sulfur.
In the present study, nonadiabatic ab initio molecular dynamics (ab initio MD)
simulation was performed to explore the photodissociation dynamics of H2 SO4 using the ZN-TSH method. The state-averaged CAS(8e, 8o)SCF/aug-cc-pVDZ isemployed for the ab initio part. The photodissociation dynamics are carried out
for the two excited states S1 and S2 for 20 fs after photoexcitation by using theon-the-fly ab initio method. If the bond length becomes twice larger than the bondlength of the optimized geometry, it is assumed that the dissociation is over andthe trajectory calculation is stopped. Fifteen and fourteen classical trajectories arepropagated for S1 and S2 excited states, respectively. It is found that nonadiabatic
transitions at potential surface crossings play crucial roles.
There are two dissociation channels, S0 → S1 → S0 and S0 → S2 → S1 (seeEqs.(i) and (ii) shown below).
{′

H2 SO4 (1 A ) + hν −→1

HSO4 (12 A′′ ) + H(2 S) · · · · · · · · · · · · (i)HSO4 (22 A′′ ) + H(2 S) · · · · · · · · · · · · (ii)

The time-profiles of the potential energy surfaces along a typical trajectory are
shown in Fig.13. It is clearly seen that there occurs a nonadiabatic transtion ataround 9.3 (9.7) fs for S1 → S0 (S2 → S1 ).
These reaction channels seem to be quite simple at first glance which provides thehydrogen atoms in the dissociation limit, but it is actually found that nonadiabatic
35

Figure 14: Bond alternation mechanism in the photolysis of H2 SO4 via the first andsecond excited states.
transitions are involved in both channels[119]. Furthermore, the motion correlatedwith the nonadiabatic transitions is SO2 antisymmetric vibration in sulfuric acid,
namely, the bond alternation between single and double SO bonds occurs as shownin Fig.14. Note that the geometrical structures at these nonadiabatic transtions arevery diﬀerent from each other.

36

7

1

9

11

6

2

I
12

10

8

13

5
3

14
4

15

NH2

Figure 15: The (12π) model system of retinal used. The numbers attached show
the positions of carbon atoms.

4.5. Nonadiabatic photoisomerization dynamics of protonated Schiﬀ base retinal in vacuo
Retinal is a core molecule responsible for the vision of vertebrate animals including
humans. The cis→trans isomerization causes electric impulse to be transmitted tothe brain and recognized as light[120, 121]. Photoisomerization of protonated Schiﬀ
base retinal in vacuo is investigated with use of the ab initio ZN-TSH method[122].
In the case of polyatomic molecules it is not recommended to use analytical functionsof potential energy surfaces, since it is a formidable task to prepare global ab initio
potential functions in advance. Instead, it is recommended to use the so called onthe-fly ab initio method[123]. The model system of retinal molecule used here is the12π system in which two methyl groups are removed from the Schiﬀ base.
First, the isomerization of the 11-cis retinal is considered. This is composed of 45
atoms and presents a 129 dimensional problem (see Fig.15). The schematic energy
diagram is shown in Fig.16. There are two conical intersections: one is in between
all-trans form and 11-cis form and the other is in between 9-cis and 11-cis. The
numbers mean the energies in kcal/mol measured from the all-trans ground state.
The blue and red lines indicate the ground and excited states.
The molecule in the 11-cis form in the ground state is photoexcited and 98
classical trajectories are run. The initial geometry in the ground state is optimized
by using the B3LYP method. The isomerization dynamics is investigated by using
the generalized ZN-TSH method explained in Section 3.1 and potential energies are

37

88.2

CI

83.7
61.2

S1 TS
69.0 78.8
S1 min(79.0)
CI

82.4

88.2
68.0
(68.0)

CI

50.6(50.6)

5.9
11-cis

0.5
9-cis

0.0
All-trans

5.9
11-cisFigure 16: Energy diagram for cis-trans photoisomerization of retinal. The blue
and red lines indicate the ground and excited states. Energy is in kcal/mol. At
the conical intersections, numbers without/in parentheses indicate energy when the
C8 − C9 = C10 − C11 torsion is twisted in the clockwise/counterclockwise direction.
CI = conical intersection, TS = transition state, min = minimum.
calculated on-the-fly by using the CASSCF(6,6) with the 6-31G basis set. Both
Molpro 2006.1[124] and Gaussian 03[125] suits of ab initio programs are used.
The dynamics clarified is summarized as follows. The average time of nonadiabatic transition after photoexcitation is 125, 107 and 118 fs for all-trans, 9-cis and11-cis state, respectively. The twisting bond is C11 − C12 double bond. The bondlength becomes longer like a single bond right after excitation, the bond twists to
60 - 100◦ , transition occurs to the ground state, then the bond goes back to doublebond. After the transition the molecule stays in a certain intermediate state for awhile and then finally stabilizes by going through a dynamical barrier which is not
an ordinary potential barrier but a barrier due to the multidimensionality of thesystem. The time evolution of the dihedral angles is shown in Fig.17 for the case of
all-trans formation. The angle ϕ(C11 − C12 ) [ϕ(C9 − C10 )] decreases (increases) inall cases of all-trans, 9-cis and 11-cis, indicating the crankshaft motion.
The twist angle |ϕ| and its velocity d|ϕ|/dt at the hopping geometries are foundto be ϕ| > 90◦ and d|ϕ|/dt > 0 in the case of trans formation and |ϕ| < 70◦ andd|ϕ|/dt < 0 in the case of cis formation. The two conical intersections are found tobe equally accessed and the overall branching ratio turns out to be ∼ 1 : 1 : 2 for
all-trans : 9-cis : 11-cis. In reality the retinal molecule is embedded in rhodopsin
protein environment and the 9-cis formation is blocked by the protein environment.
38

Dihedral angle (I) / degree

Time / fs
Figure 17: A typical time-evolution of dihedral angles of the backbone in the case of
all-trans formation. The dashed line indicates the time of transition. Reproduced
from Ref.[122] by permission of American Chemical Society (ACS).
If the 9-cis formation is counted as trans formation, the branching ratio becomes∼1 : 1 for trans : 11-cis. This agrees well with the experimental quantum yield inprotein environment which is equal to 0.67.
Next, photoisomerization of 9-cis retinal in vacuo is investigated[126]. Ab initio
quantum chemical calculations are carried out in the same way as above by using Molpro 2006.1 and Gaussian 03. The energy diagram is shown in Fig.16. A
10 kcal/mol potential barrier is found with the C8 − C9 = C10 − C11 angle twisted∼ 26◦ in the excited S1 (1 ππ ∗ ) state. This was confirmed by more sophisticatedquantum chemical computation. 71 classical trajectories are run with the time step
0.5 fs.
The overall photoisomerization process clarified by the present study is as follows.
By the photoexcitation S0 → S1 (1 ππ ∗ ) the bond length R9 (C9 = C10 ) stretches toa single bond (> 1.54 Å) and the torsion angle ϕ9 (C9 = C10 ) twists by ±30◦ . The
bond length R9 (the angle ϕ9 ) oscillates between single and double bond character(within ±30◦ ) while being trapped in the S1 potential well. When the potential
barrier is surmounted, the crankshaft motion of ϕ9 and ϕ11 is induced towards theconical intersection and finally the transition to the ground state occurs. In the
crankshaft motion ϕ9 and ϕ11 rotate clockwise and counterclockwise, respectively.
The bond length R9 , being a single bond after surmounting the barrier, becomes adouble bond (∼ 1.34 Å). The angle ϕ9 , being greater than 70◦ before transition to
all-trans form, becomes 0◦ or 180◦ finally on the ground state. Not like in the 11-cis
39

case, there is no two-step relaxation mechanism on the ground state. Because of
the potential well and barrier in the excited state the cis → trans isomerization isslow and less eﬃcient compared to the 11-cis case. This is not due to the so called
pretwist[127]. The mechanism is intrinsic property of 9-cis retinal and not triggeredby the protein environment. The branching ratio of all-trans formation is 0.13 andit becomes 0.22 if the 11-cis is counted as trans formation. This agrees well with theexperimental quantum yield[128]. The X-ray structures are also consistent with thepresent results[129, 130]. The time spent from photoexcitation to the transition to
the ground state is about 211, 292 and 284 fs for all-trans, 11-cis and 9-cis formation,
respectively. This is 2 ∼ 3 times longer than in the case of 11-cis. Human evolved
to choose the faster 11-cis retinal in rhodopsin for visual recognition.
In the above treatments the environment protein is totally neglected. In reality, the 11-cis retinal is embedded in rhodopsin and 9-cis retinal is embedded inisorhodopsin. In order to take into account the eﬀects of environment protein we
can combine the QM/MM and the ZN-TSH methods. This will be discussed in
Section 4.7.

4.6. Photo-chemical dynamics of indolylmaleimide derivatives
Compounds with highly selective fluorescence attract significant interest becauseof various applications in chemical, environmental and biological sciences. It is an
important research subject to design probes with controlled fluorescence yield andselective ON-OFF switching depending on the variations of environmental conditions. Here we investigate photo-chemical dynamics of indolylmaleimide derivatives,
molecular structures of which are shown in Fig.18 [131]. Indole derivatives which
attract much interest in biology are well known as fluorescent compounds. The fluorescence process is governed by the intra-molecular charge transfer in which indoleand maleimide play roles of electron donor and acceptor, respectively.
There are following three anionic species which exist in solution (see Fig.18), (i)′monovalent anion with a deprotonated indole NH group (IM− ), (ii) monovalent′′anion with a deprotonated maleimide NH group (IM− ) and (iii) divalent anion
with doubly deprotonated indole and maleimide NH groups (IM2− ). The photoabsorption maxima in acetonitrile solution are experimentally known to appear at
270 nm and 396 nm, but it is still unknown which species governs the photoemission
process. In order to clarify this, it is crucial to investigate the deactivation processes
due to nonadiabatic transitions by carrying out on-the-fly nonadiabatic ab initio
molecular dynamics computations.
The ab initio quantum chemical calculations are carried out by using the stateaveraged 6-electron 5-orbital complete active space self-consistent field (CAS(6e,5o)SCF) scheme with Dunning’s cc-pVDZ basis set. The package Molpro 2006.1 and2008.1 are used [132]. The initial molecular coordinates and momenta are generated
by random numbers subject to the Wigner distribution, which reflects the vibrational
40

-

IM
O

IM '

H
R2 N 12

O

11

13

10

14

H

O

R2 N

O

2

3

1
9

4

8

6

5

7 N R1

N R1

H

-

IM "

IM2O

R2 N

O

O

N R1
H

R2 N

O

N R1

Figure 18: Molecular structure of indolylmaleimides; the numbers attached to thestructures show the positions of carbons and nitrogens.
wave function at the zero-point energy. The time step of the time-evolution ofclassical trajectories is 0.25 fs and the trajectories are propagated until 60 fs. The
dynamics calculations are performed by using the ZN-TSH method. Table 3 shows
the oscillator strengths at the minimum energy geometry of the S0 and S1 states.
′′
As is clearly seen, the IM− molecule cannot be the candidate of the photoemission species, because the S1 ← S0 and S2 ← S0 photoabsorption transitions are′almost forbidden. Thus, the dynamics of IM− and IM2− are mainly investigated′here. The primary electronic configurations for both S1 and S2 states of IM− aremixed excitations of (54a)1 → (56a)1 and (55a)1 → (56a)1 . The configuration of
the S1 state of IM2− is (55a)1 → (56a)1 and that of the S2 state is (55a)1 → (57a)1 .
′The (55a)1 → (56a)1 or (54a)1 → (56a)1 excitation of IM− and (55a)1 → (56a)1 of
IM2− correspond to the electron transfer from π−MO localized on the deprotonatedindole moiety to π ∗ −MO localized on the maleimide moiety.
′
Figure 19 shows the energy diagram for S0 , S1 and S2 states of IM− and IM2− .
Photoemission from the S1 state is almost forbidden at the S1 minimum energy
structure for both IM−” and IM2− , as can be seen from Table 3. On the other hand,Fig.19 indicates that the nonadiabatic transitions through conical intersections seemto be possible for both transitions S1 → S0 and S2 → S1 . This means that it is′
not easy to determine which species, IM− or IM2− , is the one to emit photons
without carrying out dynamics calculations. Total numbers of classical trajectories′
in the case of vertical excitation to S1 for IM− and IM2− are 40 and 20, while thenumbers of trajectories in the case of excitation to S2 are 34 and 18, respectively.
41

Table 3: Oscillator strengths at the minimum energy geometry in the S0 and S1
states.

′(a) IM−
S0 -S1
S0 -S2′′(b) IM−
S0 -S1
S0 -S2(c) IM2−
S0 -S1
S0 -S2

Oscillator strength
S0 minS1 min
0.48
0.51

0.01
0.70

0.02
0.01

0.01
<0.01

0.50
0.40

<0.01
0.01

(a)

(b)

Figure 19: Relative energies for the S0 , S1 , S2 states and minimum energy conical′intersections (ME-CoIns) between S1 and S0 and between S2 and S1 for (a) IM−
and (b) IM2− where the potential energies at the most stable structures in S0 are
set to be zero.

42

(a)

(b)

(c)

(d)

Figure 20: Time-evolution of energy diﬀerence between S1 and S0 for each trajectory′′of (a)IM− and (b) IM2− after excitation to S1 and (c)IM− and (d) IM2− after
excitation to S2 , respectively.
The dynamics computations are stopped, when trajectories relax to the S0 state.
Figure 20 show time-evolutions of the energy diﬀerences between S1 and S0 states
for all trajectories after the excitation S1 ← S0 [(a) and (b)] and S2 ← S0 [(c) and′(d)] for IM− [(a) and (c)] and IM2− [(b) and (d)]. The energy diﬀerence is taken to
be positive (negative) when the trajectory is on the S1 (S0 ) state. Fig.20(a) shows′that the nonadiabatic transition from S1 to S0 in the case of IM− takes place with
60 percent probability. The average period of this nonadiabatic transition after
excitation is about 36.19 fs. Fig.20 (c) also shows that 76.5% of the excited state
trajectories are deactivated to S0 within 60 fs due to the two sequential nonadiabatic
transitions of S2 → S1 and S1 → S0 . The average period of the nonadiabatic
transition from S1 to S0 is about 26.16 fs. On the other hand, Figs.20 (b) and (d)indicate that in the case of IM2− nonadiabatic transitions between S1 and S0 do notoccur at all, although the transition from S2 to S1 occurs through the corresponding
conical intersection as shown in Fig.19 (b). This means that IM2− does not decay
to S0 non-radiatively and stays on the S1 state, irrespective of the excitation processof S1 ← S0 or S2 ← S0 . This is the most important finding of the present workand suggests that IM2− would be the species to contribute to the photoemission.
′
The two nonadiabatic coupling vectors in the case of IM− for both S2 -S1 and S1 -S0transitions have large components on the five-membered ring of the indole moiety. In
contrast in the case of transition S2 -S1 of IM2− , the vector has the main components
around the maleimide moiety. The stretching and shrinking motions of each ring43

are the driving force for the actual nonadiabatic transitions. The bond length R1 of
N(7)-C(8) and R2 of C(11)-N(12) stretch to 1.4 Å immediately after photoexcitation,
which corresponds to the positions of the S2 -S1 conical intersections. The S1 → S0′transition in the case of IM− occurs due to R1 stretching.

4.7. Photo-isomerization of retinal in opsin environment
Rhodopsin (Rh) is the photosensitive chemical in the outer segment of rod-like cells
in the light-sensing structure of the eye, i.e., retina. The 11-cis retinal chromophorein Rh changes to all-trans form upon exposure to light and bathorodopsin (bathoRh)is produced in the Rh photocycle. Isorodopsin (isoRh) is an analogue of Rh that
contains a 9-cis retinal chromophore embedded in the same opsin environment. The
9-cis retinal is also converted to all-trans form, producing bathoRh in the same way
as in the Rh case. In spite of structural similarity and the same product formation,
Rh and isoRh show quite diﬀerent rates and eﬃciencies of photoisomerization. Isomerization of 11-cis (9-cis) retinal in Rh (isoRh) to all-trans essentially completes
in about 200 fs (600 fs)[133]. The isomerization quantum yield in Rh (isoRh) is
reported to be 0.65 - 0.67 (0.22)[134, 135, 128]. A comparative study of the isomerization of Rh and isoRh is very attractive because of the big diﬀerences in reaction
time scale and quantum yield despite a small structure diﬀerence.
With use of the Landau-Zener formula Schoenlei et al. suggested that the slower
reaction time and lower quantum yield in isoRh are due to slower motion along
the reaction coordinate in the curve crossing region (see (c) of Ref.[133]). From
the X-ray structure analysis Nakamichi et al. proposed that the significant pretwistof the twist angle ϕ9 of C8 − C9 − C10 − C11 is responsible for Rh’s highly eﬃcient
isomerization[127]. In order to clarify the isomerization dynamics of Rh and isoRh
we have employed the ZN-QM/MM-TSH method, in which the QM/MM method
is combined with the ZN-TSH method as explained in Section 3.4.1[136]. The QM
part, namely, quantum chemical calculation is performed for the retinal, the sameas that in the treatment in vacuo, and the opsin environment is treated by MM(molecular mechanics) described by AMBER force field parm96.dat[137] embedded
in Gaussian09[140]. The quantum chemical calculations are carried out with thelevel of CASSCF(6,6) with the 6-31G basis set. The number of atoms taken into
account for Rh (isoRh) is 5591 (5692). A total of 162 classical trajectories arecomputed for both Rh and isoRh with 0.5 fs time step. It took typically 7 weeks for
one trajectory on Intel Core machine with a clock of 2.66 GHz. The quantum yieldor reaction probability, Φ, computed with one standard deviation for the mean is
calculated by the following formula[138],
NR (
1±
Φ=
N

√

)
1
N
(
− 1) ,
N − 1 NR

(4. 1)

where N is the total number of trajectories and NR is the number of reactive trajectories.
44

96.6

95.4
86.4

84.0
65.8
MECI

6.5
Rh

58.9
MECI

54.8
MECI

20.4
9,11-dicisRh

19.3
BathoRh
0.0
IsoRh

Figure 21: Energy diagram for the cis − trans isomerization of Rh and isoRh.
Energy is in kcalmol−1 . The blue and red lines indicate the ground and excitedstates, respectively. MECI = minimum energy conical intersection. The energy for
an MECI is the average of the energies of S1 and S0 .
In this study, two rotation centers (−C9 = C10 − and −C11 = C12 −) are closely
monitored and the isomers are classified according to these torsion angles, giving
rise to the four isomeric species: Rh, isoRh, bathoRh and 9,11-di-cis isomers. Fig.21
shows the energy diagram for the ground (S0 ) and the first excited state (S1 ) spannedby these isomers. The relative ground-state energies of the isomers are diﬀerent fromthose in vacuo[126]. The direction and the relative magnitude of the pretwist is ingood agreement with the X-ray observation[127]. The angle ϕ9 and ϕ11 deviate
rather largely from the perfect trans configuration at the bathoRh equilibrium configuration, which is also in good agreement with X-ray crystallography results[139].
On the other hand, the present excitation energies ranging from 76 to 84 kcal/mol
are larger than the experimental measurements (53 - 59 kcal/mol). Nevertheless,the relative order of the vertical excitation energies is qualitatively reproduced. Thisqualitative reproduction should be suﬃcient for the present purpose. The small potential barrier found in our previous calculation of isoRh in vacuo[126] is not found in
the present calculation. The reason for this is not clear at this stage, unfortunately.
Probably, much higher level of quantum chemical calculations would be needed toclarify this and that is currently beyond our scope. The ground and excited state
surfaces are connected by conical intersections (CI) which exist on the half way
between cis and trans forms along each twist angle. In all cases the double bondin the active torsion coordinate stretches to a single bond at minimum energy CIs.
45

In our previous calculations we found a CI that directly connects 9-cis and 11-cis,but that disappears when the opsin environment is taken into account. Overall, thecalculated results are in qualitative agreement with the experimental measurements,
and thus we may say that the present calculations well simulate the dynamics of thesystems.
The calculated quantum yields are 0.52 and 0.31 for Rh and isoRh, respectively.
The corresponding experimental values are 0.65 and 0.22, as mentioned before[135,128]. The slight overestimation of the isoRh quantum yield may be due to theunderestimation of the disturbance by the opsin environment and irreproducibilityof the S1 barrier. Photoexcitation of Rh only leads to bathoRh and the reactant, andgeneration of isoRh is totally blocked by the opsin environment. Thus, as mentioned
before, the trajectories that lead to 9-cis isomer without the opsin environment are
actually blocked away to generate all-trans bathoRh instead. It is also found thatthe opsin environment totally blocks simultaneous twisting of ϕ9 and ϕ11 for Rh
as well as isoRh, while the twist in the opposite direction of ϕ9 and ϕ11 can beseen in all cases. Photoexcitation of isoRh yields a small amount of 9,11-di-cis
product (4%) in addition to bathoRh, whereas Rh provides the bathoRh alone.
Thus, the rigorous selectivity of Rh production would be another biological reason
for living creatures to choose Rh rather than isoRh. As shown in Fig.22 for the timeevolution of classical trajectories, the population of bathoRh reaches maximum at
187 and 344 fs for Rh and isoRh, respectively. The experimentally observed periods(200 and 600 fs) are reasonably reproduced with a slight underestimation. This
underestimation would be mainly due to the overestimation of excitation energies
and the underestimation in the case of isoRh may be due to the lack of a potential
barrier in the excited state. The fact that the dynamics in the case of Rh is fastand straightforward is responsible for coherence, higher quantum yield and shorter
reaction time. In the case of isoRh, on the other hand, the back-and-forth andcomplicated dynamics leads to slower and less eﬃcient isomerization. In both Rh
and isoRh cases all trajectories exhibit ”wring-a-wet-towel” motion of ϕ11 and ϕ9 .
For all trajectories in Rh as well as most of trajectories in isoRh, the active angle(ϕ11 in the case of Rh and ϕ9 in the case of isoRh) twists in the counterclockwise
direction and the associate angle (ϕ9 in the case of Rh and ϕ11 in the case of isoRh)twists in the clockwise direction. Only exception is a small portion of isoRh (13%),
where the active angle twists in the clockwise direction and associate angle twists
in the couterclockwise direction. All the Rh trajectories go through only one CI
with ϕ11 ∼ −90◦ and slightly clockwise-twisted ϕ9 which leads to bathoRh. The
relaxation of excited isoRh, on the other hand, goes through two CIs, one leading
to bathoRh and the other leading to 9,11-di-cis. As mentioned before there is no
CI directly connecting Rh and isoRh.
In summary, the Rh protein environment ensures that the photoisomerization
proceeds in the right direction and optimizes the all-trans product; whereas thetendency to photoisomerize upon absorption of photon is intrinsic to the central
chromophore itself by Nature’s design. Also the isoRh formation from Rh is totally46

Population

1.0

(a)

BathoRh
Rh
Residual

0.8
0.6
0.4
0.2
0.0
0

100

200

300

400

Time / fs

Population

1.0

BathoRh
IsoRh
9,11-Dicis
Residual

(b)

0.8
0.6
0.4
0.2
0.0
0

100

200

300

400

Time / fs

Figure 22: Time evolution of population in the case of (a) Rh and (b) isoRh. Population is counted only after the product region is reached on the ground state.
Reproduced from Ref.[126] by permission of ACS.
blocked and vice versa. Furthermore, the opsin makes transitions to take place near
the CIs. The ZN-QM/MM-TSH scheme is a promising approach for a wide range
of applications. Wider range of QM part with higher level of quantum chemical
calculations would surely yield much better results.

4.8. Photo-isomerization of (z)-penta-2,4-dieniminium cation(benchmark system) in methanol solution
In order to demonstrate the power of the PME-ONIOM-MD method, namely the
ZN-PME-TSH method, we have investigated the title process[96]. (Z)-penta-2,4dieniminium cation (protonated Schiﬀ base, PSB3) is the minimal model based on
13-cis isomer of 11-cis-retinal protonated Schiﬀ base (PSB11), as shown in Fig.23.
First, the ordinary classical MD simulation was performed to determine the initial molecular coordinates and box size of the central layer under NPT condition
for 10 ns at 1.0 atm and 300 K. The solute molecule is embedded in 36 methanol
molecules and one chloride ion is added to neutralize the unit cell. The charges ofatoms in the unit cell are taken into account in the QM calculation by the electrostatic embedding scheme, and the external charge eﬀect is handled classically by
the PME scheme. In the MM calculation, the force field parameters are taken fromgeneral AMBER force field (GAFF) and the charges of atoms are given by restricted
electrostatic potential charge method. These parameters are generated with Leap
47

(a)

3
2

5
4
1

NH
protein

(b)

3
2

5
4

PSB3 1NH2

Figure 23: The molecular structures of (a) 13-cis isomer of PSB11 and (b) PSB3.
and ANTECHAMBER programs, which are also available in AMBER9. The size ofthe unit cell employed is 2518 Å3 , and the density for solute and solvent is 0.7606
g/cm3 , while the concentration of PSB in solution is 0.0542 g/cm3 .
For the ab initio quantum chemical calculations the state-averaged CASSCF
method with use of the MIDI4* basis set[141] in the package MOLPRO2012.1[142] isused. Other ab initio methods such as MS-MR-CASPT2/cc-pVDZ, SA-CASSCF/ccpVDZ and MRCI/cc-pVDZ are also used to check the validity of the present approach and the SA-CASSCF/MIDI4* was found to be reasonable for the present
purpose. In stead of PME the minimum-image (MI) convention is also used forcomparison. The total number of classical trajectories is 50 in each case and thetime step is 0.5 fs. All the trajectories are propagated until the isomerization isover.
It should be noted that there are two isomerization channels in solution (seeFig.24): (i) through the conical intersection with the central C=C double bondtwisting (channel 1) and (ii) through the conical intersection with the terminal
C=N double bond twisting (channel 2). The channel 2, which does not appear ingas phase[122], is revealed to appear in solution, since the solvent stabilizes therelevant conical intersection and this CoIn becomes energetically accessible[143].
Figures 25 show time-evolutions of the trajectory population for the three channels: (i)S1 state, (ii)S0 state via the C=C twisting isomerization and (iii)S0 state
via the C=N twisting isomerization.
The S1 state lifetimes mainly attributed to the channel 1 (C=C isomerization)and the channel 2 (C=N isomerization) are 3.483 ps and 81.74 fs, respectively in
48

(a)

(b)

Figure 24: The molecular structures and nonadiabatic coupling vectors at the conicalintersections of the two typical trajectories. (a) C=C twisting motion and (b) C=N
inner rotation.

50

50

40

40

number of trajectories(b)number of trajectories

(a)

30
20
10

30
20
10
0

0
0

1

2

3

4

5

6

7

8

9

0

10

50 100 150 200 250 300 350 400 450 500
time / fstime / ps

Figure 25: The population of all trajectories as a function of evolution-time; trajectories which stay in S1 are green. Those having hopped to S0 with C=C isomerization
are blue. Those with C=N inner rotation are brown. (a) and (b) show the whole
profile and the first 500 fs.
49

Table 4: Lifetime by PME-ONIOM-MD and MI-ONIOM-MD simulations
Theoretical approach
PME/SA-CASSCF/MIDI4*
MI/SA-CASSCF/MIDI4*

Channel-1 (C=C)
3.434 ps
0.4642 ps

Channel-2 (C=N)
81.74 fs
91.43 fs

the case of PME-ONIOM-MD/SA-CASSCF/MIDI4* method, as is shown in Table
4. The MI-ONIOM-MD/SA-CASSCF/MIDI4* calculations are also performed andthe corresponding life times are shown in Table 4. The life time 91.43 fs of the innerrotation (channel 2) is similar to that of PME case, but the life time 0.4642 ps of
C=C isomerization (channel 1) is much shorter than that of PME. Actually, thePME result is in accordance with the experimental finding[144] that there are tworeaction paths with 90 fs and 3.5 ps lifetimes. This is one of the most importantresults of the present calculations. The solvent charge clearly aﬀects the behaviorof the nonadiabatic phenomenon of the C=C isomerization, since this isomerizationis characterized by charge transfer. This eﬀect is nicely reproduced by the present
PME method with the long range interactions properly taken into account by theperiodic boundary condition.
A typical time-evolution of the bond lengths and its dihedral angles is analyzed
in the two channels. In the case of channel 1, the bond length gradually stretches
to 1.45 Å by about 2400 fs and the dihedral angle is twisted by about -60 degrees.
This is due to the bond character change from sp2 hybridization to sp3 by losing the
double bond character. This is in accordance with the charge transfer mechanism
mentioned above. In the case of channel 2 (C=N twisting) the bond length of C=N
remains at 1.4 Å up to 75 fs and the dihedral angle rapidly twists by more than -90
degrees.
Finally, it is noted that the life times obtained by PME-ONIOM-MD/MS-MRCASPT2/MIDI4* are 135.9 and 49.12 fs for channels 1 and 2, respectively. These
do not agree with the experimental finding, especially the time scale of channel 1 isquite diﬀerent. This is due to the large bore funnel of conical intersection createdby MS-MR-CASPT2.

50

(b)

(a)

Figure 26: The molecular structures, (a) 1,3-cyclohexadiene (CHD) and (b) 1,3,5cis-hexatriene (HT).

4.9. Photo-isomerization reaction between 1,3-cyclohexatrieneand 1,3,5-cis-hexatriene in vacuo and solution
The reversible structural transformation of photochromic molecules has fascinatedus over the 20th century. The applications are diversifying into various fields ofphotochemistry, material science and biochemistry in the 21st century. Some ofphotochromic systems are based on chemical compounds containing the Kekultype structure defined by six-carbon skeleton. The photoisomerization between
1,3-cyclohexadiene (CHD) and 1,3,5-cis-hexatriene (HT) structures is the simplest
system to represent the feature (see Fig.26) compared to diarylethenes, fulgides andvitamin D [145]-[148]. This prototype photoreaction has also attracted a great dealof interest as a simple model in both theoretical [149]-[163] and experimental [167]-[187] studies. Kosma et al. experimentally observed the CHD/HT isomerizationin gas-phase by the time-resolved photoelectron/mass spectroscopy[185]. CHD wasphotoexcited to S1 by the ultrashort third harmonic pulse at 270 nm. They reported
that molecular wavepacket passes through the conical intersection (CoIn) of 11 B/21 A(S1 /S2 ) at 56 fs after photoexcitation, and then the wavepacket is deactivated into
S0 through the CoIn of 21 A/11 A (S0 /S1 ) at 136 fs. Pullen et al. also made observations of this photoisomerization and the subsequent energy dispersion process incyclohexane solution by using two-color transient absorption kinetic measurements
and novel time-resolved absorption spectroscopy at the 260 - 300 nm[173, 174]. Theirmeasurement uncovered that the sub-picosecond photoproduct formation of CHD
51

or HT occurs on the time-scale of 250 fs, and then the isomerization to other structural isomers of HT, i.e. cZt-HT and tZt-HT, happens in the electronic ground statebecause of the energy dispersion after the ultra-fast photoproduct formation. This
energy stabilization of HT isomers in S0 was observed on the time-scale of 1 - 5 ps.
The product ratio of CHD:HT in cyclohexane solvent was reported as 60:40[175].
The CHD/HT photoisomerization in ethanol solvent was also investigated by transient absorption spectroscopy in the wavelength range between 255 and 450 nm by
Lochbrunner et al. Their study in ethanol solvent suggested that the photoproductformation of CHD or HT occurs within 300 fs[176].
On the other hand, Tamura et al. studied quantum dynamics of this reaction byusing the wavepacket propagation method on a two-dimensional potential energy
surface[160]. The branching ratio between CHD and HT was found to be 50:50 intheir model. The time duration of the first decay 11 B-21 A is approximately 15 - 25
fs, and then the molecular wavepacket goes through the next CoIn of S0 /S1 at 130
fs. Hofmann and de Vivie-Riedle also performed the wavepacket dynamics using theeﬀective Hamiltonian based on the three-dimensional model[155]. The correlation
between the theoretical CHD:HT ratios and locations of CoIn was discussed bythem. Although some other theoretical studies have also been reported so far,any simulation has not yet been carried out with all degrees of freedom taken into
account. In many cases, in order to reduce the computational cost the eﬀective
Hamiltonian in the framework of the reduced dimensional model has been employed.
Usually, the ring opening degree of CHD and the dihedral angle are taken into
account.
In the present study, nonadiabatic ab initio molecular dynamics (ab initio MD)
simulation is performed to explore the CHD/HT isomerization dynamics initiated by
the photoexcitation of CHD to the S1 and S2 states[164]. Particularly, the solvationeﬀects in hexane and ethanol solvents are investigated by the ZN-PME-TSH hybrid
method (see Fig.27)[164, 165, 166].
The MS-MR-CAS(8e,8o)PT2/cc-pVDZ is employed for the QM part. The simulation is performed for 2 ps, because the experimental time-resolved measurementof this isomerization has reported the picosecond-order lifetime. Time-profiles of
QM part (CHD) potential energy in solution phases show smaller fluctuation thanin-vacuo case (see Fig.28). This means that the molecular motions of CHD soluteare confined by the solvent molecules.
The radial distribution function (RDF) is numerically analyzed with respect
to solute-solvent interaction. As shown in Fig.29, ethanol solvent (red line) has
stronger peaks around 5 Å than hexane solvent, which means that there is a strong
cage eﬀect in ethanol solvent. The cage eﬀect would appear because of the induced
dipole moment of CHD in ethanol solvent; the dipole moment in the S2 state is 1.02
Debye [166]. Since the hydorgen-bond network is constructed by hydroxy groupof ethanol solvent, the weakly charged CHD is surrounded by the hydorgen-bondnetwork in the electronic excited state. On the other hand, hexane solvent (blackline) shows an almost uniform distribution. As Dias et al. pointed out[188], the
52

(a)

(b)

Figure 27: CHD molecule (indicated by ball and stick model) surrounded by 304
solvent molecules, (a) hexane and (b) ethanol solvents.
Table 5: The product branching ratio (CHD:HT) determined by nonadiabatic
ab initio MD simulations for photoexcitations to S1 and S2 in each solvent.
Photoexcitation
S1 ← S0S2 ← S0In hexane In ethanol
40:60
80:20
30:70
60:40

In vacuo
60:40
70:30

RDF of hexane solvent has a similar uniform distribution, albeit their simulations
are diﬀerent from our approach. This uniform distribution is generally attributed tothe high flexibility of liquid structure. The strong aﬃnity between hexane and CHD
is due to nonpolar solutions, which caused the high flexibility between solute andsolvent. The eﬃciency of the photochromic conversion of CHD to HT is evaluatedat both S1 - and S2 -excitations in two diﬀerent types of solvents, ethanol and hexane.
The CHD:HT branching ratios of the products are listed in Table 5. The solventdependent branching ratios can be explained also by the flexibility of liquid structureand demonstrate the power of the present method.

53

Potential energy / eV

10

(a)

S0
S1

8

S2

6
4
2
0

Potential energy / eV

10

(b)

S0
S1

8

S2

6
4
2
0

Potential energy / eV

10

(c)

S0
S1

8

S2

6
4
2
0
0

20

40

60

80

100

120

140

Time / fs

Figure 28: A typical time-evolution of potential energies of QM part in (a) hexane,(b) ethanol and (c) vacuum after photoexcitation to S1 .

54

2.5

hexaneethanol

g(r)

2.0
1.5
1.0
0.5
0.0
0

10

20

30

40

50

r / Angstrom

Figure 29: The radial distribution function in each solvent (ethanol and hexane)from the center of mass of solute molecule (CHD); red and black lines obtained by
the respective simulations in ethanol and hexane solvents.

55

5. Future perspectives
The ZN-TSH method was shown to be useful to clarify various nonadiabatic chemical dynamics. The method can be combined with the QM/MM method and with thePME (particle-mesh Ewald) method to investigate nonadiabatic chemical and biological dynamics embedded in environment such as protein environment and solvent.
If necessary, various phases including those induced by nonadiabatic transitions canbe incorporated. The corresponding method is called ZN-HKSCIVR, in which the
Zhu-Nakamura formulas are incorporated into the Herman-Kluk type semiclassical initial value representation method[44]. Although in this review article we havemainly reported the works carried out by the present authors’ research groups, someother groups have also employed the Zhu-Nakamura theory to clarify the nonadiabatic chemical dynamics of polyatomic molecules. Examples are nonadiabatic ab
initio molecular dynamics of photoisomerization in bridged azobenzene by Han’sgroup[189] and in azobenzene by Zhu’s group[190].
The various methods explained in this review article present promising approaches with a wide range of chemical and biological applications. Examples ofuseful future applications would be photochromism[4], molecular switches[5], molecular machines[6], molecular devices[7], photo-synthesis[23], and solar cells[191]. Another interesting field is to clarify the isotopic variations of atmospheric moleculesin geologically old age in the field of geochemistry[192]. In all these systems, nonadiabatic transitions are supposed to play crucial roles.
Quantum mechanical tunneling as another important quantum eﬀect can beincorporated into the above scheme[45]. Along classical trajectories one can detect
caustics which define the boundary between classically allowed and forbidden regions
and an optimal tunneling path can be defined starting from each caustic. By doing
this a semiclassical molecular dynamics simulation method can be developed, inwhich all the important quantum mechanical eﬀects can be taken into account.
Another important application of the above methods is laser control of chemical dynamics[193]. Conical intersections can be created by shining laser, as can beunderstood from the picture of ”dressed state” that means that potential energy
surfaces are shifted up and down by the amount of photon energy. The diabatic
coupling at this artificially created conical intersection is the laser-molecule dipole
interaction and thus by manipulating laser parameters such the frequency, strengthand polarization, not only the wave packet motion but also the nonadiabatic transitions at the conical intersections can be controlled in a desirable way. In such away, not only clarifying the dynamics but also controlling and designing reactions
would become possible.
There are some general problems to be overcome in the present methods described above. The first one is the large amount of cpu-time due to time consuming
ab initio quantum chemical computations. On-the-fly scheme saves time quite a lot
compared to the ab initio computation of global potential energy surfaces in advance. However, accurate evaluation of conical intersections requires high-level of56

quantum chemical computations, which is very much time consuming even in theon-the-fly method. One good solution is to use highly parallelized computers, sinceeach classical trajectory can be run independently. The second one is the dimensionality problem, namely, we have to figure out a good way to reduce the dimensionof the central reaction system not only in the nonadiabatic transition case but also
in the quantum mechanical tunneling case. The third problem is that we usually
don’t know in advance where conical intersections or potential barriers (transition
states) exist in the high dimensional systems. It would probably be useful to employ
a method such as the one proposed by Maeda and Ohno[194] with use of low levelof quantum chemical methods before starting real heavy computations.

Acknowledgments
First of all, we would like to thank all the collaborators of our research groups fortheir valuable contributions to all the works reported here. S.N. was supported inpart by a Grant-in-Aid for Scientific Research (S) (No. 23224013) and Grant-in-Aid
for Scientific Research (A) (No. 26248038) from the Ministry of Education, Culture,
Sports, Science and Technology (MEXT), Japan.

57

Appendix
Fortran code for the Zhu-Nakamura formulas
The attached routines written in Fortan are used for ZN-TSH computation. The
main subroutine for the ZN formulas is “hop.f” at the top of the hierarchy structure
shown below. This routine returns the hopping probability between the two potential
energy surfaces and the corresponding phase induced by the transition.

=1)RUPXODVXEURXWLQH

KRSI

17I

RXWYHFI

ILQG5B/=B17I

GHUGI

GHUGBPXUDISRWGI

RPHJDI

FDOYGYIFDOYGYI

FDOYIJDXOHJJI

FDOYGYIJHWUFI

H[WHUQDOIXQFWLRQ/=I

FJDPDI>UXQBPROSURI@

EIXQIDDEEO]IJDPPOQIFRQVWBO]I

DUJJDPI

FJDPDI>UXQBPROSURI@

Included files:

SDUDLQF

DWRPLQF

Figure 30: Hierarchy of subroutines to compute nonadiabatic transition probabilities
and phases.
Figure 30 shows all the routines and their hierarchy. The “run molpro.f” in red
represents the routine that should be supplied by the user. Here, it is assumed that
the Molpro (version 2012)[142] is employed for the on-the-fly ab initio computation.
Executing “call system(“molpro input.com”)” in Fortran code, the user can run the
ab initio program. If the user-supply main program reads the output file and gives
all data in the subroutine “hop.f”, the code starts to compute the hopping probability. The interfaces of the main subroutine (hop.f) and the user-supply subroutine(run molpro.f) are shown in the sections A and B below.
Some of the main subroutines are explained below with the information of thecorresponding section and page number in the text.
(i) findR0 LZ NT.f finds R0 as shown in Fig.1 in Section 2 on page 6.
(ii) LZ.f computes the quantities in Sections 2.1.1 and 2.2.1 on pages 6-8 and 10-11.
• the transition probability, pZN [Eqs.(2. 11) and (2. 41)]58

• the phases, ψZN [Eqs.(2. 15) and (2. 46)], σZN [Eqs.(2. 15), (2. 50) and(2. 52)] andδZN [Eqs.(2. 25), (2. 51) and (2. 53)](iii) NT.f computes the quantities in Sections 2.1.2 and 2.2.2 on pages 8-9
and 11-12.
• the transition probability, P12 [Eqs.(2. 27), (2. 54) and (2. 59)]
• the phases, ϕ¯S [Eq.(2. 29)], ϕ [Eq.(2. 67)], σZN [Eqs.(2. 28), (2. 69) and(2. 62)],
δZN [Eqs.(2. 32), (2. 70) and (2. 61)],∆12,11,22 [Eqs.(2. 36), (2. 37), (2. 38), (2. 64) and (2. 71)] andU1,2 [Eqs.(2. 40), (2. 35) and (2. 66)]
Note that the main subroutine (hop.f) returns the only hopping probability inthe argument of “prob”, but the phases listed above are shown on standard output(stdout) in computer (The default destination of stdout is the display screen on
computer.). If the phases are required to be stored in arrays, users need to improve
the arguments in “hop.f”, “LZ.f” and “NT.f.”
A. The interface of the main subroutine “hop.f”

subroutine hop(etot1d, ek1d, v1d, dir, r, v, ead, prob,&
natom, nsurf , is0, is1, keyhop, kavail)

Cartesian coordinates are employed for these arguments; the input/output is indicated in parentheses.
etot1d : energy for hopping direction. (input)
ek1d : kinetic energy for hopping direction. (input)
v1d : velocity in the hopping direction. (input)
dir : hopping direction. (input)r : current coordinate in atomic unit (a0 ). (input)
v : current velocity. (input)
ead : adiabatic potential energies in lower and upper potentials (Eh ).
(input)prob : nonadiabatic transition probability. (output)
59

natom : number of atoms. (input)
nsurf : number of potential energy surfaces taken into account. (input)
is0 : index of the relevant state. (input)
is1 : index of the adjacent state. (input)
keyhop : type for LZ transition (E ≥ EX or E < EX ) and NT transition(E ≥ Eb , Eb > E ≥ Et or E < Et ) by 1-2 and 3-5. (output)
kavail : flag indicating whether the transition is classically allowed or not.
If kavail=1, etot1d is not enough to hop vertically. (output)
The sizes of dimension of the arguments are as follows;• dir, r and v arrays are 3(x−, y − and z − components) × matom; natom ≤matom. (matom is the maximum counting number of atoms defined in the
including file, “para.inc.”)• ead array is msurf; nsurf ≤ msurf. (msurf is the maximum counting numberof potential surfaces defined in “para.inc.”)• otherwise, all arrays are single values

B. The interface of the user-supplied subroutine “run molpro.f”&

subroutine run molpro(r c, v c, dvdx c, nacme c, tdm c,
natom, nsurf , isurf1, isurf2, indgna)

Cartesian coordinates are employed for the arguments; the input/output is indicatedin parentheses. This routine controls the ab initio calculation with the argument of“indgna” by 0, 1 and 2.
r c : current coordinates (a0 ). (input)
v c : adiabatic potential energy (Eh ). (output)
dvdx c : adiabatic potential energy gradient. (output)
nacme c : nonadiabatic coupling element (nonadiabatic vector). (output)
tdm c : transition dipole moment. (output)

60

natom : number of atoms. (input)
nsurf : number of potential energy surfaces taken into account. (input)
isurf1 : index of the relevant state. (input)
isurf2 : index of the adjacent state. (input)
indgna controls the output of energy, gradients, non-adiabatic couplings
and so on. (input)• When indgna = 0, only the potential energy is stored in v c.
• When indgna = 1, potential energy and its gradient are store inv c and dvdx c.
• When indgna = 2, potential energy and nonadiabatic couplingelements are stored in v c and nacme c.
The sizes of dimension for the arguments are as follows;• r c, dvdx c, nacme c and tdm c arrays are 3 × matom; natom ≤ matom.
(matom defined in the including file, “para.inc.”)• v c array is msurf; nsurf ≤ msurf. (msurf defined in “para.inc.”)• otherwise, all arrays are single values

C. The ZN formulas software archived in Supplemental Information
The software is archived in “http://pweb.cc.sophia.ac.jp/nanbu lab/zntsh” with zipformat, “src.zip.” The archived file includes all files listed in Fig.30. The user should
modify the “para.inc” and “atom.inc” depending on the molecular system.

61

References[1] W. Domcke, D.R. Yarkony, and H. Koppel, Conical Intersections, Theory,
Computation and Experiment, (World Scientific, Singapore, 2011).
[2] H. Nakamura, Nonadiabatic Transition: Concepts, Basic Theories and Applications (World Scientific, Singapore, 2012).
[3] S. Nanbu, T. Ishida and H. Nakamura, Chem. Science 1, 663 (2010).
[4] Photochromism: Molecules and Systems ed. by H. Dürr and H. Bouas-Laurent(Elsevier, Amsterdam, 2003).
[5] Molecular Switches ed. by B. Feringa (Wiley-VCH, Weinheim, 2001).
[6] Molecular Machines and Motors (Structure and Bonding) ed. by J.P. Sauvage(Springer, Berlin, 2001).
[7] V. Balzani, M. Venturi and A. Credi, Molecular Devices and Machines: A
Journey into the Nanoworld (Wiley-VCH, Weinheim, 2003).
[8] A. Thiel, J. PPhys. G16, 867 (1990).
[9] B. Imanishi and W. von Oertzen, Phys. Rep. 155, 29 (1987).
[10] A. Yoshimori and M. Tsukada (ed.), Dynamical Processes and Ordering on
Solid Surfaces (Springer, Berlin, 1985).
[11] R. Engleman, Non-Radiative Decay of Ions and Molecules in Solids (NorthHolland, Amsterdam, 1979).
[12] G. Blatter and D.A. Browne, Phys. Rev. B37, 3856 (1988).
[13] E. Ben-Jacob and Y. Gefen, Phys. Lett. 108A, 289 (1985).
[14] Y. Gefen, E. Ben-Jacob and A.O. Caldeira, Phys. Rev. B36, 2770 (1987).
[15] J.Q. You and F. Nori, Physics Today 58, 42 (2005).
[16] S.N. Shevchenko, S. Ashhab and F. Nori, Phys. Rep. 492, 1 (2010).
[17] R. Long, N.J. English and O.V. Prezhdo, J. Phys. Chem. Lett. 5, 2941 (2014).
[18] H. Zhu, Y. Yang, K. Hyeon-Deuk, M. Califano, N. Song, Y. Wang, W. Zhang,
O.V. Prezhdo and T. Lian, Nano Lett. 14, 1263 (2014).
[19] J. Liu and O.V. Prezhdo, J. Phys. Chem. Lett. 6, 4463 (2015).
[20] B. Schwarzschild, Physics Today June 17, (1986).
62

[21] B. Schwarzschild and M.H. Pinsonneault, Rev. Mod. Phys. 64, 885 (1992).
[22] D.W. Lawlor, Photosynthesis (Springer, New York, 2001).
[23] R.E. Blankenship, Molecular Mechanisms of Photosynthesis (Blackwell Science, Oxford, 2002).
[24] T.P. Sakmar, S.T. Menon, E.P. Marin and E.S. Awad, Annu. Rev. Biophys.
Biomol. Structure 31, 443 (2002).
[25] M. Chalfie and S.R. Kain, Green Fluorescent Proteins, Properties, Applications, and Protocols (John Wiley and Sons, Hoboken, 2006).
[26] E.E. Nikitin and S.Ya. Umanskii, Theory of Slow Atomic Collisions (Springer,
Berlin, 1984).
[27] M.S. Child, Semiclassical Mechanics with Molecular Applications (Clarendon,
Oxford, 1991).
[28] E.S. Medvedev and V.I. Osherov, Radiationless Transitions in Polyatomic
Molecules, Springer Series in Chemical Physics, 57 (Springer, Berlin, 1994).
[29] W.H. Miller, J. Phys. Chem. A113, 1405 (2009).
[30] L.D. Landau, Phys. Zts. Sov. 2, 46 (1932).
[31] C. Zener, Proc. Roy. Soc. A137, 696 (1932).
[32] E.C.G. Stückelberg, Hel. Phys. Acta. 5, 369 (1932).
[33] (a) C. Zhu and H. Nakamura, J. Chem. Phys. 97, 8497 (1992); (b) ibid, 98,
6208 (1993); (c) ibid, 101, 4855 (1994); (d) ibid, 101, 10630 (1994); (e) ibid,
102, 7448 (1995); (f) ibid, 106, 2599 (1997); (g) ibid, 107, 7839 (1997);(h) ibid, 108, 7501 (1998); (i) ibid, 109, 4689 (1998); (j) Chem. Phys. Lett.
258, 342 (1996); (k) ibid, 274, 205 (1997); (l) H. Nakamura, XXI-ICPEAC
Invited Talks (Amer. Inst. Phys., Woodbury, New York, 2000) p.495; (m) H.
Nakamura, Adv. Chem. Phys. 82, 243 (1992); (n) ibid, 138, 95 (2008); (o) C.
Zhu, Y. Teranishi and H. Nakamura, Adv. Chem. Phys. 117, 127 (2001).
[34] Multidimensional Quantum Dynamics: MCTDH Theory and Applications, ed.
by H-D Meyer, F. Gatti and G.A. Worth (Wiley-VCH, Berlin, 2009).
[35] A. Bjerre and E.E. Nikitin, Chem. Phys. Lett. 1, 179 (1967).
[36] J.C. Tully and R. Preston, J. Chem. Phys. 54, 4297 (1970).
[37] J.C. Tully, J. Chem. Phys. 93, 1061 (1990).
[38] J.C. Tully, J. Chem. Phys. 137, 22A301 (2012).
63

[39] T. Nelson, S. Fernandez-Alberti, A.E. Roitberg and S. Tretiak, Accounts ofChemical Research 47, 1155 (2014).
[40] I. Tavernelli, Accounts of Chemical Research, 48, 792 (2014).
[41] L. Wang, A. Akimov and O.V. Prezhdo, J. Phys. Chem. Lett. 7, 2100 (2016).
[42] W.H. Miller, Adv. Chem. Phys. 101, 853 (1997); W.H. Miller, J. Phys. Chem.
A105, 2942 (2001); K.G. Kay, Annu. Rev. Phys. Chem. 56, 255 (2005).
[43] E.J. Heller, J. Chem. Phys. 75, 2923 (1981).
[44] M.F. Herman, Annu. Rev. Phys. Chem. 45, 83 (1994).
[45] H. Nakamura, S. Nanbu, Y. Teranishi and A. Ohta, Phys. Chem. Chem. Phys.
18, 11972 (2016).
[46] P.S. Julienne and M.J. Krauss, J. Mol. Spectr. 56, 270 (1975).
[47] J.C. Tully and R.K. Preston, J. Chem. Phys. 81, 4352 (1971).
[48] P. Oloyede, G.V. Mil’nikov and H. Nakamura, J. Chem. Phys. 124, 144110(2006).
[49] For instance, D.R. Yarkony, Rev. Mod. Phys. 68, 985 (1996).
[50] E.J. Heller, J. Chem. Phys. 68, 2066, 3891 (1978).
[51] A. Kondorskii and H. Nakamura, J. Chem. Phys. 120, 8937 (2004).
[52] E.J. Heller and R.C. Brown, J.Chem. Phys, 79, 3336 (1983).
[53] J.C. Lorquet and B. Leyh-Nihant, J. Phys. Chem. 92, 4778 (1988).
[54] A.J. Marks and D.L. Thompson, J. Chem. Phys. 96, 1911 (1992).
[55] Y. Okuno and S. Mashiko, Int. J. Quant. Chem. 102, 8 (2005).
[56] J.N. Harvey, Phys. Chem. Chem. Phys. 9, 331 (2007).
[57] R.A. Marcus and N. Sutin, Biochim. Biophys. Acta 811, 265 (1985).
[58] M. Bixon and J. Jortner, Adv. Chem. Phys. 106, 35 (1999).
[59] A.V.Z. Brazykin, P.A. Frantsuzov, K. Seki and M. Tachiya, Adv. Chem. Phys.
123, 511 (2002).

64

[60] (a) Y. Zhao, G. Mil’nikov and H. Nakamura, J. Chem. Phys. 121, 8854 (2004);(b) Y.Zhao, W. Liang and H. Nakamura, J. Phys. Chem. A110, 8204 (2006);(c) Y. Zhao and H. Nakamura, J. Theor. Comp. Chem. 5, 299 (2006); (d) Y.
Zhao, X. Li, Z. Zhen and W.Z. Liang, J. Chem. Phys. 124, 114508 (2006); (e)Y. Zhao, M.M. Han, W.Z. Liang and H. Nakamura, J. Phys. Chem. A111,
2047 (2007).
[61] Y. Teranishi, H. Nakamura and S.H. Lin, AIP Conf. Proc., Theory and Applications in Computational Chemistry: The First Decade of the second Millennium 1456, 119 (2012).
[62] S. Miertu, E. Scrocco, J. Tomasi, Chem. Phys., 55, 117 (1981).
[63] B. Mennucci, WIREs Comput. Mol. Sci., 2, 386 (2012).
[64] J. Tomasi, B. Mennucci, R. Cammi, Chem. Rev., 105, 2999 (2005).
[65] L. J. Lowden and D. Chandler, J. Chem. Phys., 59, 6587 (1973).
[66] L. J. Lowden and D. Chandler, J. Chem. Phys., 61, 5228 (1974).
[67] S. Ten-no, F. Hirata, S. Kato, Chem. Phys. Lett., 214, 391 (1993).
[68] M. Kawata, S. Tenno, S. Kato, F. Hirata, Chem Phys. Lett., 240, 199 (1995).
[69] M. Kawata, S. Tenno, S. Kato, F. Hirata Chem Phys., 203, 53 (1996).
[70] T. Miyata and F. Hirata, J. Comput. Chem., 29, 871 (2008).
[71] D. A. Case et al., AMBER9 University of California, San Francisco (2006).
[72] J. Wang, W. Wang, P. A. Kollman, D. A. Case, J. Mol. Graph. Model., 25,247 (2006).
[73] J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman, D. A. Case, J. Comput.
Chem., 25, 1157 (2004).
[74] L. Verlet, Phys. Rev., 159, 98 (1967).
[75] M. P. Allen and D. J. Tildesley, Computer Simulation in Chemical P hysics,
Kluwer Academic Publishers (1992, p.34).
[76] T. Darden, D. York, L. G. Pedersen, J. Chem. Phys., 98, 10089 (1993).
[77] U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee, L. G. Pedersen,
J. Chem. Phys., 103, 8577 (1995).
[78] L.W. Chung, X. Li, H. Hirao, K. Morokuma, WIREs Comput. Mol. Sci., 00,
1 (2011).
65

[79] T. Vreven, K.S. Byun, I. Komaromi, S. Dapprich, J.A. Montgomery, K. Morokuma, M.J. Frisch, J. Chem. Theory Comput., 2, 815, (2006).
[80] A. Heyden, H. Lin, D.G. Truhlar, J. Phys. Chem. B, 111, 2231 (2007).
[81] X. Li, L. W. Chung, H. Mizuno, A. Miyawaki, K. Morokuma, J. Phys. Chem.
B, 114, 1114 (2010).
[82] T. Kerdcharoen and K. Morokuma, J. Chem. Phys., 118, 8856 (2003).
[83] N. Takenaka, Y. Kitamura, Y. Koyano, M. Nagaoka, Chem. Phys. Lett., 524,
56 (2012).
[84] P. Sripa, A. Tongraar, T. Kerdcharoen, J. Phys. Chem. A, 117, 1826 (2013).
[85] R.S. Becker, K. Freedman, J.A. Hutchinson, L. J. Noe, J. Am. Chem. Soc.,
107, 3942 (1985).
[86] D. Huppertt, P.M. Rentzepis, J. Phys. Chem., 90, 2813 (1986).
[87] H. Kandori, Y. Katsuta, M. Ito, H. Sasabe, J. Am. Chem. Soc., 117, 2669
(1995).
[88] R.S. Becker and K. Freedman, J. Am. Chem. Soc., 107, 1477 (1985).
[89] R. Mathies and C. Brito, Science, 240, 777 (1988).
[90] K.A. Freedman and R.S. Becker, J. Am. Chem. Soc., 108, 1245 (1986).
[91] J. Tittor and D. Oesterhelt, FEBS Lett., 263, 269 (1990).
[92] H. Kandori and H. Sasabe, Chem. Phys. Lett., 216, 126 (1993).
[93] P. Hamm, M. Zurek, T. Röschinger, H. Patzelt, D. Oesterhelt, W. Zinth,
Chem. Phys. Lett., 263, 613 (1996).
[94] G. Zgrabli, K. Voı̈tchovsky, M. Kindermann, S. Haacke, M. Chergui, Biophys.
J., 88, 2779 (2005).
[95] K. Nam, J. Gao, D. M. York, J. Chem. Theo. Comput., 1, 2 (2005).
[96] O. Kobayashi and S. Nanbu, Chem. Phys. 461, 47 (2015).
[97] O. Kobayashi, “Theoretical studies of molecular electron dynamics on atto- topicosecond time scale”, Doctoral Thesis, Chaps. 2 and 3, Sophia Univ. (2016).
[98] V.G. Ushakov, K. Nobusada and V.I. Osherov, Phys. Chem. Chem. Phys. 3,
63 (2001).
[99] C. Zhu, H. Kamisaka and H. Nakamura, J. Chem. Phys. 116, 3234 (2002).
66

[100] B. Li and Ke-Li Han, J. Chem. Phys. 128, 114116 (2008).
[101] T. Murakami, “Nonadiabatic ab initio semiclassical molecular dynamics on themultidimensional photochemical reactions”, Doctoral Thesis, Chap. 5, Sophia
Univ. (2015).
[102] V. Vaida, M.I. McCarthy, P. C. Engelking, P. Rosmus, H. J. Werner, P.
Botschwina, J. Chem. Phys., 86, 6669 (1987).
[103] A. Nakajima, K. Fuke, K. Tsukamoto, Y. Yoshida, K. Kaya, J. Phys. Chem.,
95, 571 (1991).
[104] X. Zhu, J -Y. Ma, D. R. Yarkony, H. Guo, J. Chem. Phys., 136, 234301 (2012).
[105] M. Toplar and N. Makri, J. Phys. Chem. 100, 4330 (1996).
[106] S.F. Nelsen, M.N. Weaver, A.E. Konradsson, J.P. Telo and T. Clark, J. Amer.
Chem. Soc. 126, 15431 (2004).
[107] G. Nan, L. Wang, X. Yang, Z. Shuai and Y. Zhao, J. Chem. Phys. 130, 024704(2009).
[108] Y. Zhao and W.Z. Liang, Chem. Soc. Rev. 41, 1075 (2012).
[109] S.F. Nelsen, A.E. Konradsson, M.N. Weaver, R.M. Stephenson, J.V. Lockard,
J.I. Zink and Y. Zhao, J. Phys. Chem. B111, 6776 (2007).
[110] J. Farquhar, H. Bao, M. Thiemens, Science, 289, 756 (2000).
[111] J. F. Kasting, Science, 293, 819 (2001).
[112] Y. Ueno, M. S. Johnson, S. O. Danielache, C. Eskebjerg, A. Pandey, N.
Yoshida, Proc. Natl. Acad. Sci. U. S. A., 106, 14784 (2009).
[113] J. Savarino, A. Romero, J. Cole-Dai, S. Bekki, M.H. Thiemens, Geophys. Res.
Lett., 30 2131 (2003).
[114] S. Hattori, J.A. Schmidt, M.S. Johnson, S.O. Danielache, A. Yamada, Y.
Ueno, N. Yoshida, Proc. Natl. Acad. Sci. U. S. A., 110, 17656 (2013).
[115] S. Chakraborty, T.L. Jackson, M. Ahmed, M.H. Thiemens, Proc. Natl. Acad.
Sci. U. S. A. 110, 17650 (2013).
[116] H.B. Franz, S.-T. Kim, J. Farquhar, J.M.D. Day, R.C. Economos, K.D. McKeegan, A.K. Schmitt, A.J. Irving, J. Hoek, J. Dottin III, Nature, 508, 364(2014).
[117] A.R. Whitehill, C. Xie, X. Hu, D. Xie, H. Guo, S. Ono, Proc. Natl. Acad. Sci.
U. S. A., 110, 17697 (2013).
67

[118] S.O. Danielache, S. Tomoya, A. Kondorsky, I. Tokue, S. Nanbu, J. Chem.
Phys., 140, 044319 (2014).
[119] Murakami, A. Ohta, T. Suzuki, K. Ikeda, S.O. Danielache, S. Nanbu, Chem.
Phys., 452, 17 (2015).
[120] P.N. Prasad, Introduction to Biophotonics (John Wiley & Sons, 2003).
[121] L.O. Björn, The Photochemical Reactions in Biological Light Perception andRegulation In Photobiology: The Science Light and Life ed. L.O. Björn (Kluwer
Academic, Dordrecht, 2002), p.157.
[122] T. Ishida, S. Nanbu and H. Nakamura, J. Phys. Chem. A113, 4356 (2009).
[123] D. Marx and J. Hutter, Ab initio MD: Theory and Imprementation, in Modern
Methods and Algorithms of Quantum Chemistry, ed. J. Grotendorst, John von
Neumann Institute for Computing, NIC series, 2000, vol.1.
[124] H.-J. Werner and P.J. Knowles, J. Chem. Phys. 82, 5053 (1985).
[125] M.J. Frisch et al., Gaussian03, Revision D.01, Gaussian Inc. (Wallingfold, CT,
2004).
[126] W.C. Chung, S. Nanbu and T. Ishida, J. Phys. Chem. A114, 8190 (2010).
[127] H. Nakamichi, V. Buss and T. Okada, Biophys. J 92, L106 (2007).
[128] J.B. Hurley, T.G. Ebrey, B. Honig and M. Ottolenghi, Nature 270, 540 (1977).
[129] T. Okada, M. Sugihara, A.N. Bondar, M. Elstner, P. Entel and V. Buss, J.
Mol. Biol. 342, 571 (2004).
[130] S. Sekharan, M. Sugihara, O. Weingart, T. Okada and V. Buss, J. Am. Chem.
Soc. 129, 1052 (2007).
[131] T. Murakami, M. Nakazono, A. Kondorskiy, T. Ishida, S. Nanbu, Phys. Chem.
Chem. Phys., 14, 11546 (2012).
[132] H.-J. Werner et al., MOLPRO Version 2008.1; A Package of ab initio Programs, Cardiﬀ University, Cardiﬀ, UK, 2008.
[133] (a) R.W. Schoenlein, L.A. Peteanu, R.A. Mathies and C.V. Shank, Science
254, 412 (1991); (b) L.A. Peteanu, R.W. Schoenlei, Q. Wang, R.A. Mathies
and C.V. Shank, Proc. Natl. Acad. Sci. U.S.A. 90, 11762 (1993); (c) R.W.
Schoenlei, L.A. Peteanu, Q. Wang, R.A. Mathies and C.V. Shank, J. Phys.
Chem. 97, 12087 (1993).
[134] H.J. Dartnall, C.F. Goodeve and R.V. Lythgoe, Proc. Roy. Soc., A156, 158(1936).
68

[135] J.E. Kim, M.J. Tauber and R.A. Mathies, Biochemistry 40, 13774 (2001).
[136] W.C. Chung, S. Nanbu and T. Ishida, J. Phys. Chem. B116, 8009 (2012);
ibid Chem. Lett. 40, 1395 (2011).
[137] W.D. Cornell, P. Cieplak, C.I. Bayly, I.R. Gould, K.M. Merz, D.M. Ferguson,
D.C. Spellmeyer, T. Fox, J.W. Caldwell and P.A. Kollman, J. Am. Chem.
Soc. 117, 5179 (1995).
[138] H.R. May, In Dynamics of Molecules and Chemical Reactions ed. R.E. Wyatt
and J.Z.H. Zhang (Marcel Dekker, Inc., New York, 1996), p.589.
[139] H. Nakamichi and T. Okada, Angew. Chem. Int. Ed. 45, 4270 (2006).
[140] M. J. Frisch et al., Gaussian 09, Revision C.01, Gaussian, Inc., Pittsburgh,PA, (2009).
[141] Y. Sakai, H. Tatewaki and S. Huzinaga, J. Comp. Chem. 2, 108 (1981).
[142] H.-J. Werner et al., MOLPRO Version 2012.1; A Package of ab initio Programs, Cardiﬀ University, Cardiﬀ, UK, 2012.
[143] T. Mori, K. Nakano, S. Kato, J. Chem. Phys. 133, 064107 (2010).
[144] S.L. Logunov, L. Song and M.A. El-Sayed, J. Phys. Chem. 100, 18586 (1996).
[145] M. Irie, Chem. Rev., 100, 1685 (2000).
[146] M. Irie, S. Kobatake, M. Horichi, Science, 291, 1769 (2001).
[147] M. Murakami, H. Miyasaka, T. Okada, S. Kobatake, M. Irie, J. Am. Chem.
Soc., 126, 14764 (2004).
[148] K.L. Kompa, R.D. Levine, Proc. Natl. Acad. Sci. U. S. A., 98, 410 (2001).
S. Mahvidi, S. Takeuchi, S. Kusumoto, H. Sato, T. Nakagawa, Y. Yokoyama,Org. Lett., 18, 5042 (2016).
[149] P. Celani, S. Ottani, M. Olivucci, F. Bernardi, M.A. Robb, J. Am. Chem.
Soc., 116, 10141 (1994).
[150] P. Celani, F. Bernardi, M.A. Robb, M. Olivucci, J. Phys. Chem., 100, 19364(1996).
[151] M. Garavelli, P. Celani, M. Fato, M.J. Bearpark, B.R. Smith, M. Olivucci,
M.A. Robb, J. Phys. Chem. A, 101, 2023 (1997).
[152] M. Garavelli, P. Celani, F. Bernardi, M.A. Robb, M. Olivucci, J. Am. Chem.
Soc., 119, 11487 (1997).
69

[153] M. Garavelli, F. Bernardi, M. Olivucci, T. Vreven, S. Klein, P. Celani, M.A.
Robb, Faraday Discuss., 110, 51 (1998).
[154] M. Garavelli, C.S. Page, P. Celani, M. Olivucci, W.E. Schmid, S.A. Trushin,
W. Fuß, J. Phys. Chem. A, 105, 4458 (2001).
[155] A. Hofmann, R. de Vivie-Riedle, J. Chem. Phys., 112, 5054 (2000).
[156] A. Hofmann, R. de Vivie-Riedle, Chem. Phys. Lett., 346 299 (2001).
[157] L. Kurtz, A. Hofmann, R. de Vivie-Riedle, J. Chem. Phys., 114, 6151 (2001).
[158] A. Hofmann, L. Kurtz, R. de Vivie-Riedle, Appl. Phys. B, 71, 391 (2000).
[159] H. Tamura, S. Nanbu, H. Nakamura, T. Ishida, Chem. Phys. Lett., 401, 487(2005).
[160] H. Tamura, S. Nanbu, T. Ishida, H. Nakamura, J. Chem. Phys., 124, 084313(2006).
[161] C. Nonnenberg, S. Grimm, I. Frank, J. Chem. Phys., 119, 11585 (2003).
[162] M. Boggio-Pasqua, M. Ravaglia, M. J. Bearpark, M. Garavelli, M.A. Robb, J.
Phys. Chem. A, 107, 11139 (2003).
[163] J.B. Schönborn, J. Sielk, B. Hartke, J. Phys. Chem. A, 114, 4036 (2010).
[164] A. Ohta, “Theoretical study of chemical dynamics for condensed phase”, Doctoral Thesis, Chaps. 2 and 3, Sophia Univ. (2017).
[165] A. Ohta, O. Kobayashi, S.O. Danielache, S. Nanbu, Chem. Phys., 459, 45(2015).
[166] A. Ohta, O. Kobayashi, S.O. Danielache, S. Nanbu, Chem. Phys., 485, 45(2017).
[167] M.O. Trulson, G.D. Dollinger, R.A. Mathies, J. Am. Chem. Soc., 109, 586(1987).
[168] M.O. Trulson, G.D. Dollinger, R.A. Mathies, J. Chem. Phys., 90, 4274 (1989).
[169] P.J. Reid, S.J. Doig, R.A. Mathies, Chem. Phys. Lett., 156, 163 (1989).
[170] P.J. Reid, S.D. Wickham, R.A. Mathies, J.Phys.Chem., 96, 5720 (1992).
[171] P.J. Reid, S.J. Doig, S.D. Wickham, R.A Mathies, J. Am. Chem. Soc., 115,4754 (1993).
[172] M.K. Lawless, S.D. Wickham, R.A. Mathies, Acc. Chem. Res., 28, 493 (1995).
70

[173] S.H. Pullen, L.A. Walker II, B. Donovan, R.J. Sension, Chem. Phys. Lett.,
242, 415 (1995).
[174] S.H. Pullen, N.A. Anderson, L.A. Walker II, R.J. Sension, J. Chem. Phys.,107, 4985 (1997).
[175] S.H. Pullen, N.A. Anderson, L.A. Walker II, R.J. Sension, J. Chem. Phys.,
108, 556 (1998).
[176] S. Lochbrunner, W. Fuß, W.E. Schmid, K. L.Kompa, J. Phys. Chem. A, 102,
9333 (1998).
[177] W. Fuß, T. Schikarski, W.E. Schmid, S. Trushin, K.L. Kompa, Chem. Phys.
Lett., 262, 675(1996).
[178] S.A. Trushin, W. Fuß, T. Schikarski, W.E. Schmid, K.L. Kompa, J. Chem.
Phys., 106, 9386 (1997).
[179] W. Fuß, W.E. Schmid, S.A. Trushin, J. Chem. Phys., 112, 8347 (2000).
[180] N. Kuthirummal, F.M. Rudakov, C.L. Evans, P.M. Weber, J. Chem. Phys.,
125, 133307 (2006).
[181] R.C. Dudek, P.M. Weber, J. Phys. Chem. A, 105, 4167 (2001).
[182] J.D. Cardoza, R.C. Dudek, R.J. Mawhorter, P.M. Weber, Chem. Phys., 299,
307 (2004).
[183] C.Y. Ruan, V.A. Lobastov, R. Srinivasan, B.M. Goodson, H. Ihee, A.H. Zewail, Proc. Natl. Acad. Sci. U. S. A., 98, 7117 (2001).
[184] H. Ihee, V.A. Lobastov, U.M. Gomez, B.M. Goodson, R. Srinivasan, C.Y.
Ruan, A.H. Zewail, Science, 291, 458 (2001).
[185] K. Kosma, S.A. Trushin, W. Fuß, W.E. Schmid, Phys. Chem. Chem. Phys.,
11, 172 (2009).
[186] S. Adachi, M. Sato, T. Suzuki, J. Phys. Chem. Lett., 6, 343 (2015).
[187] M.P. Minitti, J.M. Budarz, A. Kirrander, J.S. Robinson, D. Ratner, T.J. Lane,
D. Zhu, J.M. Glownia, M. Kozina, H.T. Lemke, M. Sikorski, Y. Feng, S.
Nelson, K. Saita, B. Stankus, T. Northey, J.B. Hastings, P.M. Weber, Phys.
Rev. Lett., 114, 255501 (2015).
[188] A.M.A. Dias, R.P. Bonifácio, I.M. Marrucho, A.A.H. Pádua, M.F. Costa
Gomes, Phys. Chem. Chem. Phys., 5, 543 (2003).
[189] A.H. Gao, B. Li, P.Y. Zhang and K.L. Han, J. Chem. Phys. 137, 204305(2012).
71

[190] L. Yu, C. Xu and C. Zhu, Phys. Chem. Chem. Phys. 17, 17646 (2015).
[191] Organic Solar Cells: Materials, Devices, Interfaces, and Modeling ed. Q. Qiao
(CRC Press, Boca Raton, 2015).
[192] J.A. Schmidt, M.S. Johnson, S. Hattori, N. Yoshida, S. Nanbu and R. Schinke,
Atmos. Chem. and Phys. 13, 1511 (2013).
[193] A. Kondorskiy, S. Nanbu, Y. Teranishi and H. Nakamura, ȷ. Phys. Chem.
A114, 6171 (2010).
[194] S. Maeda, K. Ohno and K. Morokuma, J. Phys. Chem. A113, 1704 (2009);
ibid, J. Chem. Theory Comput. 6, 1538 (2010).

72