- Email: [email protected]
Numerical simulation on the dynamics of photoinduced cooperative phenomena in molecular crystals
Computer Physics Communications 180 (2009) 1489–1494
Contents lists available at ScienceDirect
Computer Physics Communications www.elsevier.com/locate/cpc
Numerical simulation on the dynamics of photoinduced cooperative phenomena in molecular crystals Kunio Ishida a,∗ , Keiichiro Nasu b a b
Corporate Research and Development Center, Toshiba Corporation, 1 Komukaitoshiba-cho, Saiwai-ku, Kawasaki 212-8582, Japan Institute of Materials Structure Science, KEK, Graduate University for Advanced Study, and CREST JST, 1-1 Oho, Tsukuba, Ibaraki 305-0801, Japan
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 18 October 2007 Received in revised form 20 November 2008 Accepted 25 March 2009 Available online 28 March 2009 PACS: 71.35.Lk 05.45.-a 64.70.Nd 78.20.Bh
We develop a new simulation method to study the dynamics of initial nucleation processes of photoinduced structural change of molecular crystals. In order to describe the nonadiabatic transition in each molecule, we employ a model of localized electrons coupled with a fully quantized phonon mode, and the time-dependent Schrödinger equation for the model is numerically solved. By applying a meanfield approximation in solving the Schrödinger equation, the calculation method is quite efficient on parallel computing systems. We show that coherently driven molecular distortion plays an important role in the successive conversion of electronic states which leads to photoinduced cooperative phenomena. © 2009 Elsevier B.V. All rights reserved.
Keywords: Nonequilibrium phase transition Nonadiabatic transition Parallel computing
1. Introduction Amongst various photoinduced phenomena which have been observed in solids, liquids, and/or molecules, cooperatively induced change of structural, magnetic, or ferroelectric properties [1–5] has particularly attracted our attention. These photoinduced cooperativity are considered to have a common mechanism, and many experimental and/or theoretical studies have been presented to make it clear [6–8]. We have proposed a model of molecular crystals in order to describe the initial nucleation processes of these phenomena [9–11], and showed that the nonadiabaticity of electronic transitions is a key to understand the temporal behavior of quantum-mechanical states. Dynamics of nonadiabatic transitions has been studied since the pioneering works by Landau [12] and Zener [13], and, in particular, the transition rate in general cases was analytically obtained [14]. They mainly focused on the wavefunctions before/after nonadiabatic transition, and hence the time evolution of wavefunctions itself is not within the scope of these studies. On the other hand, the dynamics of nonadiabatic processes has been considered to be important in, for example, photochemical reactions [15], and hence computational methods for nonadiabatic dynamics have been pro-
*
Corresponding author. Tel.: +81 44 549 2119. E-mail address: [email protected] (K. Ishida).
0010-4655/$ – see front matter doi:10.1016/j.cpc.2009.03.013
©
2009 Elsevier B.V. All rights reserved.
posed by many authors [16]. Since, however, these methods are based on the molecular dynamics simulations, they require that the atomic coordinates/momenta are treated as classical variables. Hence, these methods are valid after the decoherence of the atomic degrees of freedom takes place. As a result they mentioned the reaction yield or the absorption rate after various nonequilibrium processes. It has, however, been pointed out that an elementary process of photoinduced cooperative phenomena is consecutive switching of relevant potential energy surfaces (PESs) and that the coherence of the wavefunction of electrons/atoms during electronic transition should be taken into account. In other words, the wavefunction at every moment should be pursued to understand the dynamics of the whole processes, and hence, not only the bifurcation rate of the wavefunction after nonadiabatic transition but also the wavefunction as a function of time is required. This aspect of the nucleation dynamics shows that the conventional methods are not suitable for the theory of photoinduced cooperativity. Furthermore, cooperative action of electrons and phonons consists of an essential part of the photoinduced cooperativity, and thus we require a method to calculate dynamical properties of large aggregates of molecules. In this paper, we show a new computational method to study the initial dynamics of the photoinduced nucleation processes. By introducing a mean-field approximation to take into account the intermolecular interactions, we naturally derived an efficient algorithm to solve the time-dependent Schrödinger equation for a many-body Hamiltonian of localized electrons and
1490
K. Ishida, K. Nasu / Computer Physics Communications 180 (2009) 1489–1494
phonons. This method is also effective for parallel computing and thus suitable for calculations on the recent high performance computing systems. The organization of the paper is as follows: in Section 2 a model of molecular crystals is introduced and the method of calculation is described. In Section 3 the calculated results are shown, and Section 4 is devoted to discussion and conclusions. 2. Model and method 2.1. Model of molecular crystals As we mentioned in Section 1, nonadiabatic transitions between electronic states play a key role in the dynamical aspects of photoinduced cooperative phenomena. In this paper, we focus on the initial dynamics of a photoexcited state in interacting molecules, fully quantizing the relevant vibration modes. However, the dimension of the Hilbert space of the whole system drastically increases by quantizing the atomic variables, which means that the current computational resources are insufficient for numerical calculations on these systems. As a first step to overcome this difficulty, we employ a simplest model which deserves to describe the photoinduced nucleation processes. We consider molecules arrayed on a regular lattice, and all the electrons in the system are assumed to be localized in each molecule. Only two electronic levels coupled with a single vibration mode are taken into account per molecule. Each molecule has crossing diabatic PESs regarding the two electronic levels, which are relevant to the nonadiabatic transitions to drive the nucleation processes. The nonadiabaticity in the dynamics is taken into account via a “spin-flip” interaction between two electronic states, i.e., the interaction strength between two electronic levels are assumed to be a constant. We also apply the harmonic approximation to the vibration modes. We mention that this model is also known as a model to discuss the relaxation dynamics of, e.g., photoisomerization of molecules and studied by many authors [17]. As for the intermolecular interaction, we take into account three types of intermolecular interactions described below: 1. Vibrational interaction which couples the vibration modes of different molecules. It is also responsible for the dispersion of the relevant phonon mode. 2. The dipole–dipole interaction between electrons in excited molecules The interaction strength is taken up to the first order of the molecular distortions. 3. Electron-vibrational interaction which describes the distortion of molecules induced by the excited-state dipole in the adjacent ones. We have pointed out that the nucleation process is triggered only when the second and/or the third types of the above interactions in two- or higher-dimensional systems are taken into account [9–11]. This situation is different from the one-dimensional case in which the vibrational interaction is sufficient to induce the domain growth [7]. In any case, the nucleation process occurs through excitation energy transfer between molecules, and the Hamiltonian in the present study is described by:
H=
p 2 r
r
+
2
+
ω2 ur2 2
+
2h¯ ω3 sur + εh¯ ω + s2 h¯ ω nˆ r + λσxr
αω2 (ur − β nˆ r )(ur − β nˆ r )
r ,r
− V − W (ur + u r ) nˆ r nˆ r ,
(1)
Fig. 1. Schematic view of the model. Circles denote the molecules with two electronic states and a vibrational mode. Adiabatic potential energy surfaces for an individual molecule are shown in the inset.
where pr and ur are the momentum and coordinate operators for the vibration mode of a molecule at site r , respectively. The electronic states at site r are denoted by |↓r (ground state) and |↑r (excited state) and σir (i = x, y , z) are the Pauli matrices which act only on the electronic states of the molecule at site r . nˆ r denotes the density of the electron in |↑r which is rewritten as nˆ r = σzr + 1/2. The second sum which gives the intermolecular interaction is taken over all the pairs on nearest neighbor sites, where the Coulomb interaction between excited state electrons are modified by molecular distortion. The vibrational period of an individual molecule is denoted by T = 2π /ω in the rest of the paper. A schematic view of the present model in the case of square lattice is shown in Fig. 1. This figure shows that the two diabatic PESs for an individual molecule cross each other, and that the nonadiabatic coupling constant λ acts to separate them into two adiabatic PESs. The parameters for the intra/intermolecular interactions are also shown in Fig. 1. The canonical quantization of the vibrational degrees of freedom is performed. The vibronic states are denoted by |nσ r where n is the index of the vibration mode and σ = ↑, ↓, and described by [18–20]
† n a
|n↓r = √r
n!
|n↑r = e
†
|0r ⊗ |↓r ,
s(ar −ar )
(2)
† n a
√r
n!
|0r ⊗ |↑r ,
(3)
†
where ar denotes the creation operator of vibrational quanta of a molecule at site r . The coordinate of the molecule is labeled by r . We mention that the basis set for the vibronic states is composed of the Fock states shown in Ref. [19]. The phonon dispersion relation of the vibration mode is given by
Ω(k) = ω 1 + 2α (cos kx + cos k y ),
(4)
where k = (kx , k y ) denote the reciprocal lattice vector of the square lattice, and the lattice constant is taken to be unity. We note that this Ising-like model is similar to the one to study the thermodynamical properties of the Jahn–Teller effect [21], though the nonequilibrium dynamics of the excited states in the model has not been understood.
K. Ishida, K. Nasu / Computer Physics Communications 180 (2009) 1489–1494
1491
2.2. Method of calculation We obtain the numerical solutions of the time-dependent Schrödinger equation for the Hamiltonian (1) which is written as ih¯
∂
Φ(t ) = H Φ(t ) . ∂t
(5)
When we take into account up to the Mth vibronic states on each PES per molecule, the dimension of the Hilbert space of the whole system D is equal to (2M ) N where N is the total number of the molecules. Since we cannot divide the Hamiltonian into a direct sum of those for each molecule due to the intermolecular interactions, we cannot block-diagonalize the Hamiltonian, i.e. the wavefunction cannot be decomposed to a tensor product of those for each molecule. This makes it quite difficult to solve Eq. (5) numerically, since D > 1015 even for M = 16 and N = 10, for example. However, the aim of the present study is the numerical calculation on the collective phenomena induced in the system, which means that the total number of the molecules N is required to be at least 103 for quantitative discussion. To avoid this difficulty, we apply a mean-field approximation to Eq. (5), i.e., the effects of molecules in the nearest neighbor sites are substituted by their average values. This approximation allows us to decompose the wavefunction of the total system |Φ(t ) into a product of the wavefunctions at each molecule, i.e.,
Φ(t ) = φ(t ) ⊗ φ(t ) ⊗ · · · ⊗ φ(t ) . r r r 1
(6)
N
2
i h¯
Then, each component |φ(t )r obeys the following equation: i h¯
∂
φ(t ) r = ∂t
pr2
ω2 u 2
r + + 2 2
+ λσxr φ(t ) r
+
Fig. 2. A schematic view of the partition of the system for computation. All the circles denote the molecules on a square lattice. The molecular crystal is divided into blocks of molecules by the dashed lines and each block is assigned to a single processor of the computer system. Only the mean-field values (ˆnr and ur ) for the molecules at the perimeter of the blocks (filled circles) are communicated between processors.
dcnr ↑ dt
2h¯ ω3 sur + εh¯ ω + s2 h¯ ω nˆ r
+ −
αω2 (ur − β nˆ r )(ur − β nˆ r )
− V − W (ur + u r ) nˆ r nˆ r φ(t ) r ,
(7)
where the sum is taken over the nearest neighbor sites of r . After canonical quantization of pr and ur , Eq. (7) is rewritten as
† †
∂
φ(t ) r = h¯ ω ar ar + s ar + ar + ε + s2 nˆ r + λσxr φ(t ) r ∂t
+ αω2 γ ar† + ar − β nˆ r ur − βˆnr (8)
By expanding the wavefunction of each molecule by its vibronic states as (10)
n=0 σ =↑,↓
we obtained a set of differential equations for the coefficients cnσ as ih¯
dcnr ↓ dt
= nh¯ ωcnr ↓ + h¯ +
r U nm cm ↑
m
αω2 γ
√
ncnr −1↓ +
n + 1cnr +1↓
r
× u r − βˆnr ,
√
(11)
√
ncnr −1↓ +
√
√
ncnr −1↑ +
√
n + 1cnr +1↓ u r − βˆnr
V − W − u r cnr ↑
min(m,n) i =0
(9)
∞
φ(t ) = cnr σ |nσ r , r
αω2 γ
† U mn = m|e s(a −a) |n s2
where γ = h¯ /2ω . The intermolecular interaction is taken into account through the average value over the decomposed wavefunction and
Xr = r φ(t ) Xr φ(t ) r .
m
n + 1cnr +1↑
ˆnr ,
(12)
where
= e− 2
r
† − V − W γ ar + ar + u r nˆ r ˆnr φ(t ) r ,
− Wγ
r
r
i h¯
√ √
= h¯ ω n + ε + s2 cnr ↑ + s ncnr −1↑ + n + 1cnr +1↓ r + h¯ ωλ U mn cm ↓
√
m!n!
i !(m − i )!(n − i )!
sm+n−2i .
(13)
Eqs. (11) and (12) are suitable for parallel computing, since the wavefunction |φ(t )r itself is not transferred between each other during the calculation. Instead, only the average values of ur and ˆnr are communicated between CPUs, which reduces the amount of data transfer drastically. Hence, our implementation to solve the time-dependent Schrödinger equation is shown as follows: 1. We divide the whole system into blocks of molecules shown in Fig. 2. Each block contains the same number of molecules and is assigned to a CPU or a processor element (PE) of the computer system. 2. The wavefunctions cnr σ for the molecules are stored in the main memory of each processor, and their time evolution is calculated by applying the Runge–Kutta method to Eqs. (11) and (12). 3. The values of ˆnr and ur for the molecules at the perimeter of each block are transferred between processors at each time step. Thus, when a block is composed of n × n molecules, only 4(n − 1) real numbers are communicated per block. By locating these variables at a continuous block of the main memory, we can also reduce the computational cost regarding communication between processors, since the MPI coding is simpler and easier in this case (ref. filled circles in Fig. 2).
1492
K. Ishida, K. Nasu / Computer Physics Communications 180 (2009) 1489–1494
on the implementation of the MPI library and/or the system architecture, and that it varies between 0 and 1. Hence, N 0 ∝ N α where 1/3 α < 1. In any case, N 0 is an increasing function of N as shown in Fig. 3, and thus the number of molecules per block increases as N increases. This feature is also important to consider the simulation of three-dimensional systems, and will be discussed in the next section. 3.2. Dynamics of photoinduced nucleation processes We have pointed out that the population of the excited electronic state |↑r is suitable for understanding the dynamics of the initial nucleation processes [9–11]. We show in Fig. 4 the population of |↑r for 32 × 32 sites defined by
N (r , t ) = Φ(t ) nˆ r Φ(t ) , Fig. 3. Normalized execution time as a function of the number of processors. The system size is 32 × 32, 64 × 64, and 128 × 128.
In this way, we have made it possible to handle more than 104 molecules by the present method, which is sufficient to discuss the cooperative effect induced by photoexcited states. 3. Calculated results 3.1. Scalability of the computational method Although the amount of communication is reduced by the mean-field approximation, the present computation is not completely scalable as in general cases. In the rest of the paper we chose the values of the parameters as: ε = 2.3, s = 1.4, V = 1.1, W = 0.2, α = 0.1, β = 0.2, and λ = 0.2. Although those values are typical for organic molecules as for electron-vibration coupling [18] and the intermolecular Coulomb interaction [17], the other parameters are not easy to determine their values either from theoretical calculations or experimental results. We only mention that the order of magnitude for the parameters is estimated referring to those for typical organic materials. As for the time scale of the photoinduced nucleation, we mention that the vibration period of a single molecule is ∼200 fs, i.e. T ∼ 200 fs. The time step is taken to be T /36000, which is sufficient to obtain the numerical convergence of the results. Fig. 3 shows the execution time of the simulation with 6000 time steps on a 32 × 32, a 64 × 64, and a 128 × 128 lattice. The calculations were performed on the HA8000 system at the Research Center for Computational Science, National Institutes of Natural Sciences. In each plot the computation time is normalized by that for a single processor. As the number of processors P increases, the amount of communication between them increases, and thus the trade-off between the increase of CPU time and the decrease of communication burden results in the existence of the optimal value of the processors which minimizes the total computation time. On the other hand, the amount of communication is dependent on the size of the block introduced in the last section, i.e., the perimeter-to-section ratio in a square block is ∼1/n for an n × n block, and thus the relative weight of communication increases as the size of the block decreases. Hence, the optimal number of processors N 0 depends on the size of the system, and it increases as the system is larger. In the present case N 0 = 8, 16, and 32 for a 32 × 32 lattice, a 64 × 64 lattice, √ and a 128 × 128 lattice, respectively, which shows note that the total amount of communication that N 0 ∝ N. We √ is ∼ N (n − 1)/n2 ∼ N P in the present case, and that the time to Thus, if we ascalculate the time evolution in each block is ∼ N / P . √ sume that the total computation time is ∼ N / P + a( N P )x with a and x being certain √ constants, we obtain x = 2/3, which satisfies the relation N 0 ∝ N. We also note that the value of x depends
(14)
for t = 0, 5T , 10T , and 15T . Initially three molecules which form an v-tromino (see the inset of Fig. 4) are excited to the Franck– Condon state, while the others are set to the ground state. Fig. 4 shows that the number of molecules in the excited electronic state increase surrounding the initially excited molecules. Those molecules will constitute a photoinduced domain observed in many experiments [1–5] and thus the present calculation describes the initial processes of the photoinduced cooperative phenomena, i.e., photoinduced nucleation triggered by an injected excited state. The sum of the excited state population N total (t ) = r N (r , t ) indicates the measure for the growth rate of the photoinduced domain [9]. Fig. 5 shows that N total (t ) slightly decreases for 5T after photoexcitation, and turns to increase after this period. The domain size behaves superlinearly as a function of t after domain growth started. This feature is understood by the picture that the growth of the domain is predominantly driven by coherent motion of molecular distortions rather than diffusion processes. As the vibrational coherence is lost, diffusion process becomes more important and the domain growth √ will slow down to make the radius of the domain increase as ∼ t. Since vibrational coherence survives for a few ps in typical organic molecules [18], the present calculation is valid only in the time range studied in this paper, and the decoherence of the vibrational states should be taken into account to study the growth dynamics of the photoinduced domain in a longer time scale. 4. Discussion and conclusions In this paper we show a computational method to study the dynamics of the initial nucleation processes of photoinduced cooperative phenomena. By treating the intermolecular interactions by the mean-field approximation, we obtain an efficient algorithm to calculate the time-evolution of photoexcited states of more than 104 molecules on parallel computing systems. In this method we divide the crystal into blocks of molecules, and assign each block to one of the processors. When intermolecular interactions are restricted to those between adjacent molecules, the communicated data are only the average values (mean-field) of relevant physical properties at the perimeter of the blocks, and hence the amount of communication in this method is relatively smaller as the size of the block is larger. As a result we found that the optimal number of processors to minimize the computational cost depends on the size of the crystal. This feature is particularly important to consider higher-dimensional (i.e., three-dimensional) cases. In three-dimensional cases we need to take into account typically 106 (= 100 × 100 × 100) molecules or more, and hence the optimal number of processors N 0 should be much larger than that in the present case. Thus, the present method has a remarkable advantage in computation on superparallel computer systems and/or grid computing systems.
K. Ishida, K. Nasu / Computer Physics Communications 180 (2009) 1489–1494
1493
Fig. 4. Population of the excited electronic state N (r , t ) on 32 × 32 lattice for (a) t = 0, (b) t = 5T , and (c) t = 10T , and (d) t = 15T .
Fig. 5. N total (t ) as a function of time. The dotted line proportional to t 2.3 is drawn as a guide for the eyes.
We showed that the dynamics of photoinduced nucleation processes is simulated by the present method. When a cluster of molecules are excited to the Franck–Condon state, it induces distortion of adjacent molecules, and the excitation energy is transferred to the other molecules. Once the molecules start to vibrate, the electronic state conversion from |↓r to |↑r takes place and thus photoinduced domain grows. This is the basic scenario of the initial photoinduced nucleation processes where coherently driven molecular distortions play an important role. In fact, the size of the converted domain (diameter) is almost linearly increases as the nucleation proceeds, which shows that energy diffusion is subsidiary in the initial processes. However, as the coherence of molecular distortion states is lost, excitation energy propagation in the system will be dominated by diffusion processes, and hence √ the growth rate will be ∝ t after all. In calculating the time-dependence of the wavefunctions we applied a mean-field approximation in which fluctuations of ur and nˆ r are neglected. If we consider a Gaussian wavepacket on the PESs of a molecule, we understand that the fluctuation of ur is
larger for smaller value of ω . Hence, the effect of the fluctuation is large when the motion of the wavepacket is slow and the transition between PESs also slowly occurs. The validity of the meanfield approximation depends on the relation among the parameters of the model. For example, the value of λ/¯hω is a measure to estimate the validity of the approximation, i.e. as the value of λ/¯hω is large, the mean-field approximation will becomes worse. When multiple wavepackets simultaneously move on the PESs of a single molecule, the fluctuation of ur and nˆ r becomes large, which is the case when the population transfer takes place repeatedly during the nucleation processes. In the present paper, we assume that only a single relevant vibration mode exists in each molecule. However, the dynamics of wavepackets in each molecule is strongly affected by the structure of the PESs. For example, when multiple vibration modes are taken into account, the dynamics of the nucleation processes depends on the topological structure of the intersections of the PESs such as the existence of conical intersections. Hence, ab initio electronicstructure calculations of specific materials are important for a detailed discussion of such material-dependent features of the nucleation processes, and the dynamics calculation in the present paper should be extended by combining with those electronic-structure calculations in the future. We, however, stress that the present calculation method reveals the basic properties of the nucleation dynamics accompanied by the structural change in localized electron system [9–11]. Acknowledgements This work was supported by the Next Generation Super Computing Project, Nanoscience Program, MEXT, Japan, and the numerical calculations were carried out on the computers at the Research Center for Computational Science, National Institutes of Natural Sciences. References [1] S. Koshihara, Y. Takahashi, H. Sakai, Y. Tokura, T. Luty, J. Phys. Chem. B 103 (1999) 2592.
1494
K. Ishida, K. Nasu / Computer Physics Communications 180 (2009) 1489–1494
[2] S. Koshihara, Y. Tokura, K. Takeda, T. Koda, Phys. Rev. B 52 (1995) 6265. [3] A. Mino, Y. Ogawa, S. Koshihara, C. Urano, H. Takagi, Mol. Cryst. Liq. Cryst. 314 (1998) 107. [4] N.O. Moussa, G. Molnár, S. Bonhommeau, A. Zwick, S. Mouri, K. Tanaka, J.A. Real, A. Bousseksou, Phys. Rev., Lett. 94 (2005) 107205. [5] J.S. Costa, P. Guionneau, J.-F. Létard, J. Phys. Conf. Ser. 21 (2005) 67. [6] K. Nasu (Ed.), Photoinduced Phase Transitions, World Scientific, Singapore, 2004. [7] K. Koshino, T. Ogawa, Phys. Rev. B 58 (1998) 14804. [8] H. Mizouchi, K. Nasu, J. Phys. Soc. Jpn. 69 (2000) 1543. [9] K. Ishida, K. Nasu, Phys. Rev. B 76 (2007) 014302. [10] K. Ishida, K. Nasu, Phys. Rev. Lett. 100 (2008) 116403. [11] K. Ishida, K. Nasu, Phys. Rev. B 77 (2008) 214303.
[12] L.D. Landau, Phys. Zts. Sov. 2 (1932) 46. [13] C. Zener, Proc. Roy. Soc. A 137 (1932) 696. [14] C. Zhu, H. Nakamura, J. Chem. Phys. 101 (1994) 4855, and the references cited therein. [15] N. Winter, I. Chomy, J. Vieceli, I. Benjamin, J. Chem. Phys. 119 (2003) 2127. [16] J.C. Tully, in: D.L. Thompson (Ed.), Modern Methods for Multidimensional Dynamics Computations in Chemistry, World Scientific, Singapore, 1998. [17] L. Salem, Science 191 (1976) 822. [18] K. Horikoshi, K. Misawa, R. Lang, J. Chem. Phys. 127 (2007) 054104; K. Horikoshi, K. Misawa, R. Lang, K. Ishida, Opt. Commun. 259 (2006) 723. [19] K.E. Cahill, R.J. Glauber, Phys. Rev. 177 (1969) 1857. [20] K. Ishida, F. Aiga, K. Misawa, J. Chem. Phys. 127 (2007) 194304. [21] K. Boukheddaden, Prog. Theor. Phys. 112 (2004) 205.
Contents lists available at ScienceDirect
Computer Physics Communications www.elsevier.com/locate/cpc
Numerical simulation on the dynamics of photoinduced cooperative phenomena in molecular crystals Kunio Ishida a,∗ , Keiichiro Nasu b a b
Corporate Research and Development Center, Toshiba Corporation, 1 Komukaitoshiba-cho, Saiwai-ku, Kawasaki 212-8582, Japan Institute of Materials Structure Science, KEK, Graduate University for Advanced Study, and CREST JST, 1-1 Oho, Tsukuba, Ibaraki 305-0801, Japan
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 18 October 2007 Received in revised form 20 November 2008 Accepted 25 March 2009 Available online 28 March 2009 PACS: 71.35.Lk 05.45.-a 64.70.Nd 78.20.Bh
We develop a new simulation method to study the dynamics of initial nucleation processes of photoinduced structural change of molecular crystals. In order to describe the nonadiabatic transition in each molecule, we employ a model of localized electrons coupled with a fully quantized phonon mode, and the time-dependent Schrödinger equation for the model is numerically solved. By applying a meanfield approximation in solving the Schrödinger equation, the calculation method is quite efficient on parallel computing systems. We show that coherently driven molecular distortion plays an important role in the successive conversion of electronic states which leads to photoinduced cooperative phenomena. © 2009 Elsevier B.V. All rights reserved.
Keywords: Nonequilibrium phase transition Nonadiabatic transition Parallel computing
1. Introduction Amongst various photoinduced phenomena which have been observed in solids, liquids, and/or molecules, cooperatively induced change of structural, magnetic, or ferroelectric properties [1–5] has particularly attracted our attention. These photoinduced cooperativity are considered to have a common mechanism, and many experimental and/or theoretical studies have been presented to make it clear [6–8]. We have proposed a model of molecular crystals in order to describe the initial nucleation processes of these phenomena [9–11], and showed that the nonadiabaticity of electronic transitions is a key to understand the temporal behavior of quantum-mechanical states. Dynamics of nonadiabatic transitions has been studied since the pioneering works by Landau [12] and Zener [13], and, in particular, the transition rate in general cases was analytically obtained [14]. They mainly focused on the wavefunctions before/after nonadiabatic transition, and hence the time evolution of wavefunctions itself is not within the scope of these studies. On the other hand, the dynamics of nonadiabatic processes has been considered to be important in, for example, photochemical reactions [15], and hence computational methods for nonadiabatic dynamics have been pro-
*
Corresponding author. Tel.: +81 44 549 2119. E-mail address: [email protected] (K. Ishida).
0010-4655/$ – see front matter doi:10.1016/j.cpc.2009.03.013
©
2009 Elsevier B.V. All rights reserved.
posed by many authors [16]. Since, however, these methods are based on the molecular dynamics simulations, they require that the atomic coordinates/momenta are treated as classical variables. Hence, these methods are valid after the decoherence of the atomic degrees of freedom takes place. As a result they mentioned the reaction yield or the absorption rate after various nonequilibrium processes. It has, however, been pointed out that an elementary process of photoinduced cooperative phenomena is consecutive switching of relevant potential energy surfaces (PESs) and that the coherence of the wavefunction of electrons/atoms during electronic transition should be taken into account. In other words, the wavefunction at every moment should be pursued to understand the dynamics of the whole processes, and hence, not only the bifurcation rate of the wavefunction after nonadiabatic transition but also the wavefunction as a function of time is required. This aspect of the nucleation dynamics shows that the conventional methods are not suitable for the theory of photoinduced cooperativity. Furthermore, cooperative action of electrons and phonons consists of an essential part of the photoinduced cooperativity, and thus we require a method to calculate dynamical properties of large aggregates of molecules. In this paper, we show a new computational method to study the initial dynamics of the photoinduced nucleation processes. By introducing a mean-field approximation to take into account the intermolecular interactions, we naturally derived an efficient algorithm to solve the time-dependent Schrödinger equation for a many-body Hamiltonian of localized electrons and
1490
K. Ishida, K. Nasu / Computer Physics Communications 180 (2009) 1489–1494
phonons. This method is also effective for parallel computing and thus suitable for calculations on the recent high performance computing systems. The organization of the paper is as follows: in Section 2 a model of molecular crystals is introduced and the method of calculation is described. In Section 3 the calculated results are shown, and Section 4 is devoted to discussion and conclusions. 2. Model and method 2.1. Model of molecular crystals As we mentioned in Section 1, nonadiabatic transitions between electronic states play a key role in the dynamical aspects of photoinduced cooperative phenomena. In this paper, we focus on the initial dynamics of a photoexcited state in interacting molecules, fully quantizing the relevant vibration modes. However, the dimension of the Hilbert space of the whole system drastically increases by quantizing the atomic variables, which means that the current computational resources are insufficient for numerical calculations on these systems. As a first step to overcome this difficulty, we employ a simplest model which deserves to describe the photoinduced nucleation processes. We consider molecules arrayed on a regular lattice, and all the electrons in the system are assumed to be localized in each molecule. Only two electronic levels coupled with a single vibration mode are taken into account per molecule. Each molecule has crossing diabatic PESs regarding the two electronic levels, which are relevant to the nonadiabatic transitions to drive the nucleation processes. The nonadiabaticity in the dynamics is taken into account via a “spin-flip” interaction between two electronic states, i.e., the interaction strength between two electronic levels are assumed to be a constant. We also apply the harmonic approximation to the vibration modes. We mention that this model is also known as a model to discuss the relaxation dynamics of, e.g., photoisomerization of molecules and studied by many authors [17]. As for the intermolecular interaction, we take into account three types of intermolecular interactions described below: 1. Vibrational interaction which couples the vibration modes of different molecules. It is also responsible for the dispersion of the relevant phonon mode. 2. The dipole–dipole interaction between electrons in excited molecules The interaction strength is taken up to the first order of the molecular distortions. 3. Electron-vibrational interaction which describes the distortion of molecules induced by the excited-state dipole in the adjacent ones. We have pointed out that the nucleation process is triggered only when the second and/or the third types of the above interactions in two- or higher-dimensional systems are taken into account [9–11]. This situation is different from the one-dimensional case in which the vibrational interaction is sufficient to induce the domain growth [7]. In any case, the nucleation process occurs through excitation energy transfer between molecules, and the Hamiltonian in the present study is described by:
H=
p 2 r
r
+
2
+
ω2 ur2 2
+
2h¯ ω3 sur + εh¯ ω + s2 h¯ ω nˆ r + λσxr
αω2 (ur − β nˆ r )(ur − β nˆ r )
r ,r
− V − W (ur + u r ) nˆ r nˆ r ,
(1)
Fig. 1. Schematic view of the model. Circles denote the molecules with two electronic states and a vibrational mode. Adiabatic potential energy surfaces for an individual molecule are shown in the inset.
where pr and ur are the momentum and coordinate operators for the vibration mode of a molecule at site r , respectively. The electronic states at site r are denoted by |↓r (ground state) and |↑r (excited state) and σir (i = x, y , z) are the Pauli matrices which act only on the electronic states of the molecule at site r . nˆ r denotes the density of the electron in |↑r which is rewritten as nˆ r = σzr + 1/2. The second sum which gives the intermolecular interaction is taken over all the pairs on nearest neighbor sites, where the Coulomb interaction between excited state electrons are modified by molecular distortion. The vibrational period of an individual molecule is denoted by T = 2π /ω in the rest of the paper. A schematic view of the present model in the case of square lattice is shown in Fig. 1. This figure shows that the two diabatic PESs for an individual molecule cross each other, and that the nonadiabatic coupling constant λ acts to separate them into two adiabatic PESs. The parameters for the intra/intermolecular interactions are also shown in Fig. 1. The canonical quantization of the vibrational degrees of freedom is performed. The vibronic states are denoted by |nσ r where n is the index of the vibration mode and σ = ↑, ↓, and described by [18–20]
† n a
|n↓r = √r
n!
|n↑r = e
†
|0r ⊗ |↓r ,
s(ar −ar )
(2)
† n a
√r
n!
|0r ⊗ |↑r ,
(3)
†
where ar denotes the creation operator of vibrational quanta of a molecule at site r . The coordinate of the molecule is labeled by r . We mention that the basis set for the vibronic states is composed of the Fock states shown in Ref. [19]. The phonon dispersion relation of the vibration mode is given by
Ω(k) = ω 1 + 2α (cos kx + cos k y ),
(4)
where k = (kx , k y ) denote the reciprocal lattice vector of the square lattice, and the lattice constant is taken to be unity. We note that this Ising-like model is similar to the one to study the thermodynamical properties of the Jahn–Teller effect [21], though the nonequilibrium dynamics of the excited states in the model has not been understood.
K. Ishida, K. Nasu / Computer Physics Communications 180 (2009) 1489–1494
1491
2.2. Method of calculation We obtain the numerical solutions of the time-dependent Schrödinger equation for the Hamiltonian (1) which is written as ih¯
∂
Φ(t ) = H Φ(t ) . ∂t
(5)
When we take into account up to the Mth vibronic states on each PES per molecule, the dimension of the Hilbert space of the whole system D is equal to (2M ) N where N is the total number of the molecules. Since we cannot divide the Hamiltonian into a direct sum of those for each molecule due to the intermolecular interactions, we cannot block-diagonalize the Hamiltonian, i.e. the wavefunction cannot be decomposed to a tensor product of those for each molecule. This makes it quite difficult to solve Eq. (5) numerically, since D > 1015 even for M = 16 and N = 10, for example. However, the aim of the present study is the numerical calculation on the collective phenomena induced in the system, which means that the total number of the molecules N is required to be at least 103 for quantitative discussion. To avoid this difficulty, we apply a mean-field approximation to Eq. (5), i.e., the effects of molecules in the nearest neighbor sites are substituted by their average values. This approximation allows us to decompose the wavefunction of the total system |Φ(t ) into a product of the wavefunctions at each molecule, i.e.,
Φ(t ) = φ(t ) ⊗ φ(t ) ⊗ · · · ⊗ φ(t ) . r r r 1
(6)
N
2
i h¯
Then, each component |φ(t )r obeys the following equation: i h¯
∂
φ(t ) r = ∂t
pr2
ω2 u 2
r + + 2 2
+ λσxr φ(t ) r
+
Fig. 2. A schematic view of the partition of the system for computation. All the circles denote the molecules on a square lattice. The molecular crystal is divided into blocks of molecules by the dashed lines and each block is assigned to a single processor of the computer system. Only the mean-field values (ˆnr and ur ) for the molecules at the perimeter of the blocks (filled circles) are communicated between processors.
dcnr ↑ dt
2h¯ ω3 sur + εh¯ ω + s2 h¯ ω nˆ r
+ −
αω2 (ur − β nˆ r )(ur − β nˆ r )
− V − W (ur + u r ) nˆ r nˆ r φ(t ) r ,
(7)
where the sum is taken over the nearest neighbor sites of r . After canonical quantization of pr and ur , Eq. (7) is rewritten as
† †
∂
φ(t ) r = h¯ ω ar ar + s ar + ar + ε + s2 nˆ r + λσxr φ(t ) r ∂t
+ αω2 γ ar† + ar − β nˆ r ur − βˆnr (8)
By expanding the wavefunction of each molecule by its vibronic states as (10)
n=0 σ =↑,↓
we obtained a set of differential equations for the coefficients cnσ as ih¯
dcnr ↓ dt
= nh¯ ωcnr ↓ + h¯ +
r U nm cm ↑
m
αω2 γ
√
ncnr −1↓ +
n + 1cnr +1↓
r
× u r − βˆnr ,
√
(11)
√
ncnr −1↓ +
√
√
ncnr −1↑ +
√
n + 1cnr +1↓ u r − βˆnr
V − W − u r cnr ↑
min(m,n) i =0
(9)
∞
φ(t ) = cnr σ |nσ r , r
αω2 γ
† U mn = m|e s(a −a) |n s2
where γ = h¯ /2ω . The intermolecular interaction is taken into account through the average value over the decomposed wavefunction and
Xr = r φ(t ) Xr φ(t ) r .
m
n + 1cnr +1↑
ˆnr ,
(12)
where
= e− 2
r
† − V − W γ ar + ar + u r nˆ r ˆnr φ(t ) r ,
− Wγ
r
r
i h¯
√ √
= h¯ ω n + ε + s2 cnr ↑ + s ncnr −1↑ + n + 1cnr +1↓ r + h¯ ωλ U mn cm ↓
√
m!n!
i !(m − i )!(n − i )!
sm+n−2i .
(13)
Eqs. (11) and (12) are suitable for parallel computing, since the wavefunction |φ(t )r itself is not transferred between each other during the calculation. Instead, only the average values of ur and ˆnr are communicated between CPUs, which reduces the amount of data transfer drastically. Hence, our implementation to solve the time-dependent Schrödinger equation is shown as follows: 1. We divide the whole system into blocks of molecules shown in Fig. 2. Each block contains the same number of molecules and is assigned to a CPU or a processor element (PE) of the computer system. 2. The wavefunctions cnr σ for the molecules are stored in the main memory of each processor, and their time evolution is calculated by applying the Runge–Kutta method to Eqs. (11) and (12). 3. The values of ˆnr and ur for the molecules at the perimeter of each block are transferred between processors at each time step. Thus, when a block is composed of n × n molecules, only 4(n − 1) real numbers are communicated per block. By locating these variables at a continuous block of the main memory, we can also reduce the computational cost regarding communication between processors, since the MPI coding is simpler and easier in this case (ref. filled circles in Fig. 2).
1492
K. Ishida, K. Nasu / Computer Physics Communications 180 (2009) 1489–1494
on the implementation of the MPI library and/or the system architecture, and that it varies between 0 and 1. Hence, N 0 ∝ N α where 1/3 α < 1. In any case, N 0 is an increasing function of N as shown in Fig. 3, and thus the number of molecules per block increases as N increases. This feature is also important to consider the simulation of three-dimensional systems, and will be discussed in the next section. 3.2. Dynamics of photoinduced nucleation processes We have pointed out that the population of the excited electronic state |↑r is suitable for understanding the dynamics of the initial nucleation processes [9–11]. We show in Fig. 4 the population of |↑r for 32 × 32 sites defined by
N (r , t ) = Φ(t ) nˆ r Φ(t ) , Fig. 3. Normalized execution time as a function of the number of processors. The system size is 32 × 32, 64 × 64, and 128 × 128.
In this way, we have made it possible to handle more than 104 molecules by the present method, which is sufficient to discuss the cooperative effect induced by photoexcited states. 3. Calculated results 3.1. Scalability of the computational method Although the amount of communication is reduced by the mean-field approximation, the present computation is not completely scalable as in general cases. In the rest of the paper we chose the values of the parameters as: ε = 2.3, s = 1.4, V = 1.1, W = 0.2, α = 0.1, β = 0.2, and λ = 0.2. Although those values are typical for organic molecules as for electron-vibration coupling [18] and the intermolecular Coulomb interaction [17], the other parameters are not easy to determine their values either from theoretical calculations or experimental results. We only mention that the order of magnitude for the parameters is estimated referring to those for typical organic materials. As for the time scale of the photoinduced nucleation, we mention that the vibration period of a single molecule is ∼200 fs, i.e. T ∼ 200 fs. The time step is taken to be T /36000, which is sufficient to obtain the numerical convergence of the results. Fig. 3 shows the execution time of the simulation with 6000 time steps on a 32 × 32, a 64 × 64, and a 128 × 128 lattice. The calculations were performed on the HA8000 system at the Research Center for Computational Science, National Institutes of Natural Sciences. In each plot the computation time is normalized by that for a single processor. As the number of processors P increases, the amount of communication between them increases, and thus the trade-off between the increase of CPU time and the decrease of communication burden results in the existence of the optimal value of the processors which minimizes the total computation time. On the other hand, the amount of communication is dependent on the size of the block introduced in the last section, i.e., the perimeter-to-section ratio in a square block is ∼1/n for an n × n block, and thus the relative weight of communication increases as the size of the block decreases. Hence, the optimal number of processors N 0 depends on the size of the system, and it increases as the system is larger. In the present case N 0 = 8, 16, and 32 for a 32 × 32 lattice, a 64 × 64 lattice, √ and a 128 × 128 lattice, respectively, which shows note that the total amount of communication that N 0 ∝ N. We √ is ∼ N (n − 1)/n2 ∼ N P in the present case, and that the time to Thus, if we ascalculate the time evolution in each block is ∼ N / P . √ sume that the total computation time is ∼ N / P + a( N P )x with a and x being certain √ constants, we obtain x = 2/3, which satisfies the relation N 0 ∝ N. We also note that the value of x depends
(14)
for t = 0, 5T , 10T , and 15T . Initially three molecules which form an v-tromino (see the inset of Fig. 4) are excited to the Franck– Condon state, while the others are set to the ground state. Fig. 4 shows that the number of molecules in the excited electronic state increase surrounding the initially excited molecules. Those molecules will constitute a photoinduced domain observed in many experiments [1–5] and thus the present calculation describes the initial processes of the photoinduced cooperative phenomena, i.e., photoinduced nucleation triggered by an injected excited state. The sum of the excited state population N total (t ) = r N (r , t ) indicates the measure for the growth rate of the photoinduced domain [9]. Fig. 5 shows that N total (t ) slightly decreases for 5T after photoexcitation, and turns to increase after this period. The domain size behaves superlinearly as a function of t after domain growth started. This feature is understood by the picture that the growth of the domain is predominantly driven by coherent motion of molecular distortions rather than diffusion processes. As the vibrational coherence is lost, diffusion process becomes more important and the domain growth √ will slow down to make the radius of the domain increase as ∼ t. Since vibrational coherence survives for a few ps in typical organic molecules [18], the present calculation is valid only in the time range studied in this paper, and the decoherence of the vibrational states should be taken into account to study the growth dynamics of the photoinduced domain in a longer time scale. 4. Discussion and conclusions In this paper we show a computational method to study the dynamics of the initial nucleation processes of photoinduced cooperative phenomena. By treating the intermolecular interactions by the mean-field approximation, we obtain an efficient algorithm to calculate the time-evolution of photoexcited states of more than 104 molecules on parallel computing systems. In this method we divide the crystal into blocks of molecules, and assign each block to one of the processors. When intermolecular interactions are restricted to those between adjacent molecules, the communicated data are only the average values (mean-field) of relevant physical properties at the perimeter of the blocks, and hence the amount of communication in this method is relatively smaller as the size of the block is larger. As a result we found that the optimal number of processors to minimize the computational cost depends on the size of the crystal. This feature is particularly important to consider higher-dimensional (i.e., three-dimensional) cases. In three-dimensional cases we need to take into account typically 106 (= 100 × 100 × 100) molecules or more, and hence the optimal number of processors N 0 should be much larger than that in the present case. Thus, the present method has a remarkable advantage in computation on superparallel computer systems and/or grid computing systems.
K. Ishida, K. Nasu / Computer Physics Communications 180 (2009) 1489–1494
1493
Fig. 4. Population of the excited electronic state N (r , t ) on 32 × 32 lattice for (a) t = 0, (b) t = 5T , and (c) t = 10T , and (d) t = 15T .
Fig. 5. N total (t ) as a function of time. The dotted line proportional to t 2.3 is drawn as a guide for the eyes.
We showed that the dynamics of photoinduced nucleation processes is simulated by the present method. When a cluster of molecules are excited to the Franck–Condon state, it induces distortion of adjacent molecules, and the excitation energy is transferred to the other molecules. Once the molecules start to vibrate, the electronic state conversion from |↓r to |↑r takes place and thus photoinduced domain grows. This is the basic scenario of the initial photoinduced nucleation processes where coherently driven molecular distortions play an important role. In fact, the size of the converted domain (diameter) is almost linearly increases as the nucleation proceeds, which shows that energy diffusion is subsidiary in the initial processes. However, as the coherence of molecular distortion states is lost, excitation energy propagation in the system will be dominated by diffusion processes, and hence √ the growth rate will be ∝ t after all. In calculating the time-dependence of the wavefunctions we applied a mean-field approximation in which fluctuations of ur and nˆ r are neglected. If we consider a Gaussian wavepacket on the PESs of a molecule, we understand that the fluctuation of ur is
larger for smaller value of ω . Hence, the effect of the fluctuation is large when the motion of the wavepacket is slow and the transition between PESs also slowly occurs. The validity of the meanfield approximation depends on the relation among the parameters of the model. For example, the value of λ/¯hω is a measure to estimate the validity of the approximation, i.e. as the value of λ/¯hω is large, the mean-field approximation will becomes worse. When multiple wavepackets simultaneously move on the PESs of a single molecule, the fluctuation of ur and nˆ r becomes large, which is the case when the population transfer takes place repeatedly during the nucleation processes. In the present paper, we assume that only a single relevant vibration mode exists in each molecule. However, the dynamics of wavepackets in each molecule is strongly affected by the structure of the PESs. For example, when multiple vibration modes are taken into account, the dynamics of the nucleation processes depends on the topological structure of the intersections of the PESs such as the existence of conical intersections. Hence, ab initio electronicstructure calculations of specific materials are important for a detailed discussion of such material-dependent features of the nucleation processes, and the dynamics calculation in the present paper should be extended by combining with those electronic-structure calculations in the future. We, however, stress that the present calculation method reveals the basic properties of the nucleation dynamics accompanied by the structural change in localized electron system [9–11]. Acknowledgements This work was supported by the Next Generation Super Computing Project, Nanoscience Program, MEXT, Japan, and the numerical calculations were carried out on the computers at the Research Center for Computational Science, National Institutes of Natural Sciences. References [1] S. Koshihara, Y. Takahashi, H. Sakai, Y. Tokura, T. Luty, J. Phys. Chem. B 103 (1999) 2592.
1494
K. Ishida, K. Nasu / Computer Physics Communications 180 (2009) 1489–1494
[2] S. Koshihara, Y. Tokura, K. Takeda, T. Koda, Phys. Rev. B 52 (1995) 6265. [3] A. Mino, Y. Ogawa, S. Koshihara, C. Urano, H. Takagi, Mol. Cryst. Liq. Cryst. 314 (1998) 107. [4] N.O. Moussa, G. Molnár, S. Bonhommeau, A. Zwick, S. Mouri, K. Tanaka, J.A. Real, A. Bousseksou, Phys. Rev., Lett. 94 (2005) 107205. [5] J.S. Costa, P. Guionneau, J.-F. Létard, J. Phys. Conf. Ser. 21 (2005) 67. [6] K. Nasu (Ed.), Photoinduced Phase Transitions, World Scientific, Singapore, 2004. [7] K. Koshino, T. Ogawa, Phys. Rev. B 58 (1998) 14804. [8] H. Mizouchi, K. Nasu, J. Phys. Soc. Jpn. 69 (2000) 1543. [9] K. Ishida, K. Nasu, Phys. Rev. B 76 (2007) 014302. [10] K. Ishida, K. Nasu, Phys. Rev. Lett. 100 (2008) 116403. [11] K. Ishida, K. Nasu, Phys. Rev. B 77 (2008) 214303.
[12] L.D. Landau, Phys. Zts. Sov. 2 (1932) 46. [13] C. Zener, Proc. Roy. Soc. A 137 (1932) 696. [14] C. Zhu, H. Nakamura, J. Chem. Phys. 101 (1994) 4855, and the references cited therein. [15] N. Winter, I. Chomy, J. Vieceli, I. Benjamin, J. Chem. Phys. 119 (2003) 2127. [16] J.C. Tully, in: D.L. Thompson (Ed.), Modern Methods for Multidimensional Dynamics Computations in Chemistry, World Scientific, Singapore, 1998. [17] L. Salem, Science 191 (1976) 822. [18] K. Horikoshi, K. Misawa, R. Lang, J. Chem. Phys. 127 (2007) 054104; K. Horikoshi, K. Misawa, R. Lang, K. Ishida, Opt. Commun. 259 (2006) 723. [19] K.E. Cahill, R.J. Glauber, Phys. Rev. 177 (1969) 1857. [20] K. Ishida, F. Aiga, K. Misawa, J. Chem. Phys. 127 (2007) 194304. [21] K. Boukheddaden, Prog. Theor. Phys. 112 (2004) 205.